简介三角剖分是点集三角化领域的重要算法在有限元分析、计算几何与计算机图形学中常作为网格生成与空间剖分的预处理步骤。这份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地形建模从激光雷达或无人机影像提取的点云需要生成不规则三角网TINDelaunay能保证地形起伏表达不出现“尖刺”三角形。路径规划游戏AI寻路或机器人导航用的NavMesh剖分后三角形质量直接决定路径平滑度。逆向工程与3D重建从扫描点云重建曲面时高质量的三角网格是所有后续处理的基础。还有一个数学性质值得记住Delaunay三角剖分和Voronoi图互为对偶。Voronoi图的每个单元格顶点连起来就是Delaunay三角形。如果你后续要做最近邻分析、影响区域划分这两个结构经常要一起算。1.2 算法选型为什么我选了Bowyer-Watson而不是分治法Delaunay剖分的主流算法主要有三种算法平均时间复杂度实现难度适用场景逐点插入法Bowyer-WatsonO(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::setstd::pairint,int来记录所有坏三角形的边。遍历每条边时如果这条边第一次出现插入集合如果这条边已经在集合里说明它是两个坏三角形的公共边不是边界边删除它。遍历完成后集合里剩下来的边就是多边形区域的边界。写代码时要注意边的顶点顺序建议统一存储时把顶点索引按大小排序不然(a,b)和(b,a)会被当成两条边重复处理。std::setstd::pairint,int boundaryEdges; for (auto tri : badTriangles) { for (int i 0; i 3; i) { int a tri.v[i]; int b tri.v[(i1) % 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::vectorint 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::vectorPoint pts loadPoints(); removeDuplicates(pts); // 第一步去重 // 建立超级三角形 Point super[3] makeSuperTriangle(pts); pts.insert(pts.end(), super, super 3); std::vectorTriangle 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::vectorint 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::setstd::pairint,int edges; for (int ti : bad) { markDead(triangles[ti]); for (int i 0; i 3; i) { processEdge(triangles[ti].v[i], triangles[ti].v[(i1)%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 - hh是凸包顶点数所以你可以在开始时就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结合起来做最近邻图、影响区域分析。每个方向都能拆出很多实用技巧等以后有机会再单独写。本文还有配套的精品资源点击获取