拆解Botsch经典算法：手写半边结构，一步步实现Isotropic Remeshing（附C++代码）

拆解Botsch经典算法手写半边结构一步步实现Isotropic Remeshing附C代码在计算机图形学领域网格优化一直是核心挑战之一。想象一下当你从3D扫描仪获取一个粗糙的模型或者从CAD软件导出一个高多边形模型时如何将其转化为均匀、规则的三角网格这正是各向同性网格重建(Isotropic Remeshing)要解决的问题。不同于简单的网格简化这项技术能在保持几何特征的同时生成质量均匀的三角形为后续的模拟、渲染等操作奠定基础。2004年Botsch等人提出的算法因其简洁高效而成为经典。本文将带你从零开始实现这一算法不仅理解其数学原理更重要的是掌握如何用C构建核心数据结构——半边结构(Half-edge)并处理算法实现中的各种边界情况。我们不会依赖CGAL等现成库而是通过手写代码来深入理解每个细节。1. 理解半边数据结构网格操作的基石在开始算法实现前我们需要一个能够高效表示网格关系的数据结构。半边结构因其对邻接关系的完整表达而成为图形学中的标准选择。1.1 半边结构的基本组成半边结构将每条物理边拆分为两条方向相反的半边(half-edge)这种设计使得遍历操作变得直观高效。以下是我们的C定义struct Vertex { glm::vec3 position; // 顶点位置 HalfEdge* edge; // 指向任意一条以该点为起点的半边 }; struct Face { HalfEdge* edge; // 指向任意一条属于该面的半边 }; struct HalfEdge { Vertex* origin; // 半边起点 HalfEdge* twin; // 反向的半边 HalfEdge* next; // 同一面中的下一条半边 Face* face; // 所属面 };这种表示法的一个关键优势是能够高效查询各种邻接关系。例如要找到一个顶点的所有邻接点只需遍历从该顶点出发的所有半边std::vectorVertex* getAdjacentVertices(Vertex* v) { std::vectorVertex* adjVertices; HalfEdge* start v-edge; HalfEdge* he start; do { adjVertices.push_back(he-twin-origin); he he-twin-next; } while (he ! start); return adjVertices; }1.2 构建初始网格从OBJ或OFF文件加载网格后我们需要构建半边结构。这个过程需要特别注意边界情况的处理void buildHalfEdgeStructure(Mesh mesh, const std::vectorglm::vec3 vertices, const std::vectorstd::vectorint faces) { // 创建顶点 for (const auto v : vertices) { mesh.vertices.push_back(new Vertex{v, nullptr}); } // 创建面与半边 std::mapstd::pairint, int, HalfEdge* edgeMap; // 用于查找twin for (const auto faceIndices : faces) { Face* f new Face{nullptr}; mesh.faces.push_back(f); std::vectorHalfEdge* faceEdges; int n faceIndices.size(); // 创建该面的所有半边 for (int i 0; i n; i) { int j (i 1) % n; HalfEdge* he new HalfEdge{ mesh.vertices[faceIndices[i]], nullptr, nullptr, f }; // 设置顶点指向的第一条半边 if (mesh.vertices[faceIndices[i]]-edge nullptr) { mesh.vertices[faceIndices[i]]-edge he; } faceEdges.push_back(he); // 尝试查找twin auto key std::make_pair(faceIndices[j], faceIndices[i]); if (edgeMap.count(key)) { he-twin edgeMap[key]; edgeMap[key]-twin he; edgeMap.erase(key); } else { edgeMap[std::make_pair(faceIndices[i], faceIndices[j])] he; } } // 设置next指针 for (int i 0; i n; i) { faceEdges[i]-next faceEdges[(i 1) % n]; } f-edge faceEdges[0]; } // 处理边界边没有twin的半边 // ... 边界处理代码 ... }注意实际实现中需要特别处理边界边这些边只有一个半边其twin指针为nullptr。正确处理边界对后续的网格操作至关重要。2. 实现核心操作Split、Collapse和FlipBotsch算法的核心在于三个基本操作分割(Split)、塌缩(Collapse)和翻转(Flip)。这些操作需要在不破坏网格拓扑的前提下改变其连接关系。2.1 边分割(Split)操作边分割将一条过长的边分成两部分插入一个新顶点。实现时需要注意维护所有指针的正确性Vertex* splitEdge(Mesh mesh, HalfEdge* he, float splitPosition 0.5f) { // 1. 创建新顶点 Vertex* newVert new Vertex{ glm::mix(he-origin-position, he-twin-origin-position, splitPosition), nullptr }; // 2. 创建四条新半边 HalfEdge* he1 new HalfEdge{he-origin, nullptr, nullptr, he-face}; HalfEdge* he2 new HalfEdge{newVert, nullptr, nullptr, he-face}; HalfEdge* he3 new HalfEdge{newVert, nullptr, nullptr, he-twin-face}; HalfEdge* he4 new HalfEdge{he-twin-origin, nullptr, nullptr, he-twin-face}; // 3. 设置twin关系 he1-twin he4; he4-twin he1; he2-twin he3; he3-twin he2; // 4. 设置next指针 he2-next he-next; he1-next he2; he4-next he-twin-next; he3-next he4; // 5. 更新面的半边指针 he-face-edge he1; he-twin-face-edge he3; // 6. 更新顶点指向的半边 newVert-edge he2; // 7. 更新邻接半边指向 // ... 需要遍历并更新相关半边的next指针 ... // 8. 删除原来的半边 delete he; delete he-twin; mesh.vertices.push_back(newVert); return newVert; }提示splitPosition参数控制分割点的位置默认0.5表示中点分割但可以根据需要调整。2.2 边塌缩(Collapse)操作边塌缩将一条过短的边合并为一个顶点是分割的逆操作。实现时需特别注意避免产生二度点和重边bool collapseEdge(Mesh mesh, HalfEdge* he) { Vertex* v1 he-origin; Vertex* v2 he-twin-origin; // 1. 检查是否会形成二度点或重边 if (wouldCreateDegeneracy(he)) { return false; // 不能执行塌缩 } // 2. 计算新顶点位置可根据需要选择不同策略 glm::vec3 newPos 0.5f * (v1-position v2-position); // 3. 创建新顶点 Vertex* newVert new Vertex{newPos, nullptr}; // 4. 重新连接所有相关半边 std::vectorHalfEdge* incidentEdges getIncidentEdges(v1, v2); for (HalfEdge* edge : incidentEdges) { if (edge-origin v1 || edge-origin v2) { edge-origin newVert; } } // 5. 更新面的指针 // ... 需要处理所有受影响的面的edge指针 ... // 6. 删除被塌缩的边和相关顶点 deleteHalfEdgePair(he); deleteVertex(v1); deleteVertex(v2); mesh.vertices.push_back(newVert); return true; } bool wouldCreateDegeneracy(HalfEdge* he) { // 检查是否会导致任何顶点度数小于3形成二度点 std::setVertex* neighbors; // 收集v1的邻接点不包括v2 HalfEdge* start he-next; HalfEdge* current start; do { if (current-twin-origin ! he-twin-origin) { neighbors.insert(current-twin-origin); } current current-twin-next; } while (current ! start); // 收集v2的邻接点不包括v1 start he-twin-next; current start; do { if (current-twin-origin ! he-origin) { if (neighbors.count(current-twin-origin)) { // 发现重边 return true; } neighbors.insert(current-twin-origin); } current current-twin-next; } while (current ! start); // 检查是否有顶点度数会小于3 for (Vertex* v : neighbors) { if (getVertexDegree(v) - 1 3) { return true; } } return false; }2.3 边翻转(Flip)操作边翻转通过改变边的连接方式来优化顶点度数是Delaunay三角剖分的核心操作bool flipEdge(Mesh mesh, HalfEdge* he) { // 只能翻转内部边非边界边 if (he-twin nullptr || he-twin-face nullptr) { return false; } // 获取四个相关顶点 Vertex* a he-origin; Vertex* b he-twin-origin; Vertex* c he-next-twin-origin; Vertex* d he-twin-next-twin-origin; // 检查翻转是否会导致非流形或退化 if (isEdgeBetween(c, d)) { return false; } // 重新连接半边 he-origin c; he-twin-origin d; // 更新next指针 HalfEdge* he1 he-next; HalfEdge* he2 he1-next; HalfEdge* he3 he-twin-next; HalfEdge* he4 he3-next; he-next he4; he4-next he1; he1-next he; he-twin-next he2; he2-next he3; he3-next he-twin; // 更新面的半边指针 he-face-edge he; he-twin-face-edge he-twin; return true; }3. 实现各向同性网格重建算法有了上述基本操作我们现在可以组装完整的各向同性网格重建算法。算法分为四个主要步骤每个步骤都需要多次迭代以达到理想效果。3.1 计算目标边长目标边长L决定了最终网格的密度通常有以下几种计算方式平均边长法计算当前网格所有边的平均长度面积开方法根据总三角形面积和期望三角形数量推算用户指定根据应用需求直接指定float computeTargetEdgeLength(const Mesh mesh, Method method AVERAGE) { float totalLength 0.0f; int edgeCount 0; // 遍历所有边避免重复计算 std::setstd::pairVertex*, Vertex* processedEdges; for (HalfEdge* he : mesh.halfEdges) { Vertex* v1 he-origin; Vertex* v2 he-twin-origin; if (processedEdges.count(std::make_pair(v1, v2)) || processedEdges.count(std::make_pair(v2, v1))) { continue; } processedEdges.insert(std::make_pair(v1, v2)); totalLength glm::length(v1-position - v2-position); edgeCount; } switch (method) { case AVERAGE: return totalLength / edgeCount; case AREA_BASED: { float totalArea computeTotalArea(mesh); int desiredFaces mesh.faces.size(); // 可根据需要调整 return sqrtf((2.0f * totalArea) / (desiredFaces * sqrtf(3.0f))); } default: return 0.1f; // 默认值 } }3.2 主算法流程void isotropicRemeshing(Mesh mesh, int iterations 3) { for (int i 0; i iterations; i) { float L computeTargetEdgeLength(mesh); // 第一步分割过长的边4/3L splitLongEdges(mesh, 4.0f/3.0f * L); // 第二步塌缩过短的边4/5L collapseShortEdges(mesh, 4.0f/5.0f * L); // 第三步边翻转以优化顶点度数 flipEdgesForDegreeOptimization(mesh); // 第四步切平面投影 projectToTangentPlane(mesh); } } void splitLongEdges(Mesh mesh, float threshold) { std::vectorHalfEdge* edgesToSplit; // 收集所有需要分割的边 std::setstd::pairVertex*, Vertex* processedEdges; for (HalfEdge* he : mesh.halfEdges) { Vertex* v1 he-origin; Vertex* v2 he-twin-origin; if (processedEdges.count(std::make_pair(v1, v2)) || processedEdges.count(std::make_pair(v2, v1))) { continue; } processedEdges.insert(std::make_pair(v1, v2)); float length glm::length(v1-position - v2-position); if (length threshold) { edgesToSplit.push_back(he); } } // 执行分割 for (HalfEdge* he : edgesToSplit) { splitEdge(mesh, he); } } void collapseShortEdges(Mesh mesh, float threshold) { std::vectorHalfEdge* edgesToCollapse; // 收集所有需要塌缩的边 std::setstd::pairVertex*, Vertex* processedEdges; for (HalfEdge* he : mesh.halfEdges) { Vertex* v1 he-origin; Vertex* v2 he-twin-origin; if (processedEdges.count(std::make_pair(v1, v2)) || processedEdges.count(std::make_pair(v2, v1))) { continue; } processedEdges.insert(std::make_pair(v1, v2)); float length glm::length(v1-position - v2-position); if (length threshold) { edgesToCollapse.push_back(he); } } // 尝试塌缩可能因拓扑约束而失败 for (HalfEdge* he : edgesToCollapse) { collapseEdge(mesh, he); } } void flipEdgesForDegreeOptimization(Mesh mesh) { const int targetDegreeInternal 6; const int targetDegreeBoundary 4; bool changed; do { changed false; std::setstd::pairVertex*, Vertex* processedEdges; for (HalfEdge* he : mesh.halfEdges) { Vertex* v1 he-origin; Vertex* v2 he-twin-origin; if (processedEdges.count(std::make_pair(v1, v2)) || processedEdges.count(std::make_pair(v2, v1))) { continue; } processedEdges.insert(std::make_pair(v1, v2)); // 计算当前边翻转前后的顶点度数差异 int deviationBefore computeDegreeDeviation(he, targetDegreeInternal, targetDegreeBoundary); int deviationAfter computeDegreeDeviationIfFlipped(he, targetDegreeInternal, targetDegreeBoundary); if (deviationAfter deviationBefore) { if (flipEdge(mesh, he)) { changed true; } } } } while (changed); } int computeDegreeDeviation(HalfEdge* he, int targetInternal, int targetBoundary) { Vertex* v1 he-origin; Vertex* v2 he-twin-origin; Vertex* v3 he-next-twin-origin; Vertex* v4 he-twin-next-twin-origin; int dev 0; dev abs(getVertexDegree(v1) - (isBoundaryVertex(v1) ? targetBoundary : targetInternal)); dev abs(getVertexDegree(v2) - (isBoundaryVertex(v2) ? targetBoundary : targetInternal)); dev abs(getVertexDegree(v3) - (isBoundaryVertex(v3) ? targetBoundary : targetInternal)); dev abs(getVertexDegree(v4) - (isBoundaryVertex(v4) ? targetBoundary : targetInternal)); return dev; }3.3 切平面投影在分割操作后新顶点的位置需要投影到局部切平面上以保持几何形状void projectToTangentPlane(Mesh mesh) { for (Vertex* v : mesh.vertices) { if (v-isNew) { // 只处理新创建的顶点 // 计算顶点法线平均邻接面法线 glm::vec3 normal computeVertexNormal(v); // 计算切平面投影 glm::vec3 centroid computeNeighborhoodCentroid(v); glm::vec3 proj v-position - glm::dot(v-position - centroid, normal) * normal; v-position proj; v-isNew false; } } } glm::vec3 computeVertexNormal(Vertex* v) { glm::vec3 normal(0.0f); HalfEdge* start v-edge; HalfEdge* he start; do { if (he-face ! nullptr) { // 忽略边界 glm::vec3 faceNormal computeFaceNormal(he-face); normal faceNormal; } he he-twin-next; } while (he ! start); return glm::normalize(normal); } glm::vec3 computeNeighborhoodCentroid(Vertex* v) { glm::vec3 centroid(0.0f); int count 0; HalfEdge* start v-edge; HalfEdge* he start; do { centroid he-twin-origin-position; count; he he-twin-next; } while (he ! start); return centroid / float(count); }4. 性能优化与实用技巧实现基本功能后我们需要考虑算法的效率和稳定性。以下是几个关键优化方向4.1 空间索引加速对于大型网格暴力遍历所有边效率低下。使用空间索引可以显著加速邻接查询class SpatialIndex { public: void build(const std::vectorVertex* vertices) { for (Vertex* v : vertices) { glm::vec3 pos v-position; int x static_castint(pos.x / cellSize); int y static_castint(pos.y / cellSize); int z static_castint(pos.z / cellSize); grid[std::make_tuple(x, y, z)].push_back(v); } } std::vectorVertex* queryNeighbors(glm::vec3 pos, float radius) const { std::vectorVertex* result; int minX static_castint((pos.x - radius) / cellSize); int maxX static_castint((pos.x radius) / cellSize); // 类似计算minY,maxY,minZ,maxZ... for (int x minX; x maxX; x) { for (int y minY; y maxY; y) { for (int z minZ; z maxZ; z) { auto it grid.find(std::make_tuple(x, y, z)); if (it ! grid.end()) { for (Vertex* v : it-second) { if (glm::distance(v-position, pos) radius) { result.push_back(v); } } } } } } return result; } private: float cellSize 0.1f; // 根据场景调整 std::mapstd::tupleint, int, int, std::vectorVertex* grid; };4.2 并行计算算法的多个步骤可以并行化特别是边分割和塌缩操作void parallelSplitLongEdges(Mesh mesh, float threshold) { std::vectorHalfEdge* edgesToSplit; // ... 收集需要分割的边 ... // 使用OpenMP并行分割 #pragma omp parallel for for (size_t i 0; i edgesToSplit.size(); i) { splitEdge(mesh, edgesToSplit[i]); } // 注意并行操作后需要重建空间索引 }4.3 边界处理技巧边界边需要特殊处理以避免退化bool isBoundaryVertex(Vertex* v) { HalfEdge* start v-edge; HalfEdge* he start; do { if (he-twin nullptr || he-twin-face nullptr) { return true; } he he-twin-next; } while (he ! start); return false; } void handleBoundaryEdges(Mesh mesh) { // 识别所有边界边 std::vectorHalfEdge* boundaryEdges; for (HalfEdge* he : mesh.halfEdges) { if (he-twin nullptr || he-twin-face nullptr) { boundaryEdges.push_back(he); } } // 对边界边应用特殊规则 for (HalfEdge* he : boundaryEdges) { // 例如不翻转边界边 // 或者边界边的目标度数设为4 } }4.4 自适应细分策略根据不同区域的曲率自适应调整细分程度void adaptiveRemeshing(Mesh mesh) { // 计算每个面的曲率 std::unordered_mapFace*, float faceCurvature; for (Face* f : mesh.faces) { faceCurvature[f] computeFaceCurvature(f); } // 根据曲率调整目标边长 for (HalfEdge* he : mesh.halfEdges) { Face* f1 he-face; Face* f2 he-twin-face; float curvature std::max(faceCurvature[f1], f2 ? faceCurvature[f2] : 0.0f); // 曲率大的区域使用更小的目标边长 float adaptiveL baseL * (1.0f - 0.5f * curvature); // 使用adaptiveL而不是全局L来决定是否分割/塌缩 // ... } }