C++手写Delaunay三角网：Bowyer-Watson算法详解与性能优化

简介一份基于C实现的Delaunay三角网算法工程包面向计算几何初学者、GIS与有限元网格生成相关开发者目标是以完整工程示例展示Delaunay三角剖分从数学定义到代码落地的全过程。Delaunay三角网的核心特性是任一三角形外接圆内不含其他点能有效避免狭长三角形、生成均匀网格整个工程正是围绕这一特性进行数据结构与算法设计。压缩包内共42个文件以.h头文件与.cpp源文件为主体头文件用于类与接口声明源文件实现剖分逻辑同时包含Visual C工程配置、资源文件、中间编译产物以及可直接运行的DEM.exe演示程序各文件组合起来构成一个可重新编译的完整开发环境整体大小约2.66MB。已有1881人学习下载通过阅读源码可深入理解点/三角形数据结构、逐点插入或扫线构建流程、空圆检测与局部优化策略分析时可重点关注三角形邻接关系维护以及新点插入后如何通过边翻转恢复Delaunay性质。自带的MFC工程界面能直观展示三角剖分结果便于断点调试与效果验证。对于希望在C项目中实现或改造Delaunay算法的开发者这是一个信息密度较高、兼具学习与复用价值的参考资料。 Delaunay三角网算法是计算几何里绕不开的经典问题做地形建模、有限元网格生成、机器人路径规划甚至游戏里的导航网格都会碰到它。用C来实现这个算法是我觉得性价比最高的一条路线——一方面能拿到接近原生的性能另一方面实现过程中你会把三角形构造、空间关系判断这些基础功彻底打通。这篇内容不是纯理论科普更像是一个从零手写Delaunay三角网的完整记录适合刚入门C想找个硬核练手项目的同学也适合已经在项目里调过CGAL、Triangle这类库、想搞明白底层到底怎么转的开发。为了方便复现代码风格我会尽量保持简单直接主力用C11/14特性配合VS Code里配好的C/C环境就能直接编译跑通。1. 先搞清楚Delaunay三角网到底在解决什么问题1.1 从离散点到三角形网格一个经典需求我一直觉得理解一个算法的最好方式就是先弄明白它究竟在解决什么问题。假设你手头有一批平面上的离散点比如测绘回来的高程点、传感器扫描出来的点云投影第一步往往不是直接画等高线而是把这些散点连接成一张三角形网格。网格质量的好坏直接决定了后续插值、渲染或者有限元分析的精度。这里说的“好”不是随便连三个点拉倒。同一个点集存在很多种三角剖分方式有些剖分会产生特别狭长的三角形这在数值计算里非常不友好内插误差大、中间结果也不自然。Delaunay三角网就是从这个需求里长出来的它保证在所有可能的三角剖分里让三角形整体尽量饱满匀称同时具有空外接圆、最大化最小角这些漂亮的数学性质。对只想要一个稳定可靠的网格结构的人来说它就是默认答案。1.2 空外接圆规则算法的灵魂所在Delaunay三角剖分的定义非常朴素对点集P中的任意一个三角形如果它的外接圆内部不包含P中的任何其他点那么这组剖分就叫Delaunay三角剖分。换句话说每个三角形的外接圆都是“空”的。这个规则看起来简单却同时保证了全局最优的三角形质量。算法在执行过程中所有判定最后都会落回到“点是否在某个三角形的外接圆内”这一个问题上。所以写代码时判圆函数就是核心中的核心它的数值精度和计算效率会直接影响最终网格质量和整个程序的运行速度。后面我会专门说这个函数怎么实现、容差怎么取。2. 算法路线选型为什么最终选Bowyer-Watson2.1 三种主流构建方式横向对比实现Delaunay三角网业内常见的路线有三条逐点插入法Bowyer-Watson、分治法Divide and Conquer和翻边法Flip Algorithm。我一开始也纠结过用哪种实际研究下来结论很清晰。分治法的时间复杂度最优能达到O(n log n)但它对递归划分和跨区域合并的要求比较苛刻处理大量数据时内存布局也麻烦调试起来特别费劲。翻边法从任意一个合法三角剖分出发反复检查对角线是否需要翻转逻辑最简单但复杂度只有O(n²)点一多就难受。Bowyer-Watson属于增量插入的思路平均复杂度O(n log n)最坏O(n²)胜在代码直观、容易调试而且天然支持动态加点——如果以后想实现增量更新这套代码可以直接扩展。对于大多数项目的前期原型验证Bowyer-Watson几乎是最优解。2.2 Bowyer-Watson增量插入的完整流程Bowyer-Watson的核心思路可以拆成四步。第一步构造一个足够大的超三角形把所有点都包进去。第二步把点集中的点逐个插入当前三角网。第三步每插入一个点找出所有外接圆包含新点的三角形这些三角形构成一个“影响区域”删掉它们。第四步把影响区域的边界与新点相连生成新的三角形然后继续处理下一个点。这里最需要注意的是超三角形的大小和位置。它必须严格覆盖所有待插入点否则点落在三角网外面会导致边界构造错误。我习惯把超三角形的三个顶点放到点集包围盒对角线十倍距离的位置这样既保证覆盖又不会因为顶点距离过远导致后续外接圆判定时浮点误差被无限放大。别为了省事用两倍对角线实战中真会出事。3. C核心实现数据结构与关键函数3.1 点、边、三角形的数据模型设计设计数据结构时我踩过不少坑。最开始的版本用了完整的边对象每个三角形保存三条边的索引视觉上很好理解但维护成本极高——插入和删除时要去重边、更新边指针代码越写越绕。后来我换成只维护三角形列表的方案每个三角形记录三个顶点索引以及三个邻居三角形的索引。顶点索引直接指向点数组邻居索引用于快速遍历和边界提取。一条边隐式由三角形顶点对的组合表示不单独创建边对象。这样的设计有三个好处内存占用小、相邻关系清晰、删除三角形时只需要把对应邻居索引置为无效。struct Point { double x, y; }; struct Triangle { int v[3]; // 三个顶点索引 int neighbor[3]; // 与三条边相邻的三角形索引-1表示边界 };为什么用索引而不是直接存坐标因为点在插入过程中坐标不会变用索引可以避免结构体复制带来的开销同时最终输出网格时需要的就是顶点坐标数组和三角形索引数组这种布局可以直接喂给OpenGL或者写成OBJ文件非常方便。3.2 超三角形初始化与逐点插入实现初始化时我先扫描一遍所有点找到x、y方向的包围盒范围然后构造超三角形。逐点插入的主循环用伪代码描述如下对每个待插入点p先初始化一个被破坏三角形的集合badTriangles。遍历当前三角网中所有三角形t如果p位于t的外接圆内部把t加入badTriangles。从badTriangles中找出边界多边形即那些只被一个badTriangle拥有的边。删除badTriangles中的所有三角形。用边界多边形的每条边和p构建新三角形加入三角网。实际代码骨架是这样for (size_t i 0; i pts.size(); i) { std::vectorint bad; for (size_t t 0; t trigs.size(); t) { if (inCircle(pts[trigs[t].v[0]], pts[trigs[t].v[1]], pts[trigs[t].v[2]], pts[i])) { bad.push_back(static_castint(t)); } } // 提取边界边、删除坏三角形、重建新三角形 }最考验细节的地方就是“提取边界边”。我在第一版实现里直接用一个集合存所有坏三角形的边然后将出现两次的边剔除剩下就是边界边。这个方式简单且不容易出错用整数索引生成唯一边标识例如(min(v0, v1) 32) | max(v0, v1)哈希碰撞概率极低千万不要用浮点坐标拼接。3.3 非法三角形剔除与局部重建删除badTriangles之后新生成三角形的邻居关系需要重新绑定。这一步是很多教程里讲得最模糊的地方我自己在这里整整调了一下午。简单说每个新三角形由一条边界边和一个插入点构成。新三角形的邻居信息分两部分一条边和现存的非坏三角形相邻另两条边和本次重建的其他新三角形相邻。如果处理不当最终网格会残留空指针或者索引越界导致后续遍历崩溃。我的做法是重建时先把所有边界边存到一个数组并记录每条边在边界多边形中的序号然后按照顺序逐条连点成三角形最后统一遍历一遍三角形数组根据顶点对的匹配关系重新挂接邻居。关键点在于“分两步走”先集中建三角形再统一刷邻居不要边建边刷否则很容易出现中间状态不一致。4. 边界完整性与退化情况处理4.1 凸包边界为什么容易丢边我最初跑通Bowyer-Watson后迫不及待地画图结果发现生成的网格外边界不是点集的凸包而是变成了一个巨大的三角形。原因很典型超三角形的顶点也被当成真实点参与了三角剖分最后生成的边界包含超三角形的顶点。解决办法有两种。第一种是最后过滤完成所有插入后把所有包含超三角形顶点的三角形删掉。第二种是在插入过程中直接忽略超三角形顶点使用时单独标记它们的索引为负数或极大值。我在实际项目中用“过滤法”因为它更保险删除后剩下的三角形天然构成点集的凸包。如果你要的是带凹边的约束Delaunay三角网那是另一个领域的问题需要额外的约束边恢复逻辑今天先不展开。普通散点插值过滤法完全够用。4.2 四点共圆与浮点精度问题这是另一个让人头秃的坑。在判断点p是否在外接圆内时标准做法是用行列式代码可以写成这样// 返回true表示点p在三角形abc的外接圆内部 bool inCircle(const Point a, const Point b, const Point c, const Point p) { double dx1 a.x - p.x, dy1 a.y - p.y; double dx2 b.x - p.x, dy2 b.y - p.y; double dx3 c.x - p.x, dy3 c.y - p.y; double det (dx1 * dx1 dy1 * dy1) * (dx2 * dy3 - dx3 * dy2) - (dx2 * dx2 dy2 * dy2) * (dx1 * dy3 - dx3 * dy1) (dx3 * dx3 dy3 * dy3) * (dx1 * dy2 - dx2 * dy1); return det eps; // eps为极小正数通常取1e-10量级 }这个行列式涉及坐标差相乘再相减如果四个点恰好接近共圆数值误差会被放大。我在做高程点插值时一组数据里正好出现了大量近似共圆的点结果网格边反复翻转最后形成一圈锯齿状。处理办法有两个方向。一个是用精确算术比如boost::multiprecision::cpp_dec_float_50做浮点比较代价是性能下降明显。另一个是给判圆函数加一个极小容差eps把“点在圆上”统一视为“不在圆内”避免因为浮点抖动导致无意义的翻转。实际项目里我优先用容差方案只在极少数精度敏感场景才上精确算术性能和稳定性兼得。5. 性能优化与可视化验证经验5.1 用网格索引加速点的定位朴素实现里每一轮插入都要遍历当前所有三角形去计算外接圆判断整体复杂度偏高点少还行点上万就很难受。我的优化方案是引入均匀网格做索引把整个包围盒切分成若干格子每个三角形根据外接圆圆心所在的格子登记到对应桶里。插入新点时只检查该点所在格子以及周围相邻格子里的三角形。实际测试下来对1万个随机点暴力版本需要十几秒加上网格索引之后能压到一两秒。这个优化不是必须的但如果你和我一样需要处理几十万量级的点云这一处改造就是质的飞跃。另外如果数据规模继续往百万冲可以再考虑用C多线程并行遍历三角形但这时候一定要小心容器写入冲突建议先只读遍历收集结果最后再统一合并更新不要边遍历边改容器。5.2 用OpenCV画图验证结果算法写完最怕的不是不会写而是写完了不知道对不对。我的习惯是立刻可视化而不是只看坐标数据。用C配OpenCV把点画成小圆点、三角形画成线段每插入一个点刷新一张图观察整个生长过程。如果某一帧出现三角形重叠或者边界混乱立刻就能定位到对应的插入步骤。这里有个小经验可视化代码和算法核心逻辑尽量解耦分两个文件维护。我见过不少人为了让结果显示漂亮把绘制代码混进三角剖分逻辑里后面调算法时寸步难行。建议把核心算法封装成纯计算类只暴露输入点集和输出三角形索引外部想接OpenCV、Qt还是写OBJ文件都随便。6. 踩坑记录与常见问题排查表把玩这个算法快一个月我把遇到过的典型问题整理成一个速查表方便你快速定位常见问题可能原因排查与解决办法删除坏三角形后新三角形没生成邻居索引混乱或坏三角形收集不完整打印badTriangles集合核对边界边的数量是否为偶数奇数说明少收集或者多收集了三角形生成结果覆盖超三角形的顶点没过滤超三角形顶点插入结束后删除所有包含超三角形顶点的三角形程序崩溃报索引越界邻居索引没有正确置为-1重建邻居时先全部初始化为-1再逐个绑定三角形出现重叠外接圆判定的容差太小增大eps或者对共圆场景改用精确算术顶点数量大时速度骤降每轮都全局扫描三角形加网格索引或按外接圆圆心做Hash分桶边界不光滑出现锯齿近似共圆点导致反复翻转先检查eps设置再考虑消除退化点关于编译环境我在VS Code里用g编译命令一般是g -O2 -stdc17 main.cpp -o delaunay。记得开O2优化不开优化和开了优化性能差距非常明显。如果你在Windows上用的是Visual Studio的MSVC注意把结构体对齐和浮点模式调到默认值不要因为某些全局优化选项破坏浮点判断的一致性。最后再分享一个我自己很受益的小技巧除了画三角形边还可以把每个三角形的外接圆圆心画出来和原三角形叠加显示。这样能非常直观地验证“空外接圆”是否真的成立——圆心散布合理就说明剖分基本合格。我做这个项目最大的收获不是记住了某个算法流程而是养成了“写完立刻可视化验证”的习惯这东西不论放到C项目还是其他领域都非常值钱。本文还有配套的精品资源点击获取