简介这是一份以分治法实现三角剖分的完整C工程专注Delaunay三角剖分算法适合计算机图形学、几何计算与科学计算领域的初学者和开发者参考。包内包含源码、头文件、Visual Studio工程配置及可执行程序共91个文件以h/cpp源码、obj中间文件、tlog构建日志、资源脚本等类型为主压缩包约29.14MB目录结构清晰既可整体编译运行也可单独抽取核心模块研读。算法讲解覆盖种子点选取、递归划分、边界处理与局部优化等关键步骤并借助邻接表、优先队列等数据结构提升效率可帮助读者理解分治思想在复杂几何剖分中的落地方式同时熟悉一套完整的C工程组织与构建流程。已有579人学习适合用于课程设计、算法实验或项目二次开发对提升几何算法实战能力很有帮助。 最近在重构一个点云网格化的小工具又一次跟三角剖分算法面对面。说起这个领域大多数图形学背景的人第一反应都是 Bowyer-Watson 或者 Lawson 逐点插入因为代码量小、思路直白一个点一个点往图里塞不合法就翻转边。但我处理的数据经常是十几万量级的扫描点而且分布很不均匀逐点插入在点集空洞很多的区域会反复做局部修复性能波动非常大。分治法三角剖分恰好是另一个极端先排序再按 x 坐标中位线切成左右子集递归剖完最后做一次从下到上的合并扫描。它没有增量法那种“修到哪算哪”的随机性整体复杂度稳稳落在 O(n log n)。这篇不是把教科书搬过来我想记录的是从选型、递归拆分、merge 合并到调试自检的完整路径。1. 先别急着写代码三角剖分到底在求什么1.1 点集三角剖分和 Delaunay 约束三角剖分这个说法在不同领域里指的东西不完全一样。图形学里常见的多边形三角剖分是把一个多边形内部切成不重叠的三角形而计算几何里更常讨论的是点集三角剖分给定平面上 n 个点用互不相交的边把它们连起来使得所有三角形的并集恰好覆盖这些点的凸包并且三角形的顶点只能是给定点。满足这个条件的剖分方案非常多点的数量一大方案数量会爆炸式增长。Delaunay 三角剖分就是在所有合法剖分里加一条约束每个三角形的外接圆内部不允许出现其他点也就是“空圆”性质。这个约束看起来简单实际上会让整个剖分同时满足很多好性质最小角最大化、对偶到 Voronoi 图、局部修改只涉及相邻三角形。做地形网格、路径网格、有限元网格时大家都倾向用它就是因为三角形形状更“饱满”不会出现被拉得极其细长的劣质三角形。1.2 分治法的整体思路分治法的名字已经说明了一半把大问题切成两个规模接近的子问题各自解决再把结果合并。三角剖分里能做到这一点很大程度上依赖 Delaunay 三角剖分的一个局部性直觉——两个子剖分之间需要新增的边只集中在一条被称为“缝合线”的狭长区域里合并时从底部一路往上扫每条新边确定后只需要做非常局部的翻转修复。合并是整个算法最核心、也最容易写错的地方。我在第一次实现时反复在这一步翻车后来才理解merge 看似是在加边其实每一步都在回答一个问题当前这条底边已经确定下一个应当与它组成三角形的顶点究竟是左边子剖分的某个边界顶点还是右边子剖分的某个边界顶点只要这个选择逻辑是对的边翻转只是在收拾残局。2. 递归拆分与基线情形2.1 第一步永远是排序分治法能稳定拿到 O(n log n) 的前提是先把点按 x 坐标做字典序排序。注意不是只排一次就完事而是整个递归过程都用这份排好序的数组下标来切分。每次递归取mid (l r) // 2左半边是原数组里前半段右半边是后半段。因为已经按 x 排好序左半边的点天然都在右半边的点左侧至少 x 坐标不会出现大范围穿插这会让合并时的公共切线寻找变得非常简单。这里有一个工程细节递归函数尽量只传数组引用和[l, r)区间不要在每一层都copy出一个子数组。我之前为了图方便在递归入口用切片生成新数组结果数据量到十万左右时内存和拷贝开销直接让 O(n log n) 名存实亡。真正跑起来之后所有子剖分共享同一块排好序的存储索引切割就够了。2.2 一个点、两个点、三个点怎么返回基线情形通常有三种一个点返回一个不包含任何边的空结构。两个点直接返回一条边。三个点返回一个三角形但必须先做方向测试。如果三点共线就不能生成面积为零的三角形应当退化成两条边。如果不去处理共线三点后续 merge 阶段会凭空出现面积接近零的三角形然后 inCircle 测试的符号会变得极其不稳定整个剖分出乱子是迟早的事。我通常会在三个点的基线情形里加一个orient判断若面积为 0按一条折线返回若不为 0保证三个点按逆时针顺序存入邻接关系。def dc_triangulate(pts, l, r): n r - l if n 1: return empty_graph(pts[l]) if n 2: return single_edge(pts[l], pts[l 1]) if n 3: a, b, c pts[l], pts[l 1], pts[l 2] if abs(orient(a, b, c)) EPS: return two_edges(a, b, c) return triangle_in_ccw(a, b, c) mid (l r) // 2 left dc_triangulate(pts, l, mid) right dc_triangulate(pts, mid, r) return merge(left, right)2.3 用邻接结构撑起合并操作合并阶段需要频繁回答两类问题某个顶点的邻居有哪些在这些邻居里哪些落在底边的特定一侧。所以数据结构的选型会直接影响合并代码的复杂度。正式实现里当然可以上带方向的半边结构这对边翻转、遍历邻接边都很舒服。但如果你只是想把算法跑通、理解清楚用“顶点 - 有序邻居列表”的邻接表就够了。关键是有序每个顶点的邻居必须按极角从小到大维护这样才能在合并时快速判断“某个方向的区域里谁是最贴近底边的候选点”。翻转边的时候邻接表需要同步更新四条边的关系。这一步没有任何技巧就是仔细。我见过太多实现merge 逻辑写得不错最后因为翻转后邻接顺序忘记调整找候选点时漏掉关键顶点出现交叉边。3. 合并步骤是这场实现的主战场3.1 找下公共切线合并的起点两个子剖分各自是合法 Delaunay 三角剖分各自有一条凸包轮廓。合并的第一步是找到左凸包和右凸包之间的下公共切线。所谓下公共切线就是一条同时触及左右两个凸包、且所有凸包顶点都在这条线上方的直线。不能想当然地取“左子剖分最右顶点”和“右子剖分最左顶点”直接连起来那通常不是公共切线。需要从这两个初始点出发不断沿着凸包移动def lower_tangent(left, right, L, R): # L 是左剖分最右凸包顶点R 是右剖分最左凸包顶点 while True: moved False nl next_hull_vertex(left, L, clockwiseTrue) nr next_hull_vertex(right, R, counterclockwiseTrue) if orient(L, R, nl) 0: L nl moved True if orient(L, R, nr) 0: R nr moved True if not moved: break return L, R这里的orient(L, R, p) 0表示 p 落在有向线段 L-R 的下方需要把端点挪下去。旋转方向要注意左凸包沿顺时针方向推进右凸包沿逆时针方向推进因为两个子剖分的凸包朝向相对。如果方向搞反合并会进入死循环或者得到上公共切线。3.2 底边推进一侧一个候选再用外接圆定胜负找到下公共切线之后把它作为第一条底边加入最终结果。接下来整个合并循环都围绕这条底边展开。对于当前底边(L, R)左边候选是 L 的邻居中位于底边“右上方”的顶点右边候选是 R 的邻居中位于底边“左上方”的顶点。两边各自只需要保留最贴近底边的那一个不需要把所有候选都拿来两两比较。原因在于递归不变量已经保证了两个子剖分内部都是合法 Delaunay跨过缝合线的下一跳只可能在左右两个方向的“边界前沿”中产生。选出左候选lc和右候选rc后用inCircle(L, R, lc, rc)判断如果rc落在三角形(L, R, lc)的外接圆内部说明连L-lc会产生非法三角形进而应该优选rc否则选lc。这个判定方向非常容易记反我吃过亏。建议先把左右两个子剖分各放 5 个点的手算样例跑一遍再把合并结果可视化确认谁应该胜出然后回来固定这个符号约定。def choose_next(L, R, lc, rc): if lc is None and rc is None: return None if lc is None: return rc if rc is None: return lc # 注意符号与顶点顺序有关这里是按逆时针点序约定的版本 if in_circle(L, R, lc, rc): return rc return lc3.3 添加新边后的局部合法性修复新边确定后不能直接推进底边。因为新边的加入可能让缝合线周围的某些三角形不再满足空圆性质需要做一次局部合法化也就是 del 边翻转。翻转的思路是检查新边的两侧把共用边的两个三角形看成四边形如果另一个顶点落在当前外接圆内部就把公共边换成另一条对角线。这个操作在增量法里也有。分治合并里的区别是它只需要沿着新边往四周扩散检查而不是全图扫描。我习惯用一个栈来组织这些待检查边def legalize(edge, graph): stack [edge] while stack: e stack.pop() v opposite_vertex(e) # 与 e 组成四边形的另一个顶点 p opposite_vertex(e, otherTrue) if in_circle(e.p1, e.p2, v, p): flip(e) stack.append(new_edge_after_flip_1) stack.append(new_edge_after_flip_2)翻转之后原边被替换成两条对边的两条新对角边具体哪两条取决于四边形顶点顺序。所有顶点邻接关系要同步更新。3.4 merge 的完整伪代码把上面的步骤串起来merge 的骨架大概长这样def merge(left, right): L, R lower_tangent(left, right) add_edge(L, R) while True: lc best_left_candidate(L, R) rc best_right_candidate(L, R) v choose_next(L, R, lc, rc) if v is None: break if v is rc: new_edge (L, rc) else: new_edge (lc, R) add_edge(new_edge) legalize(new_edge, left, right) # 推进底边 if v is rc: R v else: L v return combine_graphs(left, right)第一次写的时候我很自然会忘记最后那两行“推进底边”。如果不更新 L 或 R循环就会一直在同一条底边上反复比较同样的候选最终栈溢出或者死循环。把这个更新动作放在 legalize 之后而不是之前是因为 legalize 需要知道新边的完整上下文推进之后再翻转容易把邻接关系弄乱。4. 复杂度分析与浮点精度两条硬约束4.1 O(n log n) 从哪来分治三角剖分的复杂度分析看起来很简单但真正理解它的来源才能在遇到诡异性能问题时定位。排序本身是 O(n log n)。递归树每一层所有子问题加起来要做多少次 merge 操作每次 merge 的循环次数与左右两个子剖分在缝合线附近扫描过的边数成正比而每一层的所有缝合线加起来只覆盖这一层的全部顶点一次因此每层的工作量是 O(n)。递归树深度是 O(log n)乘起来就是 O(n log n)。和逐点插入做对比会更明显逐点插入最坏情况下每插入一个点都可能引发一整条链的翻转修复累积起来可能到 O(n²)。实际数据里虽然很难遇到最坏情况但分布极不均匀时性能抖动很真实。分治法的优势是它把复杂度压在排序和分层合并里数据分布再难看也只是常数增大阶不会上涨。4.2 inCircle 的稳定实现合并的灵魂是 inCircle 测试它的数值稳定性直接决定剖分是否正确。我见过用atan2计算角度再比较的实现角度误差在长瘦三角形上会被放大最后得出错误的胜负判断。正确做法是用三阶行列式判断点是否在外接圆内。实现时把待测点平移到原点可以减少大数相减带来的精度损失def in_circle(a, b, c, d): ax, ay a.x - d.x, a.y - d.y bx, by b.x - d.x, b.y - d.y cx, cy c.x - d.x, c.y - d.y 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) return det 0这里的符号约定会随顶点是逆时针还是顺时针而整体翻转所以写完后务必用一组小样例标定。另一个经验是项目里如果坐标范围很大比如经纬度或者毫米级坐标建议在进入算法前先做一次归一化平移否则 x² y² 很容易把有效数字吃掉。4.3 共线点、重复点、递归深度这些边界问题共线点是 Delaunay 三角剖分里最让人头疼的输入类型。当一堆点严格落在同一条线上时Delaunay 三角剖分是没有意义的因为所有三角形面积都是零。实践中通常会先在输入阶段去重再对共线点做轻微扰动或者干脆当成特殊线段集处理不进入三角剖分主流程。重复点同样要提前去重。如果真有一模一样的点三角剖分定义的“顶点集合”就不成立了后续所有邻接关系都会错乱。我的习惯是在排序之后做一次线性扫描去重保留首次出现的点后续用到索引时会省掉很多心智负担。递归深度也要心里有数。十万量级的点在平衡切分下递归深度大约是二十层左右完全没有压力。但如果排序后点分布极度偏向某一个方向又采集中位数切分还是会有堆栈风险。更常见的隐患是每层切分复制子数组这在内存上比递归深度先爆炸。5. 调试阶段最容易翻车的三个现场5.1 合并完出现交叉边合并结果出现边相交绝大多数情况下不是几何判定的问题而是候选选择或邻接顺序的问题。我自己踩过的坑是best_left_candidate只筛选了落在右侧的邻居却没有保证“最贴近底边”。如果把一个离得很远的邻居也当成候选它会在和右候选的 inCircle 比较中获胜然后产生一条很长的新边把已有三角形区域穿过去形成交叉。排查时不要只盯着最终结果要在每次add_edge(new_edge)后立刻检查当前已有边是否相交。把这个检查写进测试断言里跑随机点集时只要第一次交叉冒出来就能精确定位到是哪一步引入了它。def has_crossing(edges): for i in range(len(edges)): for j in range(i 1, len(edges)): if segments_intersect(edges[i], edges[j]): return True return False这段 O(m²) 的自检只用于调试不要留在生产路径里。5.2 递归栈先于算法崩掉递归本身不是问题问题出在数据结构和递归之间互相放大。最典型的场景是用递归函数处理上百万个点每层还创建新的邻接表结果栈没崩内存先崩。后来我把邻接表改成预先分配的大数组顶点索引加边索引都在数组里复用只持有当前子剖分的边集合内存压力才下来。如果语言自带递归深度限制比如 Python需要在入口处设置更高的 recursionlimit或者直接用显式栈模拟递归。模拟递归会让代码难读不少但处理百万级点集时反而更可控因为你可以随时暂停、恢复、查看中间状态。5.3 用边数和 Delaunay 属性做自动自检写完 merge 之后我习惯在算法出口加一组自动断言。无重复点、无共线退化的情况下n 个点的 Delaunay 三角剖分边数 E 和三角形数 T 满足固定公式E 3n - 3 - hT 2n - 2 - h其中 h 是凸包顶点数。这个关系一旦不满足说明剖分过程中要么丢了边要么多加了边。更严格的自检是逐边验证空圆性质对每个内边找到边两侧的两个顶点检查其中一个是否落在另一个顶点形成的三角形外接圆内部。这个检查是 O(n) 量级放在测试脚本里很合适。def check_triangulation(pts, edges, triangles): n len(pts) h len(convex_hull(pts)) assert len(edges) 3 * n - 3 - h, edge count mismatch assert len(triangles) 2 * n - 2 - h, triangle count mismatch # 抽样验证空圆 for e in internal_edges(edges): a, b e c, d adjacent_vertices(a, b) assert not in_circle(a, b, c, d)我在开发阶段就是用这套断言配合随机点集一晚上抓出了三个“看起来正确但实际非法”的 bug。6. 什么时候该自己写什么时候干脆用库6.1 从分治三角剖分到受限三角剖分的扩展如果场景里需要强制保留某些边比如地形建模里的断裂线、道路边界、游戏寻路里的多边形障碍轮廓普通 Delaunay 三角剖分就不够了需要受限 Delaunay 三角剖分。分治法也能扩展出受限版本但复杂度会明显上升因为合并时公共切线、候选选择和翻转规则都要考虑约束边是否被破坏。如果只是做一次性的离线计算建议优先考虑现有库Shewchuk 的 Triangle 库是 C 语言的经典实现稳定性和鲁棒性都经过工业界多年检验Python 科研场景直接scipy.spatial.Delaunay也很省事。自己实现分治法的意义更多在于理解算法本质、定制特殊数据结构和嵌入实时渲染管线。6.2 我的选型建议我的经验是少于几千个点直接逐点插入法最舒服代码短、调试容易性能差异可以忽略。几万到几十万点不想引入复杂依赖分治法值得认真实现。上百万点、需要稳定性能和自定义内存布局时更推荐基于分治思想的成熟库或者自带迭代合并的高性能实现而不是自己从头写。6.3 收尾的调试小工具最后分享一个提升效率的习惯我把 merge 每完成一步的中间结果输出成 SVG 或 PLY 文件丢进网格查看器里回放。这种“看得到进度”的调试方式比纯打印边列表高效得多。尤其是底边推进的过程只要可视化里出现某一条新边跨过已有三角形问题几乎立刻就能定位。个人体会是分治三角剖分真正难的不是理解分治而是把合并那十几行循环的意思和邻接数据结构完全对齐。动手前先画一张三顶点、四顶点的手算图把每一步候选、inCircle 结果和翻转操作走一遍再上手写代码能少走我当初至少三天的弯路。本文还有配套的精品资源点击获取