简介一个用C实现的增量Delaunay三角剖分算法工程面向图形学、计算几何与网格处理方向的开发者和学习者。算法通过“半边”数据结构管理网格拓扑持续优化三角形的最小角度核心逻辑集中在main.cpp中该数据结构源自石溪大学顾贤峰教授的设计。压缩包共八十九个文件大小约四点一七兆字节包含C语言源文件、头文件、Visual Studio工程与解决方案、已编译的可执行文件、调试符号文件、日志及说明文档等既可直接运行也便于断点调试和阅读源码。目前已有1002人学习下载。对希望理解Delaunay三角化原理或需要复用网格数据结构的读者而言这一工程提供了从零搭建到调试运行的完整路径省去配置环境的时间适合作为算法学习与二次开发的起始模板。1. 增量 Delaunay 三角剖分这个 C 项目到底解决什么问题做网格生成、路径规划或者地形渲染的朋友大概率都撞过同一个需求给你一堆散点怎么快速连出一张“质量还不错”的三角网最直观的按坐标排序后贪心连线出来的网格细长三角形扎堆后面做插值或有限元计算时各种数值翻车。这份源码实现的是计算几何里最经典的增量式构建方案——逐点插入 Bowyer-Watson 空外接圆准则用 C 从零手写不依赖任何第三方几何库。它的核心价值在于你喂给它任意顺序的二维点集它能在 O(n log n) 到 O(n²) 的期望时间内输出一张满足 Delaunay 性质的三角网并且每个三角形的外接圆内都不含其他点。适合拿来学习算法原理、做图形学课程项目或者作为网格生成模块嵌入到自己的渲染器 / 仿真工具里。我拆这份源码时最深的感受是Delaunay 三角剖分看起来就是“逐点插入 删三角形 补三角形”三步但真正写起来坑全藏在退化情况处理和数据结构设计里。这篇笔记我会从工程视角把整个项目过一遍讲清楚它的结构、核心代码的每一处关键逻辑、跑通它的方式以及如果你要二次开发该从哪儿下手。2. 项目结构与技术选型为什么用 C数据怎么组织2.1 技术栈画像标准 C零第三方依赖代码全部基于 C11 标准编写只用了 STL 的 vector、unordered_map、queue 等基础容器没有引入 Eigen、CGAL 这类重型几何库。这意味着在任何装了现代编译器的机器上都能直接构建也方便你把剖分核心抽出来塞进别的工程。选 C 而不是 Python原因很直接三角剖分是典型的 CPU 密集 指针追踪密集的任务增量插入过程中要频繁做点定位和外接圆判断Python 的列表和对象开销会成为瓶颈。而 C 配合连续内存的 vector 存储顶点、用索引代替指针维护拓扑关系性能优势在十万级点集上会拉开一个量级。2.2 目录结构头文件、实现、示例、可视化四层分离整个工程的组织方式非常清晰属于“一个头文件对应一个实现文件”的经典布局src/核心算法目录Point.h / Point.cpp二维点定义重载了算术运算和比较运算符Edge.h边的定义用于维护三角形邻接关系Triangle.h三角形结构包含三个顶点索引、三个邻居索引、外接圆缓存DelaunayTriangulation.h / .cpp增量剖分主流程examples/demo 程序生成随机点集并调用剖分接口tests/单元测试覆盖基本剖分、退化情况、随机点集验证tools/可视化辅助工具把结果导出成 OBJ 或 SVG 格式2.3 核心数据结构不存坐标存索引这个设计是工程实现里最值得学习的一点三角形存的是三个顶点在 points 数组里的下标而不是坐标值本身。好处有两个一是内存紧凑一个 int 占 4 字节一个 double 对占 16 字节存索引能省一半以上空间二是拓扑查找时直接拿索引做 key不需要比较浮点坐标。Triangle 结构大概长这样struct Triangle { int v[3]; // 三个顶点索引逆时针排列 int neighbor[3]; // 三个邻居三角形的索引-1 表示没有邻居边界 double circumCenterX, circumCenterY, circumRadius; // 外接圆缓存 bool containsVertex(int vIdx) const { return v[0] vIdx || v[1] vIdx || v[2] vIdx; } };neighbor 数组的下标含义和 v 数组对齐neighbor[i] 是「与顶点 v[i] 相对的边」上的邻居三角形。比如 neighbor[0] 对应的是边 (v[1], v[2]) 那边的三角形。这个约定是整个算法里最容易绕晕的地方后面讲空腔修复时会反复用到。外接圆的圆心和半径在插入点后会被查询多次如果不缓存每判断一次就要算一次三点定圆n 个点插入下来就是 O(n²) 次的重复浮点运算。缓存后每个三角形只需在创建时算一次后续直接查。2.4 构建流程总览从超三角形到逐点插入剖分主流程是一个干净的循环每一轮迭代做四件事取出一个尚未插入的点 p在现有三角网中定位包含 p 的三角形 t收集所有外接圆包含 p 的三角形形成空腔删掉空腔内部三角形把 p 与空腔边界上的每条边连成新三角形这个循环不断重复直到所有点都被插入。最后一步是把一开始构造的“超级三角形”相关的三角形剔除得到最终的 Delaunay 网格。3. 核心模块拆解一个点从进入到网格更新的完整链路3.1 三角形定位随机化增量最耗时的一步每插入一个新点第一步要找到包含它的三角形。最简单粗暴的做法是遍历所有三角形做点在三角形内的判断复杂度 O(n)整体退化为 O(n²)。这个项目采用了一种带记忆的随机化策略记录上一次插入点所在的三角形从它开始往邻居方向走——如果当前三角形不包含目标点就沿着目标点所在的边跳到相邻三角形继续找。int DelaunayTriangulation::locateTriangle(const Point p, int startTri) { int current startTri; int guard 0; while (guard triangles.size()) { const Triangle t triangles[current]; // 依次判断 p 相对于三条边的方位 bool inside true; for (int i 0; i 3; i) { int a t.v[i]; int b t.v[(i 1) % 3]; if (orient2D(points[a], points[b], p) 0) { // p 在边 ab 的右侧说明不在当前三角形内 if (t.neighbor[i] ! -1) { current t.neighbor[i]; inside false; break; } else { // 走到了边界仍没找到属于浮点误差边缘情况回退全量遍历 return bruteForceLocate(p); } } } if (inside) return current; if (guard triangles.size() * 2) return bruteForceLocate(p); } return bruteForceLocate(p); }orient2D 是二维向量的叉积判断返回值大于 0 表示 p 在直线 ab 的左侧小于 0 在右侧。这里的符号约定必须谨慎它和顶点存储的绕向直接相关——这个项目保证所有三角形逆时针排列所以判断逻辑才是“左侧为内”。有一个细节值得注意guard 防死循环计数器。理论上沿邻居跳跃一定能收敛但浮点误差可能导致在边界处来回跳所以加了保险。如果超过一定步数就退化成全量遍历。这个“先快后慢”的策略保证最坏情况不会太难看。3.2 空腔收集外接圆判断的批量执行定位到了包含 p 的三角形后接下来做的是广度优先搜索从这个三角形出发检查它的邻居们如果某个邻居的外接圆包含 p就把它也拉进空腔再继续检查那个邻居的邻居。扩散过程用队列实现。std::vectorint DelaunayTriangulation::collectCavity(const Point p, int startTri) { std::vectorint cavity; std::queueint frontier; std::unordered_setint visited; frontier.push(startTri); visited.insert(startTri); while (!frontier.empty()) { int triIdx frontier.front(); frontier.pop(); cavity.push_back(triIdx); const Triangle t triangles[triIdx]; for (int i 0; i 3; i) { int nb t.neighbor[i]; if (nb ! -1 visited.find(nb) visited.end()) { visited.insert(nb); if (inCircumcircle(triangles[nb], p)) { frontier.push(nb); } } } } return cavity; }这里有个容易写错的地方visited 要在入队时标记不是在出队时标记。如果出队时才标记同一个三角形可能被多个邻居反复加入队列空腔边界会算错。inCircumcircle 判断用的是行列式判别式避免直接算距离再比较带来的平方根误差。行列式判断的几何意义是p 在外接圆内时四个点组成的行列式值为正。代码实现如下bool DelaunayTriangulation::inCircumcircle(int triIdx, const Point p) const { const Triangle t triangles[triIdx]; double ax points[t.v[0]].x, ay points[t.v[0]].y; double bx points[t.v[1]].x, by points[t.v[1]].y; double cx points[t.v[2]].x, cy points[t.v[2]].y; double matrix[4][4] { {ax - p.x, ay - p.y, (ax - p.x) * (ax - p.x) (ay - p.y) * (ay - p.y), 0}, {bx - p.x, by - p.y, (bx - p.x) * (bx - p.x) (by - p.y) * (by - p.y), 0}, {cx - p.x, cy - p.y, (cx - p.x) * (cx - p.x) (cy - p.y) * (cy - p.y), 0}, {0, 0, 0, 1} }; // 计算 3x3 行列式省去最后一列 double det matrix[0][0] * (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1]) - matrix[0][1] * (matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0]) matrix[0][2] * (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]); return det 0; }注意这里用了“平移到 p 为原点”的坐标系把 p 挪到原点后原四点共圆判定变成了三点行列式符号判定数值稳定性更好。如果 det 恰好接近 0说明 p 刚好在外接圆上——这就是退化情况后面避坑章节会单独说。3.3 空腔修复删除旧三角形创建新三角形收集完空腔后最核心的一步是重建找出空腔的所有边界边然后用 p 连接每一条边界边生成新三角形。边界边的定义是空腔内部三角形之间的边是内边要删除空腔三角形与外部三角形共享的边是边界边要保留。找出边界边的常见做法是建一个边到出现次数的映射遍历空腔内所有三角形把它们的每条边计数出现两次的说明是内部共享边出现一次的说明是边界边。std::vectorEdge DelaunayTriangulation::findBoundaryEdges( const std::vectorint cavity, const std::vectorint deletedFlag) { std::unordered_mapEdge, int, EdgeHash edgeCount; for (int triIdx : cavity) { const Triangle t triangles[triIdx]; for (int i 0; i 3; i) { Edge e(t.v[i], t.v[(i 1) % 3]); edgeCount[e]; } } std::vectorEdge boundary; for (const auto kv : edgeCount) { if (kv.second 1) { boundary.push_back(kv.first); } } return boundary; }Edge 结构需要实现 operator 和哈希函数且要保证无向性——(0, 1) 和 (1, 0) 是同一条边。这个项目用的是将两个端点索引做异或再乘大质数的方式生成哈希值。拿到边界边列表后逐条创建新三角形。每个新三角形的顶点是 (p, edge.a, edge.b)注意绕向必须保证逆时针。这里有个常见的翻车点如果边界边的方向是任意从 unordered_map 里出来的直接按 (a, b) 顺序创建三角形可能产生顺时针三角形破坏后续所有邻接关系。解决办法是创建前检查 orient2D(edge.a, edge.b, p)如果是负的就把边反过来。for (const auto e : boundaryEdges) { Triangle newTri; if (orient2D(points[e.a], points[e.b], p) 0) { newTri.v[0] e.a; newTri.v[1] e.b; newTri.v[2] pIdx; } else { newTri.v[0] e.b; newTri.v[1] e.a; newTri.v[2] pIdx; } // 计算外接圆缓存 newTri.computeCircumcircle(points); int newIdx triangles.size(); triangles.push_back(newTri); // 记录新三角形与边界外邻居的邻接关系 linkBoundaryNeighbor(newIdx, e, cavity); }linkBoundaryNeighbor 负责把新三角形和空腔外的旧三角形建立邻居关系这一步必须根据 Edge 的端点精确匹配到外部三角形的某条边。外部三角形看到了 p 这个概念但它存储的邻居下标还是旧的空腔内部三角形这些引用要在最后一步统一清理。3.4 清理超级三角形删除与超三角形顶点相连的一切超级三角形是所有增量算法的“脚手架”在点集外围构造一个足够大的三角形包住所有点让第一个点插入时也有三角形可定位、有空腔可扩散。算法结束后要把和它相关的所有三角形删除。void DelaunayTriangulation::removeSuperTriangle() { std::vectorTriangle remaining; for (const auto t : triangles) { bool connectedToSuper false; for (int i 0; i 3; i) { if (t.v[i] superA || t.v[i] superB || t.v[i] superC) { connectedToSuper super; break; } } if (!connectedToSuper) remaining.push_back(t); } triangles std::move(remaining); // 重建邻居索引 rebuildNeighborIndices(); }删除后所有三角形在 vector 中的位置都变了邻居下标全部失效。rebuildNeighborIndices 的做法是遍历所有剩余三角形对每条边在全局边表中查找对应的共享三角形重新建立 neighbor 关系。这个步骤是 O(n) 的只做一次所以无伤大雅。超级三角形的三个顶点坐标要取得足够远。常见做法是取点集包围盒然后往四周扩展一个很大的系数比如 100 倍对角线长度确保包住所有点。但太大也不行过大的超三角形会让外接圆半径变得极大导致插入点时误删过多三角形。这个“大而不巨”的尺度选择是个调参细节后面避坑章节会展开。4. 增量插入的完整流程与邻接关系维护从示例工程跑起来的每一步4.1 主流程代码组装定位、收集、修复、清理四步剖分接口对外只暴露一个函数接收点集返回三角形索引列表。内部把每一轮迭代的四个阶段串起来std::vectorTriangle DelaunayTriangulation::triangulate( const std::vectorPoint inputPoints) { points inputPoints; triangles.clear(); // 计算包围盒构造超级三角形 computeBoundingBox(); createSuperTriangle(); // 随机化插入顺序降低退化输入导致的性能抖动 std::vectorint order(points.size()); std::iota(order.begin(), order.end(), 0); std::random_shuffle(order.begin(), order.end()); int startTri 0; // 上次定位的三角形索引作为下一次定位的起点 for (int idx : order) { if (idx 3) { // 前三个点是超三角形顶点 const Point p points[idx]; int tri locateTriangle(p, startTri); std::vectorint cavity collectCavity(p, tri); std::vectorint deletedFlag(triangles.size(), 0); for (int triIdx : cavity) deletedFlag[triIdx] 1; std::vectorEdge boundary findBoundaryEdges(cavity, deletedFlag); // 标记删除可能延迟这里采用先收集后统一删除的“打标记”模式 removeMarkedTriangles(deletedFlag); createNewTriangles(p, idx, boundary); startTri triangles.size() - 1; // 新创建的三角形是最后一个 } } removeSuperTriangle(); return triangles; }把“删除”实现为“先打标记、后统一移除”而不是边遍历边删是因为 vector 的 erase 会导致迭代器失效而且在遍历过程中修改容器会让邻居索引错乱。4.2 示例工程生成随机点并输出网格examples 目录下的 demo 程序是理解整个项目怎么运作的最短路径。它做三件事生成随机点、调用 triangulate、把结果导出为 SVG。int main() { std::vectorPoint points; std::mt19937 rng(42); std::uniform_real_distributiondouble dist(0.0, 100.0); for (int i 0; i 10000; i) { points.push_back(Point(dist(rng), dist(rng))); } DelaunayTriangulation delaunay; auto triangles delaunay.triangulate(points); // 导出 SVG方便直接用浏览器查看网格 SVGDumper dumper(output.svg); dumper.dumpPoints(points); dumper.dumpTriangles(points, triangles); return 0; }SVG 导出不需要引入任何图像库就是拼文本字符串对调试几何算法非常方便。你改了算法后用同一组随机种子跑一遍直接浏览器打开 SVG 对比网格形态比看一堆数字判断断快得多。4.3 邻接关系维护的三个关键细节邻接关系是整个算法的生命线我在读代码时总结了三个必须铭记于心的细节第一新三角形的 neighbor 初始化要把三个位置都先设成 -1再逐个填充。遗漏任何一个位置后续定位阶段就可能从 -1 跳出去导致访存越界——这是最典型的崩溃来源。第二空腔边界上的外部旧三角形它的 neighbor 里还指向已经标记删除的空腔内部三角形。必须在所有新三角形创建完成后统一扫描一次所有旧三角形把指向已删除三角形的 neighbor 索引替换成新三角形的索引或 -1。第三当点 p 恰好在某个三角形的边上时collectCavity 可能只收集了一个三角形。这时候边界边有两类是共线的一类是原来的边被 p 分成两段一段保留一段删除另一类是新产生的两条边。如果不做特殊处理可能出现退化三角形。5. 避坑实录Delauany 三角剖分最容易翻车的四个场景5.1 浮点误差导致外接圆判断在边界处抖动现象点集规模超过十万后偶尔会出现“缺了一个三角形”的诡异网格明明输入没有退化情况网格却有破洞。原因inCircumcircle 的 det 值在 p 恰好接近外接圆时趋近于 0如果某个三角形被判定为“包含 p”其邻居却不包含空腔边界断裂重建后产生重叠三角形。解决在 inCircumcircle 的比较中加一个容差 epsilon把 det 0 改为 det epsilon。epsilon 的值取 1e-12 量级比较合适。注意别设太大否则会吞掉真实突出的点导致网格整体变形。5.2 超级三角形取太大导致插入第一个点时误删大量三角形现象插入前几个点时空腔收集阶段把大部分三角形全部删掉然后重建出一堆以超三角形顶点为共享点的极长三角形最后这些长条三角形又互相交叉网格质量惨不忍睹。原因超三角形三个顶点坐标设得过于极端比如超出点集范围几十万倍外接圆半径巨大覆盖了几乎所有插入点。解决超三角形取包围盒对角线长度的 20~50 倍即可。经验值是包住所有点后再往外扩 20 倍既能保证所有点在超三角形内部又不至于让外接圆半径膨胀到影响早期插入的判定。5.3 随机化插入顺序与真实坐标顺序混用现象用户给点集时如果已经按 x 坐标排好序插入过程会频繁出现长条三角形算法退化成 O(n²)一万点要跑十几秒。原因定位三角形环节用的是上次位置作为起点的启发式搜索有序插入会让这种启发式完全失效——每次都要从网格一头走到另一头。解决标准做法是插入前做一次随机洗牌。注意这不会改变最终剖分结果——Delaunay 三角剖分的结果是唯一的对一般位置的输入而言所以随机化插入顺序只影响性能不影响输出。5.4 删除三角形后索引未重建导致访问野指针现象程序在跑完一半点后随机崩溃崩溃位置的堆栈指向 triangle. neighbor 访问。原因removeMarkedTriangles 用了 erase 删除vector 中后续元素的位置前移其他 triangle 里保存的 neighbor 下标没有同步更新访问到了已经挪走或越界的位置。解决不要用 erase 逐个删先把存活的三角形拷到新 vector拷贝完成后统一重建邻居关系。这需要一次全局边表扫描但线性复杂度可以接受。6. 二次开发往这个剖分引擎里加你想加的东西6.1 从无向边到带权边把三角网格升级为半边结构如果你不只是想画网格而是要做网格简化、有限元剖分或者路径规划当前这个三顶点三邻居的 Triangle 结构不够用。常见的改造路径是引入 HalfEdge 结构把每条无向边拆成两条有向半边每个半边记录 origin、twin、next、prev 四个指针三角形的三条边对应三个半边。改造后遍历邻域、找边界、做网格简化都会方便很多。6.2 加约束把指定线段强制嵌入网格很多实际场景要求某些线段必须出现在三角网中——比如地形里的断裂线、建筑轮廓。最实用的办法是“约束 Delaunay 三角剖分”先正常剖分然后对每条约束边找出所有与它相交的三角形边把相交区域内的三角形重排。这个实现比朴素增量复杂但工程上可以直接在现成三角网上做翻转边操作代码量大约三百行以内。6.3 性能调优从 O(n²) 到近似 O(n log n)如果你点集动辄百万量级可以把定位阶段的“跳跃搜索”升级成均匀网格 index——把空间划分成 sqrt(n) × sqrt(n) 的格子每个格子维护一个三角形列表。插入点时先定位到包含它的格子再在格子内定位三角形。这一步能让整体复杂度稳定到接近 O(n log n)代价是多维护一个网格索引大约增加一百行代码。6.4 验证方法每次改完先跑三类测试这是我自己折腾这类算法时总结的黄金习惯。改完代码后固定跑三组数据一组随机均匀点、一组网格点阵会有大量四点共圆、一组包含几千个共线点的极限退化数据。随机均匀点验证网格正确性检查每个三角形外接圆内没有其他点网格点阵验证退化处理共线点验证算法不会死循环或产生零面积三角形。这套验证流程帮我在改代码时躲过了至少五次新引入的 bug。印象最深的一次是在改空腔边界查找逻辑时改了哈希函数的结构结果所有边界边的顺序变了新三角形的邻居全错位网格输出直接变成一团乱麻——就是从那次以后我每次改完邻接相关的代码一定强制走一遍那三组测试再把 SVG 网格导出来肉眼过一遍。希望这份拆解能让你少走几个我走过的弯路把算法跑通的同时也理解它每一步背后的几何直觉。本文还有配套的精品资源点击获取