简介：这是一份用C++实现增量Delaunay三角剖分算法的完整工程资源，面向计算机图形学、几何算法与网格处理方向的开发者和学生，重点演示如何通过“半边数据结构”维护动态三角网，并实时保证三角形最小角度最大化。资源包共含89个文件，包括Visual Studio解决方案（.sln、.vcxproj）、C++头文件与实现（.h、.cpp）、可直接运行的exe、GIF演示动图以及编译中间文件等，整包大小约4.17MB，目录结构包含GaussCurvature_Sample示例模块和MeshLib库，便于对照学习和二次开发。核心代码集中在GaussCurvature_Sample的main.cpp中，清晰展示了增量插入点、局部翻转等关键步骤，配合压缩包内附带的GIF动图和编译好的exe，可以直观看到三角剖分从迭代到收敛的动态过程。目前已有998人学习浏览，适合用来掌握Delaunay三角剖分的算法流程、半边数据结构的设计思路，以及Visual Studio下相关工程的组织方式；读者还可基于此扩展自己的网格生成或几何处理算法。 前阵子在做一个点云网格化工具，输入是几千个散乱坐标点，输出需要交给下游做纹理映射和碰撞检测。最初偷懒直接用了现成的开源库，功能没问题，但一遇到定制需求——比如要保留边界、要支持动态增点——就得翻源码改接口，很不顺手。于是我决定自己用 C++ 写一个增量 Delaunay 三角剖分实现，也就是标题里的 DelaunayTriangulation 项目。这篇文章把整个过程完整摊开：核心原理、数据结构选型、关键代码、性能优化，以及调试时踩过的那些坑，给正在学 C++ 或者需要把几何算法落地的同学一个可直接参考的版本。
1. 为什么最终选了"增量法"而不是 Bowyer-Watson
1.1 Delaunay 三角剖分到底在解决什么问题
先明确一下我们要做的事。给定平面上 n 个散点，把这组点用不相交的三角形连起来，覆盖所有点的凸包，这就是一个三角剖分。但三角剖分不唯一，比如四个点构成的正方形，连哪条对角线都合法。Delaunay 三角剖分在其中加入了一条准则：任何一个三角形的外接圆内部，不能包含其他顶点。这条"空外接圆"性质带来一个直观的好处——三角形会尽量饱满、接近等边，不会出现特别尖锐的狭长三角形。对后续做插值、有限元、地形网格生成来说，这个性质非常关键。
Delaunay 三角剖分的实现路线有好几条：分治法、扫描线法、逐点插入法。分治法理论复杂度最优，能达到 O(n log n)，但实现复杂度高，递归切分、合并时需要处理大量跨子集的边，调试阶段很容易崩。扫描线法对排序依赖强，代码也不短。逐点插入法（incremental insertion）虽然最坏情况 O(n^2)，但实现简单、易于理解、还天然支持动态加点，对大多数工程场景已经足够。我不会一上来就追求理论最优，能稳定跑、能改得动，比任何复杂度公式都重要。
1.2 两种逐点插入路线的取舍
同样是逐点插入，现实中还有两条分支。一条是 Bowyer-Watson 算法：每插入一个点，先找出所有外接圆包含新点的三角形，删掉这些三角形形成一个空洞，再把空洞边界上的顶点与新点逐一连接，完成重构。这条路线写起来短，但有个隐患——它需要一次"批量删除"再"批量重建"，空洞区域的边界收集和去重逻辑容易出错，而且对退化情况（多个点共圆、共线）非常敏感。
我采用的是另一条路线：Lawson 算法，也叫翻边法。它的核心思想不是先删再建，而是先找到新点落在哪个三角形里，把它拆成三个新三角形，然后对外接圆性质被破坏的邻接边逐条翻转，直到整个局部重新满足 Delaunay 条件。翻边法的好处是每次操作只影响局部几条边，拓扑变化范围小，特别适合需要动态增点的场景——球员走一步只挪一辆车，不用把整条停车场重新调度一遍。而且调试时每次只看一条边的处理结果，心理负担小很多。
2. 先把数据结构定死：索引、逆时针与三角形栈
2.1 顶点、三角形与邻接关系怎么存
写几何算法，数据结构设计比算法本身更容易决定项目成败。一开始我打算用指针互连的方式组织三角形，每个三角形持有三个指向邻居的指针，结果实现到一半就发现一个问题：vector 扩容会导致指针失效，只要三角形数量超过预分配容量，所有邻居指针都变成悬垂指针，排查起来非常痛苦。
后来我改成纯索引式结构，整个三角网只依赖两个数组：
struct Vec2 { double x, y; }; struct Vertex { Vec2 pos; int id; }; struct Triangle { int v[3]; // 三个顶点索引，约定逆时针顺序 int t[3]; // 三条邻接三角形索引，t[i] 对应边 (v[(i+1)%3], v[(i+2)%3]，-1 表示边界 };用索引代替指针，vector 扩容时只需要重新分配内存，所有索引数值不变，彻底绕开了引用失效问题。代价是取邻居时要多一次数组下标访问，但现代 CPU 对连续内存的访问远快于对散落堆对象的指针追踪，这点牺牲换来了稳定性和缓存友好性，很值。
2.2 为什么一律用"逆时针"约定
坐标结构本身不复杂，真正让代码变脆的是方向约定不一致。我吃过亏：有的函数按"外接圆包含点"判断，有的函数按"点在边的哪一侧"判断，两个函数用了不同的顶点顺序，出来的结果符号恰好相反，bug 找了一下午。后来我立了一条铁规矩：所有三角形顶点一律按逆时针存储。这样"点在边左侧"和"点在三角形内部"的判断共用同一个叉积函数，符号永远一致，不会出现"这边判断为内、那边判断为外"的灵异现象。
逆时针约定对后续核心操作还有一个关键帮助：三角形 t 的顶点是 v[0]、v[1]、v[2]（逆时针），那么边 v[0]→v[1] 的左侧就是三角形内部，而它对应的邻接三角形 t[0] 一定在右侧。翻边时我们经常需要判断"两个相邻三角形是否共享一条边且分居两侧"，有了方向约定，一句话就能写清楚。
2.3 超级三角形的引入与收尾清理
逐点插入法处理的一个前提是：新点必须落在某个现有三角形内部。但第一点插入时，整个三角网还是空的，无处可落。常规做法是准备一个大到能覆盖所有输入点的"超级三角形"（super triangle），让所有点都落在它内部。这个三角形不参与最终输出，处理完所有点之后，只需要删掉所有包含超级三角形顶点的三角形即可。
超级三角形的尺寸选择有个小讲究：不能只刚好包住输入点，最好把范围扩大个几十倍，否则靠近边界的点做外接圆判断时，可能因为超级三角形顶点参与，产生不期望的边界翻转。实测中我把包围盒扩大 20 倍，边界稳定很多。这里补充一点：如果你想让最终剖分结果只覆盖凸包，直接删掉包含超级顶点的三角形就行；如果想保留外部边界形成更大的凸包区域，需要另做边界恢复处理，这个不在本文范围内。
3. 增量插入的完整链路：点定位、拆三角、逐层翻边
3.1 点定位：别急着上高级结构
每次插入一个新点，第一步是找到它落在哪个三角形内部。最朴素的做法是遍历所有三角形，用重心坐标或方向测试判断点是否在三角形内，复杂度 O(n)，点少时完全够用。但当点量超过 2 万时，线性扫描会成为明显瓶颈。这里有一个更聪明的办法：步行法（jump-and-walk）。我们从"上一个插入点所在的三角形"出发，沿边逐步向新点方向移动，判断点在当前三角形的哪条边的外侧，就跨到那个边的邻接三角形，重复这个过程直到点落在当前三角形内部。
步行法在均匀随机点下平均移动步数很少，实测每个点定位大约只需要 5~8 次三角形跳跃，比全量扫描快一两个数量级。实现也不复杂：给定三角形三个顶点 a、b、c，判断点 p 是否在边 ab 左侧，如果在右侧，说明要继续跨越到该边的邻接三角形。循环直到三条边判断都在左侧为止。注意处理点在边上的退化情况，容差范围放小一点，不要用绝对的 0 判断，否则浮点误差会让程序陷入死循环。
3.2 拆分三角形：一次插入产生三个候选
找到新点 p 落在三角形 T 内之后，核心操作是把 T 一拆三：T 的三个顶点分别与 p 连线，生成三个子三角形。原三角形从数组里移除或标记失效，三个新三角形加入数组，同时维护好六条邻接关系：每个新三角形的三条边，两条是内部的、一条是原 T 的外边；新的外邻接继承自原 T。
这里要特别小心的是邻接表的正确性。我实现的时候用一个"待检查队列"保存新生成的三个子三角形，因为新拆出来的三角形，其外接圆性质并不一定满足 Delaunay 条件——它是从一个大三角形里"切"出来的，新点出现在边上，很可能落入相邻三角形的外接圆内。这种破坏只可能发生在共享边的两个三角形之间，所以接下来只需要检查与新点相对的那些邻接边就够了。
3.3 翻边恢复的递归逻辑（incircle 判断与算法证明直觉）
翻边是整套算法的灵魂。设新点 p 所在的某个子三角形为 T1，它有一条边与邻居 T2 共享，T2 中与 p 相对的顶点是 q。如果 q 落在 T1 的外接圆内部，说明这条共享边不满足 Delaunay 条件，需要翻转：把共享边的两个端点保持不变，删掉这条边，改连 p 和 q。
翻转前：T1(a,b,p) T2(b,a,q) 共享边 ab 翻转后：T1'(a,p,q) T2'(b,q,p) 共享边 pq翻完之后产生两个新三角形，它们各自可能又破坏了与更外层邻接三角形的关系，所以递归地对这两个新三角形的外侧边再次做判断，直到没有需要翻转的边为止。这一步对应 Lawson 算法中"legalize edge"的递归过程。它之所以保证能终止，是因为每次翻转都让两个三角形的外接圆半径平方和严格减小，整个系统存在一个单调递减的势能，不会无限翻下去。
4. 核心代码落地：从数学判断到可运行的 C++
4.1 坐标结构、外接圆行列式与浮点容差
先放下数学公式，直接给一个实用的 incircle 判断函数。判断点 p 是否落在三角形 (a,b,c) 外接圆内部，可以使用行列式，它比先求圆心再算距离效率高，而且避免了圆心坐标可能出现的极端值：
// a,b,c 为逆时针，返回 > epsilon 表示 p 在圆内 inline double incircleDet(const Vec2& a, const Vec2& b, const Vec2& c, const Vec2& p) { double ax = a.x - p.x, ay = a.y - p.y; double bx = b.x - p.x, by = b.y - p.y; double cx = c.x - p.x, cy = c.y - p.y; double ab = ax * ax + ay * ay; double bc = bx * bx + by * by; double ca = cx * cx + cy * cy; return ax * (by * ca - cy * bc) - bx * (ay * ca - cy * ab) + cx * (ay * bc - by * ab); }这个函数的精度直接决定算法稳定性。浮点数不是实数，坐标值如果量级很大（比如经纬度坐标），行列式结果可能本身很大，把微小的几何差异淹没在误差里。我这里的经验是：先对坐标做一次归一化平移缩放，让所有点落入 [0, 1] 区间，再参与计算。同时在函数外层加一个 epsilon 比较，不要直接跟 0 比较，见我后面第 6 章的踩坑内容。
4.2 插入主函数的骨架
有了数据结构、点定位函数和 incircle 判断，插入一个点的主流程写起来就非常短了：
void insertPoint(int newVertex, int startTriangle) { int tri = locatePoint(newVertex, startTriangle); // 返回包含新点的三角形 if (tri == -1) return; // 退化：点在边或已存在 int oldV[3]; // 记录原三角形三个顶点 int oldT[3]; // 记录原三角形三个邻接 int nt = splitTriangle(tri, newVertex, oldV, oldT); // 一拆三，返回新三角形起始索引 // 处理邻接关系，让三个新三角形的外侧邻居正确指向 oldT linkAfterSplit(nt, oldV, oldT); // 对与新点相对的边做翻边恢复 legalizeEdge(nt, 0, newVertex); legalizeEdge(nt, 1, newVertex); legalizeEdge(nt, 2, newVertex); }legalizeEdge 的递归实现里，比较关键的是正确找到共享边的"相对顶点"。给定当前三角形 cur，它的第 i 条边是 (v[(i+1)%3], v[(i+2)%3])，邻接三角形是 t[i]。在邻接三角形里找到与这条边相对的那个顶点，本质上是在邻居的三个顶点里找出不等于共享边两端点的那个。这个过程用顶点索引比较即可，注意别用坐标浮点比较，索引比较又快又稳。
4.3 边界条件：重复点、共线点与退化输入
代码能跑通简单用例之后，真正的考验是退化输入。我随机生成数据时遇到最多的问题有三个：
- 重复点：两个顶点完全相同，会导致三角形面积趋近于 0，外接圆判断结果不稳定。处理方式是在初始化时用哈希集合去重。
- 共线点：三个顶点共线时外接圆半径无穷大，incircle 行列式结果接近 0。工程上我把"接近 0"统一当作"不在圆内"处理，避免产生面积为零的退化三角形。
- 点在三角形边上：点定位时点恰好落在边上，拆分会得到面积为零的三角形。我的做法是当点在边上时，把点归入边一侧的较小三角形，同时额外做一次边交换避免退化。
这些边界情况每个看着都很小，但任何一处没处理好，程序跑几万个点大概率会在半路崩掉或生成破洞网格。我建议在实现阶段就准备三个独立测试集：随机均匀点、格子点（大量共线共圆）、手工特例（四点共圆、重复点混合），每完成一个阶段就全量跑一遍回归。
5. 性能实测：从 O(n^2) 到百万点可用的两次优化
5.1 理论复杂度与实测对照
先给个直观参考：我的实现最初只用线性扫描点定位，插入 1 万点耗时约 0.6 秒，5 万点飙升到约 18 秒，呈现明显二次增长。这是因为每插入一点都要扫描全部现有三角形，总复杂度 O(n^2)。虽然小数据看不出问题，但数据量翻五倍，耗时翻三十倍，这就是二次复杂度的可怕之处。
换成步行法点定位之后，均匀随机点下的跳跃次数期望是常数级别，整体复杂度接近 O(n log n)。同样的 5 万点，耗时降到 0.9 秒。再用桶网格（bucket grid）记录每个网格单元内包含的三角形列表，插入时先定位到网格单元，再从单元内的三角形出发做步行，进一步减少跳跃路径长度，10 万点也能稳定在 2 秒内完成。
5.2 桶网格点定位的工程实现
桶网格的思路特别朴素：把包围盒分成若干个格子，每个格子维护一个"与之相交的三角形"列表。插入点 p 时，先算 p 落在哪个格子里，取该格子里的三角形作为步行法起点。如果格子为空或格子里的三角形已经被删除，就往邻居格子扩散找。
实现时注意三点：一是格子数量不要贪多，大致取 sqrt(三角形数量) 即可；二是三角形可能跨越多个格子，需要在插入时把所有相交格子都加入该三角形索引；三是动态删除的三角形要惰性清理——遍历时发现无效三角形直接跳过，不要急着从 vector 里 erase，否则会破坏索引连续性。我踩过的一个坑是每次翻边后忘记更新桶索引，导致定位命中已被删除的三角形，程序看似正常但偶尔输出错误网格。所以我在代码里每次合法化边时，只更新三角形数组本身，桶索引采用"过期检测"策略，遍历时验证tris[tri].valid标记，简单可靠。
5.3 C++ 层面的细节优化（索引、预分配、缓存友好）
- vector 预分配：插入前根据输入数量
tris.reserve(initialGuess * 4)，大幅减少扩容时的拷贝。 - 索引而非对象拷贝：三角形数组元素用固定大小的结构体（三个 int 加三个 int，共 24 字节），拷贝成本很低，但频繁扩容时仍有开销，预分配能完全回避。
- 局部性优化：翻边只操作局部三角形，它们大概率位于 array 中相近位置，缓存命中率比全图遍历好得多。这也是索引式结构的隐形红利。
- 避免动态分配：递归翻边用循环加显式栈（std::vector 当栈用）实现，避免深递归带来的栈溢出和函数调用开销。
实测中这些优化叠加起来，10 万随机点的总耗时约 1.8 秒，50 万点约 11 秒。对于纯 CPU 单线程实现，这个数据已经能覆盖大多数离线网格生成场景。
6. 调试与踩坑实录：可视化是我最快的排错手段
6.1 为什么"看代码看不出 bug"
几何算法有个特点：很多 bug 不会让程序崩溃，而是让输出"看起来差不多"。比如三角形穿过了其他边、某个点的邻接三角形数量不对、网格中出现了细长缝隙。这种错误靠肉眼检查代码很难发现，因为拓扑错误往往是运行时状态累积出来的，跟某一行的逻辑对错没有直接关系。
我第一次跑通算法后，自信满满地输出坐标，用文本方式检查了几个三角形，觉得没问题。结果把数据导到可视化软件里一看，网格中间有一条贯穿始终的裂纹，还有一个点被五个三角形包围而不是应有的六个。这种问题如果只盯代码，再盯三天也未必能看出来。可视化不是可选项，是几何算法调试的必需品。
6.2 三步走：OBJ 导出—MeshLab/OpenCV 渲染—逐帧回放
我的调试流程分三步。第一步，把三角网的三角形数组导出为 OBJ 格式，文件极简：
v x y z f i j kOBJ 格式在 MeshLab、Blender、甚至 Mac 自带预览里都能打开，几行代码就能生成，十分钟之内搭好可视化通道。第二步，用 OpenCV 写一个几十行的实时渲染窗口，把三角形逐帧画出来，这样不仅能看最终结果，还能看到整个插入过程中网格的演化。
第三步是关键——逐帧回放。我会在代码里加一个调试开关，每插入一个点、每次翻边后暂停一下，把当前网格渲染出来。一旦发现某个时刻出现了非法拓扑（比如三角形面积骤变成负数、邻接关系不对称），立刻就能定位到是哪一步引入的错误。做过三次这种逐帧排查之后，我对翻边逻辑的掌控感提升了一个层级，比自己对着调试器打断点高效得多。
6.3 浮点误差导致的"幽灵翻边"与对策
最后聊一个隐蔽的坑：浮点误差引发的幽灵翻边。
incircle 计算里的行列式本质是几个大数的加减，当两个点距离非常近时，行列式的结果会很小，与浮点误差同量级。这时若直接拿结果与 0 比较，偶尔会误判成"点在圆内"，触发一次多余的翻边。多余翻边的直接后果是网格仍然合法（翻边自带局部恢复能力），但会浪费一点时间。真正危险的是反过来：当点已在圆内却因为误差判成"不在圆内"，漏掉必需翻边，最终产生违反空外接圆准则的三角形，网格质量下降。
我的对策分三层：
- 坐标归一化到 [0,1] 区间，消除量级差异；
- 在 incircle 判断外层引入
const double EPS = 1e-12，只有行列式大于 EPS 才认为"在圆内"； - 使用 long double 做中间计算，在极端退化数据下把误差再压低几个数量级。
性能测试里 long double 慢约 20%，但换来的是对怪异输入的鲁棒性。对于离线网格生成场景，这个代价完全可以接受。
还有一个小技巧：定期校验网格的欧拉公式V - E + F = 1（没有边界洞的连通三角网），作为运行时断言。任何一个三角形拓扑被破坏，这个公式立刻失衡，比我盯着渲染画面找缝隙快得多。我后来把这段校验放在每个测试用例收尾时自动执行，此后几乎没有"看起来对但实际错"的尴尬情况。
DelaunayTriangulation 做到这一步，已经从一个玩具演变成我手头可依赖的几何工具。如果你也在写类似算法，我的建议是：先把逆时针约定和索引式结构定死，再实现点定位和翻边，最后把可视化调试尽早搭起来——这个顺序能帮你避开我在前两周踩过的绝大多数坑。
本文还有配套的精品资源，点击获取
