简介本资源面向计算几何、图形学与网格处理方向的学习者与开发者提供一套用 C 实现的增量 Delaunay 三角剖分算法工程。Delaunay 剖分能最大化所有三角形的最小内角广泛用于网格生成、曲面重建与有限元前处理等场景项目采用顾贤峰教授提出的「半边」数据结构组织拓扑关系核心实现集中在 GaussCurvature_Sample 的 GaussCurvature 模块 main.cpp 中便于读者顺着主流程理解逐点插入、空外接圆检测与边翻转等关键步骤。压缩包共 89 个文件约 4.17MB包含 13 个头文件、1 个 cpp 源文件、1 个 sln 与 vcxproj 工程文件以及 exe、pdb、obj 等编译产物和 gif、htm、md 等说明素材工程可直接在 Visual Studio 中打开调试。目前已有 1002 人学习下载适合希望从源码层面掌握增量剖分实现细节、并借助半边结构搭建自有网格库的读者参考。1. 从点云到三角网为什么增量 Delaunay 值得你亲手写一遍手里有一堆散乱的点想把它变成一张干净的三角网这事在三维重建、有限元网格、地形建模里绕不开。很多人第一反应是找现成库但真到要嵌进自己的渲染管线或者做局部动态更新时才发现调库的边界比想象中窄。这份 C 实现的增量 Delaunay 三角剖分核心就是把「逐点插入 局部翻转」这套流程完整摊开让你能改、能调、能嵌。它适合已经会写 C、但没亲手搭过三角剖分的人也适合被库的黑匣子坑过、想搞清楚每一步几何判定到底在干什么的从业者。增量法的好处是逻辑线性、便于调试插一个点只动局部不像分治法那样递归展开后难以追踪中间状态。2. 增量 Delaunay 的几何底座空圆准则与翻转操作2.1 空圆准则到底在约束什么Delaunay 三角剖分的定义只有一句话任意三角形的外接圆内不包含其他点。这条约束保证了三角网的最小角最大化也就是不会出现那种又细又长的「刀片三角形」。在增量法里每插入一个新点先找到包含它的三角形把这个三角形拆成三个然后检查新生成的边是否满足空圆准则。不满足就翻转对角线直到所有受影响的边都合法。这里有个容易混淆的点空圆准则是全局性质但增量插入时只需要做局部检查。原因是插入一个点只会影响它周围一圈三角形的外接圆远处的三角形不受影响。所以算法复杂度从暴力全局检查的 O(n²) 降到期望 O(n log n)。这个「局部性」是增量法能跑起来的关键也是后面所有优化和排错的出发点。判断点是否在三角形外接圆内标准做法是用行列式。给定三角形三个顶点 A、B、C 和待测点 P构造一个 3×3 矩阵行向量分别是 (Ax-Px, Ay-Py, Ax²Ay²-Px²-Py²)B、C 同理。行列式符号决定 P 在圆内还是圆外。这个判定在浮点下会有精度问题后面避坑章节会专门讲。2.2 翻转操作把非法边掰回来翻转是增量法的核心动作。当一条边的两个相邻三角形不满足空圆准则时把这条对角线换成另一条对角线相当于把两个三角形拼成的四边形重新剖分。翻转之后要继续检查新边形成一个递归或栈驱动的过程。用栈来管理待检查边比递归更稳因为递归深度在点分布极端时可能很大。我一般会维护一个边栈每次翻转后把受影响的四条边压进去循环处理直到栈空。下面是一个简化的翻转判定骨架// 判断边 (a,b) 是否需要翻转 // 输入边两端点 a,b以及边两侧的第三点 c,d // 返回true 表示需要翻转 bool needsFlip(const Point a, const Point b, const Point c, const Point d) { // 空圆判定d 是否在三角形 abc 的外接圆内 double det orient3d(a, b, c, d); // orient3d 返回正值表示 d 在圆内需要翻转 return det 0; } // 翻转边 (a,b) 为 (c,d) void flipEdge(Edge e) { // 更新拓扑把两个三角形合并再重新剖分 // 实际实现中需要维护半边结构或邻接表 // 这里只示意几何层面的替换 e.a e.c; e.b e.d; }orient3d是几何谓词计算四个点的定向行列式。返回正值说明 d 落在 abc 的外接圆内边非法。翻转后边的端点从 (a,b) 变成 (c,d)同时要更新所有引用这条边的三角形。参数上唯一需要注意的是浮点容差纯 double 在点接近共圆时会抖动后面会讲怎么用自适应精度或者简单阈值兜住。2.3 点定位新点落在哪个三角形里插入新点之前得先找到它所在的三角形。最笨的办法是遍历所有三角形做点在三角形内测试O(n) 一次总复杂度退化成 O(n²)。常见做法是走三角形邻接关系从某个已知三角形出发朝目标点方向走类似「走迷宫」。每次判断点在哪条边的外侧就跳到相邻三角形直到点落在当前三角形内部。这个走法在点分布均匀时很快但遇到狭长三角形会绕路。更稳的方案是维护一个点定位结构比如跳表或者网格桶。不过对于这份实现走邻接已经够用代码量小、调试直观。走邻接的终止条件是点同时在三条边的内侧用三次 orient2d 判定即可。// 从 start 三角形出发定位包含点 p 的三角形 Triangle* locate(Triangle* start, const Point p) { Triangle* cur start; while (true) { // 检查 p 是否在 cur 的三条边内侧 bool inside true; for (int i 0; i 3; i) { Edge e cur-edge(i); if (orient2d(e.a, e.b, p) 0) { // p 在这条边外侧跳到相邻三角形 cur cur-neighbor(i); inside false; break; } } if (inside) return cur; } }orient2d返回正值表示 p 在边 ab 的左侧。三条边都满足左侧点就在三角形内。注意这里假设三角形顶点是逆时针排列如果顺序反了判定会全错这是新手最容易翻车的地方之一。3. 把算法落成 C 代码数据结构与插入流程3.1 半边结构还是邻接表三角网的数据结构选择直接决定代码复杂度。半边结构half-edge表达力最强每条边拆成两条有向半边能 O(1) 拿到相邻三角形和相邻边但指针多、内存管理烦。邻接表简单些每个三角形存三个顶点索引和三个邻居指针翻转时更新邻居关系稍微绕一点但够用。这份实现用的是三角形数组 邻居索引的方式。三角形存在std::vectorTriangle里邻居用整数索引而不是指针好处是扩容时不会失效坏处是删除三角形后索引会空洞。增量法只增不删所以索引方案很合适。每个三角形结构大致长这样struct Triangle { int v[3]; // 三个顶点在点数组中的索引 int n[3]; // 三条边对应的邻居三角形索引-1 表示无 bool alive; // 是否有效翻转后旧三角形标记为失效 };v[0],v[1],v[2]逆时针排列n[i]是与边(v[i], v[(i1)%3])相对的邻居。翻转时把两个旧三角形标记alivefalse新建两个三角形追加到数组末尾再修正周围邻居的索引。这种「逻辑删除 追加」的模式避免了数组中间删除的开销代价是数组会膨胀点特别多时可以考虑定期压缩。3.2 插入一个点的完整步骤插入流程分四步定位、拆分、翻转、更新。定位用上一节的走邻接。拆分是把包含新点的三角形一分为三新点作为公共顶点。翻转是检查新生成的三条边不合法就翻。更新是修正所有受影响的邻居索引。void insertPoint(const Point p) { int pi addPoint(p); // 点加入全局数组 Triangle* t locate(lastTriangle, p); // 定位 // 拆分 t 为三个三角形 int t0 newTriangle(t-v[0], t-v[1], pi); int t1 newTriangle(t-v[1], t-v[2], pi); int t2 newTriangle(t-v[2], t-v[0], pi); // 修正邻居关系 linkNeighbors(t0, t1, t2, t); // 三条新边入栈检查是否需要翻转 stack.push(Edge(t-v[0], pi)); stack.push(Edge(t-v[1], pi)); stack.push(Edge(t-v[2], pi)); // 循环翻转直到栈空 while (!stack.empty()) { Edge e stack.pop(); if (needsFlip(e)) { flipEdge(e); // 翻转后新边入栈 stack.push(newEdge1); stack.push(newEdge2); } } lastTriangle t0; // 缓存下次定位从这里出发 }addPoint把点追加到全局点数组并返回索引。newTriangle在三角形数组末尾创建新三角形初始邻居设为 -1。linkNeighbors负责把三个新三角形之间以及它们与外部邻居的索引接上这一步最容易出索引错位建议写完先拿四个点手算一遍验证。lastTriangle缓存是为了让下一次定位少走几步点分布均匀时能省不少时间。3.3 翻转的边界处理翻转不是无条件的。如果边在凸包边界上它只有一个相邻三角形没法翻。另外如果翻转后形成的对角线会与已有边重合也不能翻。判断凸包边界很简单看边的两个相邻三角形是否存在有一个是 -1 就跳过。bool canFlip(const Edge e) { // 边必须有两个相邻三角形 if (e.t1 -1 || e.t2 -1) return false; // 检查对角线是否已存在退化情况 if (diagonalExists(e)) return false; return true; }diagonalExists在标准增量法里很少触发但在点共线或重合时会出现。实际写的时候可以加一个断言如果触发说明输入点有重复应该在预处理阶段去重。点去重是插入前必做的一步重复点会让定位和翻转都进入未定义状态。4. 避坑与排查浮点精度、退化与索引错位4.1 现象三角网出现交叉边或空洞原因浮点精度不足导致空圆判定符号错误翻转了不该翻的边或者该翻的没翻。点接近共圆时行列式值接近零double 的舍入误差会让符号随机跳变。解决用自适应精度谓词比如 Shewchuk 的 robust predicates或者简单点给行列式加一个相对容差绝对值小于1e-12 * scale时视为零不翻转。容差不能太大否则该翻的不翻三角网质量下降。我一般先用固定容差跑通再根据点坐标范围调整。4.2 现象程序在插入某个点后卡死原因定位走邻接时陷入循环通常是邻居索引接错了导致三角形之间的邻接关系不是双向一致。A 认为邻居是 BB 却不认为邻居是 A走邻接就会绕圈。解决写一个validateTopology函数遍历所有存活三角形检查每条边的邻居关系是否对称。不对称就打印出来。这个检查放在每次插入后跑一遍虽然慢但能第一时间定位到出错的那次插入。定位到之后重点看linkNeighbors里索引更新的顺序常见错误是先更新了邻居再标记旧三角形失效导致邻居指向了失效三角形。4.3 现象凸包边界上的三角形被错误翻转原因canFlip没有正确判断边界边或者邻居索引 -1 被当成有效索引参与运算。解决在所有涉及邻居索引的地方加边界检查-1 直接跳过。另外凸包边界在增量法里会动态变化新点落在凸包外时需要先找到可见的边界边把凸包扩大。这一步很多简化实现会漏掉导致点插到外面后三角网不完整。处理办法是定位失败时走邻接走到边界还没找到包含三角形沿边界走一圈找可见边把新点和可见边连成新三角形。4.4 现象点集规模上万后内存暴涨原因逻辑删除的三角形没有回收数组持续膨胀。每个点插入平均产生 2 到 3 个新三角形旧三角形标记失效但不释放内存占用是有效三角形的数倍。解决定期压缩三角形数组把存活三角形紧凑排列重建索引映射。压缩频率不用太高每插入几千个点做一次即可。压缩时注意更新所有邻居索引建议压缩前后各跑一次拓扑校验。4.5 现象结果三角网与预期不符但拓扑合法原因点定位的起始三角形选择不当导致走邻接走了远路虽然最终定位正确但中间经过的三角形可能被误判。更隐蔽的情况是 orient2d 的符号约定和顶点排列顺序不一致逆时针和顺时针混用。解决统一约定所有三角形顶点逆时针排列orient2d 返回正值表示点在边左侧。写一个checkOrientation函数遍历所有三角形验证有向面积为正。面积计算用叉积符号不对就交换两个顶点。这个检查应该作为单元测试的一部分每次改完几何谓词都跑一遍。5. 进阶技巧用局部缓存和谓词过滤把插入速度再压一档增量法跑通之后性能瓶颈通常不在翻转而在点定位。走邻接的平均步数和点的空间分布强相关点越随机步数越多。一个实用技巧是缓存上一次插入的三角形下一次定位从它出发。如果点集是空间相关的比如按扫描线顺序输入这个缓存命中率很高定位步数能降到常数级。另一个技巧是谓词过滤。空圆判定涉及行列式计算浮点乘法开销不小。可以先做一个廉价的包围盒测试如果待测点明显在三角形外接圆外比如点坐标超出三角形包围盒加外接圆半径直接跳过精确判定。这个过滤在点稀疏时能省掉大量行列式运算。// 廉价过滤点是否可能在外接圆内 bool mayBeInsideCircumcircle(const Triangle t, const Point p) { // 计算三角形包围盒 double minX std::min({t.ax, t.bx, t.cx}); double maxX std::max({t.ax, t.bx, t.cx}); double minY std::min({t.ay, t.by, t.cy}); double maxY std::max({t.ay, t.by, t.cy}); // 外接圆半径上界最长边的一半 double r 0.5 * std::max({dist(t.a,t.b), dist(t.b,t.c), dist(t.c,t.a)}); // 点超出包围盒扩展 r 的范围一定不在圆内 if (p.x minX - r || p.x maxX r || p.y minY - r || p.y maxY r) { return false; } return true; // 需要精确判定 }dist是两点距离。这个过滤的几何依据是外接圆一定在包围盒扩展最长边一半的范围内超出这个范围的点点不可能在圆内。过滤掉之后再做精确行列式判定整体能省 30% 到 50% 的谓词计算。注意过滤只用于「排除」不能用于「确认」返回 true 时仍然要走精确判定。验证三角网是否正确除了拓扑校验还可以检查最小角。把所有三角形的最小内角统计出来Delaunay 剖分的最小角应该不小于某个下界理论上接近 0 但实际点集下不会太差。如果出现大量小于 1 度的角说明翻转没做干净或者精度有问题。我习惯在跑完一批点后打印最小角和平均角作为质量指标。从那以后我每次写完几何谓词都强制拿四个点手算一遍行列式符号再拿五个点跑一次插入确认翻转方向没反。这个习惯帮我省掉了无数次对着三角网发呆的时间。希望帮到你。本文还有配套的精品资源点击获取