★ FEATURED ARTICLE 简介这份资源是用C实现Delaunay三角网算法的完整工程代码面向计算机图形学、地理信息系统与有限元方向的开发者及算法学习者帮助理解空圆特性、逐点插入与空圆检测等核心构建逻辑。压缩包共42个文件约2.66MB以h头文件与cpp源文件为主体另有obj、sbr等编译中间文件以及ico、bmp、cur等界面资源并附带dsp、dsw工程文件与可执行程序便于直接编译运行与调试。资源已有1882人学习下载代码结构完整涵盖点与三角形数据结构设计、三角网构建、空圆检测、悬挂边与孤岛优化等关键环节读者可结合源码逐步跟踪算法流程验证三角剖分结果的正确性并在此基础上扩展可视化输出或移植到实际项目中适合作为算法入门与工程实践的参考范例。1. Delaunay三角网算法用C实现从点云到网格一次把边界和性能讲透做点云处理、地形建模或者有限元前处理的人迟早会撞上 Delaunay 三角网。它解决的问题很具体给你一堆散乱的点怎么连成一张“形状最合理”的三角网格。所谓合理核心就是空外接圆准则——任意一个三角形的外接圆内部不允许再落进其他点。这条约束带来的直接好处是三角形尽量接近等边不会出现又细又长的“针状”单元后续做插值、做有限元、做渲染都稳。用 C 实现它动机通常有两个一是性能点数上到十万、百万级时脚本语言的常数开销扛不住二是可控你需要自己决定数值精度、内存布局、退化点怎么处理。这篇笔记就按“先立住原理再动手复现最后讲坑”的顺序走把增量插入法、Bowyer-Watson 空腔重构、边翻转这几条主线落到能编译运行的代码上也会讲清楚哪些参数一动就翻车。适合已经会写 C、但第一次自己撸三角剖分的同学也适合想替换掉手头慢库的老手。2. 先把算法骨架立住增量插入与空外接圆判定2.1 为什么选增量插入而不是分治Delaunay 三角网的经典构造有三条路分治法、扫描线法、增量插入法。分治法理论复杂度漂亮O(n log n)但实现里递归切分、合并凸包、处理跨子集的边代码量和调试成本都高。扫描线法对点的排序和事件队列要求苛刻退化情况多。增量插入法是工程上最常被选中的逻辑直白一个点一个点插进去局部修复出问题容易定位。增量插入的核心循环只有三步定位新点落在哪个三角形里删除所有外接圆包含新点的三角形形成一个空腔把空腔边界和新点连起来重新三角化。第三步之后还要做一次合法性检查把不满足空外接圆准则的边翻转掉。听起来简单但“定位”和“翻转”这两步是性能与正确性的分水岭。我一般会先搭一个足够小的骨架用随机点跑通再逐步加规模。下面是最小可运行的结构定义点、边、三角形都用索引表达避免指针满天飞。#include vector #include cmath #include cstdint struct Point { double x, y; }; struct Triangle { int v[3]; // 三个顶点在点数组中的索引 int nbr[3]; // 三条边对应的邻接三角形索引-1 表示边界 bool alive; // 是否仍有效被删除的三角形标记为 false }; // 判断点 p 是否在三角形 abc 的外接圆内 // 用行列式符号判定避免直接算圆心和半径带来的开方误差 bool inCircumcircle(const Point a, const Point b, const Point c, const Point p) { double ax a.x - p.x, ay a.y - p.y; double bx b.x - p.x, by b.y - p.y; double cx c.x - p.x, cy c.y - p.y; double det (ax * ax ay * ay) * (bx * cy - by * cx) - (bx * bx by * by) * (ax * cy - ay * cx) (cx * cx cy * cy) * (ax * by - ay * bx); // 假设 abc 为逆时针det 0 表示 p 在圆内 return det 0.0; }这段代码的关键在inCircumcircle。它没有去求外接圆圆心而是用四点构成的 3x3 行列式符号来判断。这样做的好处是全程只有乘法和减法没有除法数值稳定性更好也更快。参数上要注意det 0这个方向依赖三角形顶点是逆时针排列的。如果你在别处把顶点顺序搞反了判定会整体翻转表现为“该删的三角形不删网格越插越乱”。所以我在插入前会统一保证三角形是逆时针用有向面积校验一次。2.2 超级三角形与边界处理增量插入有个绕不开的工程问题第一个点、前几个点怎么处理。常见做法是构造一个“超级三角形”把全部输入点包进去插入过程结束后再删掉所有与超级三角形顶点相连的三角形。超级三角形要足够大否则边界点会落在它外面导致定位失败。// 根据点集包围盒构造超级三角形 // 放大系数取 20 是经验值太小会漏点太大影响精度 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 (const 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, dy maxy - miny; double dmax std::max(dx, dy) * 20.0; double midx (minx maxx) * 0.5, midy (miny maxy) * 0.5; // 三个顶点按逆时针给出 Point a{midx - dmax, midy - dmax}; Point b{midx dmax, midy - dmax}; Point c{midx, midy dmax}; // 这里把三个点追加进点数组返回索引构成的三角形 // 实际实现中需要维护点数组与三角形数组的同步 return Triangle{{0, 1, 2}, {-1, -1, -1}, true}; }放大系数 20 不是随便写的。太小比如 2 到 5边界附近的点可能落在超级三角形外定位阶段会找不到宿主三角形太大比如 1e6会让外接圆判定里的坐标差值巨大浮点有效位被吃掉退化判定开始误判。20 到 100 之间是比较稳的区间。另外超级三角形的三个顶点要按逆时针给和前面inCircumcircle的方向假设保持一致。2.3 定位从哪个三角形开始找定位新点所在三角形最笨的办法是遍历所有存活三角形O(n) 一次整体就退化成 O(n²)。工程上常用“行走法”从上一个插入点所在的三角形出发沿着有向边朝目标点方向走平均几步就能到。行走法依赖邻接关系nbr所以维护邻接的正确性是前提。// 行走法定位从 start 出发找到包含 p 的三角形 // 返回三角形索引若 p 落在边上返回任一相邻三角形即可 int locateTriangle(int start, const Point p, const std::vectorPoint pts, const std::vectorTriangle tris) { int cur start; for (int guard 0; guard 10000; guard) { const Triangle t tris[cur]; bool moved false; for (int i 0; i 3; i) { int a t.v[i], b t.v[(i 1) % 3]; // 有向面积p 在边 ab 的右侧则朝该邻接走 double cross (pts[b].x - pts[a].x) * (p.y - pts[a].y) - (pts[b].y - pts[a].y) * (p.x - pts[a].x); if (cross 0.0) { int next t.nbr[i]; if (next 0) return cur; // 到边界认为落在当前三角形 cur next; moved true; break; } } if (!moved) return cur; } return cur; // 兜底正常不该走到这里 }guard这个计数器是后悔药。行走法在数值退化或者邻接关系被写坏时可能死循环加一个上限至少不会把程序挂死方便你回头查邻接维护的 bug。cross 0的方向同样依赖逆时针约定。定位返回后如果点恰好落在边上不同实现会返回不同侧的三角形这会影响后续空腔的起点但不影响最终正确性只要你的空腔扩展逻辑是双向的。3. Bowyer-Watson 空腔重构删哪些、连哪些、怎么翻转3.1 空腔扩展与边界收集定位到宿主三角形后从它开始做广度优先扩展凡是外接圆包含新点的三角形全部标记删除直到遇到不包含的为止。这些被删三角形的外边界就构成空腔的多边形边界。新点与这条边界上的每条边连成新三角形。#include queue #include unordered_set // 返回空腔边界边每条边用 (端点a, 端点b) 表示方向沿空腔逆时针 std::vectorstd::pairint,int buildCavity( int seed, const Point p, const std::vectorPoint pts, std::vectorTriangle tris) { std::vectorstd::pairint,int boundary; std::queueint q; std::unordered_setint visited; q.push(seed); visited.insert(seed); while (!q.empty()) { int ti q.front(); q.pop(); Triangle t tris[ti]; if (!inCircumcircle(pts[t.v[0]], pts[t.v[1]], pts[t.v[2]], p)) { continue; // 该三角形保留不作为空腔一部分 } t.alive false; for (int i 0; i 3; i) { int nb t.nbr[i]; if (nb 0 !visited.count(nb)) { visited.insert(nb); q.push(nb); } else if (nb 0) { // 边界边记录 boundary.push_back({t.v[i], t.v[(i 1) % 3]}); } } } return boundary; }这里有个容易忽略的点visited必须在入队时就插入不能等出队。否则同一个三角形可能被多个邻居重复入队空腔一大队列就爆。另外边界边的方向要沿空腔逆时针这样新三角形(a, b, p)才是逆时针和全局约定一致。如果方向反了后面翻转判定会全乱。3.2 边翻转把非法边掰回来空腔重构后新生成的三角形不一定满足空外接圆准则需要检查空腔边界上的每条边如果它和对面三角形构成的四边形是凸的且对面三角形的对角点落在当前三角形外接圆内就翻转这条边。翻转要迭代进行直到没有非法边为止。// 尝试翻转边 (a,b)成功返回 true // 需要维护邻接关系这里只给核心判定 bool tryFlip(int a, int b, std::vectorTriangle tris, const std::vectorPoint pts) { // 找到共享边 ab 的两个三角形 t1, t2 // 取 t1 中 ab 的对顶点 ct2 中 ab 的对顶点 d // 若 c 在 t2 的外接圆内则翻转 ab 为 cd // 翻转后更新四个邻接关系 // 具体索引查找略重点在判定条件 return false; // 占位实际实现按上述逻辑 }翻转的判定条件写起来琐碎但逻辑就一句四边形凸 对角点在外接圆内。两个条件缺一不可。只判外接圆不判凸性会把凹四边形也翻产生重叠三角形。翻转的终止性由 Delaunay 能量函数单调下降保证但前提是每次翻转都合法。工程上我会给翻转循环加一个迭代上限比如边数的 10 倍超过就说明有退化点或者数值问题需要单独处理。3.3 一个能跑的最小主循环把上面几块拼起来主循环就是构造超级三角形逐个插入点定位、建空腔、连边、翻转最后清理超级三角形相关的三角形。void delaunay(std::vectorPoint pts) { // 1. 追加超级三角形顶点构造初始三角形 // 2. 对每个原始点 p // - locateTriangle 找到宿主 // - buildCavity 收集边界 // - 用 p 和边界边生成新三角形维护邻接 // - 对受影响的边做 tryFlip 迭代 // 3. 删除所有含超级三角形顶点的三角形 // 4. 压缩索引输出 }主循环本身不复杂复杂的是邻接维护。每生成一个新三角形要同时更新它自己和邻居的nbr。我习惯在生成后立刻做一次一致性校验每个三角形的nbr[i]指向的三角形必须也有一条边指回来。这个校验在调试期开着能省掉大量“网格看着对、查询全错”的玄学问题。4. 避坑与排查那些让网格悄悄错掉的细节4.1 现象网格出现重叠三角形面积和为负原因通常是顶点方向约定被破坏。inCircumcircle和locateTriangle都假设三角形逆时针一旦某处生成三角形时顶点顺序反了判定方向整体翻转空腔扩展会删错三角形翻转也会翻错边。解决方法是每次生成三角形后用有向面积cross (b-a)×(c-a)校验小于 0 就交换两个顶点。这个校验开销很小但能挡住绝大多数方向类 bug。4.2 现象插入到几千个点后速度骤降原因多半是定位退化成线性扫描。行走法在点分布不均匀时会走很多步尤其是点沿一条线分布时。解决方法是给定位加一个“历史起点”缓存每次从上一个插入点的宿主开始走而不是从头。另外可以按空间分块先粗定位到块再在块内行走。实测在十万级随机点上历史起点能把定位步数从平均几百降到个位数。4.3 现象共圆点导致翻转死循环四个点共圆时空外接圆判定det接近 0浮点误差会让翻转在两个状态间反复横跳。解决方法是给判定加一个相对容差比如det eps * scale其中scale是坐标量级的平方。容差不能太大否则会破坏 Delaunay 性质一般取 1e-12 到 1e-9 之间按你的坐标范围调。另外翻转循环必须有迭代上限超限就停止并记录不要让它无限跑。4.4 现象边界点丢失或网格不闭合原因通常是超级三角形太小边界点落在它外面定位返回了错误的宿主。解决方法是把放大系数调到 20 以上并且在插入前检查每个点是否在超级三角形内部不在就报错。另一个原因是清理超级三角形时把与之相邻的合法三角形也删了需要按顶点索引精确判断只删含超级三角形顶点的三角形。4.5 现象内存暴涨程序被 OOM 杀掉原因通常是删除的三角形没有真正释放或者邻接关系里保留了悬空索引。解决方法是把alive标记和索引压缩分开插入阶段只标记最后统一压缩。压缩时重建v和nbr把所有指向已删除三角形的邻接置为 -1。这样内存峰值可控也避免悬空索引导致的随机崩溃。5. 进阶技巧用局部重建和数值技巧把稳定性再抬一档走到这里一个能跑的 Delaunay 三角网已经有了。但如果你的点集上到百万级或者坐标范围跨越好几个数量级还有两件事值得做。第一件是局部重建代替全局翻转。每次插入点后不要对整个网格做翻转检查只对空腔边界上的边做。更进一步可以把空腔边界边按“是否可能非法”排序优先处理最可能非法的减少无效判定。实测在随机点集上局部重建能把翻转次数压到全局方案的三成左右。第二件是坐标归一化。把输入点先平移到包围盒中心再缩放到 [-1, 1] 区间做完三角剖分再映射回去。这样做的好处是inCircumcircle里的行列式量级统一容差好设退化判定稳定。归一化本身是 O(n)相对于剖分可以忽略。下面是一个归一化和反归一化的最小实现。struct Normalizer { double cx, cy, scale; void fit(const std::vectorPoint pts) { double minx pts[0].x, maxx pts[0].x; double miny pts[0].y, maxy pts[0].y; for (const 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); } cx (minx maxx) * 0.5; cy (miny maxy) * 0.5; scale std::max(maxx - minx, maxy - miny); if (scale 0.0) scale 1.0; } Point apply(const Point p) const { return {(p.x - cx) / scale, (p.y - cy) / scale}; } Point revert(const Point p) const { return {p.x * scale cx, p.y * scale cy}; } };scale取包围盒长宽较大者保证归一化后坐标落在 [-0.5, 0.5] 附近。如果所有点重合scale会是 0这里兜底成 1避免除零。归一化之后前面提到的容差可以固定成 1e-12 左右不用再按数据调省心不少。还有一个我踩过的坑不要用float存坐标。float只有 24 位有效位坐标范围稍大一点inCircumcircle的行列式就会因为抵消误差给出错误符号表现为网格局部错乱而且极难复现。用double内存翻倍但正确性有保障。如果确实要省内存可以用定点数但那是另一个话题了。最后说一个验证习惯每次改完邻接维护或者翻转逻辑我都会用一组手工构造的退化点集跑一遍——正方形四角、正三角形加中心点、共线点。这些点集规模小出问题能肉眼看出哪个三角形错了。比拿十万随机点跑半天再回头查高效得多。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
