简介Delaunay三角剖分算法的C实现包面向需要进行点集三角剖分、有限元分析或几何图形处理的开发者利用C类封装了核心算法流程可直接嵌入数值分析与计算机图形学项目中。资源包共2个文件分别以.h和.cpp构成头文件负责类的接口声明源文件实现剖分细节整体仅3KB代码精简便于阅读和二次修改。已有3662人浏览学习适合希望快速掌握Delaunay构网原理并落地为代码的中高级C用户。拿到代码后可参考其点集存储结构、逐点插入与空外接圆判断等关键步骤结合Voronoi图、EMST树等衍生几何图的构建需求进行扩展能有效减少从算法描述到工程实现的转换成本。 去年做地形可视化项目时甲方给了一批散乱的高程点要求连出来的网格在地球上不能有褶皱。我第一次用最朴素的方式——每三个靠近的点连一个三角形——结果整个地形就像被揉皱的纸光照一打全是细碎的棱。换成Delaunay三角剖分之后同一批数据网格质量肉眼可见地变好狭长三角形基本绝迹。这篇文章就把Delaunay三角剖分在C里的原理、算法选型和完整实现拆一遍给正在趟这个坑的同学一份能直接抄作业的版本。1. 从一套难看的网格说起Delaunay三角剖分到底在解决什么问题1.1 随手连线为什么不行给定平面上的一组散点可以连出无数种三角形网格只要三角形互不重叠、完全覆盖所有点构成的凸包都算合法三角剖分。但合法不意味着好用。如果你随机选一种连接方式结果往往是一堆又细又长的三角形混在里面。细长三角形为什么可恶在渲染里它对光照法线计算极不友好网格表面的明暗过渡会突然撕裂在有限元分析里它会造成单元刚度矩阵病态求解器误差被放大在路径规划里它会导致寻路结果出现不自然的折线。换句话说三角形可以连接得很数学正确但工程上完全不可用。Delaunay三角剖分就是专门来解决这个问题的在无数种合法剖分里挑出一种对狭长三角形天然免疫的方案。1.2 空外接圆与最大化最小角两条等价的定义Delaunay三角剖分有两种最常用的定义方式。第一种是空外接圆属性剖分中任意一个三角形其外接圆内部不能包含点集中的任何其他点。想象一下每个三角形的外接圆画出来里面必须是干净的不允许有其他点偷偷混进来。第二种是最大化最小角属性在给定点集的所有三角剖分中Delaunay三角剖分会使得所有三角形的最小内角尽量大。换句话说它把最钝的三角形也尽量往等边三角形方向拉。这两条定义看起来完全不同但数学上被证明是等价的。空外接圆属性是算法操作时的核心判据最大化最小角属性则直观解释了为什么Delaunay剖分能避免狭长三角形。实际项目中我基本都用第一条来写代码第二条用来给同事和客户解释为什么这么连。1.3 它凭什么能避开细长三角形细长三角形本质上来自过于接近的点和过于夸张的角度差。Delaunay的空外接圆规则像一个强制约束每当某个外接圆内出现新点就必须重新划分让所有圆保持干净。这种局部调整反复发生最终使三角形趋向饱满。这套性质还带来了另一个工程上很值钱的副产品Delaunay三角剖分与Voronoi图互为对偶。Voronoi图把平面按最近邻关系划分成一个个区域Delaunay三角剖分则是把这些相邻区域中心连起来。在路径规划、最近邻搜索这类场景里一套剖分结果能同时服务两种需求。这也是我后来在多个项目里反复依赖它的原因。2. Bowyer-Watson、分治法和逐点插入C项目里的选型思路2.1 三类主流算法各自的路数实现Delaunay三角剖分有三条主流技术路线我先把它们的核心思路放在一起对比。算法平均时间复杂度核心思路实现难度逐点插入增量法O(n log n) 期望每次插入一个点局部重建受影响区域低分治法O(n log n)递归分割点集自底向上合并剖分高扫描线法O(n log n)沿扫描线动态维护已剖分区域高分治法在理论上和实现上都算正规军但合并两个子剖分时需要处理复杂的跨子区域边代码量很大容易在边界条件上翻车。扫描线法更偏向计算几何学术社区的口味生产项目里见到的频率不高。逐点插入尤其Bowyer-Watson版本则因为实现直观、判据统一成了绝大部分开源库和教材的主选。2.2 为什么我最终选择了Bowyer-Watson我最终选定Bowyer-Watson算法主要看中四点。第一核心判据只有一条InCircle测试也就是判断点是否落在三角形外接圆内。算法所有分支都围绕这一个几何测试展开调试时心智负担小很多。第二天然支持增量插入。地形项目里经常要不断追加新的高程点Bowyer-Watson可以直接在现有剖分上插入不用整体重算。分治法遇到这种需求基本等于推倒重来。第三扩展性可控。它从二维升级到三维四面体剖分时思路完全一致把外接圆换成外接球把空腔边界边换成边界三角面即可。后面想往三维走不需要重新学一套方法论。第四配合空间索引很容易加速。虽然朴素的Bowyer-Watson是O(n²)复杂度但插入点只需要检查其周边的三角形用网格索引或DAG有向无环图就能把平均复杂度压到O(n log n)。这个优化路径非常清晰。当然它也有弱点对超三角形大小、浮点精度和顶点输入顺序都比较敏感。这些坑我在第5章会专门展开。3. Bowyer-Watson原理拆解超三角形、空腔与InCircle测试3.1 整体流程一眼看懂Bowyer-Watson的完整流程可以浓缩成四个步骤。构造一个足够大的超三角形把所有点包进去。逐个插入点找到所有外接圆包含该点的三角形这些被称作坏三角形。删除坏三角形它们留下的空白区域叫空腔。把插入点与空腔边界上的每条边连起来形成新三角形。流程听起来简单但每一步都有值得深挖的细节。尤其第2步和第4步是整算法的灵魂所在。3.2 超三角形为什么不能省第一步的超三角形看似只是初始化其实决定了整个剖分能否在边界处正确收尾。它的作用体现在三点让所有输入点都落在某个初始三角形内部保证插入流程一开始就能找到坏三角形三角剖分的最终边界是输入点集的凸包超三角形在算法结束时会被切掉只剩包含真实点的区域如果没有超三角形外侧点找不到参照物剖分结果会在边界处漏掉三角形。超三角形大小的选择很考验经验。太大InCircle测试的坐标差动辄数百倍浮点误差会被放大太小万一有输入点落在超三角形外算法直接出错。我的做法是先算点集包围盒再以包围盒最大边长的4倍作为超三角形的基准尺寸把超三角形的三个顶点放在包围盒中心对称的位置上。这样既保证覆盖又不会把数值范围撑得太大。3.3 InCircle测试三角剖分唯一的裁判InCircle测试是整个算法里唯一一个几何判定函数性能和准确性全押在它身上。它判断点p是否落在由a、b、c三个顶点构成的三角形外接圆内部。最直接的想法是算出外接圆圆心和半径再量p到圆心的距离。但这样要解方程组、开平方根慢且容易丢精度。标准做法是直接用行列式判定| ax-px ay-py (ax-px)²(ay-py)² | | bx-px by-py (bx-px)²(by-py)² | 0 | cx-px cy-py (cx-px)²(cy-py)² |这个行列式展开后不需要开方只涉及加减乘运算速度更快。但要注意符号规则行列式结果大于0与否依赖三角形顶点a、b、c的绕序。实际实现里我会先判断三角形方向如果三个顶点是顺时针就交换b、c换成逆时针再做行列式。这个细节很多人会忽略导致结果在部分点集上神秘地全错。3.4 空腔重建从坏三角形到新边的关键一环找到所有坏三角形后需要把它们的12条边全部提取出来然后执行一个非常漂亮的边去重操作把每条边看成无向边如果一条边同时出现在两个坏三角形里说明它是内部共享边直接丢弃如果只出现一次说明它是空腔边界边保留。保留下来的边界边按顺序围成一个多边形空腔。由于空腔内的任何点都已经被净化而新插入点位于空腔内部所以把新点和每条边界边的两个端点分别组成三角形就能无缝填补空腔。这个拆旧房改新房的过程是Bowyer-Watson最精髓的地方坏三角形集合决定了受影响区域边界去重确定了新三角形的骨架最后一步只是简单地连接。理解了空腔重建再看任何Bowyer-Watson实现都会很通透。4. 手写C实现数据结构设计、核心函数与完整可编译Demo4.1 用索引三角形还是点数组一个容易被忽略的决定写C实现时第一个重要决定就是数据结构。很多初学者习惯把三角形直接存成三个Point对象但这会在三个地方吃亏内存里反复拷贝坐标超三角形删除时比较麻烦后续做邻居查询、边统计时完全无从下手。我采用索引三角形所有点存在一个vector 里三角形只存三个int顶点索引。这样三角形对象体积小遍历速度快删除时也只需比较索引元组。另一个要预判的需求是边界提取以索引保存边用无向边比较就能轻松实现内部边去重。4.2 核心函数排序、去重、InCircle与边界提取实现里几个核心函数各司其职。排序去重放在最前面是为了保证输入点集干净InCircle前面讲过是几何判据核心边界提取则依赖无向边的去重逻辑。我在排序上直接用了std::sort加lambda按x再按y排序。这个顺序能带来两个好处一是顶点在内存中的空间局部性更好二是能配合std::unique快速去重。去重用的是精确相等的operator不引入EPS因为浮点近似重复点应该在程序入口做专门的输入清洗而不是在算法内部偷偷抹平。4.3 一份可直接跑的完整代码下面的实现是完整可编译的C11及以上标准都能跑通。我把类封装成DelaunayTriangulation输入一堆点输出就存在triangles里。#include vector #include algorithm #include cmath struct Point { double x, y; Point(double x_ 0.0, double y_ 0.0) : x(x_), y(y_) {} bool operator(const Point p) const { return x p.x y p.y; } }; struct Edge { int a, b; }; bool sameEdge(const Edge e1, const Edge e2) { return (e1.a e2.a e1.b e2.b) || (e1.a e2.b e1.b e2.a); } struct Triangle { int v[3]; }; double cross(const Point o, const Point a, const Point b) { return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x); } bool inCircle(Point a, Point b, Point c, const Point p) { if (cross(a, b, c) 0.0) std::swap(b, c); 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; return (ax * ax ay * ay) * (bx * cy - cx * by) - (bx * bx by * by) * (ax * cy - cx * ay) (cx * cx cy * cy) * (ax * by - bx * ay) 0.0; } class DelaunayTriangulation { public: std::vectorPoint points; std::vectorTriangle triangles; void triangulate(const std::vectorPoint input) { points input; std::sort(points.begin(), points.end(), [](const Point p, const Point q) { return p.x q.x || (p.x q.x p.y q.y); }); points.erase(std::unique(points.begin(), points.end()), points.end()); triangles.clear(); int n static_castint(points.size()); if (n 3) return; double minX points[0].x, maxX points[0].x; double minY points[0].y, maxY points[0].y; for (const Point p : points) { minX std::min(minX, p.x); maxX std::max(maxX, p.x); minY std::min(minY, p.y); maxY std::max(maxY, p.y); } double dx maxX - minX, dy maxY - minY; double d std::max(dx, dy) * 4.0; double midX (minX maxX) * 0.5; double midY (minY maxY) * 0.5; Point s1(midX - d, midY - d); Point s2(midX d, midY - d); Point s3(midX, midY d); int superIdx n; points.push_back(s1); points.push_back(s2); points.push_back(s3); triangles.push_back({superIdx, superIdx 1, superIdx 2}); for (int i 0; i n; i) { insertPoint(i); } triangles.erase( std::remove_if(triangles.begin(), triangles.end(), [](const Triangle t) { return t.v[0] superIdx || t.v[1] superIdx || t.v[2] superIdx; }), triangles.end()); points.erase(points.begin() n, points.end()); } private: void insertPoint(int p) { std::vectorTriangle bad; std::vectorEdge edges; for (const Triangle t : triangles) { if (inCircle(points[t.v[0]], points[t.v[1]], points[t.v[2]], points[p])) { bad.push_back(t); edges.push_back({t.v[0], t.v[1]}); edges.push_back({t.v[1], t.v[2]}); edges.push_back({t.v[2], t.v[0]}); } } std::vectorEdge boundary; for (size_t i 0; i edges.size(); i) { bool duplicated false; for (size_t j 0; j edges.size(); j) { if (i ! j sameEdge(edges[i], edges[j])) { duplicated true; break; } } if (!duplicated) boundary.push_back(edges[i]); } triangles.erase( std::remove_if(triangles.begin(), triangles.end(), [](const Triangle t) { for (const Triangle b : bad) { if (t.v[0] b.v[0] t.v[1] b.v[1] t.v[2] b.v[2]) return true; } return false; }), triangles.end()); for (const Edge e : boundary) { triangles.push_back({e.a, e.b, p}); } } };这段代码在常规随机点集上能稳定生成正确的Delaunay三角剖分。它的目标是教学清晰不是极限性能。4.4 几个实现细节的取舍说明有几个细节我想单独拎出来说。第一边界提取用无向边重做的朴素O(m²)比较。m是坏三角形的边数量通常很小不需要在这里过度优化。真正的大性能瓶颈是逐点插入时扫描全部三角形那部分才是优化重点。第二超三角形顶点顺序。我选的三个点构成一个直角三角形覆盖范围足够且形状规则。实际使用中如果点集坐标特别极端可以把4.0这个系数再调大一些但注意别调得太大否则浮点误差会升高。第三删除三角形时用erase加remove_if。这是C标准做法比手动遍历erase更安全不会因为迭代器失效问题踩坑。5. 验证与踩坑空圆遍历、欧拉公式、共圆和浮点精度5.1 三个不依赖可视化就验证正确性的方法写完算法第一步不是跑图形界面而是用数学方法验证结果的合法性。我常用的验证手段是这三板斧。第一板斧是欧拉公式。对包含n个点、凸包边界上有h个点的三角剖分三角形数量必须等于2n - 2 - h。h可以从剖分结果里统计所有只出现一次的边就是凸包边界边数一数边数即得h。如果三角形数量对不上说明算法在某处丢了三角形或者产生了重叠。第二板斧是空圆遍历。对每个三角形遍历所有输入点检查是否有非顶点落在其外接圆内。这个做一次是O(n²)但作为验证可以接受。注意判断时要用严格在圆内而排除圆上的点否则四点共圆的合法情况会被误判。第三板斧是面积守恒。所有三角形面积之和应该等于输入点集凸包的面积。这个方法能捕捉到三角形重叠的问题而重叠在只看顶点时很容易漏掉。5.2 最常见的看起来对但实际错的几种情况我在调试这个算法时遇到过的坑归纳起来大致有四类。重复点是最隐蔽的杀手。输入里出现两个几乎重合的点时InCircle测试会退化三角形面积趋近于零后续计算全部失真。所以算法入口必须排序后去重。共线点也容易出问题。当三个点共线时外接圆半径无穷大InCircle行列式为0判据符号模糊。实际生产里需要在输入清洗时剔除共线的冗余点或者确保它们不会凑成一个三角形。四点共圆是Delaunay理论中的经典退化。空外接圆判据在点在圆上时没有明确偏向任何选择都合法但不同选择会产生不同的剖分。工程上一般用EPS将圆上点视为不在圆内让剖分结果确定下来。超三角形太小会直接导致算法失败部分输入点落在超三角形外面这些点在第一步找不到包含自己的三角形剖分结果的外边界就残缺了。5.3 浮点误差控制EPS阈值怎么给浮点误差是几何算法绕不开的话题。InCircle测试用的是行列式当点坐标很大或者三角形很扁时行列式的中间量可能达到非常大的量级精度损失随之而来。我的经验是用double类型然后给所有判断加一个相对阈值而不是固定EPS。比如在InCircle里用det 1e-9作为进入坏三角形的条件但这里的1e-9其实应该根据坐标量级缩放比如det tol * maxCoord²。写死绝对阈值在坐标范围跨越多个数量级时很容易失效。另一个实用技巧是尽量让输入点坐标归一化到相近的量级再送入算法。比如地形数据动辄几十万米的坐标可以先平移到原点附近、缩放到[0,1]区间算完再映射回去。这样能显著减小数值误差。6. 从Demo到工程性能优化、输入清洗与扩展方向6.1 空间索引和顶点排序带来的性能变化上面的实现是O(n²)级别的处理几千个点没问题到十万个点就开始吃力。关键瓶颈在insertPoint里扫描全部三角形做InCircle测试。一个简单有效的优化是网格空间索引把平面划分成均匀网格计算每个三角形外接圆的外接矩形把它挂到覆盖到的网格单元里。插入新点时只搜索新点周围若干网格单元里的三角形而不是全部三角形。这样平均复杂度能降到O(n log n)甚至更低。另一个性价比很高的优化是顶点排序。把点按x-y字典序排列或者更进一步按莫顿码排列能让逐点插入时新点与周边点在内存上更接近缓存命中率大幅提升。实测中同样十万个点莫顿码排序配合网格索引整体耗时能从十几秒降到一秒以内。6.2 重复点、共线点、超大坐标输入清洗的最佳实践我现在的习惯是任何点集进入算法之前先过一遍清洗管道。清洗管道包括三步。第一步把坐标平移到中心缩放单位避免超大坐标量级带来的精度问题。第二步对点按x-y排序用精确相等去重如果业务上允许可以做半径合并将彼此距离小于极小阈值的一簇点合并成一个。第三步检测共线三元组必要时剔除中间点或做微小的位置抖动。这套清洗逻辑看起来繁琐但能省掉后面积累下来的大量调试时间。很多时候Delaunay剖分出错不是算法问题而是输入数据本身不干净。6.3 二维到三维、约束剖分后续可以怎么扩展Bowyer-Watson从二维升到三维时思维框架完全复用InCircle变成InSphere测试空腔边界边变成空腔边界三角形面坏三角形变成坏四面体。核心流程还是找出受影响的四面体、提取边界面、连接新点生成新四面体。还有一个方向是约束Delaunay三角剖分CDT。地形建模里经常要保留道路、河流这类约束边普通的Delaunay不保证这些边出现在结果中。CDT的做法是先做无约束剖分再强行加入约束边并局部翻转相邻三角形复杂度会逐步展开。如果你和我一样需要处理几十万甚至百万级点我建议不要重复造轮子直接考虑集成CGAL或Triangle这类成熟库。手写实现的价值在于理解原理和定制化改造真到生产规模稳定性和性能还是专业库更可靠。最后说一点个人经验。这个算法的坑不在算法本身而在数据预处理和数值控制。我踩过最惨的一次是客户给的GPS轨迹点里藏了一堆重复坐标结果生成的剖分里出现了大量面积为0的三角形渲染直接花屏。从那以后我写任何几何算法都会把输入清洗放在最前面并坚持用欧拉公式和空圆遍历做回归验证。Delaunay是一个原理简洁但细节沉重的算法把细节处理好它就能成为工程里非常趁手的工具。本文还有配套的精品资源点击获取