分治三角剖分算法实战：从点集到网格，省掉两周试错

简介这份资源围绕三角剖分算法展开重点讲解如何用分治法实现 Delaunay 三角剖分面向计算机图形学、几何计算与科学计算方向的学习者和开发者适合已具备一定数据结构与算法基础、希望深入理解剖分原理与工程实现的人群。压缩包共 91 个文件约 29.14MB以 cpp 与 h 源码为核心配合 obj、pdb、ilk 等编译中间产物以及 sln、vcxproj、filters 等工程配置另有 tlog、log、idb 等调试与构建记录整体是一套可直接编译运行的完整工程。资源从种子点选择、递归划分、边界处理到优化环节逐步展开并涉及二叉堆、优先队列及邻接表等数据结构的使用有助于读者理解分治思想在几何算法中的落地方式。目前已有 590 人学习下载可作为算法课程设计、图形学实验或相关项目开发的参考素材。1. 三角剖分算法用分治法落地从点集到网格这份资源能省掉你两周试错如果你正在做有限元前处理、地形建模、散点插值或者游戏里的导航网格生成大概率绕不开三角剖分。但真正动手写的时候很多人会卡在同一个地方点集规模一上来暴力插点或者逐点搜索的复杂度直接爆炸几万个点跑几分钟甚至跑不动。这份资源给的是一个用分治法实现的三角剖分算法核心思路是把点集按空间位置递归切分先解决小规模子问题再合并成全局三角网。它适合两类人一类是正在补计算几何基础、想搞懂分治怎么落到三角剖分上的开发者另一类是有实际工程需求需要一份能读、能改、能直接嵌进自己管线里的实现参考。我拿到之后先跑了一遍随机点集又拿真实地形散点试了试下面把拆解过程、参数设置和踩过的坑一次说清。2. 分治三角剖分的核心逻辑为什么切分点集比逐点插入更稳2.1 分治法在三角剖分里的基本流程分治三角剖分不是把点随便分成两半就完事。它的标准流程是先把点集按 x 坐标排序然后递归地从中位数处切成左右两个子集分别对左右子集做三角剖分最后把两个子三角网沿着切割线合并。合并这一步是整个算法最考验实现的地方因为要找到左右子网的上切线和下切线再沿着这两条切线之间的区域逐点判断可见性把能连的边补上同时删掉跨越切割线的冗余边。常见做法是维护一个凸包或者用双向链接的边表来记录当前子网的边界合并时从最右点对开始像拉链一样往上往下走。这个过程中每次判断一个点是否在另一个三角形的外接圆内就是所谓的空圆准则保证最终结果是 Delaunay 三角剖分。如果只是要一个普通三角网可以省掉空圆判断但那样生成的三角形质量会差很多狭长三角形一多后续有限元计算或者插值都会出问题。我一般会先确认资源里有没有实现空圆检测。如果没有那它大概率是一个简化版适合做可视化或者对网格质量要求不高的场景。如果有那就要看它用的是精确谓词还是浮点容差这直接决定了大规模点集下会不会出现翻转三角形。2.2 递归切分的终止条件与合并策略递归不能无限切下去。常见的终止条件是子集点数小于等于 3直接构成一个三角形或者一条边。也有实现会设成小于等于 4 或者 5然后在小规模子集里用暴力方法建网这样能减少递归深度但合并次数会变多。资源里如果用的是 3那递归树会比较深点集几万的时候栈深度可能到十几层一般不会溢出但如果你改成 2 或者 1那就要小心了。合并策略上有两种常见写法一种是每次合并都重新计算上下切线另一种是维护一个全局的边表合并时只更新受影响的部分。前者实现简单但常数大后者写起来复杂但跑得快。我拿到资源后先看它合并时有没有重复遍历所有点如果有那大规模点集下性能会明显下降。实测下来十万点级别维护边表的版本比每次重算的版本快三到五倍这个差距在交互式应用里是能感知到的。提示如果你只是做几千个点的剖分两种合并策略体感差别不大优先选代码好读的那个。2.3 点集预处理排序和去重为什么不能省分治三角剖分的第一步是按 x 排序。如果点集里有大量重复点或者 x 坐标相同的点排序后切分时会出现左右子集不平衡极端情况下一边全是相同 x 的点递归深度退化成线性复杂度从 O(n log n) 掉到 O(n²)。资源里如果没做去重你一定要在调用前自己加一步。去重的阈值怎么设如果是浮点点集常见做法是设一个 epsilon比如 1e-9把距离小于这个值的点合并成一个。但 epsilon 不能太大否则会改变几何形状尤其是地形数据里相邻点本来就很近。我一般会先统计点集的最小间距然后取最小间距的十分之一作为去重阈值。如果最小间距是 0说明有完全重合的点那就直接去重不用犹豫。另外排序后如果发现 x 坐标相同的点很多可以考虑按 y 坐标二次排序这样切分时左右子集更均衡。这个细节很多实现会忽略但在网格状点集上效果很明显。3. 把算法跑起来从编译到出图的完整操作链3.1 环境准备与依赖确认这份资源如果是 C 写的通常只依赖标准库不需要额外装东西。如果是 Python 版本可能会用到 numpy 做数组操作但核心算法部分应该是纯 Python 或者用 C 扩展加速。我拿到之后先看目录结构确认有没有 Makefile 或者 CMakeLists.txt有的话直接按默认配置编译。# 假设资源根目录下有 CMakeLists.txt mkdir build cd build cmake .. make -j4 # 编译完成后通常会生成一个可执行文件名字可能是 triangulate 或者 demo如果编译报错说找不到某个头文件先检查是不是 C 标准设低了。分治三角剖分里经常用到std::sort和递归C11 就够了但有些实现会用到std::optional或者结构化绑定那就需要 C17。改 CMakeLists.txt 里的CMAKE_CXX_STANDARD就行。Python 版本的话先确认 Python 版本一般 3.8 以上都没问题。如果用到 numpypip install numpy即可。不建议在 Python 里跑十万点以上的剖分解释器开销太大几万点还能接受再往上建议用 C 版本或者加 Cython 封装。3.2 输入点集格式与参数配置资源一般会要求输入一个点集文件常见格式是每行两个浮点数用空格或者逗号分隔。我一般会先拿一个简单的正方形加中心点的测试用例跑通确认输出格式正确再上真实数据。# 生成一个测试点集正方形四个角加中心一个点 points [ (0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.5, 0.5), ] # 写入文件每行 x y with open(test_points.txt, w) as f: for x, y in points: f.write(f{x} {y}\n)跑完之后看输出应该得到四个三角形中心点和四个角都连上。如果输出里出现了重复边或者三角形数量不对先检查去重有没有做再检查合并时上下切线的判断逻辑。参数方面资源里可能暴露一个epsilon或者tolerance用来判断点是否重合或者是否在圆内。这个值不要随便改大默认 1e-9 或者 1e-12 就行。如果点集坐标范围很大比如地形数据里 x 从 0 到 100000那 epsilon 可以适当放大到 1e-6但再大就可能把正常点合并掉。3.3 输出格式与可视化验证输出通常是三角形列表每个三角形三个顶点索引对应输入点集的顺序。有些实现会直接输出边的列表或者输出邻接关系。我一般会先把三角形列表转成可视化能认的格式比如用 matplotlib 画出来看一眼。import matplotlib.pyplot as plt import matplotlib.tri as mtri # 假设 points 是 Nx2 的数组triangles 是 Mx3 的索引数组 triang mtri.Triangulation(points[:, 0], points[:, 1], triangles) plt.triplot(triang, bo-, lw0.5) plt.gca().set_aspect(equal) plt.show()如果画出来发现有大片空白区域或者三角形交叉那说明合并阶段有问题。常见原因是上下切线找错了或者空圆判断的容差设得不对。这时候可以先把点集缩小到 10 个点以内手动推一遍合并过程看看每一步的边表变化是否符合预期。注意可视化的时候把点也画上有时候三角形看起来没问题但点没连上说明有孤立点被漏掉了。4. 避坑与排查分治三角剖分最容易翻车的五个地方4.1 现象递归深度过大导致栈溢出原因点集按 x 排序后如果 x 坐标分布极不均匀比如大部分点集中在很小的 x 范围内少数点 x 很大切分时中位数会偏向一侧递归深度接近 n。解决在切分前先检查 x 坐标的分布如果方差很大可以先做一次坐标归一化或者改用按 x 和 y 交替切分也就是类似 KD 树的策略。资源里如果只按 x 切你可以自己改成按深度交替切 x 和 y。4.2 现象合并后出现翻转三角形法线方向不一致原因空圆判断用的是浮点运算当四个点接近共圆时误差会导致判断结果不稳定合并时连了不该连的边。解决把空圆判断改成精确谓词或者加一个容差当行列式绝对值小于容差时按固定规则处理比如优先保留短边。常见做法是用 Shewchuk 的自适应精度谓词但那个实现比较长如果资源里没有可以先用一个简单的容差顶着容差取 1e-12 左右。4.3 现象点集里有重复点输出三角形数量偏少原因重复点在排序后相邻切分时可能被分到同一边合并时又因为距离为零导致空圆判断失效。解决在输入阶段就去重用哈希表或者排序后遍历把距离小于 epsilon 的点合并。epsilon 取点集最小间距的十分之一如果没有重复点这一步开销很小。4.4 现象大规模点集跑得比预期慢很多原因合并阶段每次都在全量边表里查找或者递归时反复复制点集数组。解决把点集改成传引用或者指针避免递归时拷贝。合并时用局部边表只维护当前子网的边界不要每次遍历全局。如果资源里用的是std::vector按值传递改成const std::vector就能省不少时间。4.5 现象输出三角形有重叠或者缝隙原因上下切线合并时可见性判断漏掉了某些点或者边表更新时删边不彻底。解决在合并阶段加一个断言检查每个新生成的三角形是否与已有三角形相交。如果相交打印出当前上下切线的点和边表状态手动推一遍。常见错误是上切线和下切线的起始点搞反了或者循环终止条件写成了而不是。5. 进阶用法把分治三角剖分嵌进自己的管线5.1 用约束边处理带孔洞或边界的情况标准分治三角剖分只处理凸包内的点集如果你的应用里有孔洞或者凹边界需要加约束边。常见做法是先做一次无约束剖分然后把约束边插入删掉与约束边相交的三角形再局部重新剖分。这一步比分治本身复杂但资源里如果有边表结构插入约束边不算太难。我一般会先把约束边按长度排序从短到长插入这样每次影响的范围小重新剖分的代价低。5.2 并行化左右子集可以同时剖分分治的递归结构天然适合并行。左右子集在合并之前是完全独立的可以用两个线程分别跑。合并阶段需要同步但合并本身只涉及边界区域大部分三角形已经算好了。实测下来四核并行在十万点级别能快两倍多再往上受限于合并阶段的串行部分加速比会下降。如果你用 C可以用std::async或者 OpenMP 的parallel sectionsPython 的话因为 GIL 的存在多线程效果不好建议用多进程或者直接上 C。5.3 验证剖分质量的三个指标跑完不算完得验证结果对不对。我一般看三个指标一是三角形的最小角Delaunay 剖分的最小角应该尽可能大如果出现小于 10 度的角说明空圆判断有问题二是边的总数对于 n 个点的三角剖分边数应该是 3n - 3 - hh 是凸包上的点数如果对不上说明有重复边或者漏边三是每个点的度数内部点度数一般在 6 左右如果某个点度数为 2 或者 3可能是孤立点或者边界处理错了。# 统计三角形最小角 import numpy as np def min_angle(tri_points): # tri_points: 3x2 数组 angles [] for i in range(3): v1 tri_points[(i1)%3] - tri_points[i] v2 tri_points[(i2)%3] - tri_points[i] cos_angle np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2)) angles.append(np.degrees(np.arccos(np.clip(cos_angle, -1, 1)))) return min(angles) # 对所有三角形跑一遍看最小值如果最小角普遍偏小先检查空圆判断的容差是不是太大了再检查点集有没有共线或者共圆的情况。共线点会导致退化三角形面积为零这种要在预处理阶段就剔除。5.4 一个具体技巧用网格加速点定位合并阶段需要频繁判断点是否在三角形内如果每次都遍历所有三角形复杂度会上去。常见做法是建一个均匀网格或者 KD 树把三角形按空间位置索引起来查询时只检查附近几个格子。我一般会用一个简单的均匀网格格子大小取平均边长的两倍这样大部分查询只需要检查常数个三角形。这个优化在点集超过五万时效果很明显能省掉一半以上的合并时间。从那以后我每次拿到新的三角剖分实现都会先跑一个最小用例再跑一个带重复点和共线点的边界用例最后才上真实数据。这三个用例跑通基本就能确认实现是稳的。希望帮到你。本文还有配套的精品资源点击获取