简介：基于C++实现的增量Delaunay三角剖分算法工程包，面向计算几何学习者与需要处理三角网格的C++开发者。算法逐点插入并局部翻边，使剖分中所有三角形的最小角尽量最大化；项目采用石溪大学顾贤峰教授提出的“半边”数据结构组织网格，适合想深入理解几何拓扑实现的人群。压缩包共89个文件、4.17MB，包含13个h头文件、1个cpp主实现（位于GaussCurvature_Sample/GaussCurvature/main.cpp）、Visual Studio解决方案与工程文件、已编译的exe及构建日志，可直接打开工程运行或按需拆解调用。已有1003人浏览/学习。读者可拿到完整可编译的Delaunay三角剖分示例，厘清增量插入、局部边翻转与半边结构维护的配合关系，还可参考其网格类封装方式，迁移到其他几何算法项目中。

1. Delaunay三角剖分：为什么“最不尖锐”的网格要靠局部操作来维护

给一堆散点，Delaunay三角剖分能生成一张“尽可能匀称”的三角形网格：任何三角形的外接圆内都不允许出现第四个点。这个空圆准则听起来是全局约束——每一圈网格都要满足它。但真正让增量算法工程可用的，是一条反直觉的性质：往一张几十万个三角形的网格里插入一个新点，需要修改的只是和它相邻的一小片区域。把其中不合法的那条边翻一下，整个网格的空圆条件便自动恢复。于是插入一个点的过程可以被拆成定位、拆分、翻边三步，每一步都只围绕新点的局部邻域做文章。搞懂这三步背后的几何机制和C++实现取舍，才算真正明白这个算法为什么能work，又为什么在某些输入上会翻车。

2. 增量插入机制：一点扰动为什么不会波及整张网格

2.1 局部 Delaunay 边与空圆条件的等价关系

增量算法的正确性基石是“局部 Delaunay 边”这个概念。取一条共享边 AB，两侧分别是三角形 ABC 和 ABD。如果 D 落在三角形 ABC 的外接圆内部，边 AB 就是非法的；反之，如果四边形 ACBD 是凸的，且 D 在 ABC 外接圆外部，AB 就是合法的。这个判定与全局空圆条件之间有一个重要等价：当网格中所有内部边都满足局部 Delaunay 条件时，整张网格就是全局 Delaunay 的。

这意味着修复可以完全局部化。插入一个点之后，受影响的非法边一定与包含新点的三角形相邻，不存在隔空被破坏的边。因为空圆条件和局部判定等价，所以只要逐条处理非法边，就能恢复全局性质，不需要扫描全图。增量算法里那句口头禅可以这么记：谁碰了新点，就查谁的外接圆。

这个结论也解释了算法伪代码为什么那么短。初始化超级三角形，然后逐个插入点，每插一个点只做局部修复，最后删掉超级三角形——主循环就这么几行。真正的工作量全部藏在“如何快速定位”“如何判断一条边是否合法”和“如何翻边而不弄错邻居关系”里。

2.2 步行法定位：点与三角形的包含关系

插入新点后的第一件事，是找到它落在哪个现有三角形内部。最朴素的做法是遍历全部三角形做包含测试，每次 O(n)，在几十万点的规模下完全不可接受。工程实现里更常用的是步行法（walking algorithm）：从某个已知三角形出发，逐步向目标点移动。

每一步检查目标点 p 相对当前三角形三条边的位置。三角形按逆时针存储顶点，如果 p 在某条有向边的右侧，就说明应该从这条边跨到对面的邻居三角形去。直到 p 位于三条边的同一侧，也就是落在当前三角形内部，定位结束。

