尧图网页设计 
做三维扫描、逆向工程或者 3D 打印前处理的人大概都有过这样的经历模型渲染出来看着挺完整一转到线框模式或者丢进切片软件立马就露馅了——一个个窟窿张着嘴边缘还毛刺刺的。补洞hole filling和 fairing 算法这两个词基本就是绕不过去的坎。前者负责把缺失的面片填回去让模型重新闭合成水密实体后者负责把填完之后的补丁磨得跟周围曲面连成一片不留折痕。这不是什么锦上添花的功能而是从扫描数据到可用模型的必经环节做有限元仿真、做增材制造、做数字孪生资产入库第一步都是这个。我把这两个算法放在一起讲是因为单独做任何一个都容易出问题。只补不光顺补丁像贴了块膏药法向突变会让切片器报错、让仿真应力集中只光顺不补洞边界会越磨越缩洞反而变大。下面这套东西是我这几年在扫描件修复、模具逆向、文物数字化几个项目里反复折腾出来的代码基于 Python 的 trimesh scipy 生态原理讲清楚参数给具体值踩过的坑也都写进去了。不管你是刚接手网格处理任务的新手还是已经用过 Meshlab 自动补洞但效果不理想的老手应该都能找到能直接用的部分。1. 补洞与 fairing 到底在修什么从破洞的成因说起1.1 网格上那些洞是怎么来的很多人以为洞是软件算错了其实绝大多数洞是物理世界的信息缺失算法只是把缺失暴露出来。搞明白成因才能判断应该用哪种补法。我把常见成因归成四类每一类对应的补洞策略其实不一样。第一类是扫描遮挡。激光或者结构光扫描仪本质上是个视线传感器光打不到的地方就没数据。深腔、倒扣、叶片间的缝隙、螺纹根部这些位置永远扫不全。这类洞的特点是边界往往沿着视线方向拉长形状不规则而且洞的周围通常有噪声。第二类是表面材质问题。高反光、全透明、纯黑的表面扫描时要么过曝要么无返回。这类洞通常成片出现边界比较圆滑周围点云密度均匀。第三类是几何操作遗留。布尔运算、抽壳、镜像、网格简化每一步都可能产生退化面片。尤其是模型简化的抽取操作会把细长三角形去掉留下一个狭长裂缝。这类洞的特点是呈线状两侧边界几乎平行。第四类是文件转换丢失。STL、OBJ、PLY 之间来回倒精度截断导致顶点对不上原本共享的顶点被拆成两个边就变成了边界边。这类洞最冤明明模型是好的只是顶点没有合并。提示拿到一个有洞的模型先别急着跑补洞。做一次顶点合并merge by distance和退化面清理往往能消掉一半以上的假洞。我见过太多案例几十个所谓的洞合并一次顶点就只剩三个了。1.2 fairing 不是简单的磨皮fairing 这个词在中文里通常译作光顺但很多人把它理解成平滑这是两个层次的事。平滑是去掉高频噪声光顺是让曲面满足某种能量最小的条件——数学上就是让曲面的二阶或者三阶导数连续。放到补洞场景里为什么不能只做平滑因为补丁是凭空生成的它和周围原始曲面之间天然存在缝隙。你如果只做局部平滑边界处的法向不连续会一直在。你需要的是让补丁区域作为一个整体去承接周围曲面的曲率趋势让交界的二面角接近 180 度。打个比方补洞相当于在一块布料上剪掉一个洞然后用另一块布补上去。直接缝上去接缝处一定有褶子。fairing 做的就是让新补的布按照周围布的纹理和张力自然延展接缝摸不出来。从数学上看fairing 求解的是一个带约束的能量最小化问题。约束是补丁边界上的点必须落在原网格边界上能量是曲面的某种形变度量。不同的能量选择得到的补丁形状完全不同这也是为什么 fairing 算法有那么多种。1.3 为什么这两件事得连着做分开做的问题在于信息不对称。补洞算法如果不考虑后续 fairing它会把洞填成一个平面片——因为平面三角化最简单。但这个平面片跟周围曲面完全不搭fairing 阶段要么磨不平要么一磨就把边界也带跑了。反过来fairing 算法如果不了解补洞的结构它会把整个网格当成一个整体去光顺结果就是原始模型上那些本来该保留的锐边、特征线全被磨圆了。做机械件的逆向时这种错误是灾难性的一个倒角被磨掉装配就出问题了。所以实际流程应该是补洞阶段就引入曲率约束让初始补丁尽可能贴近周围曲面的趋势fairing 阶段只对新增的自由顶点操作把原始边界上的顶点锁死或者只允许它们在切平面内移动。这样两个阶段互相配合原始几何一点不动补丁区域又自然过渡。我在做模具修复项目的时候一开始就是分开跑的补洞用现成工具fairing 用拉普拉斯平滑结果零件上的型面分界线全糊了。后来改成一体化处理边界顶点固定只对内部顶点做双调和光顺效果立竿见影。2. 补洞算法的路线选型与背后的取舍2.1 三条主流路线直接三角化、参数化拟合、数据驱动补洞算法按复杂度大致分三档选哪档取决于洞的大小、形状、以及你对结果的要求。第一档直接三角化。把边界环当成一个多边形在三维空间里直接做三角剖分。代表算法是耳切法ear clipping和最小权重三角化minimum weight triangulation。优点是快几百个顶点的洞毫秒级出结果缺点也明显生成的补丁基本是个平面或者鞍面完全没有曲率延拓的概念。适合什么适合小洞边界顶点在 50 个以内、且洞本身近似平面的情况。螺丝孔、小裂纹、扫描小缺陷这一档够用。第二档参数化拟合。先把边界环投影到一个平面域上在二维平面里做带约束的三角化再把平面上的点映射回三维最后用周围曲面的信息去调整新增顶点的位置。这一档是工业级的主流做法。核心思想是降维求解把复杂的三维曲面问题变成二维的平面问题。代表实现就是 Liepa 的算法最小面积三角化打底然后细分加光顺。第三档数据驱动。用深度学习模型直接从残缺网格或者点云预测缺失部分。这几年在扫描补全领域有不少工作优点是能补出语义合理的形状——比如一个椅子的缺角模型知道该补出椅子腿的形状而不只是一个平面片。缺点是需要训练数据泛化能力依赖数据分布工业件上不太敢用。我的建议是批量处理同类资产比如一批同款家具扫描件可以试试单件精密零件还是老老实实走第二档。下面这张表是我对不同路线的实际体感总结供选型参考路线适用洞尺寸曲率延拓单洞耗时实现难度典型场景耳切/最小权重三角化边界点 50无毫秒级低小裂纹、螺孔、扫描缺陷平面参数化 拟合边界点 50~2000有秒级中扫描件大洞、模具型面3D 直接拟合RBF/MLS边界点 500有秒级到分钟级中高曲面件、局部缺失数据驱动补全不限语义级秒级高批量同类资产2.2 边界环提取补洞的地基工程补洞第一步永远是找洞。而找洞的本质是找那些只被一个面片使用的边。在一个封闭流形网格里每条边恰好被两个三角形共享一旦某条边只被一个三角形使用它就是边界边。道理简单实操坑很多。第一个坑是顶点没有合并。两个顶点坐标差 1e-6浮点数认为它们是两个点那条边就变成了两条边界边洞被拆得粉碎。第二个坑是非流形边一条边被三个甚至更多面共享这种网格严格来说不是二维流形边界环的串接逻辑会乱掉。我的处理顺序是先合并顶点容差取包围盒对角线的 1e-5 量级再删退化面面积接近零的再删重复面最后才提取边界。import numpy as np import trimesh def clean_mesh(path, merge_tol_ratio1e-5): mesh trimesh.load(path, processFalse) diag np.linalg.norm(mesh.bounds[1] - mesh.bounds[0]) # 合并重合顶点这一步能消掉大部分假洞 mesh.merge_vertices(digits_vertexint(np.ceil(-np.log10(merge_tol_ratio * diag)))) mesh.remove_degenerate_faces() mesh.remove_duplicate_faces() mesh.remove_unreferenced_vertices() return mesh def find_boundary_edges(mesh): # edges_sorted 是已排序的顶点索引对方便做唯一性统计 uniq, counts np.unique(mesh.edges_sorted, axis0, return_countsTrue) return uniq[counts 1]拿到边界边之后要做的是把它们串成闭环。串接的逻辑是从任意一个端点出发沿着邻接关系一路走直到回到起点。这里要处理顶点度数大于 2 的情况——数学上应该在非流形处断开工程上我通常按选择转弯角度最小的下一条边来走这样得到的环更符合直觉。from collections import defaultdict def build_loops(boundary_edges, vertices, min_len3): adj defaultdict(list) for a, b in boundary_edges: adj[a].append(b) adj[b].append(a) visited_edge set() loops [] for a, b in boundary_edges: key (min(a, b), max(a, b)) if key in visited_edge: continue loop [a, b] visited_edge.add(key) prev, cur a, b while True: cands [n for n in adj[cur] if n ! prev and (min(cur, n), max(cur, n)) not in visited_edge] if not cands: break if len(cands) 1: # 非流形分叉挑方向延续性最好的那条 d0 vertices[cur] - vertices[prev] d0 d0 / (np.linalg.norm(d0) 1e-12) best, best_dot None, -2 for n in cands: d1 vertices[n] - vertices[cur] d1 d1 / (np.linalg.norm(d1) 1e-12) if np.dot(d0, d1) best_dot: best_dot, best np.dot(d0, d1), n nxt best else: nxt cands[0] visited_edge.add((min(cur, nxt), max(cur, nxt))) if nxt loop[0]: break loop.append(nxt) prev, cur cur, nxt if len(loop) min_len: loops.append(loop) return loops注意这段代码在遇到一个顶点连着 4 条以上边界边的病态区域时仍然可能给出奇怪的结果。我的经验是先把模型丢进 Meshlab 看一眼非流形边的分布如果集中在一个小区域宁可手工切掉那块再补也不要指望算法自动处理。2.3 小洞快修最小权重三角化的实现细节对小洞我基本都用最小权重三角化。它的思路是在边界环构成的多边形上做三角剖分但每条候选对角线都有代价选总代价最小的那套剖分。代价函数一般包含两项——三角形面积以及三角形法向与边界邻域法向的偏差。为什么加法向项因为纯面积最小化会倾向于生成扁平的剖分对非平面的洞这会导致补丁凹陷。加上法向惩罚后算法会更愿意生成跟周围曲面朝向一致的三角形。Liepa 那篇 2003 年的论文里给的权重是$$\omega(f) \frac{1}{2}\left(\sum_{e \in f} |e|^2\right) \cdot \max_{e \in f}\left(\frac{1 - \cos \theta_e}{2}\right)$$其中 $|e|$ 是边长$\theta_e$ 是这条边与相邻边界边的夹角。这个式子的直观含义是面积大的三角形代价高与边界走向偏差大的三角形代价也高。两项相乘逼着算法选出既紧凑又顺滑的剖分。实现上用动态规划。设边界环上有 $n$ 个顶点定义 $M[i][j]$ 为从第 $i$ 个顶点到第 $j$ 个顶点这一段子多边形的最小代价。递推就是枚举中间的切分点 $k$取 $M[i][k] M[k][j] \omega(i,k,j)$ 的最小值。复杂度 $O(n^3)$$n$ 在 200 以内都很快。def min_weight_triangulation(pts, normals, edges_dir): n len(pts) INF float(inf) M np.full((n, n), INF) choice np.full((n, n), -1, dtypeint) def tri_cost(i, k, j): a, b, c pts[i], pts[k], pts[j] area 0.5 * np.linalg.norm(np.cross(b - a, c - a)) # 法向与边界走向的偏差惩罚 n_face np.cross(b - a, c - a) norm np.linalg.norm(n_face) if norm 1e-12: return INF n_face / norm local_n (normals[i] normals[k] normals[j]) / 3.0 local_n / (np.linalg.norm(local_n) 1e-12) cos_dev np.clip(np.dot(n_face, local_n), -1, 1) return area * (1.0 - cos_dev) area * 1e-6 for i in range(n - 1): M[i][i 1] 0.0 for span in range(2, n): for i in range(0, n - span): j i span best INF best_k -1 for k in range(i 1, j): if M[i][k] INF or M[k][j] INF: continue c M[i][k] M[k][j] tri_cost(i, k, j) if c best: best, best_k c, k M[i][j] best choice[i][j] best_k return M, choice这段代码只处理单个环实际项目里要循环处理所有环并且对每个环单独建点表。另外要注意点的顺序——边界环的顶点顺序必须是沿着边界连续排列的反了会让三角形法向整体翻转。2.4 大洞慢补参数化与曲面拟合洞的边界点超过 100 个或者洞的形状明显非平面最小权重三角化出来的补丁就开始鼓包或者塌陷了。这时候要上参数化拟合。核心流程分四步走。第一步求边界环的最优拟合平面用 PCA 就能做对边界点做中心化然后 SVD 分解最小奇异值对应的右奇异向量就是平面法向。第二步把边界点投影到这个平面上得到一组二维点。第三步在二维域内做带约束的 Delaunay 三角化注意如果边界是凹的普通 Delaunay 会跑到外面去要用带约束的版本比如triangle库或者自己实现凸分解。第四步把三角化结果映射回三维然后对内部新增顶点做位置优化。第四步是整个流程的灵魂。最朴素的做法是把平面上的点按原位置映射回三维得到的补丁是个平面。好一点的做法是用周围曲面的法向信息做外推让补丁顺着曲率走。更好的做法是求解一个径向基函数RBF拟合问题。RBF 的思路是找出一个函数 $f(x) \sum_i w_i \phi(|x - x_i|) p(x)$让它精确通过边界点及其法向约束然后在洞内采样。$\phi$ 通常取薄板样条核 $\phi(r) r^2 \log r$这个核对应的是薄板弯曲能量最小天然就是光顺的。$\phi(r) r^3$ 也常用对应的是线性化后的薄板能量。import numpy as np from scipy.spatial.distance import cdist def fit_rbf_surface(boundary_pts, interior_uv, kerneltps, reg1e-8): # boundary_pts: (m,3) 边界点 # interior_uv: (k,2) 洞内待求解的点在平面域上的坐标 # 返回洞内点的三维坐标估计 m len(boundary_pts) D_bb cdist(boundary_pts, boundary_pts) if kernel tps: def phi(r): r np.where(r 1e-12, 1e-12, r) return r * r * np.log(r) else: def phi(r): return r ** 3 K phi(D_bb) reg * np.eye(m) # 加低阶多项式项保证线性精度 P np.hstack([np.ones((m, 1)), boundary_pts]) A np.zeros((m 4, m 4)) A[:m, :m] K A[:m, m:] P A[m:, :m] P.T rhs np.zeros(m 4) rhs[:m] boundary_pts[:, 0] # 逐轴求解 w np.linalg.solve(A, rhs) return w这段代码只是骨架实际使用时要对 x、y、z 三个轴分别求解还要注意 RBF 是全局的边界点一多矩阵就巨大。我的经验是边界点在 500 以内直接稠密求解没问题超过就上紧支撑核比如 Wendland 核或者分区求解。说到分区这里有个实用技巧如果洞特别大先对边界环做一次一维的曲线光顺把边界当成一条闭合曲线用拉普拉斯算子平滑把扫描噪声引起的边界抖动消掉再做参数化。否则那些高频抖动会被 RBF 原样拟合进去补丁表面会起波纹。3. fairing 算法的数学底子与实现要点3.1 伞形算子与拉普拉斯光顺补完洞只是第一步新生成的顶点位置还很粗糙。最经典的光顺算法是伞形算子umbrella operator本质上就是离散的拉普拉斯算子。对顶点 $v_i$它的伞形算子定义为邻域顶点平均值与自身位置之差$$U(v_i) \frac{1}{\sum_j w_{ij}} \sum_{j \in N(i)} w_{ij}(v_j - v_i)$$更新规则是 $v_i \leftarrow v_i \lambda U(v_i)$$\lambda$ 一般取 0.3 到 0.6。这个操作重复若干次网格就会越来越光滑。权重 $w_{ij}$ 有两种常见取法。均匀权重 $w_{ij} 1$实现最简单但网格疏密不均的时候会出问题——顶点密的地方贡献大结果把点往疏的地方拉。余切权重cotangent weight考虑了几何形状$$w_{ij} \frac{1}{2}(\cot \alpha_{ij} \cot \beta_{ij})$$其中 $\alpha_{ij}$、$\beta_{ij}$ 是边 $ij$ 所对的两个角。这个权重有理论保证在规则三角网上用余切权重的拉普拉斯算子能精确逼近连续曲面的拉普拉斯-贝尔特拉米算子。def build_cotan_laplacian(V, F): import scipy.sparse as sp n len(V) I, J, Wdata [], [], [] for tri in F: for k in range(3): i, j, o tri[k], tri[(k 1) % 3], tri[(k 2) % 3] v1 V[i] - V[o] v2 V[j] - V[o] cross np.cross(v1, v2) area2 np.linalg.norm(cross) if area2 1e-14: cot 0.0 else: cot np.dot(v1, v2) / area2 w 0.5 * cot I [i, j] J [j, i] Wdata [w, w] W sp.coo_matrix((Wdata, (I, J)), shape(n, n)).tocsr() # 对角元取行和的相反数得到标准拉普拉斯 d np.array(W.sum(axis1)).flatten() L sp.diags(d) - W return L, W注意余切权重在钝角三角形上会出现负值。钝角大于 90 度的余切是负的这会导致拉普拉斯算子失去正定性光顺结果发散。处理方法有两个一是把负权重截断到 0二是用虚拟边界技巧在钝角顶点处补一个辅助点。我一般用截断简单有效代价是精度略降。网格质量好的话根本遇不到这个问题。3.2 从拉普拉斯到双调和薄板能量为什么更靠谱纯拉普拉斯光顺有个致命缺陷它最小化的是一阶能量 $\int |\nabla v|^2$得到的曲面是调和曲面。调和曲面有个特性——它会收缩。一个球面被反复做拉普拉斯光顺最终会缩成一个点。放在补洞场景里就是补丁会往洞中心缩补完的曲面比周围低一块。解决思路是提高能量阶数。双调和能量 $\int |\Delta v|^2$ 对应的欧拉-拉格朗日方程是 $\Delta^2 v 0$解叫双调和曲面它保证了曲率连续同时不会像调和曲面那样有强烈的收缩倾向。物理上这对应薄板弯曲能量最小——一块薄钢板被固定在边界上它自然形成的形状就是双调和曲面。离散化之后最小化 $|L x|^2$ 的问题可以写成线性系统。把顶点分成自由顶点补丁内部的和固定顶点边界上的求解$$\min |L x|^2 \quad \text{s.t. } x_{fixed} b$$对自由部分求导得正规方程 $L_{ff}^T L_{ff} , x_f -L_{ff}^T L_{fx} , b$。这是个稀疏线性系统用spsolve直接解几千个自由顶点毫无压力。import scipy.sparse as sp import scipy.sparse.linalg as spla def biharmonic_fairing(V, L, free_idx, fixed_idx, tikhonov1e-9): n len(V) A (L.T L).tocsc() if tikhonov 0: A A tikhonov * sp.eye(n, formatcsc) A_ff A[free_idx][:, free_idx] A_fx A[free_idx][:, fixed_idx] rhs -(A_fx V[fixed_idx]) V_new V.copy() sol spla.spsolve(A_ff.tocsc(), rhs) V_new[free_idx] sol return V_new这里加 Tikhonov 正则项的原因是$L$ 的零空间是常数向量理论上 $L^T L$ 的零空间也是常数向量只要自由顶点集不构成整个连通分量$A_{ff}$ 就是可逆的。但数值上接近奇异加一点点正则更保险1e-9 量级对结果几乎没有影响。双调和光顺的效果比拉普拉斯好太多。我做过一个对比测试同样一个 80 个边界点的圆洞拉普拉斯光顺 50 次后补丁中心比周围曲面低 0.3mm模型尺寸 100mm双调和一次求解后偏差只有 0.02mm。而且双调和不需要迭代一次求解直接到位。3.3 特征保持怎么磨而不圆fairing 最大的风险是把该保留的特征磨没了。具体到补洞场景有三种特征需要保护。第一种是补丁边界本身。如果洞开在一个锐边上比如零件的一个直角边缘正好缺了一块你补完洞再光顺很容易把这条边磨成圆角。处理办法是识别出边界环上属于特征边的部分把这些顶点彻底锁死。识别方法可以是计算边界两侧相邻面的二面角超过阈值比如 30 度就标记为特征边。第二种是补丁内部的形状特征。有些洞虽然大但周围曲面有明显的凸脊或者凹谷补丁应该延续这种走向。这需要在 fairing 的能量项里加约束——不是单纯最小化双调和能量而是加上采样点法向与邻域法向接近的惩罚项。实际实现中可以用 RBF 拟合代替纯能量最小化因为 RBF 天生带边界法向约束。第三种是网格本身的密度分布。双调和光顺会把顶点往曲率低的地方聚导致补丁比周围稀疏。如果后续要做有限元分析网格质量是硬指标。我的做法是 fairing 之后做一次网格重划分remeshing或者干脆在 fairing 之前就把补丁的三角形密度对齐周围。这里有个取舍需要说清楚特征保持和光顺程度是矛盾的。约束越多补丁越硬越容易出现局部应力集中。我的经验是先做无约束的双调和光顺然后检查边界二面角分布如果发现某段边界的二面角异常明显偏离 180 度再对那一段做局部约束重解。这样既保证了整体光顺又保住了关键特征。3.4 参数怎么定迭代次数、步长、权重的经验值参数没有万能值但可以给一组经过验证的起始点然后再微调。下面这张表是我在不同项目里总结的经验范围。参数含义推荐范围经验取值调整方向λ伞形算子步长0.1 ~ 0.70.45出现振荡就调小迭代次数伞形松弛轮数10 ~ 10030以位移收敛为准细分次数补丁加密轮数1 ~ 43看补丁密度是否匹配周围α双调和权重0 ~ 11.01.0 即纯双调和特征边二面角阈值判定锐边20° ~ 45°30°特征多就调低收敛判据平均位移阈值1e-5 ~ 1e-3 × 对角线1e-4 × 对角线精度要求高就调小RBF 正则拟合松弛量1e-10 ~ 1e-61e-8边界噪声大就调大关于收敛判据我用的是所有自由顶点的平均位移。计算方式是 $\frac{1}{|F|}\sum_{i \in F}|v_i^{new} - v_i^{old}|$当它小于包围盒对角线的 1e-4 倍时停止迭代。这个阈值在大多数工业件上对应 0.01mm 到 0.05mm 的精度够用。关于 λ 的取值理论上前向欧拉格式要求 $\lambda 1$ 才收敛但工程上超过 0.7 就开始振荡了。我一般从 0.5 开始试如果看到补丁表面出现橘子皮一样的波纹就降到 0.3。用余切权重的时候因为权重矩阵已经做了归一化λ 的有效范围会变需要重新标定。还有一个容易忽略的点边界顶点是否允许移动。理论上边界顶点必须固定因为它们是原始网格的一部分。但如果边界本身有扫描噪声边界上的抖动会被传导到补丁内部。折中方案是允许边界顶点在切平面内移动但法向分量锁死。这样既消掉了抖动又保证了补丁与周围的连接。4. 一套能直接抄的实操流水线4.1 环境准备与数据预处理依赖就几个trimesh处理网格 IO 和基础拓扑numpy/scipy做数值计算scipy.sparse做稀疏求解。如果要做带约束的三角化再装一个triangle库。pip install trimesh numpy scipy triangle如果你更习惯 C 生态libigl里有现成的igl::min_faces、igl::harmonic、igl::bfs_orient性能比 Python 好一个量级。用 Meshlab 做交互式处理也行它的Filters Remeshing Close Holes和Filters Smoothing Laplacian Smooth覆盖了基本需求但参数不可控批处理不方便。预处理是整条流水线里最容易被低估的一步。我的顺序是加载模型processFalse保持原始拓扑不要让它自动做任何修改。统计包围盒对角线长度作为后续所有容差参数的基准。顶点合并容差取对角线的 1e-5。删除退化面面积小于对角线平方 1e-12 的和重复面。删除孤立顶点和孤立连通分量。检查并处理非流形边。做完这六步再统计剩下的孔洞数量和每个洞的边界点数心里就有数了。4.2 孔洞分类与策略分配拿到孔洞清单后别直接一股脑全跑同一个算法。按特征分流效率和效果都好得多。分类维度有三个边界点数、平面性、边界曲率。边界点数少于 30且拟合平面的残差小于平均边长的 5%走最小权重三角化快。边界点数 30 到 200平面性一般走平面参数化 RBF 拟合。边界点数超过 200或者平面残差很大先对边界做一维曲线光顺再做分区参数化逐块补。边界上有明显的锐角转折相邻边界边夹角小于 120 度标记为特征点fairing 时锁死。平面残差的计算很简单对边界点做 PCA取最小奇异值除以边界点数开方再除以平均边长就得到归一化残差。def classify_hole(pts, avg_edge_len): c pts.mean(axis0) u, s, vt np.linalg.svd(pts - c, full_matricesFalse) planar_res np.sqrt(s[2] ** 2 / len(pts)) / (avg_edge_len 1e-12) n len(pts) if n 30 and planar_res 0.05: return planar_small elif n 200 and planar_res 0.3: return param_fit else: return large_curved分好类之后小洞批量走快路径大洞单独处理。我在一个扫描件修复项目里1200 个洞里有 1100 个是小平面洞走快路径总共花了 3 秒剩下 100 个大洞走了 RBF花了 40 秒。如果全走 RBF光小洞就要几分钟。4.3 补洞环节的代码实现把前面的零件组装起来完整的补洞流程大致是这样。def fill_all_holes(mesh, max_planar_pts30, max_fit_pts200): V np.array(mesh.vertices) F np.array(mesh.faces) diag np.linalg.norm(mesh.bounds[1] - mesh.bounds[0]) boundary_edges find_boundary_edges(mesh) loops build_loops(boundary_edges, V) new_faces [] patched_idx [] # 记录所有新生成的顶点fairing 时只动这些 for loop in loops: pts V[loop] # 边界平均边长 d np.linalg.norm(np.diff(pts, axis0, appendpts[:1]), axis1) avg_len d.mean() kind classify_hole(pts, avg_len) if kind planar_small: faces, new_pts fill_planar(pts, diag) elif kind param_fit: faces, new_pts fill_with_rbf(pts, diag) else: faces, new_pts fill_large_curved(pts, diag) base len(V) V np.vstack([V, new_pts]) if len(new_pts) else V patched_idx list(range(base, len(V))) # faces 里的索引要偏移 for f in faces: new_faces.append([loop[i] if i len(loop) else base (i - len(loop)) for i in f]) F np.vstack([F, np.array(new_faces)]) if new_faces else F return V, F, patched_idx这里我做了个简化把边界顶点和新增顶点拼在一个索引空间里边界顶点用0 ~ len(loop)-1新增顶点用len(loop) ~ ...。实际写的时候要小心索引偏移这是最容易出 bug 的地方。我一般会在补洞之后立刻做一次连通性检查——每个新面的三条边都必须能找到配对的面否则索引错了。fill_planar就是前面讲的最小权重三角化fill_with_rbf是平面投影 Delaunay RBF 映射回三维。这里不展开完整代码了思路清晰之后都是体力活。4.4 fairing 收尾与质量校验补完洞之后把所有新增顶点作为自由集原始顶点作为固定集做一次双调和光顺。def finalize(mesh, V, F, patched_idx): L, _ build_cotan_laplacian(V, F) all_idx np.arange(len(V)) free_idx np.array(patched_idx, dtypeint) fixed_idx np.setdiff1d(all_idx, free_idx) V_faired biharmonic_fairing(V, L, free_idx, fixed_idx) new_mesh trimesh.Trimesh(verticesV_faired, facesF, processFalse) new_mesh.remove_degenerate_faces() return new_mesh光顺完成后必须做质量校验不能直接交付。我固定会看这几个指标水密性。再跑一次边界边检测边界边数量必须是 0。如果还有说明有洞没补上或者补丁自身有破面。法向连续性。统计补丁面与相邻原始面的二面角看平均值和最大值。平均值应该在 175 度以上最大值不应该低于 150 度除非那里本来就有特征边。体积变化。补洞前后计算封闭体积变化率应该小于 0.5%。超过说明光顺收缩太厉害要减小 λ 或者换双调和。网格质量。统计补丁区域的最差三角形——最小角、最长边与最短边之比。最小角低于 10 度或者长宽比超过 20说明网格质量差后续做仿真会有问题需要重划分。自交检查。补丁如果跟自身或者周围面片穿插切片时会出问题。可以用trimesh的is_watertight加上碰撞检测做粗略判断精细检查要用空间哈希找相交三角形对。这套校验跑一遍几十毫秒但能拦住 90% 的交付事故。我吃过一次亏补洞后没检查自交直接送去 3D 打印切片软件报了一堆错误回头查是补丁在凹角处翻折了。从那以后校验就成了固定动作。5. 常见问题与排查技巧实录5.1 高频故障速查表下面这些是我这些年真实遇到过、并且被问得最多的问题按现象整理成表遇到问题可以直接对号入座。现象可能原因排查方法解决办法孔洞数量比预期多出几倍顶点未合并浮点误差拆散了边统计顶点距离分布看是否有大量 1e-6 量级的近邻提高 merge 容差或按空间聚类合并边界环串接不闭合存在非流形边统计每条边被使用的面数找 count 2 的边切掉非流形区域或用方向连续性强制串接补丁表面有褶皱边界噪声被拟合进去检查边界点的二阶差分看是否有高频抖动先对边界做一维光顺再补洞补丁比周围低一块拉普拉斯光顺收缩测量补丁中心与拟合曲面的偏差改用双调和光顺或加体积约束fairing 后网格出现橘子皮λ 太大迭代振荡观察连续两次迭代的位移是否反向减小 λ 到 0.3 以下或用余切权重求解时间过长RBF 矩阵稠密规模过大打印矩阵维度超过 5000 就要优化换紧支撑核或分区求解补丁与原始面片自交洞形状高度非凸平面投影重叠检查投影后的二维三角化是否有重叠分区参数化或用 3D 直接三角化特征边被磨圆fairing 未锁定特征点对比补洞前后的二面角分布识别特征边并加入固定集切片器报非水密补丁与边界顶点未焊接检查补丁边界顶点坐标与原始边界是否一致用顶点索引直接复用不要重新计算坐标补洞后体积变化超过 1%光顺能量阶数太低计算补洞前后体积升到双调和或加体积守恒约束5.2 踩过的坑与独家心得表格能覆盖的是一般性问题还有一些是只有真做过才知道的经验我单独拎出来讲。第一个坑不要在补洞之前做全局平滑。我早期做项目的时候习惯先对扫描网格做一遍全局拉普拉斯平滑去噪再做补洞。结果发现补洞效果很差边界对不齐。原因是全局平滑会把边界顶点也一起移动了边界环变了形后续的三角化就失去了参考。正确顺序是先补洞、再局部平滑而且平滑范围要限制在补丁和它周围两环内。第二个坑细分和光顺的交替节奏。Liepa 的算法是三角化 → 细分 → 光顺 → 细分 → 光顺交替进行而不是先细分到底再光顺。原因在于一次性细分出来的顶点初始位置都很差一般在边的中点或者面重心光顺一次很难收敛。交替进行的话每一轮细分都是在已经光顺好的曲面上加密初始位置就更接近最终解。我一般细分三轮每轮之后光顺 10 到 15 次。第三个坑RBF 的法向约束比位置约束更重要。只用边界点的位置做 RBF 拟合补丁表面会有一阶不连续也就是有折痕。加上边界点的法向约束之后补丁和原始曲面在边界处才真正平滑过渡。实现上就是在线性系统里多加几行方程让拟合函数的梯度在边界点上等于原始法向。第四个坑不要迷信自动参数。现在很多工具都提供一键补洞和自适应光顺参数是自动算的。这些工具在标准测试集上表现很好但真实扫描数据千奇百怪。我遇到过一个大洞自动算法判断边界是平面的用了最快的最小权重三角化结果那个洞开在一个球面上补丁塌进去一大块。我的做法是先用自动算法跑一遍看结果然后人为检查大洞边界点超过 100 的的补丁曲率是否合理不合理就手工指定策略。第五个坑网格的坐标尺度会影响参数。同一个模型用毫米做单位和对角线归一化到 1 再做光是容差参数就差了三个数量级。我的习惯是加载后立刻归一化到单位包围盒所有操作在归一化空间里做输出的时候再乘回去。这样参数的通用性好得多一套参数能适配各种尺度的模型。第六个坑图形学库的默认行为会坑你。trimesh加载 STL 时默认会做顶点合并有些情况下你不想让它合。open3d的TriangleMesh在计算法向时会自动统一朝向如果你的模型本来就是内外翻转的统一之后就全乱了。用任何库之前先看一遍它的默认参数尤其是process、merge、orient这类开关。5.3 关于收敛性的一点补充最后说一个容易被忽略的点双调和光顺求解的是一个线性系统理论上一次求解就到位。但如果你在自由顶点集上加了很强的正则项或者自由顶点集被分割成多个互不连通的区域求解结果可能不满足全局的光顺条件。这时候需要做一次迭代细化——把前一次求解的结果作为初值重新求解一次。一般两到三次就能收敛。判断是否收敛的方法很简单把两次求解的顶点位移取平均如果小于 1e-6 倍的包围盒对角线就停。这个判据比看能量值直观得多因为能量值的量级跟模型尺度强相关不好定阈值。另外如果你的补丁特别大自由顶点超过五万稀疏求解也会慢。这时候可以用多重网格multigrid的思路先在粗网格上求解再插值到细网格上做初值能提速好几倍。不过说实话工业场景里单次补洞的自由顶点超过五万的很少见真遇到了更合理的做法是先把洞分区而不是死磕求解器。 三维扫描网格修复：hole filling 与 fairing 算法实战
做三维扫描、逆向工程或者 3D 打印前处理的人大概都有过这样的经历模型渲染出来看着挺完整一转到线框模式或者丢进切片软件立马就露馅了——一个个窟窿张着嘴边缘还毛刺刺的。补洞hole filling和 fairing 算法这两个词基本就是绕不过去的坎。前者负责把缺失的面片填回去让模型重新闭合成水密实体后者负责把填完之后的补丁磨得跟周围曲面连成一片不留折痕。这不是什么锦上添花的功能而是从扫描数据到可用模型的必经环节做有限元仿真、做增材制造、做数字孪生资产入库第一步都是这个。我把这两个算法放在一起讲是因为单独做任何一个都容易出问题。只补不光顺补丁像贴了块膏药法向突变会让切片器报错、让仿真应力集中只光顺不补洞边界会越磨越缩洞反而变大。下面这套东西是我这几年在扫描件修复、模具逆向、文物数字化几个项目里反复折腾出来的代码基于 Python 的 trimesh scipy 生态原理讲清楚参数给具体值踩过的坑也都写进去了。不管你是刚接手网格处理任务的新手还是已经用过 Meshlab 自动补洞但效果不理想的老手应该都能找到能直接用的部分。1. 补洞与 fairing 到底在修什么从破洞的成因说起1.1 网格上那些洞是怎么来的很多人以为洞是软件算错了其实绝大多数洞是物理世界的信息缺失算法只是把缺失暴露出来。搞明白成因才能判断应该用哪种补法。我把常见成因归成四类每一类对应的补洞策略其实不一样。第一类是扫描遮挡。激光或者结构光扫描仪本质上是个视线传感器光打不到的地方就没数据。深腔、倒扣、叶片间的缝隙、螺纹根部这些位置永远扫不全。这类洞的特点是边界往往沿着视线方向拉长形状不规则而且洞的周围通常有噪声。第二类是表面材质问题。高反光、全透明、纯黑的表面扫描时要么过曝要么无返回。这类洞通常成片出现边界比较圆滑周围点云密度均匀。第三类是几何操作遗留。布尔运算、抽壳、镜像、网格简化每一步都可能产生退化面片。尤其是模型简化的抽取操作会把细长三角形去掉留下一个狭长裂缝。这类洞的特点是呈线状两侧边界几乎平行。第四类是文件转换丢失。STL、OBJ、PLY 之间来回倒精度截断导致顶点对不上原本共享的顶点被拆成两个边就变成了边界边。这类洞最冤明明模型是好的只是顶点没有合并。提示拿到一个有洞的模型先别急着跑补洞。做一次顶点合并merge by distance和退化面清理往往能消掉一半以上的假洞。我见过太多案例几十个所谓的洞合并一次顶点就只剩三个了。1.2 fairing 不是简单的磨皮fairing 这个词在中文里通常译作光顺但很多人把它理解成平滑这是两个层次的事。平滑是去掉高频噪声光顺是让曲面满足某种能量最小的条件——数学上就是让曲面的二阶或者三阶导数连续。放到补洞场景里为什么不能只做平滑因为补丁是凭空生成的它和周围原始曲面之间天然存在缝隙。你如果只做局部平滑边界处的法向不连续会一直在。你需要的是让补丁区域作为一个整体去承接周围曲面的曲率趋势让交界的二面角接近 180 度。打个比方补洞相当于在一块布料上剪掉一个洞然后用另一块布补上去。直接缝上去接缝处一定有褶子。fairing 做的就是让新补的布按照周围布的纹理和张力自然延展接缝摸不出来。从数学上看fairing 求解的是一个带约束的能量最小化问题。约束是补丁边界上的点必须落在原网格边界上能量是曲面的某种形变度量。不同的能量选择得到的补丁形状完全不同这也是为什么 fairing 算法有那么多种。1.3 为什么这两件事得连着做分开做的问题在于信息不对称。补洞算法如果不考虑后续 fairing它会把洞填成一个平面片——因为平面三角化最简单。但这个平面片跟周围曲面完全不搭fairing 阶段要么磨不平要么一磨就把边界也带跑了。反过来fairing 算法如果不了解补洞的结构它会把整个网格当成一个整体去光顺结果就是原始模型上那些本来该保留的锐边、特征线全被磨圆了。做机械件的逆向时这种错误是灾难性的一个倒角被磨掉装配就出问题了。所以实际流程应该是补洞阶段就引入曲率约束让初始补丁尽可能贴近周围曲面的趋势fairing 阶段只对新增的自由顶点操作把原始边界上的顶点锁死或者只允许它们在切平面内移动。这样两个阶段互相配合原始几何一点不动补丁区域又自然过渡。我在做模具修复项目的时候一开始就是分开跑的补洞用现成工具fairing 用拉普拉斯平滑结果零件上的型面分界线全糊了。后来改成一体化处理边界顶点固定只对内部顶点做双调和光顺效果立竿见影。2. 补洞算法的路线选型与背后的取舍2.1 三条主流路线直接三角化、参数化拟合、数据驱动补洞算法按复杂度大致分三档选哪档取决于洞的大小、形状、以及你对结果的要求。第一档直接三角化。把边界环当成一个多边形在三维空间里直接做三角剖分。代表算法是耳切法ear clipping和最小权重三角化minimum weight triangulation。优点是快几百个顶点的洞毫秒级出结果缺点也明显生成的补丁基本是个平面或者鞍面完全没有曲率延拓的概念。适合什么适合小洞边界顶点在 50 个以内、且洞本身近似平面的情况。螺丝孔、小裂纹、扫描小缺陷这一档够用。第二档参数化拟合。先把边界环投影到一个平面域上在二维平面里做带约束的三角化再把平面上的点映射回三维最后用周围曲面的信息去调整新增顶点的位置。这一档是工业级的主流做法。核心思想是降维求解把复杂的三维曲面问题变成二维的平面问题。代表实现就是 Liepa 的算法最小面积三角化打底然后细分加光顺。第三档数据驱动。用深度学习模型直接从残缺网格或者点云预测缺失部分。这几年在扫描补全领域有不少工作优点是能补出语义合理的形状——比如一个椅子的缺角模型知道该补出椅子腿的形状而不只是一个平面片。缺点是需要训练数据泛化能力依赖数据分布工业件上不太敢用。我的建议是批量处理同类资产比如一批同款家具扫描件可以试试单件精密零件还是老老实实走第二档。下面这张表是我对不同路线的实际体感总结供选型参考路线适用洞尺寸曲率延拓单洞耗时实现难度典型场景耳切/最小权重三角化边界点 50无毫秒级低小裂纹、螺孔、扫描缺陷平面参数化 拟合边界点 50~2000有秒级中扫描件大洞、模具型面3D 直接拟合RBF/MLS边界点 500有秒级到分钟级中高曲面件、局部缺失数据驱动补全不限语义级秒级高批量同类资产2.2 边界环提取补洞的地基工程补洞第一步永远是找洞。而找洞的本质是找那些只被一个面片使用的边。在一个封闭流形网格里每条边恰好被两个三角形共享一旦某条边只被一个三角形使用它就是边界边。道理简单实操坑很多。第一个坑是顶点没有合并。两个顶点坐标差 1e-6浮点数认为它们是两个点那条边就变成了两条边界边洞被拆得粉碎。第二个坑是非流形边一条边被三个甚至更多面共享这种网格严格来说不是二维流形边界环的串接逻辑会乱掉。我的处理顺序是先合并顶点容差取包围盒对角线的 1e-5 量级再删退化面面积接近零的再删重复面最后才提取边界。import numpy as np import trimesh def clean_mesh(path, merge_tol_ratio1e-5): mesh trimesh.load(path, processFalse) diag np.linalg.norm(mesh.bounds[1] - mesh.bounds[0]) # 合并重合顶点这一步能消掉大部分假洞 mesh.merge_vertices(digits_vertexint(np.ceil(-np.log10(merge_tol_ratio * diag)))) mesh.remove_degenerate_faces() mesh.remove_duplicate_faces() mesh.remove_unreferenced_vertices() return mesh def find_boundary_edges(mesh): # edges_sorted 是已排序的顶点索引对方便做唯一性统计 uniq, counts np.unique(mesh.edges_sorted, axis0, return_countsTrue) return uniq[counts 1]拿到边界边之后要做的是把它们串成闭环。串接的逻辑是从任意一个端点出发沿着邻接关系一路走直到回到起点。这里要处理顶点度数大于 2 的情况——数学上应该在非流形处断开工程上我通常按选择转弯角度最小的下一条边来走这样得到的环更符合直觉。from collections import defaultdict def build_loops(boundary_edges, vertices, min_len3): adj defaultdict(list) for a, b in boundary_edges: adj[a].append(b) adj[b].append(a) visited_edge set() loops [] for a, b in boundary_edges: key (min(a, b), max(a, b)) if key in visited_edge: continue loop [a, b] visited_edge.add(key) prev, cur a, b while True: cands [n for n in adj[cur] if n ! prev and (min(cur, n), max(cur, n)) not in visited_edge] if not cands: break if len(cands) 1: # 非流形分叉挑方向延续性最好的那条 d0 vertices[cur] - vertices[prev] d0 d0 / (np.linalg.norm(d0) 1e-12) best, best_dot None, -2 for n in cands: d1 vertices[n] - vertices[cur] d1 d1 / (np.linalg.norm(d1) 1e-12) if np.dot(d0, d1) best_dot: best_dot, best np.dot(d0, d1), n nxt best else: nxt cands[0] visited_edge.add((min(cur, nxt), max(cur, nxt))) if nxt loop[0]: break loop.append(nxt) prev, cur cur, nxt if len(loop) min_len: loops.append(loop) return loops注意这段代码在遇到一个顶点连着 4 条以上边界边的病态区域时仍然可能给出奇怪的结果。我的经验是先把模型丢进 Meshlab 看一眼非流形边的分布如果集中在一个小区域宁可手工切掉那块再补也不要指望算法自动处理。2.3 小洞快修最小权重三角化的实现细节对小洞我基本都用最小权重三角化。它的思路是在边界环构成的多边形上做三角剖分但每条候选对角线都有代价选总代价最小的那套剖分。代价函数一般包含两项——三角形面积以及三角形法向与边界邻域法向的偏差。为什么加法向项因为纯面积最小化会倾向于生成扁平的剖分对非平面的洞这会导致补丁凹陷。加上法向惩罚后算法会更愿意生成跟周围曲面朝向一致的三角形。Liepa 那篇 2003 年的论文里给的权重是$$\omega(f) \frac{1}{2}\left(\sum_{e \in f} |e|^2\right) \cdot \max_{e \in f}\left(\frac{1 - \cos \theta_e}{2}\right)$$其中 $|e|$ 是边长$\theta_e$ 是这条边与相邻边界边的夹角。这个式子的直观含义是面积大的三角形代价高与边界走向偏差大的三角形代价也高。两项相乘逼着算法选出既紧凑又顺滑的剖分。实现上用动态规划。设边界环上有 $n$ 个顶点定义 $M[i][j]$ 为从第 $i$ 个顶点到第 $j$ 个顶点这一段子多边形的最小代价。递推就是枚举中间的切分点 $k$取 $M[i][k] M[k][j] \omega(i,k,j)$ 的最小值。复杂度 $O(n^3)$$n$ 在 200 以内都很快。def min_weight_triangulation(pts, normals, edges_dir): n len(pts) INF float(inf) M np.full((n, n), INF) choice np.full((n, n), -1, dtypeint) def tri_cost(i, k, j): a, b, c pts[i], pts[k], pts[j] area 0.5 * np.linalg.norm(np.cross(b - a, c - a)) # 法向与边界走向的偏差惩罚 n_face np.cross(b - a, c - a) norm np.linalg.norm(n_face) if norm 1e-12: return INF n_face / norm local_n (normals[i] normals[k] normals[j]) / 3.0 local_n / (np.linalg.norm(local_n) 1e-12) cos_dev np.clip(np.dot(n_face, local_n), -1, 1) return area * (1.0 - cos_dev) area * 1e-6 for i in range(n - 1): M[i][i 1] 0.0 for span in range(2, n): for i in range(0, n - span): j i span best INF best_k -1 for k in range(i 1, j): if M[i][k] INF or M[k][j] INF: continue c M[i][k] M[k][j] tri_cost(i, k, j) if c best: best, best_k c, k M[i][j] best choice[i][j] best_k return M, choice这段代码只处理单个环实际项目里要循环处理所有环并且对每个环单独建点表。另外要注意点的顺序——边界环的顶点顺序必须是沿着边界连续排列的反了会让三角形法向整体翻转。2.4 大洞慢补参数化与曲面拟合洞的边界点超过 100 个或者洞的形状明显非平面最小权重三角化出来的补丁就开始鼓包或者塌陷了。这时候要上参数化拟合。核心流程分四步走。第一步求边界环的最优拟合平面用 PCA 就能做对边界点做中心化然后 SVD 分解最小奇异值对应的右奇异向量就是平面法向。第二步把边界点投影到这个平面上得到一组二维点。第三步在二维域内做带约束的 Delaunay 三角化注意如果边界是凹的普通 Delaunay 会跑到外面去要用带约束的版本比如triangle库或者自己实现凸分解。第四步把三角化结果映射回三维然后对内部新增顶点做位置优化。第四步是整个流程的灵魂。最朴素的做法是把平面上的点按原位置映射回三维得到的补丁是个平面。好一点的做法是用周围曲面的法向信息做外推让补丁顺着曲率走。更好的做法是求解一个径向基函数RBF拟合问题。RBF 的思路是找出一个函数 $f(x) \sum_i w_i \phi(|x - x_i|) p(x)$让它精确通过边界点及其法向约束然后在洞内采样。$\phi$ 通常取薄板样条核 $\phi(r) r^2 \log r$这个核对应的是薄板弯曲能量最小天然就是光顺的。$\phi(r) r^3$ 也常用对应的是线性化后的薄板能量。import numpy as np from scipy.spatial.distance import cdist def fit_rbf_surface(boundary_pts, interior_uv, kerneltps, reg1e-8): # boundary_pts: (m,3) 边界点 # interior_uv: (k,2) 洞内待求解的点在平面域上的坐标 # 返回洞内点的三维坐标估计 m len(boundary_pts) D_bb cdist(boundary_pts, boundary_pts) if kernel tps: def phi(r): r np.where(r 1e-12, 1e-12, r) return r * r * np.log(r) else: def phi(r): return r ** 3 K phi(D_bb) reg * np.eye(m) # 加低阶多项式项保证线性精度 P np.hstack([np.ones((m, 1)), boundary_pts]) A np.zeros((m 4, m 4)) A[:m, :m] K A[:m, m:] P A[m:, :m] P.T rhs np.zeros(m 4) rhs[:m] boundary_pts[:, 0] # 逐轴求解 w np.linalg.solve(A, rhs) return w这段代码只是骨架实际使用时要对 x、y、z 三个轴分别求解还要注意 RBF 是全局的边界点一多矩阵就巨大。我的经验是边界点在 500 以内直接稠密求解没问题超过就上紧支撑核比如 Wendland 核或者分区求解。说到分区这里有个实用技巧如果洞特别大先对边界环做一次一维的曲线光顺把边界当成一条闭合曲线用拉普拉斯算子平滑把扫描噪声引起的边界抖动消掉再做参数化。否则那些高频抖动会被 RBF 原样拟合进去补丁表面会起波纹。3. fairing 算法的数学底子与实现要点3.1 伞形算子与拉普拉斯光顺补完洞只是第一步新生成的顶点位置还很粗糙。最经典的光顺算法是伞形算子umbrella operator本质上就是离散的拉普拉斯算子。对顶点 $v_i$它的伞形算子定义为邻域顶点平均值与自身位置之差$$U(v_i) \frac{1}{\sum_j w_{ij}} \sum_{j \in N(i)} w_{ij}(v_j - v_i)$$更新规则是 $v_i \leftarrow v_i \lambda U(v_i)$$\lambda$ 一般取 0.3 到 0.6。这个操作重复若干次网格就会越来越光滑。权重 $w_{ij}$ 有两种常见取法。均匀权重 $w_{ij} 1$实现最简单但网格疏密不均的时候会出问题——顶点密的地方贡献大结果把点往疏的地方拉。余切权重cotangent weight考虑了几何形状$$w_{ij} \frac{1}{2}(\cot \alpha_{ij} \cot \beta_{ij})$$其中 $\alpha_{ij}$、$\beta_{ij}$ 是边 $ij$ 所对的两个角。这个权重有理论保证在规则三角网上用余切权重的拉普拉斯算子能精确逼近连续曲面的拉普拉斯-贝尔特拉米算子。def build_cotan_laplacian(V, F): import scipy.sparse as sp n len(V) I, J, Wdata [], [], [] for tri in F: for k in range(3): i, j, o tri[k], tri[(k 1) % 3], tri[(k 2) % 3] v1 V[i] - V[o] v2 V[j] - V[o] cross np.cross(v1, v2) area2 np.linalg.norm(cross) if area2 1e-14: cot 0.0 else: cot np.dot(v1, v2) / area2 w 0.5 * cot I [i, j] J [j, i] Wdata [w, w] W sp.coo_matrix((Wdata, (I, J)), shape(n, n)).tocsr() # 对角元取行和的相反数得到标准拉普拉斯 d np.array(W.sum(axis1)).flatten() L sp.diags(d) - W return L, W注意余切权重在钝角三角形上会出现负值。钝角大于 90 度的余切是负的这会导致拉普拉斯算子失去正定性光顺结果发散。处理方法有两个一是把负权重截断到 0二是用虚拟边界技巧在钝角顶点处补一个辅助点。我一般用截断简单有效代价是精度略降。网格质量好的话根本遇不到这个问题。3.2 从拉普拉斯到双调和薄板能量为什么更靠谱纯拉普拉斯光顺有个致命缺陷它最小化的是一阶能量 $\int |\nabla v|^2$得到的曲面是调和曲面。调和曲面有个特性——它会收缩。一个球面被反复做拉普拉斯光顺最终会缩成一个点。放在补洞场景里就是补丁会往洞中心缩补完的曲面比周围低一块。解决思路是提高能量阶数。双调和能量 $\int |\Delta v|^2$ 对应的欧拉-拉格朗日方程是 $\Delta^2 v 0$解叫双调和曲面它保证了曲率连续同时不会像调和曲面那样有强烈的收缩倾向。物理上这对应薄板弯曲能量最小——一块薄钢板被固定在边界上它自然形成的形状就是双调和曲面。离散化之后最小化 $|L x|^2$ 的问题可以写成线性系统。把顶点分成自由顶点补丁内部的和固定顶点边界上的求解$$\min |L x|^2 \quad \text{s.t. } x_{fixed} b$$对自由部分求导得正规方程 $L_{ff}^T L_{ff} , x_f -L_{ff}^T L_{fx} , b$。这是个稀疏线性系统用spsolve直接解几千个自由顶点毫无压力。import scipy.sparse as sp import scipy.sparse.linalg as spla def biharmonic_fairing(V, L, free_idx, fixed_idx, tikhonov1e-9): n len(V) A (L.T L).tocsc() if tikhonov 0: A A tikhonov * sp.eye(n, formatcsc) A_ff A[free_idx][:, free_idx] A_fx A[free_idx][:, fixed_idx] rhs -(A_fx V[fixed_idx]) V_new V.copy() sol spla.spsolve(A_ff.tocsc(), rhs) V_new[free_idx] sol return V_new这里加 Tikhonov 正则项的原因是$L$ 的零空间是常数向量理论上 $L^T L$ 的零空间也是常数向量只要自由顶点集不构成整个连通分量$A_{ff}$ 就是可逆的。但数值上接近奇异加一点点正则更保险1e-9 量级对结果几乎没有影响。双调和光顺的效果比拉普拉斯好太多。我做过一个对比测试同样一个 80 个边界点的圆洞拉普拉斯光顺 50 次后补丁中心比周围曲面低 0.3mm模型尺寸 100mm双调和一次求解后偏差只有 0.02mm。而且双调和不需要迭代一次求解直接到位。3.3 特征保持怎么磨而不圆fairing 最大的风险是把该保留的特征磨没了。具体到补洞场景有三种特征需要保护。第一种是补丁边界本身。如果洞开在一个锐边上比如零件的一个直角边缘正好缺了一块你补完洞再光顺很容易把这条边磨成圆角。处理办法是识别出边界环上属于特征边的部分把这些顶点彻底锁死。识别方法可以是计算边界两侧相邻面的二面角超过阈值比如 30 度就标记为特征边。第二种是补丁内部的形状特征。有些洞虽然大但周围曲面有明显的凸脊或者凹谷补丁应该延续这种走向。这需要在 fairing 的能量项里加约束——不是单纯最小化双调和能量而是加上采样点法向与邻域法向接近的惩罚项。实际实现中可以用 RBF 拟合代替纯能量最小化因为 RBF 天生带边界法向约束。第三种是网格本身的密度分布。双调和光顺会把顶点往曲率低的地方聚导致补丁比周围稀疏。如果后续要做有限元分析网格质量是硬指标。我的做法是 fairing 之后做一次网格重划分remeshing或者干脆在 fairing 之前就把补丁的三角形密度对齐周围。这里有个取舍需要说清楚特征保持和光顺程度是矛盾的。约束越多补丁越硬越容易出现局部应力集中。我的经验是先做无约束的双调和光顺然后检查边界二面角分布如果发现某段边界的二面角异常明显偏离 180 度再对那一段做局部约束重解。这样既保证了整体光顺又保住了关键特征。3.4 参数怎么定迭代次数、步长、权重的经验值参数没有万能值但可以给一组经过验证的起始点然后再微调。下面这张表是我在不同项目里总结的经验范围。参数含义推荐范围经验取值调整方向λ伞形算子步长0.1 ~ 0.70.45出现振荡就调小迭代次数伞形松弛轮数10 ~ 10030以位移收敛为准细分次数补丁加密轮数1 ~ 43看补丁密度是否匹配周围α双调和权重0 ~ 11.01.0 即纯双调和特征边二面角阈值判定锐边20° ~ 45°30°特征多就调低收敛判据平均位移阈值1e-5 ~ 1e-3 × 对角线1e-4 × 对角线精度要求高就调小RBF 正则拟合松弛量1e-10 ~ 1e-61e-8边界噪声大就调大关于收敛判据我用的是所有自由顶点的平均位移。计算方式是 $\frac{1}{|F|}\sum_{i \in F}|v_i^{new} - v_i^{old}|$当它小于包围盒对角线的 1e-4 倍时停止迭代。这个阈值在大多数工业件上对应 0.01mm 到 0.05mm 的精度够用。关于 λ 的取值理论上前向欧拉格式要求 $\lambda 1$ 才收敛但工程上超过 0.7 就开始振荡了。我一般从 0.5 开始试如果看到补丁表面出现橘子皮一样的波纹就降到 0.3。用余切权重的时候因为权重矩阵已经做了归一化λ 的有效范围会变需要重新标定。还有一个容易忽略的点边界顶点是否允许移动。理论上边界顶点必须固定因为它们是原始网格的一部分。但如果边界本身有扫描噪声边界上的抖动会被传导到补丁内部。折中方案是允许边界顶点在切平面内移动但法向分量锁死。这样既消掉了抖动又保证了补丁与周围的连接。4. 一套能直接抄的实操流水线4.1 环境准备与数据预处理依赖就几个trimesh处理网格 IO 和基础拓扑numpy/scipy做数值计算scipy.sparse做稀疏求解。如果要做带约束的三角化再装一个triangle库。pip install trimesh numpy scipy triangle如果你更习惯 C 生态libigl里有现成的igl::min_faces、igl::harmonic、igl::bfs_orient性能比 Python 好一个量级。用 Meshlab 做交互式处理也行它的Filters Remeshing Close Holes和Filters Smoothing Laplacian Smooth覆盖了基本需求但参数不可控批处理不方便。预处理是整条流水线里最容易被低估的一步。我的顺序是加载模型processFalse保持原始拓扑不要让它自动做任何修改。统计包围盒对角线长度作为后续所有容差参数的基准。顶点合并容差取对角线的 1e-5。删除退化面面积小于对角线平方 1e-12 的和重复面。删除孤立顶点和孤立连通分量。检查并处理非流形边。做完这六步再统计剩下的孔洞数量和每个洞的边界点数心里就有数了。4.2 孔洞分类与策略分配拿到孔洞清单后别直接一股脑全跑同一个算法。按特征分流效率和效果都好得多。分类维度有三个边界点数、平面性、边界曲率。边界点数少于 30且拟合平面的残差小于平均边长的 5%走最小权重三角化快。边界点数 30 到 200平面性一般走平面参数化 RBF 拟合。边界点数超过 200或者平面残差很大先对边界做一维曲线光顺再做分区参数化逐块补。边界上有明显的锐角转折相邻边界边夹角小于 120 度标记为特征点fairing 时锁死。平面残差的计算很简单对边界点做 PCA取最小奇异值除以边界点数开方再除以平均边长就得到归一化残差。def classify_hole(pts, avg_edge_len): c pts.mean(axis0) u, s, vt np.linalg.svd(pts - c, full_matricesFalse) planar_res np.sqrt(s[2] ** 2 / len(pts)) / (avg_edge_len 1e-12) n len(pts) if n 30 and planar_res 0.05: return planar_small elif n 200 and planar_res 0.3: return param_fit else: return large_curved分好类之后小洞批量走快路径大洞单独处理。我在一个扫描件修复项目里1200 个洞里有 1100 个是小平面洞走快路径总共花了 3 秒剩下 100 个大洞走了 RBF花了 40 秒。如果全走 RBF光小洞就要几分钟。4.3 补洞环节的代码实现把前面的零件组装起来完整的补洞流程大致是这样。def fill_all_holes(mesh, max_planar_pts30, max_fit_pts200): V np.array(mesh.vertices) F np.array(mesh.faces) diag np.linalg.norm(mesh.bounds[1] - mesh.bounds[0]) boundary_edges find_boundary_edges(mesh) loops build_loops(boundary_edges, V) new_faces [] patched_idx [] # 记录所有新生成的顶点fairing 时只动这些 for loop in loops: pts V[loop] # 边界平均边长 d np.linalg.norm(np.diff(pts, axis0, appendpts[:1]), axis1) avg_len d.mean() kind classify_hole(pts, avg_len) if kind planar_small: faces, new_pts fill_planar(pts, diag) elif kind param_fit: faces, new_pts fill_with_rbf(pts, diag) else: faces, new_pts fill_large_curved(pts, diag) base len(V) V np.vstack([V, new_pts]) if len(new_pts) else V patched_idx list(range(base, len(V))) # faces 里的索引要偏移 for f in faces: new_faces.append([loop[i] if i len(loop) else base (i - len(loop)) for i in f]) F np.vstack([F, np.array(new_faces)]) if new_faces else F return V, F, patched_idx这里我做了个简化把边界顶点和新增顶点拼在一个索引空间里边界顶点用0 ~ len(loop)-1新增顶点用len(loop) ~ ...。实际写的时候要小心索引偏移这是最容易出 bug 的地方。我一般会在补洞之后立刻做一次连通性检查——每个新面的三条边都必须能找到配对的面否则索引错了。fill_planar就是前面讲的最小权重三角化fill_with_rbf是平面投影 Delaunay RBF 映射回三维。这里不展开完整代码了思路清晰之后都是体力活。4.4 fairing 收尾与质量校验补完洞之后把所有新增顶点作为自由集原始顶点作为固定集做一次双调和光顺。def finalize(mesh, V, F, patched_idx): L, _ build_cotan_laplacian(V, F) all_idx np.arange(len(V)) free_idx np.array(patched_idx, dtypeint) fixed_idx np.setdiff1d(all_idx, free_idx) V_faired biharmonic_fairing(V, L, free_idx, fixed_idx) new_mesh trimesh.Trimesh(verticesV_faired, facesF, processFalse) new_mesh.remove_degenerate_faces() return new_mesh光顺完成后必须做质量校验不能直接交付。我固定会看这几个指标水密性。再跑一次边界边检测边界边数量必须是 0。如果还有说明有洞没补上或者补丁自身有破面。法向连续性。统计补丁面与相邻原始面的二面角看平均值和最大值。平均值应该在 175 度以上最大值不应该低于 150 度除非那里本来就有特征边。体积变化。补洞前后计算封闭体积变化率应该小于 0.5%。超过说明光顺收缩太厉害要减小 λ 或者换双调和。网格质量。统计补丁区域的最差三角形——最小角、最长边与最短边之比。最小角低于 10 度或者长宽比超过 20说明网格质量差后续做仿真会有问题需要重划分。自交检查。补丁如果跟自身或者周围面片穿插切片时会出问题。可以用trimesh的is_watertight加上碰撞检测做粗略判断精细检查要用空间哈希找相交三角形对。这套校验跑一遍几十毫秒但能拦住 90% 的交付事故。我吃过一次亏补洞后没检查自交直接送去 3D 打印切片软件报了一堆错误回头查是补丁在凹角处翻折了。从那以后校验就成了固定动作。5. 常见问题与排查技巧实录5.1 高频故障速查表下面这些是我这些年真实遇到过、并且被问得最多的问题按现象整理成表遇到问题可以直接对号入座。现象可能原因排查方法解决办法孔洞数量比预期多出几倍顶点未合并浮点误差拆散了边统计顶点距离分布看是否有大量 1e-6 量级的近邻提高 merge 容差或按空间聚类合并边界环串接不闭合存在非流形边统计每条边被使用的面数找 count 2 的边切掉非流形区域或用方向连续性强制串接补丁表面有褶皱边界噪声被拟合进去检查边界点的二阶差分看是否有高频抖动先对边界做一维光顺再补洞补丁比周围低一块拉普拉斯光顺收缩测量补丁中心与拟合曲面的偏差改用双调和光顺或加体积约束fairing 后网格出现橘子皮λ 太大迭代振荡观察连续两次迭代的位移是否反向减小 λ 到 0.3 以下或用余切权重求解时间过长RBF 矩阵稠密规模过大打印矩阵维度超过 5000 就要优化换紧支撑核或分区求解补丁与原始面片自交洞形状高度非凸平面投影重叠检查投影后的二维三角化是否有重叠分区参数化或用 3D 直接三角化特征边被磨圆fairing 未锁定特征点对比补洞前后的二面角分布识别特征边并加入固定集切片器报非水密补丁与边界顶点未焊接检查补丁边界顶点坐标与原始边界是否一致用顶点索引直接复用不要重新计算坐标补洞后体积变化超过 1%光顺能量阶数太低计算补洞前后体积升到双调和或加体积守恒约束5.2 踩过的坑与独家心得表格能覆盖的是一般性问题还有一些是只有真做过才知道的经验我单独拎出来讲。第一个坑不要在补洞之前做全局平滑。我早期做项目的时候习惯先对扫描网格做一遍全局拉普拉斯平滑去噪再做补洞。结果发现补洞效果很差边界对不齐。原因是全局平滑会把边界顶点也一起移动了边界环变了形后续的三角化就失去了参考。正确顺序是先补洞、再局部平滑而且平滑范围要限制在补丁和它周围两环内。第二个坑细分和光顺的交替节奏。Liepa 的算法是三角化 → 细分 → 光顺 → 细分 → 光顺交替进行而不是先细分到底再光顺。原因在于一次性细分出来的顶点初始位置都很差一般在边的中点或者面重心光顺一次很难收敛。交替进行的话每一轮细分都是在已经光顺好的曲面上加密初始位置就更接近最终解。我一般细分三轮每轮之后光顺 10 到 15 次。第三个坑RBF 的法向约束比位置约束更重要。只用边界点的位置做 RBF 拟合补丁表面会有一阶不连续也就是有折痕。加上边界点的法向约束之后补丁和原始曲面在边界处才真正平滑过渡。实现上就是在线性系统里多加几行方程让拟合函数的梯度在边界点上等于原始法向。第四个坑不要迷信自动参数。现在很多工具都提供一键补洞和自适应光顺参数是自动算的。这些工具在标准测试集上表现很好但真实扫描数据千奇百怪。我遇到过一个大洞自动算法判断边界是平面的用了最快的最小权重三角化结果那个洞开在一个球面上补丁塌进去一大块。我的做法是先用自动算法跑一遍看结果然后人为检查大洞边界点超过 100 的的补丁曲率是否合理不合理就手工指定策略。第五个坑网格的坐标尺度会影响参数。同一个模型用毫米做单位和对角线归一化到 1 再做光是容差参数就差了三个数量级。我的习惯是加载后立刻归一化到单位包围盒所有操作在归一化空间里做输出的时候再乘回去。这样参数的通用性好得多一套参数能适配各种尺度的模型。第六个坑图形学库的默认行为会坑你。trimesh加载 STL 时默认会做顶点合并有些情况下你不想让它合。open3d的TriangleMesh在计算法向时会自动统一朝向如果你的模型本来就是内外翻转的统一之后就全乱了。用任何库之前先看一遍它的默认参数尤其是process、merge、orient这类开关。5.3 关于收敛性的一点补充最后说一个容易被忽略的点双调和光顺求解的是一个线性系统理论上一次求解就到位。但如果你在自由顶点集上加了很强的正则项或者自由顶点集被分割成多个互不连通的区域求解结果可能不满足全局的光顺条件。这时候需要做一次迭代细化——把前一次求解的结果作为初值重新求解一次。一般两到三次就能收敛。判断是否收敛的方法很简单把两次求解的顶点位移取平均如果小于 1e-6 倍的包围盒对角线就停。这个判据比看能量值直观得多因为能量值的量级跟模型尺度强相关不好定阈值。另外如果你的补丁特别大自由顶点超过五万稀疏求解也会慢。这时候可以用多重网格multigrid的思路先在粗网格上求解再插值到细网格上做初值能提速好几倍。不过说实话工业场景里单次补洞的自由顶点超过五万的很少见真遇到了更合理的做法是先把洞分区而不是死磕求解器。