增量Delaunay三角剖分C++实现：点定位、翻边与避坑指南

简介本资源面向计算几何与图形学方向的 C 开发者提供一套增量 Delaunay 三角剖分算法的完整实现可用于最大化三角网格中所有三角形的最小内角改善网格质量。项目采用顾贤峰教授提出的「半边」数据结构组织拓扑关系核心实现集中在 GaussCurvature_Sample 的 GaussCurvature 模块 main.cpp 中适合学习算法工程化落地与网格处理。压缩包共 89 个文件约 4.17MB包含 13 个 h 头文件、1 个 cpp 源文件、1 个 sln 与 vcxproj 工程文件以及 exe、pdb、obj 等编译产物和 tlog、log 等构建日志另有 gif 演示、md 说明与 xml 配置便于直接编译运行与二次开发。目前已有 1002 人学习下载。读者可从中获取增量插入与边翻转的完整代码脉络、半边结构的组织方式以及网格曲率采样的示例工程为有限元网格生成、曲面重建等场景提供可复用的参考实现。1. 增量 Delaunay 三角剖分从点云到三角网为什么“边插边翻”才是工程正解拿到一堆散乱点想把它变成一张没有瘦长三角形、外接圆里不藏点的三角网这就是 Delaunay 三角剖分要干的事。增量法Incremental是其中最贴近工程直觉的一条路点一个一个插进去每插一个点就把它影响到的三角形局部翻边修正直到整张网重新满足空圆准则。它不像分治法那样需要递归合并、也不像扫描线那样对点集预处理挑剔天然适合流式点云、激光雷达逐帧建图、有限元网格逐步加密这类“点会陆续到齐”的场景。代价是朴素实现最坏能到 O(n²)但配上点定位 walking 和局部翻边实际跑起来对几万到几十万点完全够用。这篇就把 C 实现里那些真正卡人的地方——点定位、翻边判据、退化处理、数值精度——拆开讲透让你能照着复现一版能跑、能查、能扩的增量 Delaunay。2. 增量 Delaunay 的核心机制点定位、翻边与空圆准则2.1 为什么选增量法而不是分治或扫描线先讲选型。分治法把点集递归切成两半各自剖分再合并理论复杂度 O(n log n) 很漂亮但合并两个三角网时要处理跨网边界的可见性、公共切线代码量和调试难度陡增一旦点集带重复点或共线点合并逻辑的边界分支能写到怀疑人生。扫描线法对点按 x 排序后维护一条扫描带适合规则网格但对任意散点集的动态插入不友好。增量法的好处是“局部性”插入点 P 只影响包含 P 的三角形及其邻域翻边只在这片邻域里传播天然支持动态加点。工程上你经常遇到的是“先来 5000 个点建个粗网后面再补 2000 个点加密”增量法直接接着插就行分治法得整批重来。代价是最坏复杂度但用随机化插入顺序 walking 点定位期望复杂度能压到 O(n log n) 量级实测几万点毫秒级。我一般会这么定点集规模在百万以内、需要动态加点、或者点分布不均匀选增量法点集固定且极大千万级、追求理论最优才考虑分治。绝大多数工程场景增量法是性价比最高的。2.2 空圆准则与翻边判据的数学表达Delaunay 的核心不变量任意三角形的外接圆内不含其他点。增量插入时新点 P 落在某个三角形 T 内或边上把 T 拆成三个小三角形此时 P 是 T 三个顶点的“入侵者”需要检查 T 的三条边对面的邻接三角形看邻接三角形的外接圆是否包含 P。如果包含说明这条边不合法翻转它。翻边的判据用行列式表达最稳。给定三角形 ABC 和待测点 P判断 P 是否在 ABC 外接圆内用 4x4 行列式三维用 3x3 的 in-circle 测试| ax-Px ay-Py (ax-Px)²(ay-Py)² | | bx-Px by-Py (bx-Px)²(by-Py)² | | cx-Px cy-Py (cx-Px)²(cy-Py)² |结果 0 表示 P 在圆内假设 ABC 逆时针 0 在圆外 0 共圆。这个行列式比“算圆心半径再比距离”数值稳定得多因为避免了开方和除法。翻边时边 AB 的两个邻接三角形是 ABC 和 ABD若 D 在 ABC 外接圆内则翻转 AB 为 CD得到三角形 ACD 和 BCD然后继续检查新边。注意行列式的符号依赖顶点绕向。统一约定所有三角形逆时针存储否则判据会反。2.3 用超级三角形包住点集边界处理的工程惯例增量法第一个坑是边界。第一个点插进去没有三角形可拆怎么办惯例是造一个“超级三角形”把整个点集包进去三个顶点是虚拟点坐标取得足够远比如点集包围盒的 10 倍外扩。所有真实点先插进这个超级三角形最后再把所有含虚拟顶点的三角形删掉。超级三角形不能太小否则真实点可能落在它外面导致插入失败也不能太大否则浮点精度在超大坐标下会崩。我一般取包围盒中心为基准边长取包围盒对角线的 20 倍虚拟顶点坐标用 double 存。最后删虚拟三角形时遍历所有三角形只要含任一虚拟顶点索引就标记删除剩下的就是真实三角网。这一步的坑在于如果点集本身有共线点超级三角形可能退化。稳妥做法是给虚拟顶点加一点点随机扰动或者干脆用三个不在一条直线上的固定远点。3. 用 C 把增量 Delaunay 跑起来数据结构与插入流程3.1 三角形、边与邻接关系的数据结构设计数据结构直接决定实现难度。最朴素的是存三角形数组每次翻边时线性搜邻接三角形几万点还行几十万点就慢。工程上常用半边Half-Edge结构或三角形 邻接索引。我倾向用“三角形数组 每条边记录邻接三角形索引”的折中方案简单且够快struct Point { double x, y; int id; }; struct Triangle { int v[3]; // 顶点索引逆时针 int adj[3]; // adj[i] 是边 (v[i], v[(i1)%3]) 对面的邻接三角形索引-1 表示无 bool alive; // 是否被删除 }; class Delaunay { std::vectorPoint pts; std::vectorTriangle tris; // 点定位用 walking记录上次访问的三角形 int last_tri 0; public: void insert(int pid); int locate(const Point p); void flip(int t, int edge); bool inCircle(int a, int b, int c, int d); };adj[i]的语义要统一adj[i]是共享边(v[i], v[(i1)%3])的另一个三角形。翻边时边 AB 在三角形 T 里的局部索引是 i对面三角形是T.adj[i]翻转后要同步更新四个三角形的邻接关系这是最容易写错的地方。3.2 点定位walking 策略与退化情形点定位决定插入效率。朴素做法是遍历所有三角形找包含 P 的O(n) 每次总 O(n²)。walking 策略从上次插入的三角形出发朝 P 的方向走如果 P 在当前三角形某条边的外侧就跳到那条边的邻接三角形直到 P 落在三角形内部。int Delaunay::locate(const Point p) { int cur last_tri; while (true) { Triangle t tris[cur]; int next -1; for (int i 0; i 3; i) { int a t.v[i], b t.v[(i1)%3]; // orient(a,b,p) 0 表示 p 在边 ab 右侧需跳过去 if (orient(pts[a], pts[b], p) 0) { next t.adj[i]; break; } } if (next -1) return cur; // p 在三角形内 cur next; } }orient是二维叉积(b.x-a.x)*(p.y-a.y) - (b.y-a.y)*(p.x-a.x) 0 表示 p 在 ab 左侧逆时针三角形内部方向。walking 的坑在于如果 P 恰好落在某条边上orient返回 0循环可能卡住。处理办法是遇到 0 时按“在边上”处理直接返回当前三角形插入逻辑里单独处理边上的点拆成两个三角形而非三个。另一个坑是 walking 可能绕圈尤其在点集极度不均匀时。加一个步数上限比如 4*sqrt(n)超了就退回暴力搜索保证不死循环。3.3 插入与翻边的完整代码路径插入流程定位到三角形 T → 若 P 在 T 内部拆成三个三角形 → 对 T 的三条边依次检查翻边 → 翻边可能递归传播用栈存待检查边。void Delaunay::insert(int pid) { Point p pts[pid]; int t locate(p); // 拆分三角形 t 为三个新增两个 int t1 newTriangle(tris[t].v[0], tris[t].v[1], pid); int t2 newTriangle(tris[t].v[1], tris[t].v[2], pid); int t3 newTriangle(tris[t].v[2], tris[t].v[0], pid); // 继承邻接关系并修正 linkAdj(t1, t2, t3, t); tris[t].alive false; // 待检查边栈 std::stackstd::pairint,int stk; stk.push({t1, 0}); stk.push({t2, 0}); stk.push({t3, 0}); while (!stk.empty()) { auto [ti, ei] stk.top(); stk.pop(); if (!tris[ti].alive) continue; int opp tris[ti].adj[ei]; if (opp -1) continue; int a tris[ti].v[ei], b tris[ti].v[(ei1)%3]; int c tris[ti].v[(ei2)%3]; // 找对面三角形的对顶点 d int d -1; for (int k 0; k 3; k) if (tris[opp].v[k] ! a tris[opp].v[k] ! b) { d tris[opp].v[k]; break; } if (inCircle(a, b, c, d)) { flip(ti, ei); // 翻转后把新边压栈 } } }inCircle用前面说的行列式。flip要更新四个三角形的顶点和邻接写的时候拿张纸画一下索引映射否则邻接关系必乱。翻边后新产生的两条边要继续检查所以flip内部把新边压栈。注意翻边传播可能很长但均摊下来每条边翻转次数有限整体期望 O(n log n)。3.4 编译与最小可跑示例给一个最小 main生成随机点跑一遍int main() { Delaunay d; std::srand(42); for (int i 0; i 1000; i) d.pts.push_back({rand()%10000/100.0, rand()%10000/100.0, i}); d.initSuperTriangle(); for (int i 0; i 1000; i) d.insert(i); d.removeSuperTriangle(); d.exportOBJ(out.obj); return 0; }编译g -O2 -stdc17 delaunay.cpp -o delaunay。跑完导出 OBJ用 MeshLab 或任何三维查看器打开能看到三角网就说明基本对了。检查质量统计最小角Delaunay 的最小角应该比随机三角网大得多如果出现大量瘦长三角形八成是翻边判据符号反了或邻接更新错了。4. 避坑与排查增量 Delaunay 最容易翻车的五个地方4.1 现象插入几个点后程序卡死或死循环原因walking 点定位在退化情形下绕圈或翻边传播时邻接关系错误导致边反复翻转。常见触发是点落在已有边上orient返回 0walking 逻辑没处理。解决walking 加步数上限超限退回暴力遍历orient返回 0 时按“点在边上”处理插入逻辑单独分支。翻边时确保每次翻转后旧边标记为已处理避免同一对三角形反复翻。4.2 现象三角网出现重叠三角形或空洞原因翻边时邻接索引更新错位adj[i]的语义在四个三角形间没对齐。这是血泪经验里最高频的 bug因为索引映射纯靠脑补极易错。解决写一个assert检查每个三角形的adj[i]对面三角形的对应边确实是同一条边顶点集合相同。每次翻边后跑一遍断言错位立刻暴露。我一般会在 debug 版本里全程开这个检查release 关掉。4.3 现象共圆点导致翻边震荡原因四个点共圆时inCircle返回 0翻与不翻都合法但若判据写成 0就翻可能和邻接三角形来回翻。解决inCircle严格用 0才翻 0视为不翻。共圆点在规则网格里很常见这个约定能保证终止。代价是共圆区域的三角化不唯一但工程上无所谓。4.4 现象大坐标下精度崩坏三角网错乱原因点坐标量级到 1e6 以上时inCircle行列式里出现坐标的四次方项double 有效位不够符号判断出错。解决插入前把所有点归一化到 [0,1] 区间记录缩放和平移参数最后导出时再变换回去。或者用自适应精度谓词如 Shewchuk 的 robust predicates但实现复杂归一化是性价比最高的做法。4.5 现象删除超级三角形后边界出现缺口原因超级三角形顶点参与的真实三角形没被正确识别或者删除时把共享真实边的三角形误删。解决删除条件是“三角形含任一虚拟顶点索引”逐个检查三个顶点。删完后检查边界边每条只属于一个存活三角形的边就是边界边正常应形成凸包。如果边界凹陷说明删多了或邻接更新有残留。5. 进阶技巧用局部重建加速动态加点与质量验证增量法跑通后真正拉开差距的是动态场景下的性能。如果点是一个一个流式到达每次插入都从last_tri出发 walking点序随机时很快但点序有空间相关性比如激光雷达扫描线时walking 可能走很远。我的习惯是维护一个简单的空间哈希把三角形按质心塞进网格定位时先查哈希找候选再局部 walking实测在扫描线点序下能快 3 到 5 倍。另一个技巧是批量插入时的局部重建。如果一次要加一批点不必逐个走完整翻边流程可以先定位这批点影响的区域把区域内的三角形删掉对“区域边界 新点”重新做一次小规模增量剖分再缝合回去。这在网格加密场景里很常见能把重复翻边省掉。质量验证别只看“能不能画出三角网”。我一般统计三个指标最小角分布Delaunay 应显著优于随机、边长比最长边/最短边超过 20 的三角形要警惕、空圆违例数随机抽样三角形检查外接圆内是否有其他点正常应为 0。写个脚本跑一遍比肉眼看图靠谱得多。// 空圆违例抽样检查 int violations 0; for (auto t : tris) { if (!t.alive) continue; for (auto p : pts) { if (p.id t.v[0] || p.id t.v[1] || p.id t.v[2]) continue; if (inCircle(t.v[0], t.v[1], t.v[2], p.id)) { violations; break; } } }这个检查是 O(n·m)m 是存活三角形数几万点跑一次几秒作为回归测试足够。如果违例数不为 0先查inCircle符号再查翻边是否漏传播。最后一个习惯每次改完翻边逻辑先拿 4 个点、5 个点、共线 3 点这种最小用例跑断言全开确认邻接关系对再上大规模随机点。增量 Delaunay 的 bug 几乎都藏在边界和退化里小用例能逼出来大点集只会让你看到一团乱麻却不知从哪查。希望帮到你。本文还有配套的精品资源点击获取