简介：Delaunay三角剖分是生成高质量三角网格的经典算法，这份zip资源围绕该算法在C#场景下的工程实现打包，适合从事mesh建模、图形学或Unity开发的开发者参考。压缩包共25个文件，体积仅44KB，以6个.h头文件、5个.cpp源文件为主体，附带工程配置、资源文件与文本说明，基本覆盖了从点集处理到三角网生成的完整项目框架。已有123人学习下载。资料包含可直接编译的MFC示例工程，如testDelaunay视图与文档类、主框架代码，以及debug.txt等调试记录，能帮助读者快速理解Delaunay三角化在真实建模项目中的调用流程与数据结构组织方式。通过阅读源码可以掌握如何准备点集、调用三角剖分库并导出网格结果，对于学习计算几何和mesh生成具有直接参考价值。
1. 为什么说Delaunay是mesh建模的第一选择：一个经典三角形生成算法的完整解析
做三维重建、游戏地形、有限元仿真或者路径规划，你绕不开一个经典的三角形生成算法：Delaunay三角剖分。为什么这么多年过去，大家在做mesh建模时还是优先选它？因为它输出的是一套拓扑干净的纯三角形网格，天然适配GPU渲染管线，面片之间不重叠、不裂缝，顶点共享关系清晰；更重要的是，“空圆”这条简单约束能让三角形尽量饱满，避免出现细长到让有限元矩阵直接病态的退化单元。这篇文章不是泛泛科普，我把从性质推导、核心实现到网格导出、调试避坑的整条链路都过一遍，新手能照着跑通第一个三角网，老手可以重点看参数调整与边界条件部分。
2. Delaunay的性质与算法直觉：空圆条件、最小角最大化与Voronoi对偶
2.1 空圆条件：任意三角形外接圆内不能有第四个点
要判断一个平面点集的剖分是否是Delaunay，唯一要检查的就是空圆条件：对剖分中的每一个三角形 ABC，它的外接圆内部不允许出现点集中任意其他顶点。注意这里说的是“圆内”，不是“圆上”。如果一个点恰好落在圆边界上，并不破坏条件，这正好区分了标准位置（无四点共圆）和退化情形。很多刚接触三角形生成算法的人，第一次写校验器时把逻辑写成了“小于等于”，结果在对称数据上永远判定失败，其实就是把圆上到圆内的分界线弄混了。
实际工程里我很少逐个三角形去扫描全量点集，因为那是 O(n^3) 级别的暴力检查。更高效的做法是只检查内部边：对于任意一条内部公共边 AB，它关联两个相邻三角形 ABC 和 ABD，只要 D 不落入 ABC 的外接圆内，那么边 AB 就是合格的；对全部内部边执行一遍同样的验证，就能保证整张网格的空圆性质。内部边的数量等于三角形数量的三倍去掉边界边，线性扫描完全能接受。
def check_edge(tri_left, tri_right, eps=1e-10): """ 检查相邻三角形之间的共享边是否满足空圆条件 返回 True 表示不需要翻转 tri_left/tri_right: 顶点三元组 [(x, y), ...] """ # 找出两个三角形的公共顶点（共享边） shared = [p for p in tri_left if p in tri_right] if len(shared) != 2: return True # 各自的非共享顶点 opp_l = [p for p in tri_left if p not in shared][0] opp_r = [p for p in tri_right if p not in shared][0] # 用 in_circumcircle 判断 opp_r 是否落在左侧三角形外接圆内 return not in_circumcircle([shared[0], shared[1], opp_l], opp_r, eps)代码里没有计算距离，只比较距离平方，性能比开方快不少。in_circumcircle 内部返回的是“点是否进入外接圆并越过容差”的布尔值，这个接口设计我在多个项目里一直沿用，因为后续做边翻转判断时只需要这一个函数。如果你发现网格中出现莫名其妙的孔洞，优先怀疑这里的 eps 取值是否和坐标量级匹配：坐标只有 0.1 级别的数据集用 1e-10 和用 1e-6 判定结果差异很大，后者会把大量边界点误判为“圆内点”。
空圆条件看起来只是局部约束，但它会反向塑造网格形态：一个三角形越狭长，它的外接圆半径越大，越容易笼罩邻近的其他点。为了满足空圆条件，算法会倾向于选择较短的对角线来连接顶点，剖分结果自然出现更多的短边和饱满单元。换句话说，Delaunay 并不是“强制”三角形长得匀称，而是通过空圆这个全局不变量，把狭长三角形淘汰掉。
提示：工程里判断“点在圆内”永远比较距离平方，不要对半径开平方。这个习惯能省下大量无谓的 sqrt，也避免半径阈值调整时出现量级漂移。
2.2 最小角最大化：为什么说这套三角形数值上更稳
在所有可能的三角剖分中，Delaunay 剖分会最大化整个三角形集合里“最小内角”的数值。也就是说，给定同样的点集，不管你换哪条对角线，总会有某些三角形的角度更小，而 Delaunay 在“最小角”这一项上做到了最优。对仿真计算而言，这个性质直接决定了刚度矩阵的条件数：一个角度接近 0° 的三角形在有限元装配时会产生极小行列式，整体求解误差被放大到几十倍，渲染场景里则表现为高光走样、光影斑驳的“亮片”三角形。反过来，Delaunay 的最小角最优意味着数值表现的下限是可控的。
特别提醒一点：最小角最大化不是“所有角都最大”，而是“最差的别太差”。如果点集本身有三点近直线，任何算法都无法避免小角，Delaunay 只能保证它比别的方式略好。所以网格质量的第一道保险在输入点集，而不是三角剖分算法本身。采样地形点云、重建扫描数据时，控制最小间距和最大间距的比值，比事后调整网格参数有效得多。我见过不少团队在网格后处理上花大量精力细分、平滑，结果发现瓶颈在原始点密度分布严重不均。
| 剖分类型 | 最小角特征 | 常见后果 |
|---|---|---|
| 随机三角剖分 | 可能接近 0° | 仿真矩阵病态，渲染高光异常 |
| Delaunay 剖分 | 达到全局最大 | 网格尺寸稳定，数值误差可控 |
2.3 与Voronoi图的对偶关系：建模场景中这个视角有什么用
Delaunay 三角剖分与 Voronoi 图互为对偶：点集的每个 Voronoi 单元包含一个输入点，若两个单元共享边界，那么它们对应的输入点之间有一条 Delaunay 边。这个对偶关系看起来像纯几何定理，但在网格建模里有具体用处。做点云重建时，局部区域采样密度不均匀，Voronoi 单元格大小就是自然的内外限指标，它能直接告诉你哪个区域需要补采样，哪个区域已经过度冗余。从 Delaunay 剖分中拿到对应的 Voronoi 单元数据，不需要维护第二套数据结构，只需找到所有相邻三角形并取外接圆心，连线就是 Voronoi 边。
工程里另一个高频场景是面积权重分配。对非结构网格做节点求解时，每个节点的“控制体积”可以近似为相邻 Voronoi 单元的面积，而该面积又可以从 Delaunay 三角形的外心坐标算出来。换句话说，一套 Delaunay 网格既能用于有限元插值，又能用于有限体积法的流量计算，两种数值方法在同一个网格上切换，不需要重新生成网格。这个视角决定了数据结构设计上的取舍：三角形的外接圆心不要只当临时变量，持久化保存，后面做后处理时能省大量重复计算。
3. 从零实现Delaunay：Bowyer-Watson逐点插入法与参数解析
3.1 为什么工程实现多数选Bowyer-Watson而不是递归翻转
实现 Delaunay 剖分，常见方案有四类：递归翻边法、分治法、扫描线法、增量插入法。递归翻边从任意三角网出发，反复检查并翻转不合格的公共边，概念上最贴近性质定义，但最坏情况下要来回判断多条边，性能起伏大；分治法在二维场景下能达到 O(n log n)，但扩展到三维四面体化时需要重写大量数据结构和相交判断；扫描线法依赖顺序输入，一旦要流式插入新点，维护成本陡增。Bowyer-Watson 属于增量插入：每次插入一个点，只在外接圆覆盖到该点的局部区域做删除和重建，坏三角形删除后，边界边与新点重新连成三角形。它最坏情况 O(n²)，但真实点云数据普遍能跑到 O(n^1.2) 上下，还天然支持“先建一部分、再补新点”的流式场景，这对传感器逐步返回点位数据的建模管线特别顺手。
工程上我更看重的是 Bowyer-Watson 的局部性质：一个点的插入只影响一个联通的“空圆冲突区域”，因此删除和重建可以被限制在一个局部邻域。想用并行加速时，把点集按空间分桶，各桶内局部重建后再合并，比分治法的合并阶段简单得多。如果你只需要一次性剖分静态点集，分治法可能更快；但只要涉及增量添加、删除顶点，或者需要在局部加约束，Bowyer-Watson 大概率成为首选。
3.2 最小可运行的Bowyer-Watson实现：Python代码与逐行说明
下面这个实现是从我维护过的一个网格生成工具里剥离出来的，去掉了并行、索引缓存和法向量计算，只保留最核心的插入流程。你可以把它当成骨架，后续按自己的项目扩展。
import numpy as np class Triangle: """三角形，缓存外接圆信息，避免重复浮点计算""" def __init__(self, a, b, c): self.vertices = [tuple(a), tuple(b), tuple(c)] self.center, self.radius_sq = self._circumcalc() def _circumcalc(self): (x1, y1), (x2, y2), (x3, y3) = self.vertices # 2 倍面积：接近 0 表示三点共线 d = 2.0 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) if abs(d) < 1e-12: return (0.0, 0.0), float("inf") x1sq = x1*x1 + y1*y1 x2sq = x2*x2 + y2*y2 x3sq = x3*x3 + y3*y3 ux = (x1sq * (y2 - y3) + x2sq * (y3 - y1) + x3sq * (y1 - y2)) / d uy = (x1sq * (x3 - x2) + x2sq * (x1 - x3) + x3sq * (x2 - x1)) / d r_sq = (x1 - ux)**2 + (y1 - uy)**2 return (ux, uy), r_sq def in_circle(self, p): px, py = p return (px - self.center[0])**2 + (py - self.center[1])**2 \ < self.radius_sq - 1e-12 def bowyer_watson(points): # 1. 先做坐标去重：浮点坐标容易出现重复点 pts = list({tuple(p) for p in points}) if len(pts) < 3: return [] # 2. 构造超级三角形，要求完全覆盖全部点 xs = [p[0] for p in pts] ys = [p[1] for p in pts] min_x, max_x = min(xs), max(xs) min_y, max_y = min(ys), max(ys) dx, dy = max_x - min_x, max_y - min_y mid_x, mid_y = (min_x + max_x) / 2, (min_y + max_y) / 2 R = max(dx, dy) * 4 super_tri = Triangle( (mid_x - 2*R, mid_y - R), (mid_x, mid_y + 2*R), (mid_x + 2*R, mid_y - R) ) tris = [super_tri] # 3. 逐点插入 for p in pts: # 找出外接圆包含 p 的所有坏三角形 bad_tris = [tri for tri in tris if tri.in_circle(p)] # 提取边界边：坏三角形内部共享边出现两次，外部边界边只出现一次 edge_set = set() for tri in bad_tris: for i in range(3): edge = frozenset((tri.vertices[i], tri.vertices[(i+1) % 3])) if edge in edge_set: edge_set.remove(edge) else: edge_set.add(edge) # 删除坏三角形，把边界边与新点连成新三角形 tris = [tri for tri in tris if tri not in bad_tris] for edge in edge_set: a, b = tuple(edge) tris.append(Triangle(a, b, tuple(p))) # 4. 移除所有与超级三角形共享顶点的三角形 super_verts = set(super_tri.vertices) return [tri for tri in tris if not super_verts.intersection(tri.vertices)]这里有一个经常让新手迷惑的点：edge_set 用 frozenset 做差集。因为无向边 (a, b) 和 (b, a) 必须被视为同一条边，frozenset 恰好抹掉方向差异。如果你换成普通元组，内部公共边会残留两条记录，新三角形连接就会多出悬空边，网格随之出现裂缝。另一个关键点是超级三角形的构造：我用包围盒中心加 4 倍半径的三角形。这个倍数不是随便拍的，太小会让部分输入点落在超级三角形外部，插入过程不再正确；太大则会生成大量冗余三角形，虽然最后会被清理，但内存和计算都会白白增加。
3.3 三个直接影响结果的地方：超级三角形、浮点阈值与点集去重
首先是超级三角形尺寸。我在实测中发现，同样的代码在坐标跨度 2000 的大场景上正常，换到跨度只有 0.8 的桌面零件上，凸包边缘丢了好几个角点。原因就是超级三角形取了固定 3 倍外扩，而那次输入点里有几个点恰好落在超级三角形边界的凹陷处。解决办法是把外扩倍数提到 4，再写一段校验逻辑：临时渲染或打印所有点，确认它们都在超三角形三边内侧。代价只是多生成十几个边界三角形，后面清掉即可。
其次是浮点阈值。外接圆计算和空圆判定不要共用一个 eps，前者是单次数值计算，误差受坐标绝对量级影响；后者要跨所有三角形作比较，累积误差来自多次插入。分两个参数传，排查网格抖动时才能隔离是哪一段的问题。坐标经过平移归一化后再建网，也就是把点集中点挪到原点，能显著改善大坐标场景下的数值稳定性。我处理经纬度平面投影数据时都会先做这一步，不然外心计算里的大数相减会让半径平方失真。
第三是点集去重。浮点坐标输入时，重复点不处理会让算法陷入退化三角形，外接圆半径变成无穷大，in_circle 直接失真。我见过有人在这个坑里调了一周，最后发现点云里有两千多个完全相同的坐标。去重时用 tuple 而不是数组作为集合元素，因为 numpy 数组的哈希行为不稳定，tuple 才能保证正确去重。
4. 从散点到Mesh：把Delaunay三角剖分接进建模管线的三种落地方式
4.1 地形网格TIN：按高度场投影构建三角网
最直接的落地场景是地形网格，也就是把高程点云转成 TIN（不规则三角网）。做法是先把三维点 (x, y, z) 投影到二维平面，对 (x, y) 做 Delaunay 剖分，得到连通关系后再把 z 值粘回去。
# 假设 points_xyz 是 N x 3 的 numpy 数组 pts_2d = points_xyz[:, :2] # 只保留 x, y 用于剖分 tris = bowyer_watson(pts_2d) # 得到二维三角形集合 # 后续重建三维三角形时，通过顶点坐标映射回原始 Z 值这段流程的关键是：二维剖分只关心平面位置，不感知地形起伏。坡度大的悬崖处，二维 Delaunay 会在投影后产生很窄的面片，导致三维网格出现“拉链式”的交叉面。我在处理山地地形时，会在预处理阶段做高差过滤：相邻点的高差除以平面距离超过某个阈值（比如 45°），就把这对点中的一个保留、一个剔除，避免剖分连接了本来不该连接的山体两侧。阈值设置取决于地形用途，飞行仿真地形网格可以放宽到 60°，有限元分析则建议收紧到 30° 左右。
4.2 从三角形数组到OBJ/STL：Mesh导出必经之路
拿到 Triangle 数组还不算真正完成建模，渲染工具和仿真软件需要的是顶点索引和面列表。OBJ 格式是最常见的中转格式，结构简单、可读性好，而且对第三方工具兼容性也最好。导出时最需要注意的坑是索引从 1 开始计数，很多第一次写导出器的人在这里翻车，导出的文件面全部指向错误的顶点。
def export_obj(verts_3d, triangles, filepath): """ verts_3d: [(x, y, z), ...] 唯一顶点列表 triangles: [(idx_a, idx_b, idx_c), ...] 三角形顶点索引，0 起始 """ with open(filepath, "w") as f: for v in verts_3d: f.write(f"v {v[0]:.6f} {v[1]:.6f} {v[2]:.6f}\n") for t in triangles: # OBJ 索引从 1 开始，所以每个下标 +1 f.write(f"f {t[0]+1} {t[1]+1} {t[2]+1}\n")除了索引起始，还有一个隐性问题：三角形绕序。如果相邻三角形绕序不一致，导入三维软件后会出现部分表面发黑的现象，因为法向量指向背离了观察方向。可以在导出前加一遍统一绕序的检查，用右手定则计算每个三角形的法向量，再对比它与网格中心的连线方向，不一致就交换两个顶点。这个小习惯省掉了渲染阶段大量排查时间。
STL 格式则是二进制的实体网格格式，适合 3D 打印和简单的碰撞体表示。它不保存顶点索引，只保存每个三角形的独立坐标和法向量，所以文件体积会比 OBJ 大不少。如果项目目标是有限元仿真，我建议直接导出含拓扑连接的格式，如 OBJ 或带索引的 VTK，STL 缺少拓扑信息，后续做网格细化、重划分时需要重建邻接关系，是一笔额外的开销。
4.3 在真实建模管线中嵌入Delaunay：预处理与后处理不要乱序
实际业务系统里，Delaunay 很少独立运行，通常夹在预处理和后处理之间。我习惯固定这样的顺序：先去噪、降采样，再剖分，然后补洞、平滑。这个顺序一旦颠倒，会出现大量无效顶点和碎片三角形被烘焙进网格里的情况。
# 伪命令流，工具名按实际项目的可执行文件替换 geometry_preprocess input.xyz --radius 0.02 --out clean.xyz geometry_mesh --algo delaunay clean.xyz --out init.obj geometry_postprocess init.obj --fill-holes --smooth --out final.obj预处理阶段我通常做两步：第一步基于体素或格网做均匀降采样，保证点之间的最小间距接近；第二步基于曲率做自适应采样，曲率大的区域保留更多点。这两步组合使用，比任何单一策略都稳定。后处理阶段，孔洞填补和网格平滑的顺序也不要乱改：平滑前补洞可以让新填充的面片和周围顶点共享同一套邻接信息，避免平滑过程中拉出新裂缝。
网格平滑我常用拉普拉斯平滑，简单实现就是把每个顶点移动到相邻顶点的平均值位置。它的参数只有一个迭代次数，通常 3 到 5 次即可，过多次数会把细节全部磨掉。迭代 10 次以上的地形网格，山头会被明显削平，这是我在对比多次实验结果后得出的经验值。如果担心整体收缩，可以采用限制位移量的版本，让每次移动距离不超过边长的 1/4。
5. 避坑指南：Delaunay三角剖分最常见的5个工程翻车点
5.1 共圆点：对称数据导致网格输出不一致
现象：四个点恰好共圆时，例如一个正方形的四个顶点，每次运行剖分输出的对角线连接方向可能不同。最后一次插入顺序变了，对角线就跟着翻转。 原因：这四个点同时满足空圆条件，Delaunay 剖分在这个局部不唯一，“圆内”无法区分谁更合法。工程数据里，规则格网的采样点最容易出现这种退化配置。 解决：输入阶段给坐标叠加上一个非常小的扰动，例如网格最小间距的千分之一，把退化配置转成一般位置。扰动会略微改变网格大小，但对宏观建模几乎无感知。注意扰动值本身要一致，同一个点在所有处理阶段都使用同一份扰动坐标，不然会出现同一个点在两套坐标间跳变。
5.2 浮点误差：外接圆边界判定引发“抖动”
现象：某个点明明位于三角形外接圆外，判定结果却是“圆内”，插入流程删掉了不该删除的三角形，网格出现重叠面或拓扑缝隙。 原因：外接圆心坐标通过除法求解，当三角形接近退化时圆心会被推到很远，距离平方从 1e10 和 1e9 这种量级比较时，浮点精度不够稳定。 解决：统一使用dist_sq < r_sq - eps的形式，给比较留出容差。另外可以在剖分前对点集做去中心化，把平均值平移到原点，把坐标值从百万级降到个位数级，外心计算精度显著提升。我在处理像素坐标时经常这样做，效果好到直接省去了后续补洞步骤。
5.3 超级三角形设置过小：凸包边缘丢点
现象：剖分结果凸包少了一个角，边缘变成一条直线；或者边界出现悬挂的孤立线段。 原因：超级三角形没有完全覆盖所有输入点，导致外部点插入时，它的空圆冲突区域没有被正确识别，算法漏掉了和该点相连的三角形。 解决：把超级三角形的外扩倍数从 3 提高到 4，并用可视化方式检查全部点确实落在超三角形内部。调试时可以临时把超三角形设为半透明色输出到预览窗口，扫一眼就能发现问题。有些实现会使用“动态超级三角形”，在插入点超出当前覆盖范围时自动扩容，但这样会让代码复杂不少，我一般建议静态超三角形加足够大的余量。
5.4 大数据量下性能塌方：从 O(n log n) 到 O(n²)
现象：1 万点跑得很流畅，点集涨到 10 万后构建时间从毫秒级跳到秒级，内存占用同时飙升。 原因：每插入一个新点，代码都遍历所有现存三角形来寻找坏三角形，整体复杂度退化到 O(n²)。Python 里还在反复创建对象，常数因子被放大，问题雪上加霜。 解决：给三角形集合建一个均匀分桶，按外接圆圆心或包围盒的格网索引存储。插入新点时只去当前格和邻近格寻找候选三角形，不做全局扫描。这个优化把绝大多数点集的实际构建时间拉回了 O(n log n) 附近。
# 利用外接圆圆心做空间分桶，缩小候选三角形范围 cell_size = max(dx, dy) / max(1, int(np.sqrt(len(pts)))) bucket = {} for idx, tri in enumerate(tris): cx, cy = tri.center key = (int(cx // cell_size), int(cy // cell_size)) bucket.setdefault(key, []).append(idx) # 插入新点 p 时，只需检查 p 所在桶和 8 个邻桶 px, py = int(p[0] // cell_size), int(p[1] // cell_size) candidates = [] for kx in range(px - 1, px + 2): for ky in range(py - 1, py + 2): candidates.extend(bucket.get((kx, ky), []))代码里的 cell_size 取的是包围盒边长除以点数的平方根，这个值要让每个桶平均容纳 5 到 10 个三角形。桶太大，优化效果不明显；桶太小，邻桶列表本身变成瓶颈。针对十亿级场景服务器，可以进一步用八叉树或 R-tree，但对绝大多数桌面级点云工具，格网方案已经足够了。
5.5 边界约束缺失：网格在孔洞处出现“跨海大桥”
现象：点集本身有缺口或者多个区域，剖分结果却连成了一张完整的三角形面，不该连接的地方生成了一条跨越空白的细长三角形。 原因：Delaunay 剖分的作用域是点集的凸包，它不感知边界约束和孔洞。只要两个区域在凸包内，算法就会“热心”地把它们连起来。 解决：换用约束 Delaunay 剖分，把边界边和孔洞轮廓强制设置为约束边。实现上的改动集中在边对象上：增加 is_constraint 标志，在 Bowyer-Watson 删除坏三角形的阶段跳过含约束边的三角形，新生成的边如果穿过约束边则禁止连接。
class Edge: def __init__(self, a, b, is_constraint=False): self.a, self.b = a, b self.is_constraint = is_constraint约束边的判断核心是“不可破坏”，做二维地图建模、建筑平面图网格划分时几乎必用。如果没有约束边支持，最简单的替代方案是在孔洞区域先放一些“虚拟点”，剖分后删除关联三角形，但这种方式容易让孔洞边界变成锯齿状，后续还得花精力做边界裁剪。
6. 进阶：约束Delaunay、局部加密与网格质量维护
6.1 让三角网尊重边界：CDT的最小改造
三维建模里常见的河流边界、断层线、地物轮廓都不允许三角网跨越。给 Bowyer-Watson 增加约束边支持，核心是两处改动：三角形对象记录自己包含哪些边，边对象携带约束标志；删除坏三角形时，若某条约束边的两侧都包含在待删除集合里，则强制保留其一侧为新三角形的一部分。这个逻辑可以保证约束边最终出现在输出网格中。
约束边的初始输入通常来自二维多边形轮廓。排序处理后把轮廓点直接插入点集，同时给相邻轮廓点之间的边加上约束标志，剩下的流程与普通 Delaunay 一致。要处理多个彼此不相邻的区域，可以先各自构建局部三角网，再在缝合边界上补充带约束边的过渡网格，比一次性全局处理容易控制质量。
6.2 外接圆圆心插入法：低成本自适应加密
做有限元前处理时，经常要对指定区域加密网格。我不建议全局重新剖分，那样会把已满足质量要求的部分打回原形，成本也高。更精准的做法是局部加密：找出质量最差的三角形（比如面积和周长的比值低于阈值），在其外接圆圆心处插入一个新点，然后只对这个局部邻域重新执行 Bowyer-Watson 更新。插入点选外心而不是重心，因为外心正是当前空圆判定的薄弱中心，一次性让该区域的三角形自动重组，往往比在重心插入更有效地消掉最坏三角形。这个操作可以迭代执行，每轮只处理最差的几个三角形，网格质量稳步提升，整体规模也不会爆炸式增长。
网格质量维护方面，我习惯在每次更新后跑两条统计：所有三角形的最小角分布，以及共享边一致性。最小角低于 15° 的三角形占比超过 2%，就说明输入点集或参数有问题；共享边一致性的检查则是遍历每条内部边，验证两端各只有一个相邻三角形，出现两个以上相邻三角形意味着拓扑已经损坏。这两项检查加上约束边逻辑，Delaunay 剖分在建模管线里的稳定性就有了基本保障。
我自己做网格工具时，最后的习惯是任何网格交付前都跑一遍这两项校验，数据达标再谈性能优化、并行加速。毕竟三角形生成算法再好，也扛不住在拓扑错误的基础上反复调试。希望这个习惯能帮到你的项目。
本文还有配套的精品资源，点击获取
