C++增量Delaunay三角剖分原理与实现避坑指南

简介面向希望深入理解计算几何与C工程实现的开发者这份资源提供了一套基于“半边”数据结构的增量Delaunay三角剖分算法实现。算法采用逐点插入的增量策略以最大化所有三角形的最小角度为目标适合作为网格生成与几何处理的学习样例。核心代码集中在GaussCurvature_Sample目录下的main.cpp中配合MeshLib基础库和Visual Studio工程文件可直接编译运行便于单步调试和结果可视化。压缩包共89个文件、约4.17MB以h头文件与cpp源文件为主另有sln/vcxproj构建配置、生成的exe与pdb调试文件以及gif演示和说明文档等完整覆盖从源码阅读到运行验证的链路。读者可借助工程中的注释与演示程序快速掌握半边数据结构的遍历方式及其在网格处理中的应用。已有1002人学习下载适合计算几何入门者、图形学学习者及需要处理三角网格的开发者参考。1. 增量 Delaunay 三角剖分这份 C 代码包到底值不值得下做网格生成、路径规划或者点云重建的同行大概率都跟 Delaunay 三角剖分打过照面。市面上的库不少但要么是 CGAL 那种重依赖的庞然大物要么是只给二进制的黑匣子。我拿到这份用 C 实现的增量 Delaunay 三角剖分源码时第一反应是看两个东西顶点插入流程是不是标准 Bowyer-Watson 思路以及翻边之后有没有处理退化三角形。这两个点直接决定能不能移植到自己的项目里。这份源码的亮点在于不依赖第三方库核心剖分逻辑集中在单个实现文件里适合想彻底搞懂算法细节、或者需要把三角剖分嵌进自研引擎的从业者。它解决的是“从零实现”的问题而不是“调一个接口”的问题。2. 核心原理与数据结构空外接圆判据、超三角形与半边思想2.1 Delaunay 判据剖分质量的第一道门槛增量 Delaunay 三角剖分的理论基础只有一个任意一个三角形的外接圆内部不能包含点集中的其他顶点。这个判据决定了最终网格的“株连”程度——三角形尽可能接近等边最小角最大化。实现的时候最朴素的做法是逐点插入每来一个新点找到包含它的三角形连接新点与该三角形的三个顶点然后对生成的边做翻转检查。翻转检查的本质是判断两个相邻三角形共享的那条边是否需要换成另一条对角线。判断依据就是对边上的第四个点做空外接圆测试。如果第四个点落在共享边对面的三角形外接圆内部说明当前边不与 Delaunay 判据兼容要翻转否则保留。这个操作每插入一个点都要对周围若干三角形做局部更新增量这个名字就是这么来的。判定函数是核心中的核心它不涉及任何三角函数。直接利用行列式符号判断点与圆的位置关系代码写出来大概四五行却是整个算法最容易因为符号弄反而出错的地方。我见过不少实现把“内部”和“边界”搞反导致翻边后网格出现裂缝。2.2 为什么用 C内存控制与删除操作效率这份源码选择 C 而不是 Python 或 Java背后有很实际的考量。Delaunay 三角剖分在顶点数量大的时候性能瓶颈在顶点定位和翻边操作这两类操作都强烈依赖数组的随机访问效率。Python 的 list 虽然也能存三角形索引但每次翻边要新建对象、回收旧对象GC 压力会让性能急剧下降。C 的优势在于可以预先分配顶点数组三角形用固定大小结构体存储三个顶点索引和三个邻居索引。翻边操作本质上是修改索引不涉及内存重排。另外删除超三角形、清理无效三角形时直接用一个布尔标记数组加延迟清除避免频繁erase导致的内存搬移。这些细节在几千个点的时候体现不出来点集上万以后差距就非常明显了。2.3 数据结构设计索引三角形而非指针三角形源码里的核心数据结构是一个三角形结构体包含三个顶点索引和三个邻居索引。顶点坐标单独放在一个结构体数组里。这里有个设计选择值得注意不用指针而用整数索引。原因有两点。第一整数索引天然支持输出到文件——写 OBJ 格式或者自定义格式时直接把索引值写出去就行指针得先换算成下标。第二翻边操作需要交换邻居关系整数索引在 vector 扩容后不会失效但指针会。某开发者没想清楚这点一开始用指针存邻居结果顶点数组一旦reserve不够触发扩容所有邻居指针全部悬垂排查了整整一个下午。三角形数据布局的示意代码如下// 顶点坐标加索引标记 struct Vertex { double x, y; }; // 三角形三个顶点索引逆时针 三个邻居索引 struct Triangle { int v[3]; int n[3]; // n[i] 是与 v[i] 相对的邻居三角形索引,-1 表示边界 bool valid; // 延迟删除标记 };顶点数组按输入顺序存三角形数组动态追加。每插入一个新点至少产生一个新三角形而翻边操作会修改已有三角形的索引不会改变数组大小。需要删除的三角形比如包含超三角形顶点的统一标记为无效最后一次性清理。这里把邻居索引和顶点索引交错存放翻边的时候要同时维护两端三角形的顶点和邻居逻辑容易绕晕。一个实用的习惯是定义半边数组edge(va, vb)作为查询键记录该边属于哪个三角形就能快速定位共享边的两侧三角形。3. 增量插入实现顶点定位、翻边修复与超三角形清理3.1 顶点定位从蛮力扫描到局部寻路增量插入的第一步是找到新点落在哪个现有三角形内部。最简单的方式是遍历所有有效三角形逐个用重心坐标判断。这个写法在三五百个点的时候完全够用但点集到五千以上就会明显卡顿。因为每插入一个点都要扫全表复杂度退化到 O(n²)。源码里用的是有向寻路walking策略从上一个插入点所在的三角形出发沿着新点相对当前三角形的位置跨边进入相邻三角形。平均只需要走根号 n 步就能到达目标三角形。实现要点是维护“上一个操作三角形”作为缓存因为新点和旧点在空间上往往相邻从缓存出发比从头扫快得多。判断点是否在三角形内部可以用重心坐标也可以简化成同向法对三条边依次做叉积如果符号一致就说明点在三角形内部或者边界上。不过浮点情况下点落在边上的情形需要单独定义。源码在这地方给了一个容差参数叉积绝对值小于某个阈值就认为在边上直接复用边所在三角形的处理逻辑。3.2 翻边修复维护 Delaunay 性质的循环检查找到包含新点的三角形之后把新点与三角形的三个顶点连成三个新三角形。此时局部网格并不一定满足 Delaunay 判据因为新点可能落在邻居三角形的外接圆内部。修复手段是翻边如果两个相邻三角形共享的边不符合判据就把共享边换成另一条对角线。翻边操作会引入新的边新的边又可能破坏更外圈的 Delaunay 性质所以必须用一个栈来记录待检查的边。每次翻完边把受影响的新边压入栈循环处理直到栈空。这个循环终止性依赖一个事实每次翻转都会让最小角严格增大而最小角的取值是离散的所以最终会停下来。这里会碰到一个典型问题共享边两侧的三角形可能是无效的比如包含超三角形顶点的旧三角形。翻边时要把这些无效邻居跳过去否则访问到validfalse的三角形索引关系就是错的。3.3 关键代码插入一个顶点的完整流程把上面这些组织到一起单次插入的代码骨架长这样// 插入一个点,返回创建的新三角形数量 int insertVertex(double x, double y) { Vertex pt {x, y}; int idx vertices.size(); vertices.push_back(pt); // 1. 定位:从 walkCache 出发找包含 pt 的三角形 int tri locateTriangle(pt); // 2. 把 pt 与 tri 的三个顶点相连,分裂成三个新三角形 int t0 splitTriangle(tri, idx); // 3. 三个新边压入修复栈 stackint edgeStack; pushEdgesOfTri(t0, edgeStack); // ... 另外两个新三角形同理 // 4. 循环翻边 while (!edgeStack.empty()) { auto [va, vb] edgeStack.top(); edgeStack.pop(); if (shouldFlip(va, vb)) { flipEdge(va, vb); // 把新生成的边压栈 pushAdjacentEdges(va, vb, edgeStack); } } walkCache tri; // 缓存定位起点 return validTriangles; }splitTriangle负责把旧三角形的valid标记改成false创建三个新三角形并修正邻居关系。这里的坑点在于邻居索引修正原三角形的三个邻居要分别挂到三个新三角形对应的边上次序必须和顶点排列一致一错就是整片网格的连锁错误。shouldFlip就是空外接圆测试返回值是布尔量共享边若不符合 Delaunay 判据则返回真。这个测试的输入是四个点——共享边的两个端点和两个对顶点。行列式符号判断内部点的函数要特别小心坐标差值的计算次序。我一般会把四个点的坐标按共享边端点对齐避免多次计算时符号不自洽。flipEdge是另外一个复杂度来源。两个三角形共享一条边时它们构成一个四边形——如果这个四边形是凸的才能翻边如果是凹的翻边会产生退化网格。源码里有个在shouldFlip之前做的凸性检查遗漏它就会出现三角形重叠。4. 边界处理与退化情形超三角形、共圆容差与输出组织4.1 超三角形选择范围计算的“玄学”与可靠做法增量 Delaunay 实现必须从某个初始结构开始。代码包的做法是先构造一个覆盖所有输入点的大三角形——超三角形。超三角形的三个顶点不是真实点插入完所有点之后要连带着所有包含这些顶点的三角形一起删除。超三角形的尺寸选取有讲究。取得太小部分真实点落在超三角形外面破坏“所有点都位于某个三角形内部”的假设取得太大浮点计算时矩离很大的点做叉积容易丢失精度。源码里的做法是先扫一遍点集求包围盒然后取包围盒对角线的若干倍数作为超三角形的边长。我习惯再用一个安全系数通常是 10 倍以上尽量避免边界点恰好落在超三角形边上的情况。边界点落入超三角形边所在的直线附近时会引发三角形面积趋近于零的退化问题哪怕概率低也影响确定性。4.2 共圆四点与退化三角形容差、抖动与确定性真正常见的翻车场景是四个点共圆。矩形网格的四个顶点就是典型——外接圆测试的行列式算出来刚好接近零此时翻转与不翻转在数学上都合法但不同的判断会导致不同的网格结构。源码的做法是给行列式计算加一个很小的 epsilon判断符号时分三种情况正、负、近似零。近似零时不翻转保持现有边。这个方法简单有效但要注意 epsilon 的值不能设太大否则会把本应翻转的边也过滤掉。经验值是1e-12到1e-10这个量级具体要结合输入坐标的尺度来调。另外一个更稳妥的做法是在做测试前对整数型输入做小扰动。比如坐标都是整数的输入点可以给每个点加一个基于索引的微小平移。扰动需要保证不改变点的相对顺序否则会引入真正的几何错误。源码没有内置扰动逻辑但保留了扩展接口——想用抖动方案时只需要在顶点插入前做一个坐标微调即可。4.3 输出组织从内部结构到可用网格三角形剖分完成之后内部结构是带邻居索引的。用户通常只需要输出顶点坐标列表和三角形顶点索引列表。源码里提供了简单的序列化函数支持输出成带索引的文本格式。输出时要额外做的一件事是过滤无效三角形。延迟删除策略留下的数组里全是validfalse的空壳输出前压缩数组会把三角形索引重排一遍。这里容易踩坑邻居索引指向的是旧三角形数组的下标压缩后必须同步映射否则跳过这一步直接写文件写出来的三角形网格会有大量“飞边”和重复面。这一步的逻辑如下// 压缩三角形数组,重建索引映射 std::vectorint oldToNew(trigs.size(), -1); std::vectorTriangle compact; for (int i 0; i trigs.size(); i) { if (!trigs[i].valid) continue; oldToNew[i] compact.size(); compact.push_back(trigs[i]); } // 遍历 compact,把 n[0..2] 替换成 oldToNew 映射后的值如果是写 OBJ 格式顶点从 1 开始编号输出三角形顶点索引时全部加 1。这个细节看着不起眼拿 MeshLab 打开的时候少一行的错位直接导致整个网格变成一团乱麻。5. 避坑与常见问题排查增量 Delaunay 的五个高频翻车现场5.1 翻边死循环现象插入某个点之后程序进入无限循环CPU 占用拉满加日志后发现shouldFlip一直返回真。原因flipEdge修改了三角形顶点索引之后边缘栈里压入的新边索引没有正确剔除旧边导致同一条边反复被翻转。更隐蔽的原因是凸性检查缺失四边形是凹的情况下翻边会让三角形面积变成负值进而影响后续判据。解决在shouldFlip里先验四边形凸性——两条对角线的叉积符号必须一致。另外在flipEdge成功后把共享边附近的新边逐个压栈并且压栈前检查边上两个三角形是否都有效。加上这两个防护当前点集的死循环基本消失。5.2 删除超三角形后残留“飞边”现象最终输出的网格里有一圈细长的三角形顶点坐标极大或者极小明显不属于输入点集。原因超三角形的三个顶点参与了很多三角形的构造但删掉“包含超三角形顶点”的三角形时只判断了三角形的顶点索引是否等于超三角形顶点索引。没考虑到某些三角形虽然不直接包含超三角形顶点但与这些三角形共享边。其邻居索引指向的三角形索引已经失效却没被清理。解决删除时做一轮传播清理。先标记所有包含超三角形顶点的三角形再遍历所有有效三角形如果三个邻居中有被标记的那自己也标记为无效。循环几次直到没有新标记产生。这个方案不会误删正常三角形但胜在简单可靠。5.3 定位函数越界现象程序在插入第几个点之后崩溃断点停在locateTriangle里访问trigs[tri]tri的值是 -1。原因walking定位算法在跨边时走到了边界外——也就是说新点超出了当前剖分的凸包范围。正常情况不会发生因为初始超三角形覆盖了所有点但如果超三角形的范围取得不够大或者插入点存在大量重复坐标导致某个三角形面积趋近于零寻路会从面积很小的缝隙跨出去。解决给locateTriangle加一个护栏跨边时判断目标三角形是否为 -1如果是就退回到蛮力遍历。这个 fallback 路径虽然慢但只会触发极少数点。之后再把超三角形范围从对角线 10 倍提到 100 倍基本不再触发。5.4 浮点误差导致翻边不彻底现象输出的网格里存在大量瘦长三角形局部看起来没完全剖分直观感受是“有些地方该翻没翻”。原因空外接圆测试用的行列式计算在坐标尺度接近 1e8 时精度丢失严重。比如点坐标是经纬度180 的量级两个巨大数值相减后的小数位根本保不住。此时epsilon按 1e-12 设置完全失效。解决把坐标归一化到单位包围盒再做三角剖分输出时再映射回去。归一化操作在内存里不占地方但对浮点精度帮助极大。源码的接口没直接提供归一化但坐标是double在输入层之前做一个平移缩放的预处理即可简单有效。5.5 大规模点集的性能瓶颈现象点集从几千涨到五万之后插入速度显著下降实测耗时从秒级变成分钟级。原因walking定位的平均复杂度是 O(√n)但翻边修复的复杂度跟局部三角形的数量相关。五万点规模时局部凸壳周围的三角形数量增多边缘栈的处理次数也跟着涨。更关键的是每次插入都从缓存三角形出发如果输入点的空间分布是随机均匀的缓存命中率表现不错但如果输入点是按扫描线顺序排列的缓存与新点距离常常很远定位复杂度退化。解决在插入前做一次空间排序。可以先按 y 坐标分桶再按 x 排保证相邻插入点在空间上靠近。这个预处理不改变点集内容只是打乱插入顺序能显著提升缓存命中率。实测同一份点集排序前后耗时差距接近 40%。6. 顶点数量自检、欧拉公式验证与性能剖析拿到源码跑通之后需要做一个系统的验证来确认实现正确。第一件必做的事是欧拉公式自检。对于没有洞的三角网格V - E F 1这里是平面剖分包含外边界。V 是输入点数F 是三角形数E 是边数可以由三角形数算出来每条边被两个三角形共享边界边只属于一个三角形所以 E (3F B) / 2B 是凸包边界边数。代入 V - E F 1 化简后能得到一个跟 B 相关的关系式。自己写个小脚本随机生成几千组均匀分布的点验证这个恒等式是否对每个输出都成立比肉眼检查靠谱得多。第二件事是验证 Delaunay 性质本身。写一个独立于剖分算法的验证器对输出的每个三角形计算其外接圆圆心和半径再验证所有其他顶点都不在圆内。这个验证绝对不能复用实现里的shouldFlip函数否则等于自己检查自己。独立验证器用最直接的距离比较法不追求性能只追求正确。做到这一步之后再随机抽几十组小点集把输出网格画出来看有没有裂缝、重叠或者明显的瘦长条。性能剖析方面把顶点定位、翻边操作和清理阶段分开计时。常见做法是在insertVertex的三个阶段分别加chrono计时累计到三个独立的计数器里。点集从一万递增到十万记录三阶段占比。我发现定位时间一般占 30%翻边占 60%清理占 10%——翻边比例过高就说明边缘栈里压入了太多无效边应该回溯凸性检查逻辑。我自己的习惯是任何一次对共享边的修改都要重新验证一遍这条边两侧的三角形是否仍然满足 Delaunay 判据。从那以后每次跑完剖分我都会强制走一遍“欧拉公式 独立判据验证”的流程。算法没有玄学每一步都能被数学验证代码就不会翻车。这份 C 增量 Delaunay 实现把核心流程拆得比较干净顺着源码走一遍能彻底搞清楚 Bowyer-Watson 思路的每个细节。希望帮到你。本文还有配套的精品资源点击获取