简介：这是一份面向C++开发者的Delaunay三角剖分算法学习资源，项目采用增量式构建思路，使用石溪大学顾贤峰教授开发的半边数据结构维护网格拓扑，目标是最大化剖分所得三角形的最小角，适合计算机图形学、几何建模与网格生成方向的读者参考。压缩包共89个文件，包含头文件与C++源码、Visual Studio工程文件（sln/vcxproj）、可执行程序与调试符号（exe/pdb），另有构建日志、资源定义文件与说明文档等，整体约4.17MB；目录按GaussCurvature_Sample、MeshLib等模块组织，核心实现集中在主程序源码中，便于对照阅读。资源已有1000人学习，适合希望快速理解Delaunay增量算法代码组织与半边结构用法的进阶读者。通过案例工程与源码阅读，可掌握增量点插入、局部合法化及网格拓扑更新等关键思路，为自定义网格处理工具打下基础。代码分层清晰，附带的GaussCurvature示例工程还能帮助读者在真实网格数据上观察剖分结果，加深对算法效果的理解。
1. 为什么我用 C++ 重写增量 Delaunay 三角剖分，而不是调库
如果你处理过点云重建、地形网格生成或者有限元前处理，大概率遇到过同一个尴尬：调完现成的三角剖分库，结果一换数据就出破洞、出现细长三角形、甚至直接崩掉。我最初用 Python 的 scipy.spatial.Delaunay 处理十万个点还算顺利，但一旦进入需要动态加点的场景——比如实时扫描点云逐步接入、或者交互式编辑网格——重建整个三角网的成本，会从「可接受」变成「肉眼可见的卡顿」。这个项目标题里的两个关键词：C++ 和增量，恰恰是冲着这个痛点来的。
所谓增量 Delaunay 三角剖分，意思是点不是一次性全部给齐，而是一个一个（或一小批一小批）插入到已有的三角网中，每插入一个点，只修改局部受影响的三角形，而不是全量重算。配合 C++ 的手动内存管理和直接控制底层数据结构，能把单点插入的均摊成本压到很低，这也是我选择自己实现而不是继续调库的原因。本文会把我的完整实现思路、核心数据结构、局部优化（Lawson 翻转）的判定条件、以及我踩过的几个让人头大的坑，全部摊开来讲。适合已经学过计算几何基础、想自己动手写一个能跑在真实数据上的三角剖分程序的人，也适合想搞清楚 Delaunay 三角剖分内部到底怎么转的人。我不打算给你贴一个几百行的成品代码然后说「拿去用」，而是按我实际调试的顺序，把每一步为什么这么做、参数怎么设、报错时看哪里讲清楚。
如果你只是想在某个工具链里快速出个网格，那我得先泼盆冷水：手写增量 Delaunay 工作量不小，边界情况和数值稳定性问题能磨掉你一周时间。但如果你需要动态加点、需要把剖分结果和自家渲染或仿真流程深度整合、或者干脆是想把这块硬骨头啃下来作为算法功底，那这条路值得走。下面开始，从我在内存布局和点定位上吃的亏讲起，因为这决定了整个实现的性能和稳定性。
2. 数据结构与点定位策略：先让三角网“转”起来，再谈效率
2.1 为什么带头部的三角形链表比直接存顶点索引好用
一开始我照着教科书上的伪代码实现，用一个std::vector<Triangle>存所有三角形，每个 Triangle 里存三个顶点索引。插入新点时，遍历所有三角形找包含新点的那个。这样在几百个点上没问题，但一旦到一万个点，每次插入都做全量遍历，总复杂度逼近 O(n²)，和全量重算几乎没有区别。真正的增量算法需要一个机制，让「找点落在哪个三角形」这个操作变快，同时让「删除旧三角形、创建新三角形」不引发大量的内存拷贝或索引失效。
我先说存储结构。我最终用的是带头部的双向链表，每个节点是一个三角形。之所以不用 vector，是因为增量插入会频繁删除和新增三角形，vector 的中间删除是 O(n) 的，而且迭代器失效规则会让记录「当前三角形」的指针变成野指针。链表正好相反，删除一个节点只动前后两个指针，新增节点就是改两个指针，代价是访问某个三角形需要遍历过去，但这恰恰是增量算法不需要的——我们从来不需要随机访问第 k 个三角形，我们只需要从当前三角形走到邻居，或者从某个点出发遍历它的周围一圈。
每个三角形节点里存的不只是三个顶点索引，还有三个邻居指针。邻居的序号用 0、1、2 表示，约定邻居neigh[i]是与顶点i相对的边对面的三角形。举个例子，三角形(a, b, c)，其中顶点 0 是 a，1 是 b，2 是 c，那么neigh[0]就是与边bc共边的三角形，neigh[1]与边ca共边，neigh[2]与边ab共边。这个约定必须全程序一致，否则后面做局部翻转时找共享边会乱套。我给每个三角形起一个顶点索引成员v[3]和一个邻居指针成员neigh[3]，另外留一个isDeleted标志位。删除三角形时先标记，等批量清理时再释放内存，避免频繁 new/delete 导致内存碎片。
存储方案定下来后，点定位就成了下一个关键问题。我第一版用的是遍历所有三角形找包含点，后来改成从「上次插入点所在的三角形」开始，沿三角形邻居链向目标点「走」过去。这个策略叫 walking 或 jumping-and-walking，在均匀分布的点上表现很好，均摊接近 O(n^1.5) 到 O(n log n) 之间，具体取决于点分布和起始位置。每个三角形有一个重心坐标判定函数，判断目标点是否在当前三角形内部或者边上，如果在边上还要处理特殊逻辑。
2.2 点定位的三种做法，以及我为什么选了 walking
教科书上常见的点定位方案有三种：第一种是全量遍历，简单但是 O(n) 每次；第二种是构建一个网格索引或者 bounding volume hierarchy，把空间分桶，每个桶里存三角形列表，查询时先定位桶再找三角形，这种方法对静态数据很好，但增量插入时维护索引本身也是一笔开销；第三种就是 walking，从某个已知三角形开始，根据目标点相对于当前三角形三条边的有向位置，决定下一步走到哪个邻居。
我最终选了 walking 并做了一点改进。核心代码大致是这个样子，这是我在项目中实际调通的第一版：
// 返回值是包含 p 的三角形指针 Tri* locatePoint(Point* p, Tri* start) { Tri* cur = start; // 防止极端情况下死循环，设置最大步数 int maxSteps = 100000; while (maxSteps-- > 0) { // 检查 p 是否在 cur 内部 int outEdge = isInsideOrOnEdge(p, cur); if (outEdge == -1) return cur; // 否则越过与该边相对的邻居 Tri* next = cur->neigh[outEdge]; // 边界情况：穿过边界边时需要特殊处理，见避坑章节 if (next == nullptr) { // 这是边界，说明点可能在凸包外，或者数值误差导致越界 return handleBoundaryCase(p, cur); } cur = next; } throw std::runtime_error("locatePoint: too many steps"); }这段代码的逻辑是：isInsideOrOnEdge返回 -1 表示点在三角形内部或边上（此时可以直接用），返回 0/1/2 表示点在相对于顶点 0/1/2 的外侧。所谓外侧，就是与neigh[i]所在的那一侧相反。每执行一步，我们从当前三角形跨到它与该边相邻的邻居，重复直到找到包含点为止。maxSteps 是防止数值误差导致路径反复震荡、无限循环的兜底，正常情况下均匀分布的点走到目标三角形的步数大概是点密度的立方根量级。
需要提醒的是，isInsideOrOnEdge的实现不能直接用面积坐标判断，因为浮点数精度问题会导致点在边上时被判定为内部或外部，进而引发后面三角形拆分的错误。我用的方法是分别计算点与三条边的叉积符号，设定一个很小的 epsilon（比如 1e-9 量级，具体取决于点坐标的量纲），只要叉积绝对值小于 epsilon 就认为是「在边上」，这样能减少由于浮点误差引起的误判。这个 epsilon 的取值是个玄学，我后面专门写一节讲讲怎么定。
为什么选 walking 而不是网格索引？主要因为增量场景下，每次插入点后，三角网局部变化，网格索引里的三角形列表需要同步更新，这个维护成本不小；而 walking 不需要额外空间，只需要一个可靠的起始三角形。对于点集按空间顺序输入的情况（比如扫描线逐行扫描得到的点），上一次插入点所在的三角形往往就在新点附近，walking 几乎一步就能到达，这才是真正的增量优势。如果点完全是随机的，walking 的性能会退化，但那也已经比全量遍历快了。
2.3 从点定位到拆分三角形：把新增点塞进网格的关键一步
找到包含新点的三角形后，增量剖分的第一步是把这个三角形拆成三个小三角形。假设原三角形是(a, b, c)，新点是 p，拆出来的三个三角形分别是(a, b, p)、(b, c, p)、(c, a, p)。这时候要更新原三角形的邻居关系，并设置新三角形之间的邻居。这一步看起来简单，但邻居指针非常容易接错。
我的实现顺序是：先创建一个新三角形 T0（存储 a、b、p），再创建 T1（存储 b、c、p），再创建 T2（存储 c、a、p），原来的三角形 T 标记为删除。然后手动调整邻居指针。这里有个容易翻车的细节：T 原来的邻居需要正确地分配给相应的新三角形。比如 T 的neigh[2]是与边ab共边的外部三角形，而 T0 的边ab就是原来 T 的边ab，所以这个邻居应该给 T0。同理 T 的neigh[0]对应于边bc，要交给 T1；T 的neigh[1]对应于边ca，要交给 T2。
下面是拆分时的邻居接线代码片段，我专门写了注释：
// T 是原来的三角形，p 是新插入点 // 建立 T0 = (a,b,p), T1 = (b,c,p), T2 = (c,a,p) Tri* T0 = new Tri(a, b, p); Tri* T1 = new Tri(b, c, p); Tri* T2 = new Tri(c, a, p); // 设置三个新三角形之间的邻居关系 // T0 的边 bp 与 T1 的边 pb 共享，这里用索引匹配 T0->neigh[2] = T1; // T0 的边 a-b 外侧是原来的邻居，稍后处理 T0->neigh[1] = T2; T1->neigh[0] = T2; T1->neigh[2] = T0; T2->neigh[0] = T1; T2->neigh[2] = T0; // 原 T 的外部邻居分配到新三角形 T0->neigh[2] = T->neigh[2]; // 边 ab T1->neigh[0] = T->neigh[0]; // 边 bc T2->neigh[1] = T->neigh[1]; // 边 ca正如注释里体现的，这一步最需要小心的是「哪个新三角形的哪条边对应原三角形的哪条边」。我不建议靠记忆，而是把三角剖分在纸上画出来，标好顶点顺序再写代码。实际调试时，我吃过一次大亏：把 T0 和 T1 的邻居接反了，导致后续的空外接圆检查把原本不需要翻转的三角形翻来翻去，最终产生破洞。所以这一步值得停下来，用三个实际点坐标画个图核对一遍。
拆分完成后，还不能收工。新插入的点 p 出现在三个新三角形里，但这三个新三角形未必满足 Delaunay 性质——它们的邻居三角形可能有空外接圆冲突。这就是下一章要讲的 Lawson 翻转，也是整个算法最核心的环节。
3. 局部优化与 Lawson 翻转：每插一个点，如何让网格重新满足 Delaunay 性质
3.1 空外接圆判定，用行列式而不是角度
增量 Delaunay 的第二大步是「修复」。新点 p 插入后，需要检查与新三角形相邻的所有三角形，看看它们的外接圆是否包含了 p。如果包含，说明共享边不合法，需要将两个三角形的公共边翻转（即删除这条公共边，连接两个非公共顶点，形成新的两个三角形），然后继续检查新生成的三角形，直到所有相邻三角形都满足空外接圆性质。这一过程称作 Lawson 翻转。
判断一个点 d 是否在三角形(a,b,c)的外接圆内，最经典的方式是计算四个点构成的矩阵的行列式。我一开始用了角度法——比较 ∠acb 和 ∠adb 的大小，结果因为三角函数和反三角函数的精度问题，在共圆或接近共圆时不断误判，翻转算法出现振荡，甚至死循环。后来换了行列式法，稳定性好很多。这个四行列式的定义是：
// 返回 true 表示 d 在三角形 (a,b,c) 的外接圆内（或圆上，按严格排除处理） bool inCircle(Point* a, Point* b, Point* c, Point* d) { double ax = a->x - d->x, ay = a->y - d->y; double bx = b->x - d->x, by = b->y - d->y; double cx = c->x - d->x, cy = c->y - d->y; double det = (ax*ax + ay*ay) * (bx*cy - by*cx) - (bx*bx + by*by) * (ax*cy - ay*cx) + (cx*cx + cy*cy) * (ax*by - ay*bx); return det > 0; // 逆时针三点时，det > 0 表示 d 在内 }这个函数的参数顺序和方向性极其重要。这里的(a,b,c)是某个三角形的三个顶点，并且要求(a,b,c)是逆时针顺序。如果原始数据给的点是顺时针，需要先做一次叉积检测并交换任意两个顶点把它们改成逆时针。否则这个行列式的符号会完全反转，翻转判定全反，最终生成的不是 Delaunay 剖分而是某种奇怪的剖分。我因为坐标系不同（Y 轴向上还是向下），在这个方向上翻过车，后面避坑章节会细说。
det 的数值大小和点坐标的数量级直接相关。如果点坐标都是小数（比如 0.001 量级），det 可能很小，和误差混在一起；如果点坐标是几万的大数，det 又可能溢出。因此我通常把点坐标归一化到一个合适的范围，比如将所有点平移到原点附近，再缩放到 [-1, 1] 区间。这个归一化不影响 Delaunay 结构（因为空外接圆性质在相似变换下保持不变），但能显著改善数值稳定性。我一般会在读入点集后做一次平移缩放，不然det的计算很容易出现精度灾难，尤其在共圆点很多的情况下。
3.2 递归翻转的过程与终止条件：如何避免死循环
拆完三角形后，我把所有新生成的内部边（即包含 p 作为顶点的边，或者与 p 相对的外边）放入一个待检查队列。常见做法是用一个std::vector<Tri*>当栈，每次从栈里弹出一对共享边的三角形，判断是否需要翻转。如果需要，则删除这两个旧三角形，创建两个新三角形，并把新产生的边（新三角形的两条公共边和与 p 相关的边）重新压入栈中。
我实际用的伪代码下面这段，是围绕队列展开的：
std::vector<Tri*> queue; // 初始时将包含 p 的三个三角形的三条内部边对应的邻居对加入队列 // 这里用边（一个面与它的邻居）表示 queue.push_back(T0); queue.push_back(T1); queue.push_back(T2); while (!queue.empty()) { Tri* t1 = queue.back(); queue.pop_back(); // 如果 t1 被删了，跳过 if (t1->isDeleted) continue; for (int i = 0; i < 3; i++) { Tri* t2 = t1->neigh[i]; if (t2 == nullptr || t2->isDeleted) continue; // 找共享边的两个端点 v1, v2 int v1, v2; getSharedEdge(t1, t2, v1, v2); // 找各自的第三个顶点 int p1 = thirdVertex(t1, v1, v2); int p2 = thirdVertex(t2, v1, v2); // 判断 p2 是否在 t1 的外接圆内（或者 p1 是否在 t2 的外接圆内，两者等价） if (inCircle(t1->v[0], t1->v[1], t1->v[2], pointOfIndex(p2))) { // 执行翻转 flip(t1, t2); // 翻转后新生成的三角形需要继续检查 queue.push_back(newTri1); queue.push_back(newTri2); break; // 重新处理新的邻居集合 } } }这段代码有几处细节需要说明。第一，getSharedEdge要根据两个三角形各自的顶点索引找出共享边的两个端点，这个可以直接通过顶点坐标比较得到，也可以借助预先建立的顶点到三角形邻接表。我选择直接比较坐标值，因为三角形数量在后期的局部区域内很少，比较成本可忽略。第二，flip(t1, t2)并不是简单地交换顶点，而是要正确更新四个三角形的邻居关系，比拆分还要繁琐一点。第三，break很重要：因为翻转后 t1 和 t2 已经不存在了，原来的 for 循环里的邻居信息全部失效，所以必须跳出来重新从栈中取新的三角形做检查。如果不 break，继续循环会使用已经被删除的三角形指针，轻则读到垃圾数据，重则直接段错误。
翻转的终止条件理论上是有限的——Delaunay 翻转的局部序列必然在有限步内结束，这一点有数学保证。但在浮点数环境下，如果inCircle判断接近 0 时忽正忽负，可能出现两个三角形之间反复翻转的死循环。我在flip函数里加了一个计数器：如果某个共享边在同一轮里被翻转超过一定次数（比如 10 次），就认为检测到数值振荡，强制使用共圆判定（即det >= -epsilon）来打破循环。这个方法不优雅但有效，实际运行中很少触发，触发的时候通常说明 epsilon 设得不对或者点坐标有重复。
3.3 翻转时邻居指针的重新接线：把图画出来再写循环
flip的接线逻辑是整个项目中最容易写错的地方，我在这里花的时间比前面所有部分都多。设两个共享边的三角形是 T1 = (a, b, c) 和 T2 = (a, b, d)，共享边是 ab，c 和 d 分别在 ab 的两侧。翻转操作就是删掉边 ab，连接 cd，形成 T3 = (a, c, d) 和 T4 = (b, c, d)。这个过程会改动以 ab 为公共边的两个三角形，同时 T1 的另外两个邻居（一边为 ac，一边为 bc）和 T2 的另外两个邻居（一边为 ad，一边为 bd）都需要重新指派。
我实现的时候并不是直接修改 T1 和 T2，而是创建两个新三角形 T3、T4，标记 T1、T2 为删除，然后依次设置 T3 的三个邻居：其中一个应该是 T4，另外两个分别是 T1 和 T2 的对应邻居。这里的关键在于「对应」二字。画图才能看清：
- T3 = (a, c, d)：它的边 ac 的原邻居是 T1 中边 ac 的邻居；边 ad 的原邻居是 T2 中边 ad 的邻居；边 cd 的邻居是 T4。
- T4 = (b, c, d)：它的边 bc 的原邻居是 T1 中边 bc 的邻居；边 bd 的原邻居是 T2 中边 bd 的邻居；边 cd 的邻居是 T3。
除了邻居指针，T3 和 T4 内部的顶点顺序也必须保持逆时针，否则后续的inCircle会出错。所以我再强调一遍：所有三角形的顶点顺序统一为逆时针，这是一个全局约定。一旦出现顺时针的三角形，inCircle的符号判断就会反过来，翻转逻辑会完全崩坏。
void flip(Tri* t1, Tri* t2) { // 找出共享边的两个顶点和各自的第三个顶点 int v1, v2, third1, third2; getSharedEdgeAndThirds(t1, t2, v1, v2, third1, third2); // 创建新三角形，注意顶点顺序保持逆时针 // 这里假设 v1, v2, third1, third2 已经按逆时针顺序排列好 Tri* t3 = new Tri(v1, third1, third2); Tri* t4 = new Tri(v2, third1, third2); // 接线 t3 与 t4 t3->neigh[?] = t4; t4->neigh[?] = t3; // 从 t1, t2 的邻居接过来 // 具体索引以逆时针顶点顺序来定 // 标记旧三角形已删除 t1->isDeleted = true; t2->isDeleted = true; // 把新三角形加入全局列表 mesh->insertTri(t3); mesh->insertTri(t4); }上面的?不是真的留空，而是表示你必须在实现时根据你定义的顶点顺序填写。我给不出统一答案，因为每个人定义的三角形顶点索引顺序不同。这也是我为什么建议你拿一张纸，把翻转前后的图形画出来，在每个三角形内部标号 0/1/2，然后对着图写代码，比凭空硬记要可靠得多。
4. 从点集合到凸包边界：超三角形与边界维护的细节
4.1 超三角形的选法与大小：不能太小，也不能大到让数值崩溃
增量 Delaunay 算法的一个重要前提是，初始时需要一个足够大的三角形，把所有输入点都包在里面。这个三角形叫做超三角形（super-triangle）。它本身不是最终网格的一部分，等所有点插入完毕，需要删除所有包含超三角形顶点的三角形，剩下的就是点集的 Delaunay 三角剖分（凸包内部）。
超三角形怎么选，直接影响后续数值稳定性。教科书上只说「取一个足够大的三角形」，但这个「足够大」在实践中有讲究。如果超三角形太小，导致有些点落在它外面，算法会直接把点丢掉或者产生未定义行为；如果太大，比如取 (0,0), (1e9,0), (0,1e9)，那么在计算空外接圆行列式时，超三角形顶点坐标的平方会达到 1e18，和其他点（坐标在 1 附近）的数值相差悬殊，浮点数精度会严重丢失，导致本该判定为在圆内的点被判定为在圆外，最终得到错误的三角剖分。
我一般这样选取超三角形：先遍历所有点求得包围盒 (minX, maxX, minY, maxY)，然后把包围盒中心平移到原点，再把包围盒尺寸放大一定倍数。常见的做法是取包围盒对角线长度的 20 倍作为超三角形的边长，使超三角形的三个顶点分别位于：
// 假设包围盒中心在 (cx, cy)，半径 r = 20 * 包围盒对角线长度 double r = 20.0 * sqrt((maxX-minX)*(maxX-minX) + (maxY-minY)*(maxY-minY)); Point* s1 = new Point(cx - 2*r, cy - r); Point* s2 = new Point(cx + 2*r, cy - r); Point* s3 = new Point(cx, cy + 2*r);为什么用 20 倍而不是 2 倍或 100 倍？这是我多次实验后的折中：20 倍足够大，能覆盖所有可能的外部点（即使点分布在凸包外很远，也不至于超出超三角形）；同时它的坐标数值在r不超过 1e6 量级时，平方后是 1e12 量级，还在 double 的安全精度范围内。如果你的点坐标本身已经很大（比如地理坐标），那么先做归一化再选超三角形是最稳妥的。我通常会写一个normalizePoints()函数，把所有点平移到原点、等比缩放到 [-1, 1] 范围，剖分完成后再把三角形顶点坐标映射回去。后续如果读者需要做地形网格，这种归一化处理是基本功。
超三角形选定后，初始三角网就只有一个三角形，即超三角形本身。所有三个顶点都标记为「虚拟点」，在最后清理阶段，只要三角形的任何一个顶点是虚拟点，就删除这个三角形。清理之后剩下的三角形构成凸包内部的三角剖分。注意：如果点集本身有共线点或者有若干个点严格在凸包边上，清理后可能会留下一些退化三角形，这些需要在后处理中根据边长阈值或面积阈值去除。
4.2 插入点在凸包外部时怎么办：重新定位到超三角形边界
增量插入过程中，点不一定总落在当前三角网内部。因为当前三角网可能只是部分点集的凸包，新点可能落在凸包外部，但仍然在超三角形内部。此时 walking 定位会走到边界三角形附近，发现目标点在当前三角网的外侧。还好超三角形足够大，新点依然处于超三角形内部，所以不会找不到位置，只是需要越过当前凸包的边，进入更外层的三角形直到找到包含点。
被我标记为handleBoundaryCase的地方要做的事情，就是当 walking 走到一个三角形，且目标点在这个三角形的外侧边所对应的邻居是nullptr（即这条边是当前三角网的凸包边界）时，不能返回错误，而是要继续往更外侧找。但此时当前三角网的边界不是超三角形的边，而是普通三角形的边，所以在边界上不存在邻居。这该怎么办？
我的方案是：如果next == nullptr且当前三角形不是超三角形自身，就把这个三角形从网格中删除，然后用一个「虚拟三角形」继续定位。虚拟三角形的顶点由超三角形顶点和当前凸包边界上的某些点组成，这些虚拟三角形只用于定位，不参与最终的 Delaunay 翻转。具体实现略复杂，但有一个更简单的替代方案：在所有点插入之前，先找到凸包外围的几个点，将超三角形扩充为超四边形或更多边形的套壳，这样点定位时永远有可走的邻居，不会碰到 nullptr。不过这个方案需要额外维护多边形，代码量更大。
我实际采用的是一种折中：在 walking 遇到 nullptr 时，以该边界边为基准，把它向超三角形方向「扩出去」生成一个临时三角形，这个临时三角形只用于定位，等找到包含目标点的三角形后，再把这个临时三角形删除。这种临时三角形的创建和删除虽然频繁，但每次插入点最多遇到一两次，性能代价可忽略。关键是代码里要区分「正常的三角形」和「临时定位三角形」，我通过一个bool isGhost标志来标记。凡是isGhost的三角形不参与 Delaunay 翻转，也不进入最终的网格输出。这样一来，插入外部点时的逻辑就和插入内部点完全一致了：定位、拆分、翻转，只是翻转过程中 ghost 三角形不参与inCircle检查，也不需要翻转。
4.3 超三角形顶点的删除顺序：先清 ghost，再清虚拟点，别搞反
最后清理阶段，我踩过一个很隐蔽的坑：如果先删除了包含超三角形顶点的三角形，再去清理 ghost 三角形，可能会出现 ghost 三角形引用已删除顶点导致野指针问题。正确的顺序是：先遍历所有三角形，把isGhost == true的全部标记为删除；然后遍历一遍剩余三角形，删除任何包含超三角形顶点索引的三角形；最后再做一次碎片整理，把链表中未被删除的三角形搬到一个连续的 vector 里，输出给下游使用。
清理代码大致如下：
for (Tri* t : allTris) { if (t->isGhost) t->isDeleted = true; } for (Tri* t : allTris) { if (t->isDeleted) continue; for (int i = 0; i < 3; i++) { if (t->v[i] >= numRealPoints) { // 超三角形顶点的索引排在实际点之后 t->isDeleted = true; break; } } } // 最后输出未删除的三角形 std::vector<Tri*> finalTris; for (Tri* t : allTris) { if (!t->isDeleted) finalTris.push_back(t); }这里有个假设：我给真实点和虚拟点分配的顶点索引是不同的范围，虚拟点索引从numRealPoints开始。这样做的好处是删除判断只需要比较索引，不用存储指针。如果你用 Point 指针直接作为顶点标识，就需要额外维护一个虚拟点指针集合。两种方式我都试过，索引方式更简洁，特别是后续做网格拓扑分析时，整型索引比指针更友好。当然，索引方式要求你最初读入点集时就知道总点数，如果点数是动态增长的，则需要一个动态索引分配器，这也不难，但会让代码不那么直白。
5. 常见问题与避坑清单：我花了三天才搞定的五个坑
5.1 三角形顶点顺序不统一：逆时针约定被破坏后，inCircle判定全反
现象：算法跑起来不会崩，但生成的三角网出现大量「交叉三角形」——三角形的边的方向错乱，网格看起来像一团打结的线。用输出 OBJ 文件后，在渲染器里看到的是黑色破面。
原因：插入三角形或翻转时，新三角形的顶点顺序没有保持逆时针。inCircle行列式的符号依赖于顶点绕序，一旦某个三角形变成了顺时针，它的外接圆判定符号就反向，该翻转的不翻转，不该翻转的翻个不停。更麻烦的是，这个错误不会立即暴露，而是在后续多次插入后累积成拓扑混乱。
解决：我在所有创建三角形的地方统一调用一个makeCCW(Tri*)函数，传入三个顶点后计算叉积，如果为负则交换前两个顶点索引。这样无论上层怎么传参，最终存储的三角形都是逆时针。同时我写了一个调试断言：每次翻转或拆分后遍历所有三角形的顶点叉积，如果发现负值，立即中断并打印所在函数和报错信息。虽然断言在 Release 模式会被禁用，但 Debug 模式下能快速定位是哪一步产生了顺时针三角形。这个习惯帮我节省了大量排查时间。
5.2 点在边上被误判为内部或外部：一个 epsilon 引发的破洞
现象：插入点正好落在某条边上时，偶尔会出现点被插入两次的情况，导致网格中有一条边被重复使用，最终形成破洞或 T 形裂缝。
原因：isInsideOrOnEdge判断点在边上的条件用了绝对 epsilon，比如fabs(cross) < 1e-9。但当点坐标被放大到 1e4 量级时，叉积的正常大小也在 1e4 量级，1e-9 的 epsilon 相对太小，导致原本在边上的点被判定为在边的一侧；而当点坐标缩小到 1e-4 量级时，叉积本身只有 1e-4，1e-9 的 epsilon 又太大，导致离边还有一定距离的点被误判为在边上。
解决：不要用绝对 epsilon，改用相对 epsilon：fabs(cross) < epsilon * (|AB| * |AP|)，其中 AB 是边向量，AP 是从 A 到目标点的向量。或者更简单，先把所有点归一化到 [-1, 1]，然后使用一个固定的相对 epsilon，比如 1e-12。我在归一化后，叉积的尺度大约在 1 附近，此时 epsilon 设为 1e-12 就足够应付大多数情况。如果你有大量共线或近共线点，建议再降低到 1e-13，但再低就要小心 double 本身的限制。
5.3 相邻关系更新时遇到已删除三角形：野指针与死循环
现象：程序偶尔会崩溃，报错在访问t->neigh[i]->v[j]时出现段错误。或者在翻转循环中，两个三角形反复相互翻转，CPU 飙到 100%。
原因：删除三角形时只做了标记，但它的邻居指针还在被其他三角形引用。如果处理队列时没有判断isDeleted，就会拿到一个已标记删除的三角形，它的顶点数据可能已被释放（如果用的是 new/delete），或者数据虽然还在但状态已经不属于当前网格，导致逻辑混乱。
解决：所有从队列或邻居遍历中取到的三角形，先检查isDeleted，是则跳过。这是一个万能的安稳解法。另外，我还在删除三角形时将所有邻居指针置为 nullptr，这样即使访问到已删除三角形，顶多是空指针崩溃，而不是读到随机数据，这能更快定位问题。翻转死循环的问题，除了加计数器，还可以在 debug 模式下记录每次翻转的边端点，如果同一组端点反复出现 10 次以上，就打印警告并停止该条边的翻转，把网格保存出来供分析。
5.4 凸包边界上的点被忽略或错误连接：超三角形太小或 ghost 处理不当
现象：输入点在凸包边界上有较多分布时，剖分结果缺少边界上的某些点，或者边界三角形被错误地翻转成内部形状，导致凸包形状变凹。
原因：最常见的是超三角形太小，边界点落在了超三角形之外，导致这些点根本没有被插入。另一种可能是 ghost 三角形参与了 Delaunay 翻转，ghost 三角形的外接圆范围极大，使得边界边的翻转条件被错误触发，局部优化越改越乱。
解决：检查超三角形的半径是否真的是包围盒对角线长度的 20 倍，并确认所有点都在超三角形内部。对于 ghost 三角形，我单独维护一个std::set<Tri*> ghostSet，所有 ghost 三角形的inCircle判定直接跳过，同时在翻转时的邻居配对里，如果遇到 ghost 邻居，就不单独翻转该边，而是强制接受当前状态，等后续点插入时再自然修正。这里没有完美的数学解，工程上这样做能得到稳定的结果。
5.5 输出三角形方向与预期相反：坐标系的 Y 轴朝向背锅
现象：最终三角形导入到 OpenGL 或 DirectX 后，背面剔除开启时有一半三角形看不见。换到别的软件又是好的。
原因：Delaunay 剖分本身只关心三角形集合，不关心法线方向。如果你使用 Y 轴向上（数学坐标系），生成三角形的顶点顺序是逆时针，但在 Y 轴向下（屏幕坐标系）的应用里显示时，顶点顺序在屏幕空间看起来变成顺时针，导致渲染背面剔除失效。
解决：输出时根据目标平台做一次顶点顺序调整。通用的做法是输出 OBJ 文件时把索引改为交换前两个：index[0], index[2], index[1]，这样就能在 Y 轴向下的坐标系里保持逆时针。如果你的下游是渲染引擎，建议在 mesh 类里提供一个flipOrientation()方法，调用后遍历所有三角形交换前两个顶点索引并重新计算法线，这样在剖分层保持标准逆时针，在渲染层按需翻转，职责分明。
6. 用 Debug 绘制和随机点测试验证实现：一个让我跪过的「共圆点」教训
如果你按前面的章节把代码写完，接下来最重要的事情不是急着接真实数据，而是先做正确性验证。我强烈建议你花一个下午做一个可视化的调试器：把当前三角网的三角形边画到窗口上，每个插入点用一个醒目的颜色标出，每执行一次翻转刷新一次画面。这种直观反馈比任何单元测试都能更快暴露问题。常见做法是用 OpenGL 或 SDL 画线框，也可以把每一帧导出成 SVG，然后用浏览器逐帧查看。我偏向于导出 SVG，因为可以记录完整的插入和翻转过程，出问题时不用重新复现。
下面是我用来做随机点测试的一段核心代码，它生成 1000 个随机点，执行剖分后检查三条不变量：
// 验证1：每个三角形顶点都是逆时针 for (Tri* t : finalTris) { double det = cross2D(t->v[0], t->v[1], t->v[2]); if (det <= 0) { std::cerr << "Clockwise triangle detected!\n"; return false; } } // 验证2：每条内部边被两个三角形共享，边界边只被一个三角形拥有 std::map<std::pair<int,int>, int> edgeCount; for (Tri* t : finalTris) { for (int i = 0; i < 3; i++) { int a = t->v[i], b = t->v[(i+1)%3]; if (a > b) std::swap(a, b); edgeCount[{a, b}]++; } } for (auto& kv : edgeCount) { if (kv.second > 2) { std::cerr << "Edge shared by >2 triangles: " << kv.first.first << " " << kv.first.second << "\n"; return false; } } // 验证3：空外接圆性质在数值容差内成立 for (Tri* t : finalTris) { for (Point* p : points) { if (p == t->v[0] || p == t->v[1] || p == t->v[2]) continue; if (isInsideCircumcircle(t, p, 1e-9)) { std::cerr << "Empty circle violated!\n"; return false; } } }第三个验证的复杂度是 O(n * t)，对于 1000 个点还可以接受，超过一万点就要改用相邻三角形验证而不是全点遍历。验证时isInsideCircumcircle要用一个宽松的 epsilon（比如 1e-9），因为我们允许数值误差导致极少数退化情况。如果验证失败，不要急着改 epsilon，先定位到具体是哪一步插入或翻转导致的不变量破坏，大多数情况下是你的邻居指针接错了。
我最后想说一个让我跪了很久的教训：共圆点。当输入点集中有四个或更多点严格共圆时，Delaunay 三角剖分不是唯一的——不同的插入顺序会产生不同的合法剖分。我的实现没有做特殊处理，结果在一组正方形四个角加中心点的测试数据上，生成的三角形总有一条对角线是「随机」的，用 A 点作为起始点是一种结果，换 B 点作为起始点又是一种结果。这不是 bug，但如果你下游的网格生成需要形状一致，就必须对共圆情况做打破平局的规则：常见做法是给inCircle判定增加一个微小的扰动值，或者按点编号作为 tie-breaker，在 det 的绝对值小于某个epsilon 时，固定认为点不在圆内。这样能保证同一个点集无论插入顺序如何，最终生成的三角剖分在拓扑上一致。
写到这里，我的建议是千万别跳过验证环节直接上生产。我自己第一次写完觉得没问题，拿着真实地形点去跑，结果有 0.1% 的三角形出现了自交，后来用可视化调试器一看，是翻转时一个邻居指针没接对，在十几次插入后才引发连锁错误。现在我的习惯是每次改完数据结构或翻转逻辑，先跑一遍随机点测试，再跑一遍共圆点测试，最后才放真实数据。如果这套流程也能帮你在调试时少走弯路，我会很高兴，希望帮到你。
本文还有配套的精品资源，点击获取