int locatePoint(const vector<Point>& pts, const vector<Tri>& tris, const Point& p, int start) { int t = start; for (int guard = 0; guard < (int)tris.size() * 4; ++guard) { bool inside = true; for (int i = 0; i < 3; ++i) { int a = tris[t].v[i]; int b = tris[t].v[(i + 1) % 3]; double o = orient2d(p, pts[a], pts[b]); if (o < 0) { // p 在边 (a,b) 的右侧 int next = tris[t].n[i]; if (next < 0) return t; // 理论上不会走到边界 t = next; inside = false; break; } } if (inside) return t; } return t; // 保护性退出，正常不会执行到这里 }

这里的orient2d(p, a, b)是计算几何里的方向测试：返回正值表示 p 在从 a 到 b 的有向边左侧，负值表示右侧。三角形顶点按逆时针存储时，如果 p 在三条边的左侧，就说明它在三角形内部。每跨一步，离目标三角形都更近一点；只要网格拓扑正确，这个过程必然收敛。

代码里的 guard 循环上限是必须的。步行法理论上最坏走 O(n) 步，但浮点误差可能让判断在两条边的临界处来回振荡，不加保护会死循环。生产实现里还有一种常见优化：把 for+break 展开成三次顺序判断，省掉分支跳转。在几十万点的规模下，这个细节能带来可感知的定位速度差异。

2.3 三角形拆分与级联翻边

定位完成后，新点 p 落在三角形 ABC 内。先把 ABC 从网格中摘除，再用 p 与 A、B、C 生成三个新三角形 ABp、BCp、CAp。几何上这一步等于把旧三角形劈成三块，开销很小。但劈完之后，新三角形可能与它的邻居构成共边四边形：比如三角形 ABp 与邻居 ABD 共享边 AB，如果 D 落在三角形 ABp 的外接圆内，AB 就不再是合法边。

这时需要翻边（flip）：把共享边 AB 换成 Dp，两个三角形变成 ADp 和 DBp。翻边之后，新的 ADp 边或 DBp 边可能又让更外圈的三角形不合法，于是翻边必须递归进行。这个递归过程像涟漪一样向外扩散，但因为非法边只能沿着含新点的区域传播，扩散范围被限制在邻域内。

翻边必然会终止，原因在于它等价于“选择正确的对角线”。对一个凸四边形，两条对角线中恰好有一条满足局部 Delaunay 条件。每翻一次边，四边形内的最小角严格变大，而网格三角形总数不变，所以递归深度存在上界。随机顺序插入时，这个深度的期望值接近常数，这也正是增量算法平均效率的关键所在。

2.4 一个完整插入序列的动线

把定位、拆分、翻边串起来，一次插入完整流程可以概括为：

int t = locatePoint(pts, tris, newPt, lastTriangle); int a = tris[t].v[0], b = tris[t].v[1], c = tris[t].v[2]; // 创建三个新三角形，手工维护邻居关系 // 删除旧三角形 t // 对三个新三角形分别执行 legalize legalizeTriangle(ta, 1); legalizeTriangle(tb, 1); legalizeTriangle(tc, 1);

每个legalizeTriangle处理一块共享边：如果是内部边且不满足局部 Delaunay 条件，就翻边，然后对新形成的两条边继续递归。每次插入的“动线”都完全一致：先定位到包含新点的三角形，把它拆开，再把新点带来的非法边一条条修好。整个增量算法的骨架到这里就清楚了，剩下的难点全在数据结构——邻居关系怎么存、翻边时如何精准地改写四个三角形的连接信息。这部分放到下一章展开。

3. C++ 实现：数据结构、几何谓词与翻边的源码级拆解

3.1 三角形的数据表示：顶点数组加邻居数组

增量算法里最常用的三角形结构是顶点数组加邻居数组：

struct Tri { int v[3]; // 顶点索引，按逆时针排列 int n[3]; // 邻居索引，n[i] 是边 (v[i], v[(i+1)%3]) 对面的三角形 };

例如n[0]表示与边(v[0], v[1])共享的邻居三角形。如果某条边没有邻居，说明它是当前网格的边界边，用 -1 表示。

为什么用 int 索引而不是指针？第一个原因是 vector 扩容时元素地址会移动，裸指针全部失效，索引不会。第二个原因是内存局部性：一个 Tri 占 24 字节，存放在连续数组里，几十万个三角形遍历时缓存命中率远高于用 new 逐个分配的散落对象。第三个原因是翻边时要频繁改写邻居指向，索引赋值就是一次 int 拷贝，比指针更简单也更安全。

这个结构有一个隐藏约定：v 数组的排列顺序必须始终是逆时针。所有几何谓词和邻居关系都依赖这个约定。一旦某个三角形变成顺时针，后面的翻边条件会整体取反，产生难以排查的镜像网格。

3.2 orient2d 与 inCircle：两个决定正确性的几何谓词

方向测试orient2d是全部几何判断的地基。它计算的是向量 AB 和 AC 的叉积：

double orient2d(const Point& a, const Point& b, const Point& c) { return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x); }

返回值大于零表示 c 在 b 的左侧，即三角形 (a,b,c) 为逆时针。定位、翻边前的凸四边形检查、超三角形顶点选择，都依赖这个函数的符号。

外接圆判断inCircle是另一个核心谓词。它判断第四个点 d 落在三角形 abc 的外接圆内还是圆外。直接计算圆心半径再比较距离当然可行，但更高效的做法是构造一个四阶行列式：

double inCircle(const Point& a, const Point& b, const Point& c, const Point& d) { double adx = a.x - d.x, ady = a.y - d.y; double bdx = b.x - d.x, bdy = b.y - d.y; double cdx = c.x - d.x, cdy = c.y - d.y; double abd = adx * bdy - ady * bdx; double bcd = bdx * cdy - bdy * cdx; double cad = cdx * ady - cdy * adx; double alift = adx * adx + ady * ady; double blift = bdx * bdx + bdy * bdy; double clift = cdx * cdx + cdy * cdy; return alift * bcd + blift * cad + clift * abd; }

先做坐标平移，把 d 作为原点，能有效减小坐标值的绝对大小，从而降低浮点误差。符号约定是：a、b、c 逆时针排列时，返回正表示 d 在圆外，返回负表示在圆内。这两个谓词的实现都不复杂，但它们共同决定了算法会不会在某个角落静默地生成非法三角形。

3.3 拆分三角形时的邻居接线

插入新点 p 到三角形 ABC 内部时，旧三角形被替换成三个新三角形。最清晰的实现方式是先算出旧三角形三个顶点的索引，再按固定模式创建新三角形：

void insertPoint(int pi) { int t0 = locatePoint(pts, tris, pts[pi], lastTriangle); const Tri& old = tris[t0]; int a = old.v[0], b = old.v[1], c = old.v[2]; // 新三角形 tA 的顶点顺序为 (a, b, pi) // 新三角形 tB 的顶点顺序为 (b, c, pi) // 新三角形 tC 的顶点顺序为 (c, a, pi) // 每个新三角形的邻居需要同时连接旧邻居和另外两个新三角形 }

创建时最需要注意的是内部接线：三个新三角形之间互连，形成一个小环。例如 tA 的边 (b, pi) 对面是 tB，边 (pi, a) 对面是 tC。手工写这三个三角形的邻居数组极易出错，一个常见的做法是先统一置 -1，再用一个辅助函数linkTri把两个三角形在指定边上互相指向对方：

void linkTri(int t1, int e1, int t2, int e2) { tris[t1].n[e1] = t2; tris[t2].n[e2] = t1; }

然后把拆分时的所有邻居关系写成一行行的linkTri调用。代码虽然长，但每行只做一件事，出错了也容易定位。这类“接线错误”不会立刻报错，而是会在后续翻边或输出网格时以拓扑断裂的形式出现，排查起来非常耗时。

3.4 legalize 与 flip 的源码实现

legalize负责检查和修复一条共享边。它的输入是一个三角形索引和一条边的索引：

int thirdVertex(const Tri& tri, int u, int v) { for (int i = 0; i < 3; ++i) { if (tri.v[i] != u && tri.v[i] != v) return tri.v[i]; } return -1; } void legalize(int t, int e) { int nb = tris[t].n[e]; if (nb < 0) return; // 边界边，直接返回 int u = tris[t].v[e]; int v = tris[t].v[(e + 1) % 3]; int w = tris[t].v[(e + 2) % 3]; int z = thirdVertex(tris[nb], u, v); if (inCircle(pts[u], pts[v], pts[w], pts[z]) < 0) { flip(t, e, nb, u, v, w, z); // 需要翻边 } }

flip把共享边为 uv 的两个三角形（第三顶点分别是 w 和 z）替换成共享边为 wz 的两个新三角形。实现时要特别注意新三角形的顶点必须重新排成逆时针，并且四个邻居关系要整体重接。另一种更稳妥的写法是“先断后连”：先把四个相关三角形的邻居索引清零，再按新拓扑重新linkTri。少做一步，残留的旧邻居索引就会让下一次legalize读到错误的三角形。

翻边之后，新形成的两条边可能又和它们的邻居发生冲突，所以flip内部还需要对两条新边递归调用legalize。这与 2.3 节说的级联翻边是对应的。递归深度由最小角增大性质保证有界，不会栈溢出。

3.5 超级三角形与收尾清理

增量算法需要一块初始地盘作为定位起点，这就是超级三角形。它要覆盖所有输入点，最常用的构造方式是先算输入点集的包围盒，把中心作为三角形重心，再向外扩到对角线长度的两倍左右。

Point minP = pts[0], maxP = pts[0]; for (auto& p : pts) { minP.x = std::min(minP.x, p.x); minP.y = std::min(minP.y, p.y); maxP.x = std::max(maxP.x, p.x); maxP.y = std::max(maxP.y, p.y); } double d = std::sqrt((maxP.x - minP.x) * (maxP.x - minP.x) + (maxP.y - minP.y) * (maxP.y - minP.y)); Point c = { (minP.x + maxP.x) / 2.0, (minP.y + maxP.y) / 2.0 }; // 三个顶点围绕 c，半径取 2d 左右

超级三角形的顶点索引一般从 n 开始保存，n 是原始输入点数。当所有点都插入完毕后，任何包含超级三角形顶点的三角形都必须删除，剩下的就是最终 Delaunay 网格。

int n = pts.size(); for (auto it = tris.begin(); it != tris.end(); ) { if (it->v[0] >= n || it->v[1] >= n || it->v[2] >= n) { it = tris.erase(it); } else { ++it; } }

这一步看似简单，但它隐藏了一个细节：删除超级三角形后，原来与超级三角形相接的边界边会暴露出来，成为凸包边界。如果插入顺序或翻边过程有过任何轻微错误，这里往往会出现外挂的碎片三角形，一眼就能看出网格不对。

4. 避坑指南：增量 Delaunay 最容易翻车的五个细节

4.1 inCircle 符号约定与顶点绕序不一致

现象：网格生成过程几乎正常，但输出结果像是被镜像了，部分三角形重叠，或者某些内部边缺失。

原因：inCircle返回值的符号依赖三角形顶点的绕序。三角形按逆时针存储时，返回负表示第四点在圆内；如果某个三角形的顶点被存成顺时针，符号恰好互换，翻边条件整体取反，等于一直在维护“局部非 Delaunay”网格。

解决：在每次创建三角形后加一个 debug 断言，检查orient2d(pts[v0], pts[v1], pts[v2]) > 0。只要有一个三角形违反逆时针约定，立刻终止程序。这种断言只在排查时开，release 里可以关闭，但它能拦截掉一大类隐蔽的接线错误。

4.2 共圆点：inCircle 恰好等于零

现象：某个局部区域在两种对角线之间反复横跳，legalize递归不停止，或者同一份输入在不同优化等级下得到不同结果。

原因：四个点精确共圆时，inCircle 行列式值为零，翻边逻辑没有定义。真实数据中完全共圆很少见，但对称的栅格点、人工生成的测试点很容易触发。零值本身不可怕，可怕的是它在不同编译器或优化设置下的浮点表现不一致。

解决：最常见的是引入符号扰动（symbolic perturbation），给每个点一个与坐标相关的微小虚拟偏移，打破共圆对称性。工程上更省事的做法是对输入做轻量去对称化，并在inCircle返回绝对值小于某个极小阈值时默认不翻边。对正确性要求高的场景，应该换用精确自适应精度谓词，而不是依赖 double 的运气。

4.3 重复点与共线点没预处理

现象：剖分过程中出现 NaN，定位死循环，甚至数组越界访问。

原因：重复点会让两个顶点在几何上完全重合，inCircle 距离项退化，行列式信息丢失；共线点则让“点在三角形内部”的判定退化为“点在线上”，orient 值反复为零，步行法在相邻三角形之间不停摆动无法收敛。

解决：插入前先按坐标排序，再去重。排序可以用std::sort配合自定义比较器，去重后用唯一序列重新建点集。共线点属于另一种退化，不能硬塞进同一个剖分流程，需要根据业务规则做降维处理或单独分段。数据清洗是这一步的硬要求，侥幸心理基本都会以神秘崩溃收场。

4.4 超级三角形选太大

现象：靠近超级三角形的三角形面积异常大，整体质量统计被拉低；某些翻边操作产生极细长三角形。

原因：超级三角形顶点坐标比输入点大几个数量级时，大三角形上的外接圆半径也成倍放大。inCircle 里坐标差的平方项变成“大数减大数”，浮点有效位数被大量消耗，翻边判断变得不可靠。

解决：超三角形的半径取包围盒对角线长度的两倍左右，而不是拍脑袋取一个 1000 或 10000。另一个有效的补充手段是在插入前把所有点平移到以包围盒中心为原点，降低绝对坐标数值。几何变换不改变拓扑结构，平移对整个剖分结果零影响，却能让浮点误差明显变小。

4.5 步行法在边界点上的振荡

现象：定位函数触发 guard 退出，返回一个错误的三角形索引，后续翻边全部跑偏。

原因：目标点恰好在两个三角形的公共边上时，orient 结果在两个微小正负之间摇摆，两次判断分别选中不同的跨越方向，于是步行法在相邻三角形之间来回跳动，无法收敛。

解决：对 orient 的结果做钳制，当绝对值小于极小阈值时直接置为零，并把“点在边上”按“点在三角形内”处理。阈值不能定得太大，否则会把本应通过跨越离开的邻居点误判成内部点。一个实用的经验值是按输入坐标量级的 1e-12 到 1e-10 来做钳制。

5. 设计取舍：增量翻边、Bowyer-Watson 与分治法的路线选择

增量翻边并不是唯一的 Delaunay 构建路线。另一种思路 Bowyer-Watson 算法会先找出所有外接圆包含新点的三角形，把它们批量删除，再把留下的空洞与新点连成新三角形。它比增量翻边更直观，但每次插入都需要维护一个候选集合，删除和重建的数量更多，常数更大。分治法能在最坏情况下稳定做到 O(n log n)，但合并两个子剖分时要处理复杂的凸包边界缝合，实现量远大于增量法，而且完全无法支持动态加点。

增量法最大的工程优势在于局部性和动态性。局部性体现在每次插入只需处理与新点相邻的少量三角形；动态性体现在它天然支持流式添加点，后续任意时刻都能继续插入。代价是复杂度对输入顺序敏感：随机顺序插入时，步行法定位和翻边次数的期望都很友好，整体期望 O(n log n)；如果输入点按某种特定顺序排列，定位可能要穿越大半个网格，退回 O(n^2)。所以很多实现会强制对点集做随机打乱，牺牲输出顺序的确定性来换取稳定的性能。

有一次我处理某公司提供的逐行扫描地形数据，数十万个点按栅格顺序排列，直接喂进增量实现，生成时间从预期的几秒膨胀到两分多钟。排查了很久才意识到，不是算法写错了，而是输入顺序太规整，规整到每次插入的定位都要从上一个点所在的三角形跨过大半个网格。从那以后，我每次跑增量 Delaunay 之前都会先检查输入点的排列习惯，必要时加一次 shuffle；同时把 inCircle 和 orient2d 的符号约定用断言锁死，杜绝绕序问题在半夜偷袭。三角剖分本身不难，难的是让每一次局部修复始终聪明地工作。希望这篇拆解能帮到你，把这套增量机制的底层逻辑理顺。

本文还有配套的精品资源，点击获取