简介一份基于C的Delaunay三角网算法实现压缩包面向计算机图形学、地理信息系统与有限元分析的开发者和学习者解决二维点集最优三角剖分问题。包内共42个文件以C头文件与实现文件、Visual C工程配置、可执行程序及编译产物为主另含少量图标和文本资源整体仅2.66MB下载与本地编译成本低。实现覆盖数据结构定义、增量式构网、空圆检测与局部优化等核心步骤并带有可直接运行的exe和说明文档方便对照源码观察算法效果。已有1882人学习下载适合需要动手实践几何算法、提升C工程能力的读者。1. Delaunay三角网算法为什么C实现值得自己造轮子拿到一堆散乱坐标点要做地形网格、有限元剖分或者路径规划的可视化底图第一步往往逃不开三角网构建。Delaunay三角网的核心特性就一句话任何三角形的外接圆内不存在其他点同时最大化三角剖分中的最小角避免出现瘦长三角形。网上现成的几何库能调用但工作里一旦要嵌入自定义渲染管线、要控制顶点索引布局、要对接自研数据格式黑匣子就非常难受。这份C实现的资源直接把Bowyer-Watson逐点插入算法的数据结构、外接圆判断、边界重建、退化处理完整拆开适合需要亲手掌控网格生成的开发者。它解决的问题很具体给定一组二维点坐标输出符合Delaunay特性的三角形索引数组。2. Bowyer-Watson算法原理逐点插入与三个核心数据结构2.1 为什么选Bowyer-Watson而不是分治或扫描线主流的Delauany构网算法有三条路线。分治法递归分割点集再合并最坏时间复杂度O(n log n)很漂亮但合并时处理交叉边、维护递归状态代码量轻松上千行边界情况极其刁钻。扫描线法按x坐标顺序扫入点同时维护前沿边集状态机复杂碰到多点共线时容易崩。Bowyer-Watson算法走的是逐点插入路线每插入一个点找出所有外接圆包含该点的三角形删掉它们形成空洞再把空洞边界与新点连接。这个逻辑可以浓缩在一个函数里配合合理的点索引排序平均表现接近O(n log n)最坏退化到O(n²)但真实数据里很少触发。我一般在项目里首选Bowyer-Watson理由很朴素每一段逻辑都能单步调试出现异常时可以逐点回放不会陷入递归合并的泥潭。分治算法适合对性能极致敏感的生产级场景但如果目标是理解算法、快速落地并验证正确性逐点插入是性价比最高的起点。这份资源就是沿着这条路线实现的。2.2 数据结构定义点、三角形、边实现第一步是把三个基础结构定清楚。点结构很简单两个double存坐标。三角形结构需要三个顶点索引同时缓存外接圆圆心和半径平方插入过程中每个新点都要和大量三角形做外接圆包含判断如果每次现算圆心时间复杂度会翻好几倍所以建三角形时就把圆心和半径平方算好存进去。边结构用于重建空洞边界只需要两个端点索引。这里有个关键设计顶点用int索引而不是指针。因为插入过程中triangle数组会不断push_backvector扩容会导致指针失效而索引是稳定的。三角形里顶点顺序也必须固定为逆时针否则外接圆圆心公式的符号会错后面章节会专门说这个坑。struct Point2D { double x, y; }; struct Triangle { int v[3]; // 三个顶点索引v[0]-v[1]-v[2] 保持逆时针 double cx, cy; // 外接圆圆心坐标 double r2; // 外接圆半径平方 }; using Edge std::pairint, int; // 边存储时保证 first second三角形缓存外接圆参数的意义在于inCircle判断从“计算圆心计算距离”变成“一次距离比较”省掉大量重复的叉积运算。参数说明r2存的是半径平方而不是半径因为后续比较距离平方与半径平方就能判断包含关系无需开根号这在十万级点集上能明显缩短耗时。2.3 外接圆判断只用乘法不开根号外接圆圆心计算用的是标准叉积公式。给定三个逆时针点a(x_a, y_a), b(x_b, y_b), c(x_c, y_c)先算分母d再算圆心坐标。d接近零说明三点共线此时外接圆退化成半径无穷大的情况需要单独标记。void computeCircumcircle(const Point2D a, const Point2D b, const Point2D c, double cx, double cy, double r2) { double d 2.0 * (a.x * (b.y - c.y) b.x * (c.y - a.y) c.x * (a.y - b.y)); if (std::fabs(d) 1e-12) { // 三点近似共线返回一个“无穷大”的外接圆 cx cy 0.0; r2 1e18; return; } double a2 a.x * a.x a.y * a.y; double b2 b.x * b.x b.y * b.y; double c2 c.x * c.x c.y * c.y; cx (a2 * (b.y - c.y) b2 * (c.y - a.y) c2 * (a.y - b.y)) / d; cy (a2 * (c.x - b.x) b2 * (a.x - c.x) c2 * (b.x - a.x)) / d; r2 (a.x - cx) * (a.x - cx) (a.y - cy) * (a.y - cy); } bool inCircle(const Point2D p, const Triangle t) { double dx p.x - t.cx; double dy p.y - t.cy; return dx * dx dy * dy t.r2; }注意d的符号依赖顶点逆时针排列。如果三角形顶点是顺时针d会反号但圆心坐标公式里分子分母同时反号结果仍然正确所以这个公式其实对顺逆时针都能给出正确圆心。真正要命的是理论上共线时d0浮点下会得到一个随机方向的大圆心所以先判d的绝对值是否小于阈值。参数说明阈值1e-12是经验值适合点坐标在1e-6到1e6量级的场景坐标特别大或特别小时要按包围盒比例调整。3. C实现主流程从点集到三角网的核心代码3.1 超三角形怎么选边界才不会被误删Bowyer-Watson的第一步是构造一个能包住所有输入点的“超三角形”。这个三角形足够大使得任何输入点都不会落在它外面但同时也不能大到让浮点误差失控。常见做法是在点集的包围盒基础上外扩很大倍数。我一般取包围盒最大边长的20倍作为外扩距离太小会导致后续插入时真实边界点被误删太大会让外接圆半径巨大前几个点的插入几乎要遍历所有已有三角形造成性能浪费极端情况还会让圆心坐标超出double精度。std::vectorPoint2D points; // 输入点集 // 求包围盒 double min_x 1e18, min_y 1e18, max_x -1e18, max_y -1e18; for (const auto p : points) { min_x std::min(min_x, p.x); min_y std::min(min_y, p.y); max_x std::max(max_x, p.x); max_y std::max(max_y, p.y); } double dx max_x - min_x, dy max_y - min_y; double margin std::max(dx, dy) * 20.0; std::vectorPoint2D allPoints points; // 把超三角形顶点追加到点集末尾 int superA allPoints.size(); allPoints.push_back({min_x - margin, min_y - margin}); int superB allPoints.size(); allPoints.push_back({min_x 2.0 * margin, min_y - margin}); int superC allPoints.size(); allPoints.push_back({min_x, min_y 2.0 * margin}); std::vectorTriangle tris; Triangle superTri; superTri.v[0] superA; superTri.v[1] superB; superTri.v[2] superC; computeCircumcircle(allPoints[superA], allPoints[superB], allPoints[superC], superTri.cx, superTri.cy, superTri.r2); tris.push_back(superTri);这段代码的关键是margin的选取。用max(dx, dy)*20而不是固定值是为了让不同量级的坐标都能获得合适的尺度。超三角形顶点必须按逆时针顺序排列这里把superA放在左下、superB放在右下、superC放在正上方天然满足逆时针。参数说明20倍这个系数是经验值10倍在点集分布极端时可能不够50倍以上在坐标值很大时会损失精度20是一个安全区间。3.2 逐点插入与坏三角形删除核心插入逻辑分四步标记坏三角形、收集边、去重提取边界、重建三角形。坏三角形指外接圆包含新插入点的三角形。把这些三角形所有边收集起来如果一条边被两个坏三角形共享它是内部边在空洞重建时要丢弃如果只出现一次它是边界边要用来和新点连成新三角形。void insertPoint(int pi, std::vectorTriangle tris, const std::vectorPoint2D allPoints) { std::vectorchar badFlag(tris.size(), 0); std::vectorEdge edges; for (size_t i 0; i tris.size(); i) { if (inCircle(allPoints[pi], tris[i])) { badFlag[i] 1; for (int j 0; j 3; j) { int a tris[i].v[j], b tris[i].v[(j 1) % 3]; if (a b) std::swap(a, b); edges.push_back({a, b}); } } } // 边去重出现两次的是内部边只出现一次的是边界边 std::sort(edges.begin(), edges.end()); std::vectorEdge boundary; for (size_t i 0; i edges.size();) { size_t j i; while (j edges.size() edges[j] edges[i]) j; if ((j - i) % 2 1) boundary.push_back(edges[i]); i j; } // 删除坏三角形用 remove_if 配合 badFlag std::vectorTriangle keep; keep.reserve(tris.size()); for (size_t i 0; i tris.size(); i) { if (!badFlag[i]) keep.push_back(tris[i]); } tris.swap(keep); // 用边界边构造新三角形 for (const auto e : boundary) { Triangle nt; nt.v[0] e.first; nt.v[1] e.second; nt.v[2] pi; computeCircumcircle(allPoints[nt.v[0]], allPoints[nt.v[1]], allPoints[nt.v[2]], nt.cx, nt.cy, nt.r2); tris.push_back(nt); } }这段代码最容易被忽略的是边去重的方式。排序后统计每组相同边的数量奇数则保留偶数则丢弃这样内部边自然抵消、边界边留下不需要哈希表。参数说明边界边的顺序不影响最终结果但新三角形的顶点顺序不能直接沿用e.firste.second, pi的排列因为这样可能产生顺时针三角形需要在computeCircumcircle前修正。3.3 清理与输出移除超三角形相关三角形全部点插入完成后凡是包含超三角形三个顶点中任意一个的三角形都必须删除。这些三角形是外围的“假三角形”不包含真实数据点。void removeSuperTriangle(std::vectorTriangle tris, int superA, int superB, int superC) { std::vectorTriangle keep; keep.reserve(tris.size()); for (const auto t : tris) { bool hit false; for (int i 0; i 3; i) { if (t.v[i] superA || t.v[i] superB || t.v[i] superC) { hit true; break; } } if (!hit) keep.push_back(t); } tris.swap(keep); }输出格式建议用最朴素的文本第一行两个整数点数N和三角形数M接着N行是顶点坐标再M行是每行的三个索引。这样下游无论是接渲染管线还是写离线检查工具都很方便。清理时机必须是在所有点插入完成后而不是每插一个点就清理一次。原因在于后续插入的新点可能仍然需要与后台的“过渡三角形”发生连接提前删掉会让空洞边界不完整。另外输出前最好验证每个三角形的r2是否非负如果是负数说明外接圆计算出了问题这往往是顶点顺序导致的。4. 边界与退化场景四点共圆、重复点、共线点的处理4.1 重复点与极近点插入前的清洗输入的原始点集里经常混入重复点或距离极近的点。Bowyer-Watson碰到两个完全重合的点时第二个点会落在已有三角形的外接圆圆周上inCircle返回true但生成的三角网会出现零面积三角形。更隐蔽的问题是极近点会让外接圆圆心计算中的d值接近零得到一个数值不稳定的巨大外接圆。好的做法是在插入前按x坐标排序然后线性扫描把距离小于阈值的相邻点合并阈值取包围盒对角线长度的1e-9倍比较稳妥。注意这里的排序结果还可以复用Bowyer-Watson按x递增顺序插入点可以让每个新点影响的三角形集中在局部区域减少遍历范围。std::vectorPoint2D deduplicate(std::vectorPoint2D pts) { std::sort(pts.begin(), pts.end(), [](const Point2D a, const Point2D b) { return a.x b.x || (a.x b.x a.y b.y); }); std::vectorPoint2D out; double eps 1e-9; // 实际应按包围盒对角线比例动态计算 for (const auto p : pts) { if (out.empty() || std::fabs(p.x - out.back().x) eps || std::fabs(p.y - out.back().y) eps) { out.push_back(p); } } return out; }注意用fabs(p.x - out.back().x) eps而不是直接计算欧氏距离因为排序后x接近的点才有可能重复先比较x可以快速跳过大部分点只有x相近时才需要继续比较y。这里eps用固定值1e-9只是示意实际应该根据点集包围盒动态计算为max(dx, dy) * 1e-10避免坐标量级不同导致误删。4.2 共圆点外接圆判断的EPSILON策略四个点共圆是Delaunay构网里最经典的玄学问题。从几何定义看共圆时两种连接方式都满足空外接圆特性算法理论上任选其一即可。但浮点运算下“共圆”这个状态几乎不可能严格出现内点是否在外接圆内的判断会抖动导致插入不同点时选到不同的连接方式最终结果可能是一个“非局部Delaunay”的网格。处理方式是给inCircle判断加一个相对容忍度。判断时不是r2 距离平方就返回true而是距离平方 r2 EPS其中EPS要按坐标尺度动态计算。如果全用一个绝对值比如1e-12点坐标在1e4和1e8时表现完全不同。bool inCircleEps(const Point2D p, const Triangle t, double eps2) { double dx p.x - t.cx; double dy p.y - t.cy; return dx * dx dy * dy t.r2 eps2; }调用时传入的eps2取全点集包围盒对角线长度的平方乘以一个极小系数比如1e-20。这样做的目的是把浮点误差吸收掉让“理论上共圆”的点被一致性处理。但提醒一点eps2过大会把本不应该包含的点也判定为坏导致生成的三角形穿透真实点集。这个参数需要根据实际数据微调我一般从1e-18起步每个量级测试一遍直到输出网格的min角不再抖动。4.3 三角形顶点顺序逆时针存储的陷阱新建三角形时如果顶点顺序是顺时针外接圆圆心计算式里的d值会反号虽然圆心坐标偶然正确但在后续的LOP翻转和面积计算中会持续产生错误。修正顶点顺序用叉积判断按(v0, v1, v2)计算有向面积若为负则交换v1和v2。void ensureCCW(Triangle t, const std::vectorPoint2D pts) { const Point2D a pts[t.v[0]]; const Point2D b pts[t.v[1]]; const Point2D c pts[t.v[2]]; double cross (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x); if (cross 0) std::swap(t.v[1], t.v[2]); }这里会踩的坑是超三角形的顶点顺序本身已经是逆时针但新插入点构建的三角形可能因为边界边的搜索顺序而变成顺时针。所以每次构造新三角形后必须强制调用ensureCCW且必须在computeCircumcircle之前调用因为圆心公式对顶点顺序敏感。建议把这两步封装成一个makeTriangle函数避免漏调。5. 常见问题与避坑五个容易翻车的位置5.1 现象十万点插入时程序卡到怀疑人生原因有两个一是每个点插入时都线性扫描全部三角形复杂度退化到O(n²)二是inCircle判断里写了sqrt每次比较都做一次浮点开根号代价是乘法比较的五到十倍。解决先用按x排序优化插入顺序再把r2直接存成半径平方比较时全部用平方量级。排序的收益在于新插入点只需检查x坐标接近的三角形配合一个简单的网格索引可以把实际扫描量降到常数级。5.2 现象边界三角形被莫名其妙删除原因是超三角形选小了真实点落在超三角形边上甚至外面导致外围三角形的外接圆被这些“外部点”影响清理时误删。解决检查margin系数不要低于20。另一个隐蔽原因是超三角形坐标取值过大比如直接用了int的极限值1e9导致外接圆计算时x²项溢出变成inf所有点都被判定为坏三角形。超三角形坐标应该基于点集包围盒的有限外扩而不是一个拍脑袋的大数。5.3 现象空洞重建后出现裂缝或重叠边根源在边去重这一步。如果收集坏三角形的边时不把端点归一化让(a,b)和(b,a)同时存在于列表里排序后它们不会相邻无法正确配对抵消边界边集合里就会混入内部边。解决在push边之前强制交换保证first second然后再排序统计。这一步是漏掉最多的地方写完之后用五个点的简单用例手工验证一遍空洞边界是否闭合。5.4 现象四点共圆时三角网出现非局部翻转这是最玄学的一个同样的输入点改变插入顺序输出的对角线不同。原因不是算法错了而是浮点下外接圆包含判断在共圆点附近不稳定。解决给inCircle加动态EPSILON值让共圆点被一致处理如果对稳定性要求极高可以在这四个点周围做一次局部LOPLawson翻转强制把不满足Delaunay条件的对角线换过来。注意EPSILON太小没用、太大会把好三角形也翻坏需要实际数据里测试。5.5 现象程序崩溃在删除坏三角形的erase循环里用索引遍历并原地erase时vector元素前移导致后续索引全部失效。解决像前面代码那样先标记badFlag再两阶段删除——收集保留的三角形到keep向量再swap。这个坑在C里特别容易犯因为用erase配合迭代器看起来天经地义但边收集和边删除交错时迭代器就悬空了。尽早改成标记回收模式后这个问题就再也没出现过。6. 验证与优化让三角网从能跑到好用6.1 空外接圆验证器自己检查算法对不对算法写完先别急着上真实数据写一个O(n²)的验证器对每个三角形、每个非顶点数据点用inCircle判断是否有外点落入外接圆内部。只要有一个违规打印这个三角形的三个索引和违规点的坐标。这个验证器只在小数据量时用我一般随机生成100到1000个点反复跑二十轮确认零违规后再接入生产管线。bool verifyDelaunay(const std::vectorPoint2D pts, const std::vectorTriangle tris) { for (const auto t : tris) { for (int i 0; i (int)pts.size(); i) { bool isVertex (i t.v[0] || i t.v[1] || i t.v[2]); if (isVertex) continue; double dx pts[i].x - t.cx; double dy pts[i].y - t.cy; if (dx * dx dy * dy t.r2 - 1e-12) { return false; } } } return true; }参数说明这里用t.r2 - 1e-12而不是直接t.r2是给验证器一个极小容差避免浮点边缘情况误报如果生产数据里出现连续失败加大EPSILON再做一轮。6.2 局部翻转LOP修复浮点撕裂Bowyer-Watson生成的网格理论上是Delaunay但浮点误差可能让某些共边三角形对不满足条件。LOP的思路是遍历每对相邻三角形检查它们共享的边是否应该翻转如果其中一个三角形的第四个点落入另一个三角形的外接圆内把公共边换成另一条对角线。这个过程能局部修复共圆和误差问题。void legalizeEdge(Triangle t1, Triangle t2, int sharedA, int sharedB, const std::vectorPoint2D pts) { int p3 /* t1 中非 sharedA/sharedB 的顶点 */; int p4 /* t2 中非 sharedA/sharedB 的顶点 */; if (inCircle(pts[p4], /* 由 sharedA, sharedB, p3 构造的临时三角形 */)) { // 翻转新三角形 (sharedA, p3, p4) 和 (sharedB, p3, p4) } }实际编码时LOP需要维护三角形邻接关系这要求数据结构里存邻居索引。这份资源如果没带邻接信息可以只在最后的验证阶段跑特定局部修复而不做全量LOP。多数场景下Bowyer-Watson加EPSILON已经够用。6.3 实践顺序与最后的习惯完整跑通的流程我一般固定成这样先用五个手算点验证中心插入和边界清理逻辑再上两百个随机点配合验证器检查空外接圆特性然后跑一万点看看性能和稳定性最后才接真实坐标。数据量超过十万时把点按x排序并加速扫描逻辑这时候Bowyer-Watson才算真正落地。这套流程帮我排掉过无数个隐蔽问题。以前我图省事写完直接扔进去十万个地形点结果边界漏了一圈三角形排查了一整天才发现是超三角形margin系数太小导致外围点误删。从那以后我每次跑真实数据之前都强制自己先走一遍小样本验证器再检查一遍超三角形覆盖范围最后才放数据进去。这个习惯救过我好几次希望帮到你。本文还有配套的精品资源点击获取