简介这是一份用C实现增量Delaunay三角剖分算法的完整工程资源面向计算机图形学、几何算法与网格处理方向的开发者和学生重点演示如何通过“半边数据结构”维护动态三角网并实时保证三角形最小角度最大化。资源包共含89个文件包括Visual Studio解决方案.sln、.vcxproj、C头文件与实现.h、.cpp、可直接运行的exe、GIF演示动图以及编译中间文件等整包大小约4.17MB目录结构包含GaussCurvature_Sample示例模块和MeshLib库便于对照学习和二次开发。核心代码集中在GaussCurvature_Sample的main.cpp中清晰展示了增量插入点、局部翻转等关键步骤配合压缩包内附带的GIF动图和编译好的exe可以直观看到三角剖分从迭代到收敛的动态过程。目前已有998人学习浏览适合用来掌握Delaunay三角剖分的算法流程、半边数据结构的设计思路以及Visual Studio下相关工程的组织方式读者还可基于此扩展自己的网格生成或几何处理算法。 前阵子在做一个点云网格化工具输入是几千个散乱坐标点输出需要交给下游做纹理映射和碰撞检测。最初偷懒直接用了现成的开源库功能没问题但一遇到定制需求——比如要保留边界、要支持动态增点——就得翻源码改接口很不顺手。于是我决定自己用 C 写一个增量 Delaunay 三角剖分实现也就是标题里的 DelaunayTriangulation 项目。这篇文章把整个过程完整摊开核心原理、数据结构选型、关键代码、性能优化以及调试时踩过的那些坑给正在学 C 或者需要把几何算法落地的同学一个可直接参考的版本。1. 为什么最终选了增量法而不是 Bowyer-Watson1.1 Delaunay 三角剖分到底在解决什么问题先明确一下我们要做的事。给定平面上 n 个散点把这组点用不相交的三角形连起来覆盖所有点的凸包这就是一个三角剖分。但三角剖分不唯一比如四个点构成的正方形连哪条对角线都合法。Delaunay 三角剖分在其中加入了一条准则任何一个三角形的外接圆内部不能包含其他顶点。这条空外接圆性质带来一个直观的好处——三角形会尽量饱满、接近等边不会出现特别尖锐的狭长三角形。对后续做插值、有限元、地形网格生成来说这个性质非常关键。Delaunay 三角剖分的实现路线有好几条分治法、扫描线法、逐点插入法。分治法理论复杂度最优能达到 O(n log n)但实现复杂度高递归切分、合并时需要处理大量跨子集的边调试阶段很容易崩。扫描线法对排序依赖强代码也不短。逐点插入法incremental insertion虽然最坏情况 O(n^2)但实现简单、易于理解、还天然支持动态加点对大多数工程场景已经足够。我不会一上来就追求理论最优能稳定跑、能改得动比任何复杂度公式都重要。1.2 两种逐点插入路线的取舍同样是逐点插入现实中还有两条分支。一条是 Bowyer-Watson 算法每插入一个点先找出所有外接圆包含新点的三角形删掉这些三角形形成一个空洞再把空洞边界上的顶点与新点逐一连接完成重构。这条路线写起来短但有个隐患——它需要一次批量删除再批量重建空洞区域的边界收集和去重逻辑容易出错而且对退化情况多个点共圆、共线非常敏感。我采用的是另一条路线Lawson 算法也叫翻边法。它的核心思想不是先删再建而是先找到新点落在哪个三角形里把它拆成三个新三角形然后对外接圆性质被破坏的邻接边逐条翻转直到整个局部重新满足 Delaunay 条件。翻边法的好处是每次操作只影响局部几条边拓扑变化范围小特别适合需要动态增点的场景——球员走一步只挪一辆车不用把整条停车场重新调度一遍。而且调试时每次只看一条边的处理结果心理负担小很多。2. 先把数据结构定死索引、逆时针与三角形栈2.1 顶点、三角形与邻接关系怎么存写几何算法数据结构设计比算法本身更容易决定项目成败。一开始我打算用指针互连的方式组织三角形每个三角形持有三个指向邻居的指针结果实现到一半就发现一个问题vector 扩容会导致指针失效只要三角形数量超过预分配容量所有邻居指针都变成悬垂指针排查起来非常痛苦。后来我改成纯索引式结构整个三角网只依赖两个数组struct Vec2 { double x, y; }; struct Vertex { Vec2 pos; int id; }; struct Triangle { int v[3]; // 三个顶点索引约定逆时针顺序 int t[3]; // 三条邻接三角形索引t[i] 对应边 (v[(i1)%3], v[(i2)%3]-1 表示边界 };用索引代替指针vector 扩容时只需要重新分配内存所有索引数值不变彻底绕开了引用失效问题。代价是取邻居时要多一次数组下标访问但现代 CPU 对连续内存的访问远快于对散落堆对象的指针追踪这点牺牲换来了稳定性和缓存友好性很值。2.2 为什么一律用逆时针约定坐标结构本身不复杂真正让代码变脆的是方向约定不一致。我吃过亏有的函数按外接圆包含点判断有的函数按点在边的哪一侧判断两个函数用了不同的顶点顺序出来的结果符号恰好相反bug 找了一下午。后来我立了一条铁规矩所有三角形顶点一律按逆时针存储。这样点在边左侧和点在三角形内部的判断共用同一个叉积函数符号永远一致不会出现这边判断为内、那边判断为外的灵异现象。逆时针约定对后续核心操作还有一个关键帮助三角形 t 的顶点是 v[0]、v[1]、v[2]逆时针那么边 v[0]→v[1] 的左侧就是三角形内部而它对应的邻接三角形 t[0] 一定在右侧。翻边时我们经常需要判断两个相邻三角形是否共享一条边且分居两侧有了方向约定一句话就能写清楚。2.3 超级三角形的引入与收尾清理逐点插入法处理的一个前提是新点必须落在某个现有三角形内部。但第一点插入时整个三角网还是空的无处可落。常规做法是准备一个大到能覆盖所有输入点的超级三角形super triangle让所有点都落在它内部。这个三角形不参与最终输出处理完所有点之后只需要删掉所有包含超级三角形顶点的三角形即可。超级三角形的尺寸选择有个小讲究不能只刚好包住输入点最好把范围扩大个几十倍否则靠近边界的点做外接圆判断时可能因为超级三角形顶点参与产生不期望的边界翻转。实测中我把包围盒扩大 20 倍边界稳定很多。这里补充一点如果你想让最终剖分结果只覆盖凸包直接删掉包含超级顶点的三角形就行如果想保留外部边界形成更大的凸包区域需要另做边界恢复处理这个不在本文范围内。3. 增量插入的完整链路点定位、拆三角、逐层翻边3.1 点定位别急着上高级结构每次插入一个新点第一步是找到它落在哪个三角形内部。最朴素的做法是遍历所有三角形用重心坐标或方向测试判断点是否在三角形内复杂度 O(n)点少时完全够用。但当点量超过 2 万时线性扫描会成为明显瓶颈。这里有一个更聪明的办法步行法jump-and-walk。我们从上一个插入点所在的三角形出发沿边逐步向新点方向移动判断点在当前三角形的哪条边的外侧就跨到那个边的邻接三角形重复这个过程直到点落在当前三角形内部。步行法在均匀随机点下平均移动步数很少实测每个点定位大约只需要 5~8 次三角形跳跃比全量扫描快一两个数量级。实现也不复杂给定三角形三个顶点 a、b、c判断点 p 是否在边 ab 左侧如果在右侧说明要继续跨越到该边的邻接三角形。循环直到三条边判断都在左侧为止。注意处理点在边上的退化情况容差范围放小一点不要用绝对的 0 判断否则浮点误差会让程序陷入死循环。3.2 拆分三角形一次插入产生三个候选找到新点 p 落在三角形 T 内之后核心操作是把 T 一拆三T 的三个顶点分别与 p 连线生成三个子三角形。原三角形从数组里移除或标记失效三个新三角形加入数组同时维护好六条邻接关系每个新三角形的三条边两条是内部的、一条是原 T 的外边新的外邻接继承自原 T。这里要特别小心的是邻接表的正确性。我实现的时候用一个待检查队列保存新生成的三个子三角形因为新拆出来的三角形其外接圆性质并不一定满足 Delaunay 条件——它是从一个大三角形里切出来的新点出现在边上很可能落入相邻三角形的外接圆内。这种破坏只可能发生在共享边的两个三角形之间所以接下来只需要检查与新点相对的那些邻接边就够了。3.3 翻边恢复的递归逻辑incircle 判断与算法证明直觉翻边是整套算法的灵魂。设新点 p 所在的某个子三角形为 T1它有一条边与邻居 T2 共享T2 中与 p 相对的顶点是 q。如果 q 落在 T1 的外接圆内部说明这条共享边不满足 Delaunay 条件需要翻转把共享边的两个端点保持不变删掉这条边改连 p 和 q。翻转前T1(a,b,p) T2(b,a,q) 共享边 ab 翻转后T1(a,p,q) T2(b,q,p) 共享边 pq翻完之后产生两个新三角形它们各自可能又破坏了与更外层邻接三角形的关系所以递归地对这两个新三角形的外侧边再次做判断直到没有需要翻转的边为止。这一步对应 Lawson 算法中legalize edge的递归过程。它之所以保证能终止是因为每次翻转都让两个三角形的外接圆半径平方和严格减小整个系统存在一个单调递减的势能不会无限翻下去。4. 核心代码落地从数学判断到可运行的 C4.1 坐标结构、外接圆行列式与浮点容差先放下数学公式直接给一个实用的 incircle 判断函数。判断点 p 是否落在三角形 (a,b,c) 外接圆内部可以使用行列式它比先求圆心再算距离效率高而且避免了圆心坐标可能出现的极端值// a,b,c 为逆时针返回 epsilon 表示 p 在圆内 inline double incircleDet(const Vec2 a, const Vec2 b, const Vec2 c, const Vec2 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 ab ax * ax ay * ay; double bc bx * bx by * by; double ca cx * cx cy * cy; return ax * (by * ca - cy * bc) - bx * (ay * ca - cy * ab) cx * (ay * bc - by * ab); }这个函数的精度直接决定算法稳定性。浮点数不是实数坐标值如果量级很大比如经纬度坐标行列式结果可能本身很大把微小的几何差异淹没在误差里。我这里的经验是先对坐标做一次归一化平移缩放让所有点落入 [0, 1] 区间再参与计算。同时在函数外层加一个 epsilon 比较不要直接跟 0 比较见我后面第 6 章的踩坑内容。4.2 插入主函数的骨架有了数据结构、点定位函数和 incircle 判断插入一个点的主流程写起来就非常短了void insertPoint(int newVertex, int startTriangle) { int tri locatePoint(newVertex, startTriangle); // 返回包含新点的三角形 if (tri -1) return; // 退化点在边或已存在 int oldV[3]; // 记录原三角形三个顶点 int oldT[3]; // 记录原三角形三个邻接 int nt splitTriangle(tri, newVertex, oldV, oldT); // 一拆三返回新三角形起始索引 // 处理邻接关系让三个新三角形的外侧邻居正确指向 oldT linkAfterSplit(nt, oldV, oldT); // 对与新点相对的边做翻边恢复 legalizeEdge(nt, 0, newVertex); legalizeEdge(nt, 1, newVertex); legalizeEdge(nt, 2, newVertex); }legalizeEdge 的递归实现里比较关键的是正确找到共享边的相对顶点。给定当前三角形 cur它的第 i 条边是 (v[(i1)%3], v[(i2)%3])邻接三角形是 t[i]。在邻接三角形里找到与这条边相对的那个顶点本质上是在邻居的三个顶点里找出不等于共享边两端点的那个。这个过程用顶点索引比较即可注意别用坐标浮点比较索引比较又快又稳。4.3 边界条件重复点、共线点与退化输入代码能跑通简单用例之后真正的考验是退化输入。我随机生成数据时遇到最多的问题有三个重复点两个顶点完全相同会导致三角形面积趋近于 0外接圆判断结果不稳定。处理方式是在初始化时用哈希集合去重。共线点三个顶点共线时外接圆半径无穷大incircle 行列式结果接近 0。工程上我把接近 0统一当作不在圆内处理避免产生面积为零的退化三角形。点在三角形边上点定位时点恰好落在边上拆分会得到面积为零的三角形。我的做法是当点在边上时把点归入边一侧的较小三角形同时额外做一次边交换避免退化。这些边界情况每个看着都很小但任何一处没处理好程序跑几万个点大概率会在半路崩掉或生成破洞网格。我建议在实现阶段就准备三个独立测试集随机均匀点、格子点大量共线共圆、手工特例四点共圆、重复点混合每完成一个阶段就全量跑一遍回归。5. 性能实测从 O(n^2) 到百万点可用的两次优化5.1 理论复杂度与实测对照先给个直观参考我的实现最初只用线性扫描点定位插入 1 万点耗时约 0.6 秒5 万点飙升到约 18 秒呈现明显二次增长。这是因为每插入一点都要扫描全部现有三角形总复杂度 O(n^2)。虽然小数据看不出问题但数据量翻五倍耗时翻三十倍这就是二次复杂度的可怕之处。换成步行法点定位之后均匀随机点下的跳跃次数期望是常数级别整体复杂度接近 O(n log n)。同样的 5 万点耗时降到 0.9 秒。再用桶网格bucket grid记录每个网格单元内包含的三角形列表插入时先定位到网格单元再从单元内的三角形出发做步行进一步减少跳跃路径长度10 万点也能稳定在 2 秒内完成。5.2 桶网格点定位的工程实现桶网格的思路特别朴素把包围盒分成若干个格子每个格子维护一个与之相交的三角形列表。插入点 p 时先算 p 落在哪个格子里取该格子里的三角形作为步行法起点。如果格子为空或格子里的三角形已经被删除就往邻居格子扩散找。实现时注意三点一是格子数量不要贪多大致取 sqrt(三角形数量) 即可二是三角形可能跨越多个格子需要在插入时把所有相交格子都加入该三角形索引三是动态删除的三角形要惰性清理——遍历时发现无效三角形直接跳过不要急着从 vector 里 erase否则会破坏索引连续性。我踩过的一个坑是每次翻边后忘记更新桶索引导致定位命中已被删除的三角形程序看似正常但偶尔输出错误网格。所以我在代码里每次合法化边时只更新三角形数组本身桶索引采用过期检测策略遍历时验证tris[tri].valid标记简单可靠。5.3 C 层面的细节优化索引、预分配、缓存友好vector 预分配插入前根据输入数量tris.reserve(initialGuess * 4)大幅减少扩容时的拷贝。索引而非对象拷贝三角形数组元素用固定大小的结构体三个 int 加三个 int共 24 字节拷贝成本很低但频繁扩容时仍有开销预分配能完全回避。局部性优化翻边只操作局部三角形它们大概率位于 array 中相近位置缓存命中率比全图遍历好得多。这也是索引式结构的隐形红利。避免动态分配递归翻边用循环加显式栈std::vector 当栈用实现避免深递归带来的栈溢出和函数调用开销。实测中这些优化叠加起来10 万随机点的总耗时约 1.8 秒50 万点约 11 秒。对于纯 CPU 单线程实现这个数据已经能覆盖大多数离线网格生成场景。6. 调试与踩坑实录可视化是我最快的排错手段6.1 为什么看代码看不出 bug几何算法有个特点很多 bug 不会让程序崩溃而是让输出看起来差不多。比如三角形穿过了其他边、某个点的邻接三角形数量不对、网格中出现了细长缝隙。这种错误靠肉眼检查代码很难发现因为拓扑错误往往是运行时状态累积出来的跟某一行的逻辑对错没有直接关系。我第一次跑通算法后自信满满地输出坐标用文本方式检查了几个三角形觉得没问题。结果把数据导到可视化软件里一看网格中间有一条贯穿始终的裂纹还有一个点被五个三角形包围而不是应有的六个。这种问题如果只盯代码再盯三天也未必能看出来。可视化不是可选项是几何算法调试的必需品。6.2 三步走OBJ 导出—MeshLab/OpenCV 渲染—逐帧回放我的调试流程分三步。第一步把三角网的三角形数组导出为 OBJ 格式文件极简v x y z f i j kOBJ 格式在 MeshLab、Blender、甚至 Mac 自带预览里都能打开几行代码就能生成十分钟之内搭好可视化通道。第二步用 OpenCV 写一个几十行的实时渲染窗口把三角形逐帧画出来这样不仅能看最终结果还能看到整个插入过程中网格的演化。第三步是关键——逐帧回放。我会在代码里加一个调试开关每插入一个点、每次翻边后暂停一下把当前网格渲染出来。一旦发现某个时刻出现了非法拓扑比如三角形面积骤变成负数、邻接关系不对称立刻就能定位到是哪一步引入的错误。做过三次这种逐帧排查之后我对翻边逻辑的掌控感提升了一个层级比自己对着调试器打断点高效得多。6.3 浮点误差导致的幽灵翻边与对策最后聊一个隐蔽的坑浮点误差引发的幽灵翻边。incircle 计算里的行列式本质是几个大数的加减当两个点距离非常近时行列式的结果会很小与浮点误差同量级。这时若直接拿结果与 0 比较偶尔会误判成点在圆内触发一次多余的翻边。多余翻边的直接后果是网格仍然合法翻边自带局部恢复能力但会浪费一点时间。真正危险的是反过来当点已在圆内却因为误差判成不在圆内漏掉必需翻边最终产生违反空外接圆准则的三角形网格质量下降。我的对策分三层坐标归一化到 [0,1] 区间消除量级差异在 incircle 判断外层引入const double EPS 1e-12只有行列式大于 EPS 才认为在圆内使用 long double 做中间计算在极端退化数据下把误差再压低几个数量级。性能测试里 long double 慢约 20%但换来的是对怪异输入的鲁棒性。对于离线网格生成场景这个代价完全可以接受。还有一个小技巧定期校验网格的欧拉公式V - E F 1没有边界洞的连通三角网作为运行时断言。任何一个三角形拓扑被破坏这个公式立刻失衡比我盯着渲染画面找缝隙快得多。我后来把这段校验放在每个测试用例收尾时自动执行此后几乎没有看起来对但实际错的尴尬情况。DelaunayTriangulation 做到这一步已经从一个玩具演变成我手头可依赖的几何工具。如果你也在写类似算法我的建议是先把逆时针约定和索引式结构定死再实现点定位和翻边最后把可视化调试尽早搭起来——这个顺序能帮你避开我在前两周踩过的绝大多数坑。本文还有配套的精品资源点击获取