用C++手写增量Delaunay三角剖分：从原理到工程实现

简介面向计算几何方向学习者和 C 开发者的 Delaunay 三角剖分增量实现主程序位于 GaussCurvature 示例工程的 main.cpp。项目基于石溪大学顾贤峰教授的半边数据结构以逐点插入方式构建三角网格使所有三角形最小角最大化适用于网格生成、曲面重建与数值模拟。压缩包共 89 个文件包含以 h/cpp 为主的源码文件与 VS 工程文件sln/vcxproj同时带有编译产物exe/pdb/obj、mesh 数据、gif 演示和 README 文档整体约 4.17MB。已有 1001 人学习浏览。读者可以借此深入理解半边网格的邻接关系维护与增量插入流程也可直接编译运行查看三角化效果作为几何算法学习与二次开发的基础参考。1. 用 C 手写增量 Delaunay 三角剖分从暴力 O(n²) 到可用的工程实现如果你处理过散点网格生成、地形 TIN 构建或路网骨架提取一定绕不开 Delaunay 三角剖分——它保证三角形的外接圆内不含其他点让网格质量在几何上最优。市面上的库CGAL、Triangle都很重想魔改局部分支逻辑时你大概率会陷入黑匣子。这份资源是我拆过的一个 C 实现核心是增量算法Bowyer-Watson 思路不依赖任何第三方几何库从数据结构到翻边修复全部手写骨架清晰适合那些需要把剖分逻辑嵌进自己渲染或仿真管线的场景。下文从原理选型开始沿着数据结构、插入流程、边界处理和性能调优一路拆到能直接跑的代码。2. 为什么选增量法三种主流剖分算法的取舍实现 Delaunay 三角剖分业界常见的路线是分治法、扫描线法和增量法。分治法理论复杂度最优O(n log n)但递归分割与合并时处理跨子集的空外接圆判定极其麻烦边界的拓扑缝合容易写崩扫描线法的状态树维护复杂而增量法虽然最坏情况是 O(n²)但实现直观、内存布局可控实际随机点集表现接近 O(n log n)并且天然支持动态加点——如果你的场景需要边插入点边重建局部网格增量法是唯一不用全量重算的选项。2.1 增量法的核心不变量空外接圆与 Lawson 翻边整个算法的正确性依赖一个几何判定给定三角形 T 的三顶点 (a,b,c) 和点 p若 p 落在 T 的外接圆内部则 (a,b,c,p) 四点构成的四边形对角线需要翻转。实现上这个判定表现为一个返回值符号——使用行列式判定四点共圆方向注意浮点误差会在这里咬人后文避坑章节会单独说。// 判断点 d 是否在 a,b,c 的外接圆内 // 返回值 0 表示 d 在圆内 0 在圆外 0 四点共圆 double in_circle(const Point2d a, const Point2d b, const Point2d c, const Point2d d) { // 构造 3x3 行列式叉积展开 double adx a.x - d.x, ady a.y - d.y; double bdx b.x - d.x, bdy b.y - d.y; double cdx c.x - d.x, cdy c.y - d.y; double bdxcdy bdx * cdy, cdxbdy cdx * bdy; double cdxady cdx * ady, adxcdy adx * cdy; double adxbdy adx * bdy, bdxady bdx * ady; double det adx * (bdxcdy - cdxbdy) - ady * (bdx * cdxady - cdx * adxbdy) // 注意展开式要逐项验算 (bdx * cdxady - cdx * adxbdy); return det; }这段代码有三个容易踩的地方。第一坐标相减的顺序必须统一以 d 为原点否则行列式的符号含义反转第二det 的展开式我在这份资源里做了逐项符号验证直接抄网上的版本有可能带着错误符号第三当 det 接近于零时共圆点你需要一个 epsilon 阈值做容差判断工程上一般取 1e-12 到 1e-10 之间的值具体看你的坐标量级。2.2 超级三角形用人为外包络规避边界特判增量插入点之前必须有一个初始三角形能覆盖所有点。最省心的做法是构造一个超级三角形super-triangle把全部点包在其内部。这个三角形会在算法末尾连同与它相连的所有边一起删除。超级三角形的顶点选择有个细节不是随便取三个很大的数就行顶点坐标太大会导致外接圆半径极大in_circle 判定的浮点误差被放大。常见做法是取点集包围盒的对角线长度的 20~40 倍而不是直接填 DBL_MAX。Triangle build_super_triangle(const std::vectorPoint2d pts) { // 计算点集包围盒 double min_x pts[0].x, max_x pts[0].x; double min_y pts[0].y, max_y pts[0].y; for (const auto p : pts) { min_x std::min(min_x, p.x); max_x std::max(max_x, p.x); min_y std::min(min_y, p.y); max_y std::max(max_y, p.y); } double dx max_x - min_x, dy max_y - min_y; double dmax std::max(dx, dy); double mid_x (min_x max_x) * 0.5; double mid_y (min_y max_y) * 0.5; // 放大系数取 20~40 倍顶点坐标保持在远离点集的半平面 double scale dmax * 20.0; // 这里不是 1000 倍太大会爆精度 return Triangle( Point2d{ mid_x - scale, mid_y - scale }, Point2d{ mid_x scale, mid_y - scale }, Point2d{ mid_x, mid_y scale } ); }选择等腰直角三角形作为超级三角形是经过考量的——直角顶点朝上左右两个底角是 45°这样外接圆半径在三个方向上相对均衡不会出现极端扁长的外接圆。放大系数从 20 倍起步如果点集分布不均匀比如很多点集中在某个狭长带再加大到 40 倍。这里最怕的是放大系数太小导致某个真实点落在三角形某条边的外侧初始化就先崩了这点在自测里务必覆盖。3. 核心数据布局半边结构让翻边操作变成指针交换增量 Delaunay 的实现难点不在几何判定而在三角形邻接关系的维护。每次插入一个点你可能要删除 3~8 个三角形再新建 6~10 个三角形中间还要做几次翻边如果邻接关系用「三角形存三个顶点」的朴素结构删除和重建时的邻接查找成本是 O(三角形总数)点越多越慢。这份资源用了半边结构Half-Edge把每个三角形拆成三条有向边邻接查找退化成数组索引访问。3.1 半边结构的内存布局与索引语义半边结构的最基本单元是一条有向边它记录了起点顶点索引、所属三角形、下一条半边和对面的半边。一个物理三角形包含三条半边三条半边的循环构成了三角形的边界。关键的语义是一对相邻三角形共享的那条边在两个三角形中表现为方向相反的两条半边通过 twin 指针互指。struct HalfEdge { int origin; // 起点顶点索引指向顶点数组 int twin; // 对偶半边索引反向边-1 表示边界边 int next; // 同三角形内下一条半边逆时针 int face; // 所属三角形索引 }; struct TriangleFace { int edge[3]; // 三条半边索引按逆时针排列 bool alive; // 删除标记懒删除避免数组重排 };使用半边结构后翻边flip操作不再是「记下四顶点坐标再重建两个三角形」而是直接交换两条半边的 twin 指针和它们的 next 链接时间复杂度 O(1)。这里有一个工程取舍老练的实现会做「懒删除」——被移除的三角形不真正从数组里抹掉而是打上 alivefalse 标记新三角形继续追加到数组尾部等整个剖分完成后再统一压缩数组。这样做的好处是三角形索引在插入过程中保持稳定后续要回溯局部修改历史时非常方便。3.2 点定位从暴力遍历到有向边行走增量插入一个点时首先得找到它落在哪个三角形里。最朴素的做法是遍历全部三角形做点在三角形内的测试总复杂度 O(n²)点少还行点超过几万个就肉眼可见地卡。这份资源实现了有向边行走walk策略从上一个插入点的位置开始沿着三角形边不断跨越到相邻三角形直到找到包含目标点的那个。因为点集是随机插入的平均跨步数很少实测十万个点内效率接近线性。int locate_point(Point2d p, int start_tri, const std::vectorTriangleFace faces, const std::vectorHalfEdge edges, const std::vectorPoint2d verts) { int cur start_tri; int guard 0; // 防死循环 while (guard faces.size() * 4) { // 理论最多跨越面数次 // 对当前三角形检查 p 是否在某条边的外侧 bool found true; for (int i 0; i 3; i) { const HalfEdge e edges[faces[cur].edge[i]]; const Point2d a verts[e.origin]; const Point2d b verts[edges[e.next].origin]; // 判断 p 是否在边 ab 的外侧即位于三角形外部 if (orient2d(a, b, p) 0) { cur edges[e.twin].face; // 跨到邻居 found false; break; } } if (found) return cur; } return -1; // 定位失败返回错误标记 }这段代码的细节在 orient2d 的方向约定上。这里约定顶点按逆时针存储因此如果 p 在边 ab 左侧orient 为正说明 p 在三角形外侧需要翻到这条边的对面邻居。若是你的顶点顺序是顺时针所有正负号都要反过来否则定位永远找不到三角形。另外 guard 变量的边界条件要设为 faces.size() 的倍数而不是固定大数因为三角形总数在插入过程中是增长的固定常数可能提前触发误报。3.3 插入点的完整流程拆除与重建的边界条件点定位找到目标三角形后接下来的流程是将目标三角形从网格中拆除以新点为顶点连接它的三条边的端点生成三个新三角形然后对每条新边递归检查是否需要翻边。这个过程最需要小心的地方在拆除阶段——不只是删一个三角形而是要把所有外接圆包含新点的三角形团块全部找出来并一起删除这个团块叫 Cavity。漏删一个三角形最终网格就会出现重叠边。std::vectorint find_cavity(Point2d p, int start_tri, const std::vectorTriangleFace faces, const std::vectorHalfEdge edges, const std::vectorPoint2d verts) { std::vectorint cavity; std::vectorint stack { start_tri }; std::vectorbool visited(faces.size(), false); while (!stack.empty()) { int t stack.back(); stack.pop_back(); if (visited[t]) continue; visited[t] true; // 检查三角形外接圆是否包含 p if (!in_circle(verts[faces[t].v0], verts[faces[t].v1], verts[faces[t].v2], p)) continue; cavity.push_back(t); // 将三个邻居压栈 for (int i 0; i 3; i) { int twin_edge edges[faces[t].edge[i]].twin; if (twin_edge ! -1) { // 非边界边才压栈 int nf edges[twin_edge].face; if (!visited[nf]) stack.push_back(nf); } } } return cavity; }find_cavity 是整个算法正确性的命门它用 BFS/DFS 从起始三角形向外扩散把所有外接圆包含新点的三角形找出来。注意压栈条件是「外接圆包含 p」不是「三角形边被 p 穿过」——这俩不一样用后者会导致 Cavity 边界不是凸多边形后续重建时产生裂缝。Cavity 的边界应该是简单多边形重建时以新点为中心连接边界上的每条边形成新三角形。4. 翻边修复策略递归传播不变量与退化处理Cavity 拆除重建后生成的初始三角形只保证包含新点不保证整个区域重新满足 Delaunay 空外接圆性质必须通过翻边操作把不变量传播出去。这里有个初学者常犯的错误以为只需检查新建的三角形就够了。实际上新建三角形破坏的是局部空圆性质但受影响的邻居可能离新点有一圈距离所以翻边过程必须递归向外层传播直到整片区域收敛。4.1 翻边的判定条件与邻接关系更新对一条共享边 AB 和它的两个对顶点C 和 D若 D 落在三角形 ABC 的外接圆内则需要把 AB 这条共享边翻转为 CD。翻转后原来的两个三角形变成两个新三角形它们的顶点一样但连接顺序完全不同。这个操作在几何上等价于更换四边形的对角线。void flip_edge(int edge_idx, std::vectorHalfEdge edges, std::vectorTriangleFace faces) { // 取共享边的两个三角形 int e0 edge_idx; int e1 edges[e0].twin; int f0 edges[e0].face; int f1 edges[e1].face; // 获取四个顶点a-b 是共享边c 属于 f0d 属于 f1 int a edges[e0].origin; int b edges[e1].origin; int c /* 从 faces[f0] 中取非 a、b 的顶点 */; int d /* 从 faces[f1] 中取非 a、b 的顶点 */; // 重建两个三角形新的边是 c-d faces[f0] TriangleFace({/* c,d 和两者间的边 */}, true); faces[f1] TriangleFace({/* d,c 和两者间的边 */}, true); // 更新相关半边的 twin 指针 // 关键与原 c、d 相连的外部半边不变只改内部连接 }翻边操作的指针更新顺序极容易出错。我的习惯是先用纸笔画一个四边形加对角线标注每条半边的 origin 和前驱后继然后写代码时对照着逐一更新 twin 和 next而不是凭记忆乱指。这份资源里给了完整的指针交换顺序表——先改 faces 的 edge 数组再改四个端点半边的 twin最后修正 next 链顺序颠倒会出现半边指向错误三角形的现象而且极难排查。4.2 退化情况四点共圆与重复点工程上最烦的输入数据是存在共圆点集——比如正方形网格的四个顶点。它们的外接圆判定 in_circle 返回值接近零由于浮点误差可能一会儿判定为圆内、一会儿判定为外。后果是翻边决策不稳定最终生成的剖分不一定错但每次运行结果可能不同。处理手段分为两层。第一层是预处理在插入前对点集做去重去重阈值一般取坐标最小间距的 1e-9 倍把真正重合的点合并。第二层是容差in_circle 返回值的绝对值小于 epsilon 时统一按「不在圆内」处理这样共圆四点会选择一条固定的对角线连接方向网格拓扑稳定。epsilon 的选取直接影响稳定性量级在 1e-12 左右通常够用但如果你用 float 存坐标epsilon 需求量级要放大到 1e-6这时候更建议把所有坐标转成 double 计算。bool in_circle_eps(const Point2d a, const Point2d b, const Point2d c, const Point2d d, double eps 1e-12) { double det in_circle(a, b, c, d); return det eps; // 注意不是 0而是 eps }这里有个反直觉的细节共圆时 det 可能在正负间抖动如果你用 abs(det) eps 做归零处理再与符号判断组合使用逻辑会复杂化。我最终采用的方案是直接把判定从「det 0」改成「det eps」负方向的抖动自动被吸收逻辑链最短而且测试表明对重复点输入也稳定。4.3 边界边的保留避免删除超级三角形时撕毁网格剖分完成后的清理阶段需要把所有连接到超级三角形顶点的三角形删掉。 naive 做法是遍历全部三角形检查其顶点是否包含超级三角形顶点包含则删除。但这里有个隐藏的坑如果 Cavity 爆炸到网格边界时边界边上的半边 twin 指向超级三角形此时如果直接删掉超级三角形连接的三角形边界边上的某些三角形可能被误伤导致最终网格出现孔洞。std::vectorint remove_super_triangle_connected( int super_v0, int super_v1, int super_v2, std::vectorTriangleFace faces, std::vectorHalfEdge edges) { std::vectorint survivors; for (size_t i 0; i faces.size(); i) { if (!faces[i].alive) continue; if (faces[i].v0 super_v0 || faces[i].v0 super_v1 || faces[i].v0 super_v2 || faces[i].v1 super_v0 || faces[i].v1 super_v1 || faces[i].v1 super_v2 || faces[i].v2 super_v0 || faces[i].v2 super_v1 || faces[i].v2 super_v2) faces[i].alive false; // 标记删除 else survivors.push_back((int)i); } return survivors; }5. 避坑排查增量 Delaunay 最常见的五个翻车现场这章把我在实现过程中实际踩过的五个坑记下来每条都是「现象 → 原因 → 解决」结构排查周期最短的一次花了我三个小时希望你看完能少走弯路。5.1 点定位死循环现象插入到点集中段时locate_point 函数进入死循环guard 数始终增长程序卡死。原因超级三角形设置过小某个真实点落在了超级三角形的边外侧。此时有向边行走到边界时沿着边界边跨不出去又因为点不在任何三角形内部找不到目标循环回绕。解决把超级三角形的放大系数从 20 倍调大到 40 倍同时给 locate_point 增加一个断言——连续跨越次数超过三角形总数的两倍时直接触发重新定位从最后一个三角形开始重新走一遍。5.2 删除超级三角形后出现孔洞现象剖分结束后网格内部出现几个孤立的空心区域三角形缺失。原因清理阶段虽然删除了所有与超级三角形相连的三角形但某些应该恢复的边界边没有被正确标记为「保留」从 Cavity 重建的旧三角形在翻边时被误删。解决我的修法是在清理超级三角形之前先遍历所有与超级三角形相邻的边界边把这些边标记为 no-flip。这样翻边时碰到标记边就直接跳过确保边界拓扑不被破坏。5.3 结果依赖插入顺序现象同样一批点洗牌后剖分结果不同某些对角线方向翻来覆去。原因四点共圆的情况没有做容差处理det 在正负边界抖动导致不同的翻边决策。解决统一换成 in_circle_eps 判定epsilon 取 1e-12。同时把输入点集在插入前做一次去重和随机洗牌——去重是为了避免退化输入洗牌是为了降低最坏情况复杂度。5.4 半边指针错乱导致遍历越界现象程序跑完后遍历三角形某些三角形的边索引指向了尚未初始化的内存读出的邻接三角形完全不对。原因C 中半边数组用 vector 存储push_back 添加新半边导致旧索引失效——不对vector 扩容时底层数组迁移你以为的索引指向的内存已经被重新分配。解决所有半边、三角形的索引不用裸指针而是用整数索引并在访问前做范围断言。或者提前预分配足够大的容量比如点数的 3 倍避免扩容迁移问题。5.5 性能劣化到 O(n²) 以上现象五千个点就要跑几十秒明显不正常。原因点定位的起始三角形始终固定为 0 号三角形每插入一个新点都要从 0 号开始走一遍远路跨步数随点数线性增长。解决记住上一个插入点所在的三角形作为下一次定位的起点。这样相邻插入点之间通常只有几步跨越实测千万级别点集也能在可接受时间内跑完。6. 进阶验证与性能调优如何确认你的网格真的满足 Delaunay 性质剖分写完第一件事不是看渲染效果而是做不变量验证。写一个独立函数遍历所有三角形对检查共享边的对面点是否落在对面三角形的外接圆内——只要有一对违反说明剖分有 bug。我每次改动几何判定或指针交换逻辑后都会强制走一遍这个验证它能在几何错误演变成大面积崩溃之前把问题揪出来。性能调优层面除了本文已经提到的点定位记忆化记住上一个点所在三角形还有三个实用技巧值得尝试。第一输入点集排序按 x 坐标分桶后每个桶内局部插入减少跨越距离你可以在插入前做一次桶排序分成 16 个桶效果显著。第二OpenMP 并行初始剖分把点集分成多个子集各自独立做增量剖分最后再用增量法把子集的凸包边连接融合多核下能拿到接近线性的加速比。第三内存预分配给半边数组预留 3 倍点数容量给三角形数组预留 2 倍点数容量这两个压测值来自大量实验——插入过程中的平均三角形总量会涨到点数的 1.6~2.0 倍预留太少触发扩容预留太多浪费内存取 2~3 倍最均衡。网格质量评估方面我一般算两个指标最小角比剖分三角形的最小角与理论最优 60° 的比值和边长比分布。随机点集的 Delaunay 剖分最小角比应该在 0.3 以上如果你跑出来大量三角形的某个角小于 20°说明算法或者输入数据有问题。当然这个阈值不是绝对的极度不均匀的点集可能出现细长三角形这本质上是最优剖分下的几何必然不代表算法错误。最后说一个我保留的习惯在这份资源的测试套件里我用三类数据集压测——均匀随机点、格点扰动、共圆同心圆环。均匀随机点验证整体正确性格点扰动验证退化输入稳定性同心圆环验证大量共圆点集中出现时网格会不会崩。从那以后我每改一次几何内核都会强制走一遍这三组测试确认不变量验证通过再合入。希望这篇拆解能帮你在 C 里写出属于自己的那份 Delaunay有问题的地方多跑几组数据坑踩过了就不会再摔。本文还有配套的精品资源点击获取