简介这是一份基于C的Delaunay三角网算法完整实现工程面向计算机图形学、地理信息系统、有限元网格生成领域的学习者与开发者重点解决二维点集的最优三角剖分构建问题。代码中体现了空圆特性、逐点插入与局部优化等核心思路适合用于算法理解、课程设计或工程二次开发。资源共42个文件以.h头文件与.cpp源文件为主另含Visual C 6.0工程文件.dsp/.dsw、资源脚本.rc、编译中间文件.obj/.sbr、图标及可执行程序等压缩包约2.66MB目录结构清晰可打开工程直接查看、调试和运行。已有1881人学习下载。通过阅读源码可以理清点、三角形及邻接关系的数据结构设计追踪从初始化、逐点插入到空圆检测的完整流程并结合附带的可执行程序直观验证三角网生成效果是掌握Delaunay算法原理与C实现技巧的实用参考。 写 C 版的 Delaunay 三角网说实话不是个轻松的活儿。我最早接触这个算法是在做点云降面与地形网格生成的时候当时用第三方库确实方便但一遇到需要自定义顶点类型、动态插入约束边或做局部网格加密就发现还是自己维护一套核心算法最顺手。网上关于 Bowyer-Watson 算法的讲解很多但大多数停留在理论层面真正能把 C 实现细节、数值稳定性和性能优化讲透的很少。这篇文章我会从工程落地角度完整拆解如何用 C 从零构建一套可用的 Delaunay 三角网算法覆盖数据结构选型、核心循环流程、外接圆判定的数值陷阱以及实测中容易踩的坑。1. 为什么在众多三角化算法里首选 Bowyer-Watson1.1 Delaunay 三角网到底解决了什么问题先简单交代一下背景给定平面上的一组散点我们希望生成一个三角形网格让这些三角形尽量“饱满”避免出现特别狭长的三角形。Delaunay 三角网的核心准则叫空外接圆准则也就是任何一个三角形的外接圆内部不包含其他顶点。这个准则带来的直接好处是三角形最小内角最大化网格质量整体最优。正因为这个特性它在有限元网格划分、地形建模、路径规划、图像配准等领域都被广泛应用。理论上有好几种构造 Delaunay 三角网的思路包括分治法、扫描线法和逐点插入法。分治法时间复杂度最优能达到 O(n log n)但实现复杂度极高处理退化情况要格外小心扫描线法适合流式数据但代码组织比较绕。Bowyer-Watson 算法属于逐点插入这一类平均复杂度 O(n log n)最坏 O(n^2)但它代码量最少、最容易正确实现也最容易扩展成带约束边的版本。对于绝大多数工程场景100 万个点以内 Bowyer-Watson 的性能完全够用。1.2 为什么我坚持用 C 而不是 Python 或第三方库如果你只是跑一次离线建网格那 Python 的 scipy.spatial.Delaunay 很方便几分钟就能出结果。但在我实际做过的几个项目里Python 暴露的问题很明显一是处理千万级点云时内存占用高得吓人二是无法精细控制网格局部密度三是生产环境里要嵌入到 C 服务或实时管线中跨语言调用本身就增加复杂度。自己用 C 实现一遍可以精确管理内存布局、按需分配、甚至针对特定硬件做 SIMD 优化。另一个重要原因是学习价值Delaunay 三角网是计算几何里少有的“看起来简单、写起来全是细节”的算法完整实现一遍对指针、容器、浮点舍入的理解都会上升一个层次。2. 动手前先定数据结构顶点、边与三角形怎么设计最省心2.1 从零定义三个核心结构体在设计数据结构时我的原则是“先简单、后优化”第一版跑通了再去考虑节省内存。顶点结构最简单两个浮点数加一个索引号就够struct Point { double x, y; int id; // 原始点索引方便回溯结果 };三角形结构需要存三个顶点索引以及三个相邻三角形索引。这里有个工程细节很多入门实现会直接把三个顶点坐标复制进 Triangle 结构体这样后续在 LOPLocal Optimization Procedure局部优化过程交换对角线时更新坐标会非常别扭而且浪费内存。只存索引顶点统一放在外部数组里这样所有操作都基于索引完成代码逻辑清晰很多struct Triangle { int v[3]; // 三个顶点索引按逆时针顺序存储 int neighbor[3]; // 邻接三角形索引neighbor[i] 对应 v[i] 的对边 };这里我把“邻接关系”和“顶点索引”放在同一个结构体里是为了后续做局部优化时能快速找到共享边。neighbor[i] 的定义要刻意记住它存的是与 v[i] 对边相邻的那个三角形。这个约定初看有点绕但写代码时非常顺手。2.2 管理三角形列表vector 与 free list 的取舍Bowyer-Watson 算法在插入点时涉及删除与新增三角形的高频操作。我的做法是使用 std::vector 配合一个 std::vector 的空闲索引栈free list。删除三角形时不真正从 vector 中 erase而是把它的索引压入栈新增三角形时优先复用空闲索引。这样做的好处是避免了频繁内存移动也保证了所有三角形索引在迭代过程中相对稳定便于追踪调试。如果你用可变的 std::list 或反复 erase vector代码会简单一点但实测在 10 万点规模下性能相差了接近 5 倍而且 list 的缓存局部性很差遍历时 CPU 缓存命中率低。用 free list 配合紧凑的 vector 布局是性能和实现复杂度之间的最优平衡。真正要拥抱这些工程细节不能只做“能跑就行的 Demo”。2.3 关于空间索引的提前规划如果不做任何空间索引每次插入新点都要遍历当前所有三角形来查找“受影响三角形”效率是 O(n^2)跑 10 万点会明显卡顿。常见优化手段是网格哈希uniform grid或者四叉树。我在第一版里先不做索引保证逻辑正确后再加。原因很简单空间索引把问题复杂化了如果核心算法本身有 bug你会分不清是索引的错还是三角剖分的错。先把无索引版本调试到正确再引入网格哈希每次只检查落在邻近网格单元里的三角形可以显著缩小搜索范围。3. 核心主循环拆解超级三角形、定位与 LOP 局部优化3.1 超级三角形一个技巧解决边界问题Bowyer-Watson 算法对点集范围之外的三角形处理一直是个麻烦。核心技巧是构建一个足够大的“超级三角形”把所有点包进去使得算法运行过程中不会出现没有任何三角形的空区域。这个超级三角形会在最后被移除。超级三角形要多大为了数值安全顶点坐标范围在 [min, max] 的情况下我会让超级三角形覆盖到三倍范围。因为 Bowyer-Watson 的“空外接圆”检测在某些精度不足的边界点附近可能因为浮点误差出现误判留出足够大的余量可以极大减少这类问题。实际操作时直接取包围盒中心点再向周围扩张三倍半径即可。// 构造超级三角形覆盖所有输入点 Triangle makeSuperTriangle(const std::vectorPoint pts) { double minX pts[0].x, maxX pts[0].x; double minY pts[0].y, maxY pts[0].y; for (auto p : pts) { minX std::min(minX, p.x); maxX std::max(maxX, p.x); minY std::min(minY, p.y); maxY std::max(maxY, p.y); } double dx maxX - minX; double dy maxY - minY; double midX (minX maxX) / 2; double midY (minY maxY) / 2; double r std::sqrt(dx * dx dy * dy) * 3.0; // 三倍半径 Triangle sup; sup.v[0] addPoint(midX - r, midY - r); sup.v[1] addPoint(midX r, midY - r); sup.v[2] addPoint(midX, midY r); return sup; }3.2 逐点插入寻找受影响三角形集合核心流程可以用一个简洁的循环概括对每个输入点找到外接圆包含该点的所有三角形把它们统一删除形成一个“空腔”然后连接新点与空腔边界上的每一条边生成新的三角形。这里有个容易写错的点你在循环中不能一边遍历一边删除三角形。我的做法是先遍历一遍三角形列表将所有外接圆包含当前点的三角形索引存入一个临时集合 badTriangles然后再统一处理。因为删除操作会影响索引如果边遍历边删很可能跳过某些三角形或访问到无效索引。std::vectorint badTriangles; for (int i 0; i (int)triangles.size(); i) { if (inCircumcircle(triangles[i], pts, newPoint)) { badTriangles.push_back(i); } }定位“从哪开始搜索”也很关键。无索引版本从索引 0 扫到末尾复杂度高加入网格哈希后只需要从新点所在网格单元及其周围的单元取候选三角形这里需要额外维护“三角形中心落在哪个单元”之类的元数据。3.3 提取空腔边界用边计数取代集合运算删除受影响三角形后需要找出空腔的边界边。最直观的思路是用两组集合分别记录“出现一次的边”和“出现两次的边”以此确定哪些边是边界。初次实现我也这么干但后来发现这里有个更高效的方法用一个边的计数 map初始化时把每个受影响三角形的三条边的计数加一最后计数为 1 的边就是边界边。实际用 std::map 或 std::unordered_map 来做边计数都能跑通。但要注意边表示的一致性我用的是有序对比如 (minIndex, maxIndex) 来避免 (a,b) 和 (b,a) 被当成两条不同边。struct Edge { int p1, p2; // p1 p2 }; std::mapEdge, int edgeCount; for (int tIdx : badTriangles) { for (int k 0; k 3; k) { int a tri.v[k]; int b tri.v[(k 1) % 3]; edgeCount[makeEdge(a, b)]; } } // 边界边是出现次数为 1 的边 for (auto kv : edgeCount) { if (kv.second 1) { int a kv.first.p1, b kv.first.p2; // 连接新点与 a, b 形成新三角形 } }这种做法比集合差运算要好不需要维护额外的端点集合代码也容易验证正确性。3.4 LOP 局部优化交换对角线保证 Delaunay 性质连接边界边生成新三角形后还需要对这些三角形与其邻接三角形做局部优化LOP。LOP 是这算法里最容易写错但又最关键的部分。基本原理是如果一个新三角形与它的某个邻接三角形组成的凸四边形违反了空外接圆准则就把这对三角形的公共对角线交换。LOP 的实现我用了一个栈来迭代传播因为一次交换可能会引发相邻区域再次不满足 Delaunay 性质。这里我用 BFS 式的队列保证所有受影响的边都被处理void legalizeEdge(Triangle tri, int edgeIdx, std::vectorTriangle triangles, std::vectorint freeList, std::stackint checkStack) { int neighborIdx tri.neighbor[edgeIdx]; if (neighborIdx -1) return; Triangle neighbor triangles[neighborIdx]; // 找到邻接三角形中与 tri 共享边的那个顶点 int oppositePt ...; // 遍历 neighbor.v 找出不在共享边上的那个点 if (inCircumcircle(tri, pts, pts[oppositePt])) { // 交换对角线 flipEdge(tri, neighborIdx, edgeIdx, oppositePt, ...); checkStack.push(tri.id); checkStack.push(neighborIdx); } }这里有几点值得强调交换对角线时要同时更新三角形的顶点索引和三个邻居关系最容易漏掉的是“邻居的邻居”也要重定向。交换完成后无论当前 tri 还是 neighbor都需要重新检查它们的另外几条边是否依然合法。我用栈来存放需要再检查的三角形索引。很多实现省略了这一步或者只检查新产生的三角形而不检查老三角形这样输出的网格会存在个别的非 Delaunay 三角形肉眼难见但计算误差会累积。我自己第一版就吃过“忘了传播检查”的亏导致 80% 的区域网格正确但一些局部区域出现细长的非 Delaunay 三角形排查了很久才定位到是 LOP 迭代不完整造成的。4. 外接圆判定的数值稳定性与退化情况处理4.1 行列式法取代斜率的直觉计算判断点是否在三角形外接圆内最容易想到的方法是用圆的方程解出圆心坐标再去比较距离。这个方法有个隐藏风险当三点接近共线或点距很小时圆心坐标的数值会非常大距离比较时精度不够甚至因为两个超大数相减导致灾难性抵消。实践中我采用的是计算几何里经典的行列式判定式。给定三个点 A(x1,y1)、B(x2,y2)、C(x3,y3) 和待测点 P(x,y)计算以下行列式D | x1 y1 x1^2y1^2 1 | | x2 y2 x2^2y2^2 1 | | x3 y3 x3^2y3^2 1 | | x y x^2y^2 1 |如果三角形按逆时针排列那么 D 0 表示 P 在外接圆内D 0 表示在外侧D 0 表示在圆上。为了减少溢出风险我会在计算前先把所有坐标平移到以 A 为原点的局部坐标系里再算行列式。平移不只是微优化而是能显著降低超大坐标下平方项溢出概率的操作。4.2 共圆与三点共线的边界策略工程数据里“四个点恰好共圆”并不罕见例如整数网格点上的正方形顶点就会导致共圆情况。面对 D 0 的情况处理策略要分场景如果只是普通散点三角化任意决定即可但必须保证程序不崩溃如果做地形网格共圆时保留对角线还是交换对角线会影响网格的形态。我在代码里给 inCircumcircle 加了一个极小阈值 epsilon当 |D| epsilon 时视为“在圆上”按“不在圆内”处理。这个阈值的选择也要因地制宜对于坐标在几千数量级的点云我一般取 1e-10 量级但对于归一化到 [0,1] 区间内的点集阈值可以放宽到 1e-14 左右。不要试图用一个固定的全局阈值应对所有数据范围否则要么误判太多要么容错太差。4.3 重复点和退化输入先清洗再计算另一个很容易被忽视的问题是输入点里可能包含重复点或几乎重合的点。两点距离小于 1e-12 时外接圆计算会严重退化产生 NaN 或无穷大。因此算法主循环开始前我写了一个简单的点清洗函数按坐标排序后去重距离够近的点直接丢弃。这个预处理看起来“不优雅”但对稳定性至关重要。std::vectorPoint deduplicatePoints(std::vectorPoint pts) { std::sort(pts.begin(), pts.end(), [](const Point a, const Point b) { if (a.x ! b.x) return a.x b.x; return a.y b.y; }); std::vectorPoint out; for (auto p : pts) { if (out.empty() || hypot(p.x - out.back().x, p.y - out.back().y) 1e-10) { out.push_back(p); } } return out; }这步清洗的代价很低但对后续所有计算的稳定性提升巨大。我在处理 LiDAR 点云数据时经常发现同一位置被多次采集如果不做去重生成的网格在重复点附近会有明显的裂缝或重叠三角形。5. 完整实现的主循环把上面的零件组装起来5.1 整体代码骨架有了前面所有的铺垫主循环的代码就比较清爽了。我把整个流程封装在一个类里核心接口就两个一个传入点集一个输出三角形列表。class DelaunayTriangulation { public: std::vectorTriangle triangulate(const std::vectorPoint inputPts) { pts deduplicatePoints(inputPts); triangles.clear(); freeList.clear(); // 1. 构造超级三角形 Triangle super makeSuperTriangle(pts); triangles.push_back(super); // 2. 逐点插入 for (auto p : pts) { insertPoint(p); } // 3. 删除所有包含超级三角形顶点的三角形 removeSuperTriangle(super); return triangles; } private: std::vectorPoint pts; std::vectorTriangle triangles; std::vectorint freeList; void insertPoint(const Point p) { std::vectorint badTriangles; for (int i 0; i (int)triangles.size(); i) { if (inCircumcircle(triangles[i], pts, p)) { badTriangles.push_back(i); } } std::mapEdge, int edgeCount; for (int tIdx : badTriangles) { for (int k 0; k 3; k) { int a triangles[tIdx].v[k]; int b triangles[tIdx].v[(k 1) % 3]; edgeCount[makeEdge(a, b)]; } } // 删除受影响三角形 for (int tIdx : badTriangles) { removeTriangle(tIdx); } // 连接新点与边界边生成新三角形 std::vectorTriangle newTriangles; for (auto kv : edgeCount) { if (kv.second 1) { Triangle t; t.v[0] kv.first.p1; t.v[1] kv.first.p2; t.v[2] addPoint(p.x, p.y); // 这里的索引会被后续复用 newTriangles.push_back(t); } } // 添加新三角形并初始化邻居关系 for (auto t : newTriangles) { addTriangleWithNeighbors(t); } // 对新三角形做 LOP 优化 std::stackint checkStack; for (auto t : newTriangles) checkStack.push(t.id); while (!checkStack.empty()) { int tIdx checkStack.top(); checkStack.pop(); for (int k 0; k 3; k) { legalizeEdge(triangles[tIdx], k, checkStack); } } } };这个版本的代码足够跑通中小规模数据。后面的流程就不复杂了但每一步的实现都有细节尤其是邻居索引的维护需要多花点心思。5.2 一个容易被忽略的细节addPoint 与 addTriangle 时的索引管理如果你按上面的思路写会发现 addPoint 时如果去重新点的索引可能不是简单地递增。这里我在pts里保存的是清洗后的去重点数组插入超级三角形时会把超级三角形的顶点也追加到 pts 里这样在处理三角形时顶点索引都指向 pts 数组。最后 removeSuperTriangle 时只需检查三角形里是否包含超级三角形的三个顶点索引把包含的剔除即可。5.3 邻居关系的初始化策略如果图省事新生成的三角形可以先完全不用设置邻居等所有三角形都生成后再做一遍“邻居构建”遍历。这个方法代码简单但多了一次全量扫描用时较长。另一种做法是每次插入新点后利用 badTriangles 原有的邻接关系来初始化新三角形的邻居。第二种效率高但代码逻辑较复杂。为了可靠我第一版用的还是“全量重建邻居”的方式正确跑通后再优化成增量方式。全量重建邻居的代码非常直观void rebuildNeighbors() { for (auto t : triangles) { t.neighbor[0] t.neighbor[1] t.neighbor[2] -1; } std::unordered_mapEdge, std::pairint,int edgeToTri; for (int i 0; i (int)triangles.size(); i) { for (int k 0; k 3; k) { Edge e makeEdge(triangles[i].v[k], triangles[i].v[(k1)%3]); if (edgeToTri.count(e)) { int j edgeToTri[e].first; int jk edgeToTri[e].second; triangles[i].neighbor[k] j; triangles[j].neighbor[jk] i; } else { edgeToTri[e] {i, k}; } } } }这里的 unordered_map 需要为 Edge 类型提供哈希函数可以直接用有序对的组合乘以一个大质数比较简单。6. 性能实测与优化方向从能用走向好用6.1 未优化版本的性能基准我在一台普通 i7 台式机上用随机生成的 10 万点做了测试坐标范围 [0, 10000] x [0, 10000]。未加任何空间索引的 Bowyer-Watson 实现耗时约 8.2 秒。这个结果比纯 O(n^2) 的预期好一些因为随机分布下每插入一个点实际影响的三角形数量比较有限平均可能在十几到几十个三角形之间。但如果数据分布极不均匀比如大量点密集在很小区域受影响三角形数量会显著增加耗时也会明显上涨。加了网格哈希索引后同一份数据降到约 1.1 秒。网格大小的选择对性能影响非常明显网格过小每个单元点太少需要检查的邻域单元多网格过大单元内三角形数量多退化成近似全量遍历。我的经验是让每个网格单元平均包含 4 到 8 个三角形时效果最好。这个经验值可以根据数据规模微调。6.2 内存上的进一步优化对于超大点集原来的 Triangle 结构体里每个三角形包含 3 个 int 和 3 个 int共 24 字节。10 万三角形大概是 2.4 MB看起来不大但如果是 5000 万点网格三角形数量轻松上亿内存会变得很紧张。这时候可以考虑用 int32 存储索引而不是 size_t。64 位系统下 size_t 是 8 字节int 是 4 字节省一半内存。把 neighbor 数组缩短成 3 个 int但需要用特殊值表示边界。如果不需要频繁访问邻居只在输出前重建邻居关系可以节省大量内存开销。6.3 并行化思路分块-合并策略Delaunay 三角化的并行化一直是计算几何里的热点。一个可行的策略是分治把点集按空间划分为多个子块各子块独立三角化最后在边界处做缝合。缝合过程比较复杂要处理跨子块的空外接圆检查和 LOP 优化。另一种思路是“波前并行插入”不过实现起来比较繁琐我自己暂时没有在生产环境里做大规模并行只做了分块收效一般。这块的改进空间还很大后续计划尝试 GPU 版本的增量式 Delaunay。7. 实测中踩过的三个大坑7.1 超级三角形不够大导致的边界断层第一版实现里我把超级三角形的半径设成包围盒直径的 1.5 倍以为足够了。结果在跑一个分布范围跨度大的点集时反复出现边界断裂有些点没有连入网格。排查很久才发现问题不是算法逻辑错了而是某些原本应该属于边界的点在计算外接圆时其外接圆超出了超级三角形的覆盖范围导致受影响三角形集合漏选。把半径改成 3 倍后问题彻底消失了。7.2 浮点误差导致的“错误翻转”在做 LOP 对角线交换时因为行列式阈值设置得不合理会出现一对反复交换的三角形第一次检查认为该翻转翻转后检查新的四边形又认为该翻转回去形成死循环。这是典型的浮点精度与阈值设计问题。我的解决方法是在 legalizeEdge 中如果检测到“两个三角形共用一个顶点且该顶点到圆心的距离与半径几乎相等”就保守地放弃交换用容错阈值给判断留出缓冲空间。这个处理让我避免了很多诡异的死循环和抖动。7.3 三角形遍历顺序引发的调试痛苦如果三角形的顶点顺序不统一一部分是顺时针、一部分是逆时针计算外接圆行列式时会导致符号混乱inCircumcircle 的结果时对时错非常难排查。我的经验是在构造函数里就强制所有顶点按逆时针顺序排列。每次创建或交换三角形后都做一次 orientation 检查必要时交换 v[1] 和 v[2]。虽然多了一点计算但整个后续逻辑的稳定性大幅提升。8. 实战衍生从三角网到 Voronoi 图与网格质量评估Delaunay 三角网和 Voronoi 图是一对对偶结构连好三角网后要生成 Voronoi 图只需要找每个三角形的外心然后连接相邻三角形的外心即可。这个延伸对于做最近邻搜索分析和计算几何可视化特别有用。我在项目里就用这套 C 三角网代码直接生成了 Voronoi 图用于蜂窝网格的形态分析效果很好。网格质量评估可以直接基于三角网计算每个三角形的最小角和最大角。通过统计最小角小于一定阈值的三角形比例可以衡量这套三角化实现在实际数据上的表现。比如随机均匀点云的三角网中最小角小于 30 度的三角形占比一般在 1% 到 3% 之间如果超过 5%基本说明 LOP 没做完整或数值判定有偏差。这个指标可以作为实现是否正确的快速检验。提示我建议每写完一版算法都用随机点与恶意点共线、共圆、大量退化分别测试把能力边界搞清楚。只看闪亮亮的渲染图并不能证明算法正确边界情形才是真正见功力的时候。实现一套 C 的 Delaunay 三角网算法虽然过程曲折但带来的收益远超过代码本身你会对整个增量式构造的逻辑、数值稳定性、邻接表维护、性能优化方向都建立直觉。这些经验在以后处理网格算法、路径规划、几何建模时都能复用。最后再说个小技巧调试时打开控制台模拟 20 个点的插入过程逐步对比每个阶段三角形的数量与邻居关系能让你快速发现逻辑漏洞。比起白白盯着代码看这种调试方式省力太多了。本文还有配套的精品资源点击获取