简介:一份面向几何计算与计算机图形学中高级学习者的三角剖分算法实现资源,基于分治策略完成Delaunay三角剖分,覆盖种子点选择、递归划分、边界处理与局部优化等经典步骤,展示复杂多边形如何逐步分解为互不重叠的三角形。压缩包内共91个文件,包括C++源码、头文件、Visual Studio工程配置、可执行程序、调试记录与界面资源等,整体约29.14MB,工程结构完整,可直接编译运行观察算法行为。已有591人学习浏览。通过对比源码与调试信息,读者可掌握分治法的递归拆分与合并思路,理解内切圆判断、邻接表维护、优先队列管理等关键细节;工程同时保留Visual Studio升级报告和调试目录,便于追踪构建过程,为渲染、物理模拟、网格生成等下游应用提供扎实基础。
1. 三角剖分用分治法实现:为什么切两半再缝合比逐点插入快
拿到几千个散点,第一反应是把它们连成三角形网格——有限元前处理、点云表面重建、地形建模都要走这一步。最常见的做法是逐点插入法,简单但点一多就慢,几万个点交互时就明显卡顿,插完还要反复做局部优化。这份资源用的是分治法三角剖分:先把点集按坐标排序,递归切成小块,小块内部直接构造三角形,再自底向上缝合,最终得到满足Delaunay条件的完整三角网。它把“切”和“缝”两个阶段拆得很干净,适合做网格前处理、计算几何入门验证,也适合给逐点插入法做对比评测。推荐给经常处理点集、被三角化性能卡住的人。
2. 分治三角剖分的拆解逻辑:排序、切分、叶子阈值怎么定
分治法听着高深,核心思想其实就两句话:把点集分成两半,各自生成完整的三角剖分,再用一条一条边把两个子三角网缝合起来。但“切”和“缝合”各有各的讲究,切不好会交叉,缝不好会翻车。这一章先把三条核心设计决策讲透,后面写代码时你才知道每行在干什么。
2.1 分治 vs 逐点插入:数据量、初始形态与最坏复杂度
逐点插入是多数人最早接触的Delaunay三角剖分实现。它一次插入一个点,找到包含该点的三角形,拆成三个或四个三角形,然后对邻居边做翻转恢复Delaunay性质。这个流程对几千点很友好,但最坏时间复杂度是O(N^2),平均也在O(N^1.5)左右。点分布越均匀,逐点插入表现越好;点分布越差,定位包含三角形的开销越大。
分治法走的是另一条路:排序后按中点对切,递归求解子问题,最后合并。期望时间复杂度是O(N log N),而且不受插入顺序影响,结果完全确定。代价是合并阶段比逐点插入的局部优化复杂得多,需要处理两个子凸壳之间的缝合。
选哪个,数据规模是第一判据。少于五千点,逐点插入可能更快,因为代码短、开销低、不容易出边界问题;超过五万点或者需要反复三角化不同点集,分治法的优势就很明显。如果你要处理的是百万点云,那基本不用犹豫,直接上分治法。
2.2 预处理排序:x坐标主序、y坐标次之,重复点必须去重
分治法的第一步永远是排序。排序规则直接决定切分是否均匀,也决定合并阶段能不能用简单的索引区间来表示子集。
常见做法是先按x坐标排序,x相同再按y坐标排序,形成字典序。排序之后再从中间切一刀,左半部分和右半部分的x坐标天然分开,缝合时不会有交叉。这里的坑是:如果只按x排序,x重复率高的数据会让切割线附近点集混乱,递归两边的点在x上重叠,合并时你会发现有三角形越过分割线,怎么调都调不干净。
去重同样不能省。两个坐标相差微小的点会被当成两个顶点,生成面积接近零的三角形,合并时又会被翻转逻辑删掉,浪费大量时间。我一般会在排序后做一遍线性去重:相邻两个点距离小于EPS就合并,保留先出现的点。EPS取值根据坐标量级来,屏幕坐标用1e-6到1e-7,浮点大地坐标用1e-9甚至更小,不要盲目用一个固定值。
2.3 递归终止阈值:为什么常见实现把LEAF_THRESHOLD设成32
递归不能一路切到底。切到只剩三个点时虽然可以直接生成一个三角形,但合并次数会非常多——每N/3个叶子都要和相邻叶子做一次凸壳缝合,缝合本身又有指针维护、切线查找、边翻转等开销,总体反而更慢。
所以工程实现里都会设一个叶子阈值LEAF_THRESHOLD,点数小于等于这个值就停止递归,改用暴力三角化。常见取32或64,也有实现取16的。块越小,每个叶子内部的构建越快,但合并次数越多;块越大,叶子内部构建越慢,但合并次数减少。我曾在一组十万点数据上测过,从16调到32,总耗时下降了百分之十五左右,再从32调到64又反而回升了,因为叶子内部暴力三角化的O(k^2)开始拖后腿。
叶子内部三角化怎么打底?做法是:先对这小块点集求凸壳,把凸壳顶点按顺序连成边,然后逐个插入内部点,用局部优化恢复Delaunay性质。这一步不用分治,直接暴力就行,因为块内点数少,怎么折腾都快。
3. 数据结构与递归主流程:边表、LUT查找表与divide_conquer实现
有了第2章的三个决策,这一章进入代码层面。数据结构的设计会直接影响合并代码的复杂度和内存占用,尤其是边的表示和查找效率。下面用Python给出核心结构,语言不同但思路通用。
3.1 点、边、三角形的对象设计:索引比坐标更省内存
点在排序时会被反复交换位置,如果点对象里存完整坐标,交换时拷贝面积会比较大。常见做法是点对象里放坐标和原始索引,排序时按坐标比较,但保留idx字段用于唯一标识。
class Point: __slots__ = ("x", "y", "idx") def __init__(self, x, y, idx): self.x = x self.y = y self.idx = idx边的表示要注意方向。一条边在左右两侧各有一个邻接三角形,如果存成有向边,每次查邻接都要处理方向问题;如果存成无向边,需要保证Edge(a, b)和Edge(b, a)哈希结果相同。最省心的写法是构造时把idx小的放前面。
class Edge: __slots__ = ("a", "b") def __init__(self, a, b): self.a = a if a.idx < b.idx else b self.b = b if a.idx < b.idx else a def __hash__(self): return hash((self.a.idx, self.b.idx)) def __eq__(self, other): return (self.a.idx, self.b.idx) == (other.a.idx, other.b.idx)逻辑说明:__slots__让每个对象减少约一半内存占用,处理十万点集时差别非常明显,这一点值得习惯性地用起来。Edge构造时把idx小的存为a,保证两个方向创建的边是同一个对象,合并阶段不需要反复确认“边的方向反没反”。
三角形则保存三个顶点和三边。顶点顺序必须逆时针,这是后面所有翻转和inCircle判断的前提,如果某处生成了顺时针三角形,Delaunay性质校验会全部错乱。
3.2 divide_conquer主流程:切分、递归与返回出口
排序完成后,整个分治结构就是一个深度为log2(N)的递归树。递归函数只需要接收全局点集数组和左右索引区间,不要传切片,避免每层都复制一份点集。
def divide_and_conquer(points, lo, hi): """points 是已按 (x, y) 字典序排序的全局数组 lo, hi 是左闭右开区间,表示当前要处理 points[lo:hi]""" if hi - lo <= LEAF_THRESHOLD: return triangulate_leaf(points, lo, hi) mid = (lo + hi) // 2 left_mesh = divide_and_conquer(points, lo, mid) right_mesh = divide_and_conquer(points, mid, hi) return merge_two_meshes(left_mesh, right_mesh, points)逻辑说明:mid由区间长度二分得到,左右子问题点数差距不超过1,递归深度稳定在 log2(N) 附近。triangulate_leaf负责在第2章说的叶子内部构建三角网,返回的是一个边列表和一个“边到邻接三角形”的映射。merge_two_meshes是分治法的重头戏,下一章专门展开。
参数说明:LEAF_THRESHOLD是模块级常量,初始建议32,调优时按第2章的经验搜索。返回的边列表是最终的Delaunay边集合,渲染或者做有限元输出时直接消费这个列表即可。整个递归过程不允许修改points的排序顺序,否则上层递归拿到的索引区间就失效了,这一点很容易被忽略。
3.3 LUT查找表:合并时的边邻接查询
合并两个子三角网时需要频繁回答一个问题:某条边已经被哪个三角形使用了?每次还去遍历整个边列表,时间复杂度直接爆炸。常见做法是维护一个边的查找表(Lookup Table,LUT),把边映射到它的邻接三角形集合。
from collections import defaultdict def build_lut(triangles): lut = defaultdict(set) for tri in triangles: for e in (tri.e0, tri.e1, tri.e2): lut[e].add(tri) return lut逻辑说明:一条边在完整的三角网里最多被两个三角形共享。合并过程中查LUT就能快速知道某条边是内部边还是边界边,以及左右邻居分别是谁。某条边在LUT里只有一个邻居,说明它在子网格的凸壳边界上,是缝合时要关心的候选边。
参数说明:lut的键是Edge对象,值是一个集合。集合里的三角形对象要支持哈希,所以三角形类要定义__hash__,或者用三角形内部顶点的元组做标识。如果不做哈希,合并时频繁判断“这个三角形是否已存在”会退化成线性匹配,整个缝合过程会被拖成O(N^2)。
4. 合并两个子三角网:LR-EDGE扫描与inCircle判定的配合
合并阶段是分治法的核心,也是所有实现里最容易写错的地方。它的目标是:给定左右两个各自满足Delaunay条件的三角网,把它们缝合在一起,缝合过程中产生的所有新三角形也必须满足Delaunay条件。这一章按实际推进顺序拆开讲。
4.1 先从两个凸壳找下切线和上切线
两个子三角网分别有自己的凸壳边界。合并的第一步,是把这两个凸壳的“下公切线”和“上公切线”找出来。下公切线是同时与左右两个凸壳相切、且位于两个凸壳下方的那条线段,上公切线同理。
下公切线就是合并后的第一条LR-EDGE。LR-EDGE这个名字直观描述了它的两个端点:L端在左子网格,R端在右子网格。缝合过程从下公切线开始,逐条向上推进,最终结束于上公切线。
找切线时有几个细节要处理。左右凸壳都以逆时针顺序存储顶点,下切线的算法是:先随便取左凸壳最右点、右凸壳最左点,然后分别判断该点是否能继续“向下”移动,直到两边都无法移动,得到的线段就是下公切线。这个过程的时间复杂度是O(h)的,h是凸壳顶点数。
4.2 LR-EDGE推进:左右候选点只保留一个
有了第一条LR-EDGE之后,问题就变成:以这条边为底边,找一个第三个顶点,构造一个新的Delaunay三角形。候选点来自这条边上方的左右两侧,每侧各选一个最优候选,然后通过空外接圆判定决定谁胜出。
下面是一个简化但结构完整的合并循环:
def merge_two_meshes(left, right, points): left_hull = convex_hull(left) right_hull = convex_hull(right) base = lower_common_tangent(left_hull, right_hull) upper = upper_common_tangent(left_hull, right_hull) result = left.triangles + right.triangles lut = build_lut(result) while (base.a.idx, base.b.idx) != (upper.a.idx, upper.b.idx): cand_l = best_candidate_left(base, left, lut) cand_r = best_candidate_right(base, right, lut) if (cand_l is not None and cand_r is not None and incircle(base.a, base.b, cand_l, cand_r) > 0): winner = cand_r elif cand_l is not None: winner = cand_l else: winner = cand_r tri = make_triangle(base.a, base.b, winner) add_to_result(result, lut, tri) base = next_base_edge(tri, base) return result逻辑说明:best_candidate_left不是把所有点遍历一遍,而是从LUT里找与base.a相连且位于基线上方的点,再从中选一个使角度最合适的候选。右侧同理。当左右两侧都有候选时,用incircle判定:如果左侧候选点落在右侧候选点与基线构成的外接圆内,说明右侧候选点更优,选右侧,反之选左侧。被选中的候选点与基线构成新三角形,然后基线更新为新三角形靠上的那条边,循环继续。
参数说明:base.a必须来自左子网格,base.b来自右子网格。方向约定颠倒,切线和推进方向会全部反转,轻则死循环,重则生成交错三角形。我在初版实现里犯过这个错,调试了很久才发现切线方向写反了,从那以后我都在代码注释里强制标注“base.a=left,base.b=right”。
4.3 inCircle判定的精确实现:方向与容差缺一不可
inCircle判定是整个合并过程的“法官”。两个候选点竞争时,谁能在当前三角形外接圆内把另一个“挤走”,靠的就是它。inCircle的标准数学形式是4×4行列式,但对二维点集可以降维成3×3行列式,计算量更小。
def orient2d(a, b, c): return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x) def incircle(a, b, c, d): # a, b, c 必须逆时针;若顺时针则交换前两点 if orient2d(a, b, c) < 0: a, b = b, a ax = a.x - d.x; ay = a.y - d.y bx = b.x - d.x; by = b.y - d.y cx = c.x - d.x; cy = c.y - d.y det = ( ax * (by * (cx * cx + cy * cy) - (bx * bx + by * by) * cy) - ay * (bx * (cx * cx + cy * cy) - (bx * bx + by * by) * cx) + (ax * ax + ay * ay) * (bx * cy - by * cx) ) return det > 1e-12逻辑说明:det > 0表示点d在a、b、c的外接圆内部,det < 0表示在外部,det接近0表示四点近似共圆。共圆情况下谁胜出都可以,但结果会不稳定,所以要通过容差将共圆归到某一侧,保证算法有确定行为。
参数说明:1e-12是一个经验容差,适用于坐标量级在1e3到1e5之间的普通数据。如果你的点坐标是米制经纬度或者超大坐标,这个固定容差会失效,需要改成相对容差,比如先求所有点的最大坐标范围scale,再比较det > 1e-12 * scale ** 4。另外,orient2d的方向统一不能省,逆时针是inCircle判定的前置条件,三角形类构造时就该保证,但为了保险,inCircle内部再做一次判断也不亏。
5. 避坑与常见问题:分治三角剖分最容易翻车的四个环节
分治三角剖分写起来比逐点插入难,难点全在边界条件和极端数据上。这一章整理四个高频踩坑点,都是我实际调试中遇到过的,按“现象、原因、解决”三步说清。
5.1 排序规则不统一:两侧三角网在分割面上交叉
现象:分治递归跑完后,合并阶段生成的三角形横跨左右子网格的分割线,从整体上看三角形相互交叉,面积校验直接失败。
原因:递归排序时用了按x坐标排序然后二分,但在叶子三角化或者其他步骤里,代码又按某种别的规则把点重新组织过,比如先按y再按x,导致左子集的点出现在右子集的索引区间里。索引区间和空间位置对不上,合并时LR-EDGE的L端和R端实际位置就没法保证左右关系。
解决:全流程只允许一次排序,排序键统一为 (x, y)。排序之后所有子问题都用索引区间表示,任何步骤都不要再对点数组做重排。如果数据本身带有严格的空间分层,可以考虑按莫顿码排序,但前提是整个递归都使用同一个排序结果,不要中途换规则。
5.2 共线点导致面积为零的退化三角形
现象:叶子三角化时出现面积刚好为零的三角形,边表里出现重复边,合并时next_base_edge跳来跳去回到旧基线,程序卡死。
原因:叶子块内有三点近似共线,暴力三角化把这三点强行连成了三角形,但它的面积为零,在后续翻转和缝合判定里被当成“非法三角形”剔除,边关系就断裂了。
解决:叶子三角化之前做一次退化检测。三点共线时只保留最长边和两条短边,不生成三角形。如果是重复点,在排序后去重时直接去掉。关键是检测标准不要只用orient2d == 0判断,要用绝对值小于容差判断,否则浮点坐标下永远检测不到共线。
5.3 inCircle浮点判定翻车,Delaunay性质被破坏
现象:合并完成后检查所有内部边,部分边对侧的顶点落在外接圆内,理论上不该出现但确实出现了,且每次都出现在坐标数值很大的区域。
原因:inCircle行列式的计算涉及坐标平方差,坐标值达到1e6量级时,行列式结果可能放大到1e12量级,固定容差1e-12完全失效,导致接近共圆的点被错误判定。
解决:把固定容差换成相对容差。常见做法是先算数据包围盒的尺度scale,判定条件改成det > 1e-12 * scale ** 4;或者干脆引入自适应精度判定,先算float结果,结果接近零时再用高精度重算。后者是生产级库的标准做法,如果你只是教学演示,相对容差已经够稳。
5.4 递归深度失控:十万点级数据集直接StackOverflow
现象:几千点跑得好好的,一换成十万点就栈溢出,或者十几万点时要跑几十秒,明显不是O(N log N)的表现。
原因:递归没有真正按中点二分。有些实现为了“按坐标空间分割”,先计算点集包围盒,再按包围盒的x方向中值切分,这在点分布极度不均时会让一侧子问题包含大部分点,递归深度变成O(N)而不是O(log2(N))。
解决:排序之后直接按索引mid = (lo + hi) // 2切分,不要按几何中值切分。即使点集中有大量相同x坐标,按字典序排序后索引区间依然是均衡的。如果担心递归栈本身的开销,可以设置最大深度计数,超过log2(N) + 8强制把当前子区间按暴力方式处理,避免极端数据导致栈溢出。
6. 验证与进阶:面积守恒、Voronoi对偶和约束扩展
代码写完只是第一步,验证分治结果是否正确才是真正花时间的地方。这一章给三个实用方案:第一个用欧拉公式和面积守恒做快速自检,第二个用Voronoi对偶做结构性验证,第三个说清楚怎么往分治框架里加入约束边。
6.1 用欧拉公式与面积守恒做自动化体检
对于N个点的Delaunay三角网,设凸壳点数为h,三角形数量T和边数量E有关系式:
T = 2N - 2 - h,E = 3N - 3 - h
这个关系在每次合并完成后都成立。检查脚本里直接数一数三角形和边数,对比公式是否一致。然后再把每个三角形面积累加,与凸壳多边形面积做差,相对误差应小于1e-9。这组检查不依赖具体坐标方向和容差,是最稳定的“体检项”。
6.2 用对偶Voronoi图做第二轮交叉验证
Delaunay三角剖分与Voronoi图互为对偶。每个三角形的外心对应Voronoi顶点,每两个相邻三角形的外心连线对应Voronoi边。若三角剖分正确,则这条外心连线应垂直平分其公共Delaunay边。
检查方法:遍历每条内部Delaunay边,取出左右两个三角形外心,计算外心连线与公共边的点积是否为0,同时检查外心连线中点到公共边两端点距离是否相等。误差控制在1e-6以内。我一般用这个做最终验收,它比面积守恒更能暴露局部缝合错误。
6.3 向带约束的Delaunay三角剖分扩展
分治框架天然可以扩展为带约束的版本。约束边(比如地形断裂线)在叶子三角化时就要被标记为“不可翻转”。合并阶段,如果LR-EDGE与约束边相交,不能直接缝合,需要把相交区域拆开,以约束边为新的边界重新构造三角形。具体做法是:在merge循环里增加一个约束边相交判断,相交时暂停LR-EDGE推进,先处理约束边两侧的三角形,再继续推进。
处理约束边之后,结果就从标准Delaunay变成了约束Delaunay三角剖分(CDT),可以直接接有限元网格生成或道路建模流程。这一步工作量不小,但好处是复用当前的排序、叶子三角化和LUT结构,只需要在merge里加一个分支。
从那以后我每次拿到三角剖分包,都强制先跑一遍欧拉公式、面积守恒和Voronoi对偶三重校验,跑过了再谈性能调优。分治法的坑都在缝合阶段,校验脚本就是我的后悔药。希望帮到你。
本文还有配套的精品资源,点击获取