简介：三角剖分是点集三角化领域的重要算法，在有限元分析、计算几何与计算机图形学中常作为网格生成与空间剖分的预处理步骤。这份C++实现围绕Delaunay三角剖分的基本原则展开，适合需要了解或集成该算法的开发者，尤其适合数值分析或图形学方向的入门与二次开发。压缩包内共2个文件，包含1个cpp实现文件与1个头文件，整体体积约3KB，代码轻量、结构清晰，便于阅读、修改与移植。目前已有3662人学习/下载。阅读源码时可以看到三角剖分的函数接口与实现路径，理解算法如何处理点集输入，并生成满足“最小角最大化”特性和任意四点不共圆唯一性条件的三角网，参考实现还有助于进一步理解Voronoi图、EMST、Gabriel图等与Delaunay相关的几何结构，可直接作为自制三角剖分模块的参考骨架，整体思路清晰，适合快速上手实践。 前阵子帮一个做激光点云处理的朋友调Delaunay三角剖分代码，他用的C++，自己照着论文写了个Bowyer-Watson增量法，结果一跑起来就出各种诡异三角形：有的点明明在外接圆内却判不出来，有的区域出现重叠三角形，有的三角形长成“针形”直接把后续网格生成搞崩。我帮他Debug了两天，最后定位到是数据结构和几何精度两个方向的问题。今天就把这套完整思路和实现细节写出来，从原理到C++代码骨架，再到性能优化和踩坑实录，一次性说透。
这个内容适合什么人看？你如果要做地形TIN生成、游戏地图导航网格、点云网格化、逆向工程或任何需要把散乱点集变成“良好形状三角形网格”的工作，这篇文章都值得看完。不需要你有计算几何基础，但至少要熟练C++的基本语法和STL容器操作。
1. 原理先行：Delaunay三角剖分到底在干什么
1.1 一个容易被忽略的前提：空外接圆与最大化最小角
Delaunay三角剖分的定义听起来简单：给定平面上一组点，用不相交的三角形把这些点连起来，且保证任意一个三角形的外接圆内不包含点集中其他点。这个“空外接圆”性质是整个算法的灵魂。
但很多人没意识到，这个性质背后还藏着一个更实用的推论：在所有可能的三角剖分中，Delaunay剖分能够最大化所有三角形的最小内角。换句话说，它能尽量避免那种“瘦长条”三角形，让网格整体形状更均匀。这一点在数值计算里至关重要，比如有限元网格、流体模拟，如果网格里有极端细长的三角形，求解器很容易数值发散。
我见过的实际项目中，至少有三类场景会优先考虑Delaunay：
- 地形建模：从激光雷达或无人机影像提取的点云，需要生成不规则三角网（TIN），Delaunay能保证地形起伏表达不出现“尖刺”三角形。
- 路径规划：游戏AI寻路或机器人导航用的NavMesh，剖分后三角形质量直接决定路径平滑度。
- 逆向工程与3D重建：从扫描点云重建曲面时，高质量的三角网格是所有后续处理的基础。
还有一个数学性质值得记住：Delaunay三角剖分和Voronoi图互为对偶。Voronoi图的每个单元格顶点连起来就是Delaunay三角形。如果你后续要做最近邻分析、影响区域划分，这两个结构经常要一起算。
1.2 算法选型：为什么我选了Bowyer-Watson而不是分治法
Delaunay剖分的主流算法主要有三种：
| 算法 | 平均时间复杂度 | 实现难度 | 适用场景 |
|---|---|---|---|
| 逐点插入法（Bowyer-Watson） | O(n^2) | 低 | 十万点以内，动态加点场景 |
| 分治法 | O(n log n) | 高 | 百万级点集，一次性构建 |
| Sweep-line算法 | O(n log n) | 高 | 静态数据，对性能有硬要求 |
我个人的建议是：如果你的数据量在十万点以内，无脑选Bowyer-Watson。它的最大优势是逻辑清晰、实现简单，而且天然支持动态增量插入。什么意思？就是你可以随时往已有的三角网格中加一个新点，不需要重新剖分所有点，这在实际项目里太常用了。
分治法虽然渐进复杂度更低，但它的实现涉及递归划分、跨子问题的合并、基准边搜索，代码量至少是Bowyer-Watson的三倍，而且边角情况特别多。我见过很多朋友折腾分治法最后在合并阶段怎么调都调不对，基本都弃坑了。除非你是百万级点集的离线处理，否则增量法完全够用。即便要处理百万点，也建议先了解数据分布，很多场景可以用分块加去噪把点集缩小到十万级别。
2. 动手前的关键设计：数据结构与几何判断
2.1 点、边、三角形怎么组织最顺手
写C++实现时，第一步不是写算法，是选数据结构。数据结构设计得好，后面所有环节都能少踩很多坑。我推荐的方案：
struct Point { double x, y; int id; // 便于追踪调试 }; struct Triangle { int v[3]; // 三个顶点索引 int neighbors[3]; // 三个邻接三角形索引，-1表示边界 bool alive; // 标记是否已经被删除 };点用索引而不是直接存坐标值，这是关键设计。因为算法中删除和新建三角形非常频繁，如果每次复制坐标，性能会很差，而且容易因浮点比较产生不一致问题。用索引的好处是数据只存一份，三角形只存引用，内存访问也连续。
邻接关系neighbors[3]是可选的，但强烈建议从一开始就维护好。它对于遍历边界、寻找相邻三角形、后续做等值线或灰度可视化都很有用。维护成本也很低，更新三角形时顺手赋值即可。
2.2 inCircle测试：整个算法的灵魂
Bowyer-Watson算法的核心判断就是“一个点是否在某个三角形的外接圆内”。这个判断决定了哪些三角形要被删除，它错了，后面全部跟着错。
判断方法最常用的是行列式法。给定三角形三个顶点 A(ax, ay)、B(bx, by)、C(cx, cy) 和待插入点 D(dx, dy)，计算如下行列式：
double incircle(double ax, double ay, double bx, double by, double cx, double cy, double dx, double dy) { double adx = ax - dx; double ady = ay - dy; double bdx = bx - dx; double bdy = by - dy; double cdx = cx - dx; double cdy = cy - dy; double det = (adx*adx + ady*ady) * (bdx*cdy - cdx*bdy) - (bdx*bdx + bdy*bdy) * (adx*cdy - cdx*ady) + (cdx*cdx + cdy*cdy) * (adx*bdy - bdx*ady); return det; // < 0 则在外接圆内 }你可能注意到这个公式里全是平方项，如果坐标数值很大，比如经纬度坐标或毫米级坐标，行列式结果很容易溢出或者因为浮点精度不足产生错误判断。解决方法是先对点做平移和缩放归一化。我在实际工程里通常是先把所有点平移到以某个中心为原点，然后缩放到[0,1]区间，这样不仅计算稳定，调试时看着也直观。
这里有一个更隐蔽的问题：当D正好落在外接圆上时，行列式等于0。理论上这是“共圆”退化情况，实际中由于浮点误差，结果会在0附近微小摆动。处理方式我会在后面的踩坑实录里详细说，这里先记住一个词：隧道式容差，别用> 0这种硬判断。
2.3 边界提取与“半边思维”
在增量插入过程中，找到所有坏三角形之后，需要提取这些三角形组成的多边形区域的边界边，然后与插入点构造新三角形。
提取边界边的经典做法很巧妙，用一个std::set<std::pair<int,int>>来记录所有坏三角形的边。遍历每条边时：
- 如果这条边第一次出现，插入集合；
- 如果这条边已经在集合里，说明它是两个坏三角形的公共边，不是边界边，删除它。
遍历完成后，集合里剩下来的边就是多边形区域的边界。写代码时要注意边的顶点顺序，建议统一存储时把顶点索引按大小排序，不然(a,b)和(b,a)会被当成两条边重复处理。
std::set<std::pair<int,int>> boundaryEdges; for (auto& tri : badTriangles) { for (int i = 0; i < 3; ++i) { int a = tri.v[i]; int b = tri.v[(i+1) % 3]; if (a > b) std::swap(a, b); auto edge = std::make_pair(a, b); auto it = boundaryEdges.find(edge); if (it == boundaryEdges.end()) { boundaryEdges.insert(edge); } else { boundaryEdges.erase(it); } } }这段逻辑用“如果出现两次就是内部边”的思路，相当于用集合做了一次边计数。实际跑起来效率不错，代码也简洁。如果追求极致性能，可以换成哈希容器，但点量在十万以内时 set 足够。
3. 核心实现：增量插入的全流程
3.1 先搭一个超级三角形
Bowyer-Watson的第一步是构造一个足够大的三角形，把点集中所有点都包住，这个三角形叫做超级三角形。
很多人在这里踩坑。超级三角形太小了，会有一批点落在外面，导致剖分结果缺了一块；超级三角形太大了，外接圆计算会产生数值问题，而且会引入大量最终会被删除的“多余三角形”，影响性能。
我的做法是：先算出点集的包围盒，取包围盒中心作为 c，包围盒对角线长度为 d。然后构造三个顶点为：
Point super[3]; super[0] = {c.x - 20*d, c.y - d, -1}; super[1] = {c.x + 20*d, c.y - d, -1}; super[2] = {c.x, c.y + 20*d, -1};注意20倍这个经验系数。太小了有风险，太远了会溢出。我实测下来20倍在对角线长度在1e3量级的数据上是稳定的，如果你的坐标范围特别大，建议先归一化再建超级三角形。
3.2 插入点的坏三角形检查
每插入一个新点，要做的是遍历当前所有“存活”的三角形，找出外接圆包含新点的三角形。这里有一个性能优化点：你不需要从头到尾遍历所有三角形，可以用上一轮构建的邻接关系从某个种子三角形BFS扩散查找，但实现复杂度会上升。
std::vector<int> badTriIndices; for (int i = 0; i < (int)triangles.size(); ++i) { if (!triangles[i].alive) continue; auto& t = triangles[i]; double det = incircle(points[t.v[0]].x, points[t.v[0]].y, points[t.v[1]].x, points[t.v[1]].y, points[t.v[2]].x, points[t.v[2]].y, pt.x, pt.y); if (det < 0) { // 注意容差，见第5节 badTriIndices.push_back(i); triangles[i].alive = false; } }这段代码里有个细节：标记alive = false时不要真的从容器中删除元素，否则所有索引都会失效。正确做法是先标记，等全部处理完再统一清理，或者用墓碑标记法复用空间。
3.3 重建网格与去重
找到坏三角形集合后，提取边界边，对于每条边界边和插入点构造新三角形：
int newIdx = points.size(); points.push_back(pt); for (auto& edge : boundaryEdges) { Triangle newTri; newTri.v[0] = edge.first; newTri.v[1] = edge.second; newTri.v[2] = newIdx; newTri.alive = true; triangles.push_back(newTri); }注意一个问题：如果点集中有两个点坐标相同，或者三个点共线，那么这里构造新三角形时v[0]、v[1]、v[2]就有可能是共线的，导致退化三角形。所以在算法开始前一定要做预处理去重，把重复点剔除，这能省掉后面大量排查时间。
3.4 完整代码骨架
为了更直观，我把整个增量插入的核心流程用C++伪代码串一遍：
std::vector<Point> pts = loadPoints(); removeDuplicates(pts); // 第一步：去重 // 建立超级三角形 Point super[3] = makeSuperTriangle(pts); pts.insert(pts.end(), super, super + 3); std::vector<Triangle> triangles; triangles.push_back({superIdx[0], superIdx[1], superIdx[2], -1, -1, -1, true}); for (int pi = 0; pi < (int)pts.size(); ++pi) { // 1. 找坏三角形 std::vector<int> bad; for (int ti = 0; ti < (int)triangles.size(); ++ti) { if (!triangles[ti].alive) continue; if (inCircleTest(triangles[ti], pts[pi]) < 0) bad.push_back(ti); } // 2. 提取边界边 std::set<std::pair<int,int>> edges; for (int ti : bad) { markDead(triangles[ti]); for (int i = 0; i < 3; ++i) { processEdge(triangles[ti].v[i], triangles[ti].v[(i+1)%3], edges); } } // 3. 重建三角形 for (auto& e : edges) { triangles.push_back({e.first, e.second, pi, -1, -1, -1, true}); } } // 最后删除所有包含超级三角形顶点的三角形 std::erase(std::remove_if(triangles.begin(), triangles.end(), [&](const Triangle& t) { return !t.alive || t.v[0] >= n || t.v[1] >= n || t.v[2] >= n; }), triangles.end());最后一步判断v[i] >= n用的是插入超级三角形之前的点数量 n，这样能把所有与超级三角形相关的三角形清理干净。注意这里使用的是索引判断而不是坐标判断，可以避免坐标数值比较产生的精度问题。
4. 性能优化与工程化落地
4.1 空间哈希索引加速定位
如果只是10万点以内，上面写的基础版Bowyer-Watson是够用的，时间复杂度大约 O(n^2)，跑起来大概几秒。但如果点集超过20万，几秒会变成几十秒，这时候就需要优化。
最实用的手段是网格空间哈希。做法是：把整个点集所在的包围盒划分成均匀网格，每格记录包含哪些三角形的外接圆信息。插入新点时，先计算新点落在哪个格子，只需要检查该格子以及相邻格子里的三角形，而不是遍历全部。
我实现过一个版本，用均匀网格索引后，20万点从23秒优化到3秒左右，提升接近8倍。网格尺寸怎么选？实践中取点集平均密度对应的直径的2倍比较均衡。太小了格子太多，维护开销大；太大了退化回全量遍历。
4.2 内存和容器的real-world经验
C++实现里一个隐形性能杀手是容器频繁扩容。三角形数量可以达到点数的2倍左右（准确说是2n - 2 - h，h是凸包顶点数），所以你可以在开始时就triangles.reserve(2 * n + 100)，避免vector反复rehash。
还有一个小技巧：用对象池管理三角形而不是直接push_back。因为增量算法中三角形会被“删除”和“新建”，频繁的push和erase会导致内存碎片。对象池配合空闲列表可以显著降低分配开销。对于追求极致性能的场景，这是值得做的。
4.3 与OpenCV等库的集成思路
很多找我调这个东西的朋友，最终目的不只是得到三角形网格，还要展示结果或做后处理。这里给几个实用的集成方向：
- 可视化：使用OpenCV的
polylines画出每个三角形的三条边。对每个三角形取三个顶点，连线画上去即可。 - 填充效果：如果需要看剖分密度，可以用
fillPoly逐三角形填充，颜色根据三角形面积或外接圆半径映射。 - 后续处理：如果要从网格中提取边界轮廓，
findContours可以和三角形邻接关系配合使用，先找到所有边界边，再拼接成多边形。
我之前把Delaunay三角剖分结果和OpenCV的fillPoly结合起来做点云密度热力图，效果非常直观。具体做法就是先剖分，然后计算每个三角形的面积，面积越小说明点越密集，再映射到颜色值填充。
5. 踩坑实录：浮点误差、退化输入与调试技巧
5.1 共圆与退化点：inCircle接近零怎么办
这是所有实现Delaunay的人都会遇到的问题。理论上一组点可以完美共圆，但计算机里浮点数导致incircle的结果在0附近抖动。如果你判断条件是det < 0，有时会得到正，有时负，同一个几何场景在不同插入顺序下会产生完全不同的拓扑结构。
解决办法常用的有三种：
- 轻微扰动：给每个点坐标加上e-9量级的随机噪声，打破退化。简单但会破坏精确性，不适合科研场景。
- 符号容差：设置一个极小阈值，比如
if (det < -1e-12)才判定在外接圆内。这个阈值怎么定？取决于坐标量级，归一化后一般取1e-12比较安全。 - 精确几何谓词：用Shewchuk的
robust predicates，这是计算几何圈标准的鲁棒性方案。代价是引入外部代码，但能彻底解决浮点问题。
我个人的建议是：先用容差法，如果后续发现拓扑错误，直接切换到Shewchuk的精确谓词。网上能下到源码，单独封装成一个头文件即可。
5.2 调试三角剖分的三个常用手段
这算法写出来跑一遍发现不对，怎么定位问题？我总结了一套调试流程：
- 极简用例单步跟踪：构造5-10个点的小点集，每个点的坐标都是整数，然后逐点插入，每步打印当前所有存活三角形的顶点索引。对照手绘草图，很容易看出哪一步开始出错。
- 图形化输出：把每步结果输出成SVG或者PNG。OpenCV几行代码就能画出来，然后用
imwrite保存。我有一次找了两个小时没找到的bug，画成SVG后一眼看出是大三角形边界上的边被错误删除。 - 断言验证：在每个点插入后检查所有三角形是否满足空外接圆性质，复杂度高但仅用于Debug构建。对每个三角形遍历所有其他点做inCircle检查，一旦发现违反就打印当前点索引和三角形索引。
5.3 超级三角形参数的坑
上面提到了超级三角形取包围盒对角线长度的20倍。这个参数不是随便拍的，我最早用2倍，结果有条边离点集太近，点集凸包外包了一圈超长的“贴边三角形”。用5倍在中型数据集上没问题，但点集范围特别不均匀时有风险，后来统一改成20倍才稳定。
另一个实践心得是：判断一个三角形是否与超级三角形相关，一定要用索引判断，不要用坐标值。因为超级三角形的顶点索引可以预先设定，比如superIdx0 = n，那么只要检查t.v[i] >= n就足够了，干净利落。用坐标值判断极容易因为浮点误差把普通三角形误删。
5.4 数据预处理不能省
最后再说一个实战中极其常见的坑：输入数据里带了重复点、极近点或者三个以上点近似共线的退化结构。
我处理过一份地形点云，里面有一个区域扫描了三次，重复点高达3000多个。如果不先去重，剖分结果会出来一堆面积几乎为零的异常三角形，最直接的后果是后续算面积总和时数值爆炸。
去重的方法很简单，先对所有点按坐标排序，把距离小于阈值的点合并或删除即可。阈值取点集平均间距的1/100通常比较合理。这一步一定要在算法开始前做。
写在最后的工程体会
从“能跑”到“跑得稳”这段路，我在这算法上花了挺长时间。最有价值的一个体会是：几何算法的坑大部分不在算法本身，而在数据输入的病理特征和浮点数精度。写C++实现时多做防御性检查，把归一化、去重、边界处理做扎实，比靠算法调参要省心得多。
如果接下来你还想往深了走，可以看三个方向：约束Delaunay三角剖分（CDT），处理带边界约束的多边形内部剖分；3D Delaunay / 表面重建，用三维点云生成四面体网格；以及把Delaunay和Voronoi结合起来做最近邻图、影响区域分析。每个方向都能拆出很多实用技巧，等以后有机会再单独写。
本文还有配套的精品资源，点击获取
