C++实现Delaunay三角网：Bowyer-Watson增量法核心解析

简介C实现的Delaunay三角网算法项目面向计算机图形学、地理信息系统和有限元分析等领域的开发者与学习者用于解决二维点集的高质量三角剖分问题。资源包共42个文件压缩包大小2.66MB包含6个头文件、5个C源文件、工程配置与资源文件.dsp/.dsw/.rc以及编译生成的exe、obj、pdb等既适合阅读源码理解算法也可直接运行程序观察结果。目前已有1882人学习下载。项目围绕Delaunay三角网的空圆特性展开详细实现了点与三角形数据结构设计、逐点插入构建流程、空圆检测以及三角网优化等核心环节代码中体现了分治策略和扫线算法的应用思路能够帮助读者深入掌握增量式构网原理与C实现细节。对照源码调试还可以直观观察三角形动态生成与调整过程验证算法正确性并进一步应用于地形建模、网格生成等实际场景。1. Delaunay三角网算法为什么三维建模和路径规划都绕不开它如果你处理过点云重建、地形网格生成或者写过有限元前处理一定听过这个名字Delaunay三角网算法。用C实现它不是简单地“把点连成三角形”就完事而是要让每个三角形的外接圆内不包含其他任何点这个性质被称为空圆特性。正是这个特性让Delaunay三角网在同类三角剖分中长得最“匀称”不会出现又长又扁的狭长三角形而这恰恰是网格质量的生命线。这篇文章面向的是已经能写C、但没深入接触过计算几何的开发者。我会从算法选型讲到数据结构再到核心实现和边角处理把每一步的参数选择和翻车原因都摊开讲。读完你会有一套能跑的Bowyer-Watson增量法实现也清楚它为什么在动态点集场景下比分治法更顺手。先给结论Delaunay三角网不是某一种唯一算法而是一族满足空圆特性的三角剖分总称选对实现方式比背下算法名更重要。2. 先选对算法家族Bowyer-Watson增量法到底赢在哪2.1 三种主流Delaunay实现路线对比常见的Delaunay三角网算法有三条路线分治法、扫描线法和增量插入法。分治法Divide and Conquer由Guibas和Stolfi在1985年提出把点集按x坐标排序后递归分割再合并两个子三角网理论时间复杂度O(n log n)是最快的确定性算法。扫描线法则适合处理带约束边的情况但在通用点集上实现复杂度偏高边界缝合逻辑容易出错。增量插入法里最有名的是Bowyer-Watson算法思路极其直白逐个把点插入已有的三角网中找到所有外接圆包含新点的三角形删掉它们形成一个“空洞”然后把新点与空洞边界的每条边连接成新三角形。这个过程摊还时间复杂度接近O(n log n)虽然平均表现略逊于分治法但胜在实现直观、支持动态加点——这在实时场景里是刚需。三种路线我实际都用C写过说说感受分治法第一次写时合并步骤的“两条基线的选择与递归”几乎让我怀疑是玄学调试两天都找不到是左右子网哪条边没连对扫描线法在带洞多边形内剖分时洞的环方向要自己维护很容易翻车。Bowyer-Watson从超三角形开始逐点插入每步只涉及局部三角形出错时打印当前三角形集合就能定位这个对开发调试的友好度是文档里看不出来的。2.2 为什么实际项目默认选Bowyer-Watson接到的项目需求大概是三类一次性给全部点生成网格点分批到达需要增量更新给定三角形后还要动态删点。第一种场景分治法稍快但差距在十万点级别才明显第二种场景只有增量法能处理分治法每次都要全量重算第三种场景需要在Delaunay基础上再维护拓扑结构增量法的局部修改天然契合。另一个现实因素是内存访问模式。分治法递归合并时需要频繁读取左右子网的半边结构缓存命中率低Bowyer-Watson的操作集中在“当前插入点附近的三角形”这一小片区域用vector存三角形、用索引访问对CPU缓存很友好。我负责过的某个商业网格生成模块最初用的分治法后来改成Bowyer-Watson配合空间哈希生成百万点网格的时间反而降了20%原因就是内存局部性更好。2.3 确定实现目标最小可行产物的功能边界在动手写第一行代码前我建议先想清楚两个问题要支持的维度是二维还是三维曲面是否需要输出拓扑信息每个三角形的邻居、每条边的左右三角形如果只是生成STL网格用于渲染只需要顶点坐标和三角形索引就够但如果要做有限元网格或路径规划寻路网格邻居关系必不可少。我一般从二维Bowyer-Watson入手输出三个数组顶点数组、三角形数组每个三角形存三个顶点索引、邻居数组每个三角形存三个邻居索引-1表示边界边。这样的结构既是后面加约束边的基础也能直接用于可视化调试。没错优先做二维是因为三维曲面的Delaunay三角剖分要处理投影平面选择、法向量一致性等问题复杂度是二维的好几倍先跑通二维再拓展是性价比最高的路线。3. 数据结构先行C里怎么设计Point、Triangle和邻接关系3.1 从最小结构开始顶点、边、三角形的取舍定义数据结构前想清楚是存“边”还是存“三角形上的顶点组合”。最直观的做法是Triangle结构里存三个顶点索引和三个邻居索引边的信息通过三角形的边隐式表示。下面是我常用的一套定义// point.h #pragma once struct Point { double x, y; Point(double x_ 0.0, double y_ 0.0) : x(x_), y(y_) {} }; // triangle.h #pragma once #include point.h struct Triangle { int v[3]; // 三个顶点索引逆时针排列 int neighbors[3]; // 邻居三角形索引neighbors[i] 对应 v[i] 的对边 bool deleted; // 标记是否被删除Bowyer-Watson中临时标记 Triangle() { v[0] v[1] v[2] -1; neighbors[0] neighbors[1] neighbors[2] -1; deleted false; } };这里的关键决策是不单独定义Edge结构而是用“顶点对 所属三角形”来表示边。边的查找操作统一为“给定两个顶点索引找到共享这条边的两个三角形”这个操作会频繁出现写成单独函数。邻居数组的下标规则是neighbors[i]存储的是与顶点v[i]相对的那条边另一侧的三角形也就是不包含v[i]的那个邻居。规则不统一边翻转和空洞边界查找全乱套这是踩坑后得出的血泪经验。3.2 为什么用索引而不是指针内存复用和调试友好性一开始我习惯用指针Triangle*互相指向写起来自然。但增量插入会不断创建和删除三角形用指针管理需要手动delete而且三角网被切割时旧三角形会被新三角形替代指针悬空问题防不胜防。改用整数索引后所有三角形统一存在一个vector 里删除用标记deleted跳过空间换时间避免频繁vector扩容导致索引失效。还有一个实际好处调试时能直接打印“三角形12的邻居是7和34”这类信息用指针打印出来是一串十六进制地址看多了头痛。索引方案配合断点能看到完整的拓扑演变过程这对理解Bowyer-Watson每一轮插入至关重要。如果你需要极致性能预先reserve足够的容量避免扩容即可索引方案的内存开销只是每多一个三角形多了三个int完全可以接受。3.3 邻接关系的建立初始化与增量维护增量插入新点时原来包含新点的旧三角形被删除新三角形被创建。新旧交接的关键是正确构建新三角形的邻居关系。具体来说找出空洞的边界边每条边连接两个相邻的边界三角形新三角形与它们互为邻居同时新三角形之间共用一个顶点也需要互相设置邻居。void buildNeighborRelations(std::vectorTriangle triangles, const std::vectorint boundaryEdges, int newPointIndex) { // boundaryEdges: 每对相邻顶点代表一条空洞边界边 // 每条边界边创建一个新三角形并与边界两侧的旧邻居互设关系 }这段代码的核心是两步第一遍历空洞边界边为每条边创建新三角形记为新三角形与空洞外部那个旧三角形互为邻居第二对相邻两条边界边创建的新三角形它们共享一个顶点和一条公共边需手动设置邻居关系这步容易漏漏掉一条边就会导致后续空圆检测访问越界或死循环。我的习惯是每创建完一个三角形先打印它的三个邻居排查可见的缺口比最后一次性检查高效得多。4. 手写Bowyer-Watson核心代码超三角形、逐点插入、边翻转4.1 构造超三角形覆盖全部点的初始外包三角形Bowyer-Watson需要一个初始三角形完全覆盖所有点常见做法是把点集的包围盒放大后取外接三角形。超三角形必须足够大否则边缘点的插入会破坏它又不能太大因为超大三角形会引入数值精度问题。经验法则是把包围盒的两个方向尺寸都放大100倍然后以包围盒对角点构造等腰直角三角形。计算超三角形用下面代码#include vector #include limits #include algorithm struct Triangle { int v[3]; int neighbors[3]; bool deleted; }; // 输入点集返回超三角形的三个顶点索引 std::vectorint computeSuperTriangle(const std::vectorPoint points) { double minX std::numeric_limitsdouble::max(); double maxX std::numeric_limitsdouble::lowest(); double minY maxX, maxY minX; for (const auto p : points) { minX std::min(minX, p.x); maxX std::max(maxX, p.x); minY std::min(minY, p.y); maxY std::max(maxY, p.y); } double dx (maxX - minX) * 100.0; double dy (maxY - minY) * 100.0; double midX (minX maxX) / 2.0; double midY (minY maxY) / 2.0; std::vectorint super(3); // 三个顶点正下方、左上、右上形成覆盖整个点集的三角形 super[0] addPoint(Point(midX - dx, midY - dy * 2.0)); super[1] addPoint(Point(midX, midY dy)); super[2] addPoint(Point(midX dx * 2.0, midY - dy)); return super; }放大100倍这个倍数不是拍脑袋放大倍数太小点集边缘点可能落在超三角形外导致外接圆检测漏判放大倍数太大超三角形三个顶点坐标差异过大在计算外接圆圆心和半径时浮点误差会被放大。100倍在常见坐标尺度几百到几万下足够安全但如果你的点集本身跨尺度很大比如坐标在0.001到1000000之间先做归一化再放大比单纯调倍数更靠谱。4.2 逐点插入主循环定位目标三角形与空洞扩展插入每个点时最核心的步骤是找到所有“外接圆包含该点”的三角形。最朴素的做法是遍历所有三角形但每次插入都全量扫描时间复杂度O(n²)点多了就卡死。我一般用一个简单的空间网格加速定位把当前三角网所在区域划分成固定大小的格子维护每个格子覆盖的三角形列表插入新点时先定位它所在的格子只检测格子内的三角形命中后再递归扩展相邻格子从而缩小搜索范围。空洞扩展的终止条件是关键——当某个三角形的外接圆包含新点但其所有邻居的外接圆都不包含新点时空洞就停止扩展了。这个边界条件写错会导致空洞扩大或缩小生成的三角网不满足Delaunay性质。void insertPoint(std::vectorPoint points, std::vectorTriangle triangles, int pointIndex, SpatialGrid grid) { // 1. 找到外接圆包含新点的初始三角形 int startTri grid.findContainingTriangle(points[pointIndex]); // 2. 用栈进行深度优先搜索收集所有需要删除的三角形 std::vectorint toDelete; std::vectorint boundaryEdges; std::stackint stack; stack.push(startTri); while (!stack.empty()) { int tIdx stack.top(); stack.pop(); Triangle t triangles[tIdx]; if (t.deleted) continue; if (circumcircleContains(t, points[pointIndex])) { t.deleted true; // 延迟删除标记后统一清理 toDelete.push_back(tIdx); // 收集边界边该三角形与其他未删除三角形的公共边 for (int i 0; i 3; i) { int nb t.neighbors[i]; if (nb -1 || triangles[nb].deleted) continue; if (circumcircleContains(triangles[nb], points[pointIndex])) { stack.push(nb); } else { // 这条边是空洞边界 int vA t.v[i]; int vB t.v[(i 1) % 3]; boundaryEdges.push_back(vA); boundaryEdges.push_back(vB); } } } } // 3. 根据边界边创建新三角形并建立邻居关系 buildTrianglesFromBoundary(boundaryEdges, pointIndex, points, triangles); }这里用“栈标记”代替递归防止点集规模大时递归深度过深导致栈溢出。边界边的收集用的是“相对邻居引理”一个三角形外接圆包含新点它的邻居不包含则这两个三角形共享的那条边属于最终空洞的边界。边界边收集完成后新三角形的邻居设置逻辑就比较清晰了每个新三角形与空洞边界的另一侧三角形互为邻居新三角形之间共享两个顶点则互为邻居注意此时新旧三角形的邻居数组要双向设置否则后续翻转会找不到对面的三角形。4.3 边翻转让三角网满足Delaunay性质的局部修正Bowyer-Watson算法本身保证插入过程完成后三角网是Delaunay的但前提是外接圆检测和空洞扩展逻辑严格正确。如果点集里有四点共圆等退化情况插入后可能会残留违反空圆特性的边。边翻转是全局修正手段遍历三角形每条内部边检查共享这条边的两个三角形若其中一个外接圆包含另一个三角形的第三个顶点则交换四边形的对角线。翻转后新形成的边依然需要继续检查所以用栈保存所有可能受影响的边直到找不到需要翻转的边为止通常在退化点集上迭代次数有限不会死循环。void edgeFlip(std::vectorPoint points, std::vectorTriangle triangles) { std::stackstd::pairint, int edges; // 边用 (三角形索引, 顶点索引) 表示 for (size_t i 0; i triangles.size(); i) { if (triangles[i].deleted) continue; for (int j 0; j 3; j) edges.push({static_castint(i), j}); } while (!edges.empty()) { auto [triIdx, edgeIdx] edges.top(); edges.pop(); Triangle t triangles[triIdx]; if (t.deleted) continue; int nbIdx t.neighbors[edgeIdx]; if (nbIdx -1) continue; // 边界边不翻转 Triangle nb triangles[nbIdx]; if (nb.deleted) continue; // 获取共享边的两个端点以及两侧的第三顶点 int vA t.v[edgeIdx]; int vB t.v[(edgeIdx 1) % 3]; int vC t.v[(edgeIdx 2) % 3]; // t 的第三顶点 // 找到 nb 中不同于 vA、vB 的顶点 int vD -1; for (int k 0; k 3; k) { if (nb.v[k] ! vA nb.v[k] ! vB) { vD nb.v[k]; break; } } if (circumcircleContains(points[vA], points[vB], points[vC], points[vD])) { // 需要翻转将公共边 (vA, vB) 换为 (vC, vD) // 更新 t 和 nb 的顶点以及邻居关系 flipDiagonal(points, triangles, triIdx, nbIdx, vA, vB, vC, vD); // 受影响的边重新进栈检查 for (int k 0; k 3; k) { edges.push({triIdx, k}); edges.push({nbIdx, k}); } } } }这段逻辑的注意点翻转操作会同时修改两个三角形的顶点顺序和邻居关系顺序错了会破坏逆时针定向影响后续外接圆检测的符号判断。我通常在翻转后调用一个verifyOrientation函数检查每个三角形的有向面积符号是否一致不一致就直接断言失败。翻转本身是局部操作收敛速度很快实测十万点集翻转次数在几千到几万之间性能开销可接受。4.4 删除超三角形清理辅助边的最后一步所有点插入完成后把所有包含超三角形顶点的三角形标记为deleted然后统一从vector中移除或压缩。这里最容易踩坑的是“包含超三角形顶点的三角形”和“与超三角形相邻的三角形”混淆。包含超三角形顶点的三角形一定在三角网最外层全部删除后剩下的就是有效三角网但如果只删除“相邻”的三角形会留下一些跨越超三角形边界的畸形三角形。void removeSuperTriangle(std::vectorTriangle triangles, const std::vectorint superVertices) { // 先标记再压缩避免遍历时vector元素移动 for (auto t : triangles) { for (int v : t.v) { for (int sv : superVertices) { if (v sv) { t.deleted true; break; } } } } std::erase_if(triangles, [](const Triangle t) { // C20后可用erase_ifC11用vector::erase remove_if 配合 return t.deleted; }); }注意删除后邻居数组里可能还有指向已删除三角形的索引直接访问会拿到无效数据。解决方法是删除前先对每个剩余三角形做一次邻居校验把指向已删除三角形的邻居索引置为-1。如果不做这步后续任何遍历邻居的算法都会读到悬空索引轻则越界重则默默算错结果。这也是为什么我在最初设计数据结构时坚持deleted标记而不是物理删除压缩在最后统一做中间过程所有算法都跳过deleted标记的三角形即可。5. 踩坑实录浮点精度、退化点集、索引混乱的四个翻车现场5.1 四点共圆导致的外接圆检测“薛定谔结果”现象对一组规则网格点比如10×10的整数坐标点生成三角网结果偶尔出现对角线方向不统一同一个点集生成两次得到不同的剖分而且每次都有少量狭长三角形。原因规则网格里四个点经常精确共圆此时外接圆检测的数值结果对浮点误差极其敏感。inCircle判断差一个ulp就可能从true变false导致空洞边界判断不一致最终剖分结果不稳定。解决在circumcircleContains里引入显式的容差项eps 1e-12 * (dx*dx dy*dy)判断“外接圆半径距离减去新点距圆心距离”与eps的关系不再用严格大于零的浮点比较。bool circumcircleContains(const Triangle tri, const Point p, double eps 1e-12) { // 用叉积和行列式计算外接圆然后直接判断p是否在圆内 // 注意这里返回true表示p在圆内或圆上用于Bowyer-Watson的空洞扩展 }这个eps不能拍脑袋取太大否则会把原本不共圆的点也包含进去导致空洞过度扩展生成多余三角形。经验值是取点集坐标尺度的1e-12如果你的点集坐标单位是毫米且有大量共圆情况调到1e-10也不会对结果质量产生明显影响。5.2 超三角形放大倍数不够边界点被“挤出”有效区域现象点集的凸包边缘处生成的三角形形状诡异有些三角形顶点组合明显不对甚至出现负面积三角形。原因超三角形的顶点坐标与真实点集坐标尺度相差不够大某些边缘点的外接圆超出了超三角形的包围范围导致这些点被错误地包含到空洞扩展中生成了跨越超三角形边界的新三角形。删除超三角形时这些跨边三角形因为包含超三角形顶点而被移除但它们的邻居——那些本应该保留的三角形也被误标记为deleted。解决把超三角形的放大倍数从100倍提高到200倍只是缓解根治方案是先归一化点集到单位正方形超三角形取归一化后的外包三角形删除超三角形后再把剩下的三角形顶点坐标反向映射回原始坐标。归一化也能解决切片数据坐标尺度过大的问题。5.3 邻居数组没有双向更新导致边翻转直接死循环现象程序在边翻转阶段陷入死循环同一个四边形反复翻转每次结果都一样CPU占用100%但就是结束不了。原因在buildTrianglesFromBoundary中只设置了新三角形的邻居索引忘记反向设置旧三角形的邻居索引导致邻接关系不对称。边翻转时遍历到的“邻居”信息是错的翻转后四边形被改坏检查“是否需要翻转”的判断永远为真死循环出现。解决建立一条原则——每次修改邻居关系无论创建、删除还是翻转都必须同时更新两边的邻居数组不存在单向设置。我后来在Debug构建里加了一个一致性校验函数遍历所有三角形检查triangles[t.neighbors[i]].neighbors[j] t是否成立不成立直接抛异常。这个校验在Release构建里关闭不影响性能但Debug阶段能立刻抓到所有索引混乱的问题。5.4 大量删除三角形时vector迭代器失效现象删除超三角形阶段程序崩溃报错信息指向vector::erase但代码看起来没有逻辑问题。原因在遍历vector的同时调用erase迭代器指向的元素被移动遍历越界。更隐蔽的是std::remove_if配合erase时自定义的谓词函数里访问了已经被移动的元素导致读野指针。解决固定采用“mark-delete-compact”三步走先统一标记再统一压缩。C20后我用std::erase_ifC17必须手动remove_iferase两段式。任何情况下不在遍历循环体里直接erase还是更稳的删除逻辑全部延迟到循环结束后统一处理代码可读性也会变好。5.5 空间网格尺寸对性能的影响格子个数是玄学还是科学现象点集规模到10万级插入变得异常慢定位三角形花费的时间占总耗时的60%以上。原因空间网格的格子尺寸太大每个格子覆盖大量三角形起不到“缩小搜索范围”的作用格子尺寸太小则一个点会穿过好几个格子查找时会重复检测大量空格子甚至因为跨格子边界而漏检。解决格子尺寸取点集平均最近邻距离的2~3倍是一个经验值或者在点集数量已知时格子边长取sqrt(包围盒面积 / 点数量) * 1.8。后者更适合预处理阶段的点数已知场景。实测10万点集单次插入的检测三角形数量从数千降到几十性能提升显著。6. 进阶实践从静态网格到动态更新的最后一个技巧如果你已经把上面的代码跑通恭喜你已经能处理静态点集的Delaunay三角网生成。但实际项目很少只做一次剖分——点云配准后要反复加点、滤波后要删除离群点、路径规划里障碍物在移动。这时候纯粹靠重新执行整个Bowyer-Watson是没有性能优势的。最后一个技巧如何把增量插入扩展到点删除和局部重建。点删除的常见做法是删除目标点后收集所有与该点相连的三角形把它们从三角网中移除形成一个“多边形空洞”。空洞的边界顶点构成一个多边形对这个多边形重新Delaunay三角化。这里的关键是——重新三角化的结果必须与原来的整体三角网能无缝缝合即新的三角形边界与原空洞边界完全重合。void deletePoint(std::vectorPoint points, std::vectorTriangle triangles, int pointIndex) { // 1. 收集与该点关联的所有三角形 std::vectorint affectedTriangles; for (size_t i 0; i triangles.size(); i) { if (triangles[i].deleted) continue; for (int j 0; j 3; j) if (triangles[i].v[j] pointIndex) { affectedTriangles.push_back(static_castint(i)); break; } } if (affectedTriangles.empty()) return; // 2. 收集空洞边界所有受影响三角形的外边界不与这些三角形共享的边 std::vectorint boundaryVertices; collectBoundary(boundaryVertices, affectedTriangles, triangles); // 3. 删除受影响三角形 for (int idx : affectedTriangles) triangles[idx].deleted true; // 4. 对空洞多边形做三角形扇或递归剖分并缝合到外网 reTriangulateHole(boundaryVertices, pointIndex, triangles); }点删除的坑在于空洞多边形可能是非凸的直接三角形扇会生成非法三角形。正确处理方法是先对空洞多边形做三角剖分剖分线必须完全在多边形内部否则会穿到外部与已有三角形重叠。实在不想提前实现这个逻辑可以在空洞较小时对所有边界点跑一遍局部遍历把多边形看作独立点集重新做一次Bowyer-Watson剖分然后把剖分结果合并回主三角网代价是稍微多一点冗余计算但实现简单不易错。做完增删改后强烈建议加一个验证函数对最终每个三角形计算外接圆遍历点集中所有点检查是否有点落在圆内O(n²)的复杂度不适合大点集但作为Debug阶段的校验彻底可靠。同时检查整个三角网的拓扑一致性每条内部边恰好被两个三角形共享所有三角形的顶点索引都在有效范围内。这套实现我在几个项目里反复用过。某次做地形网格实时更新点云以每秒5000点的速度流入我用上述增量插入配合格子索引整体延迟保持在几十毫秒内远优于每帧全量重建。回头想想最开始踩过的浮点坑和索引混乱坑其实都是对“算法的不变性”理解不够深。Delaunay三角网的精髓不在那几行代码而在于对空圆特性和拓扑结构的敬畏——任何时候改变网格结构先验证再上生产环境这个习惯帮我省了不少后悔药。希望这篇笔记对你有点用动手跑起来遇到问题欢迎再拿出来对一遍。本文还有配套的精品资源点击获取