牙齿STL网格分割实战：投影栅格化与牙龈外轮廓提取

简介面向牙科数字化诊断与三维建模开发者的牙齿STL网格模型分割算法资料以投影算法曲面栅格化为核心解决牙齿与牙龈分离、牙龈外轮廓计算等问题。压缩包共54个文件大小18.37MB主体为40个MATLAB脚本.m涵盖STL读取、范围图像生成、边缘检测、区域分割、轮廓拟合等环节另含4个Java辅助程序、4个MAT数据文件以及项目配置和说明文档便于理解工程结构。已有594人学习。资料包含teeth_segmentation-master全部源码、示例数据与运行主程序提供inner_outer_face_filter、calc_edge_main、fitContour等关键函数实现可帮助研究者复现投影分割流程并针对牙科CT/扫描模型进行二次开发与算法优化。通过阅读代码可掌握栅格化预处理、阈值分割、空洞填充及外轮廓追踪等具体做法适合具备一定MATLAB基础的图形学或医疗影像学习者参考。1. 牙齿STL网格模型分割投影栅格化为什么比几何切割更适合牙颌模型做数字化正畸、隐形牙套设计或者种植导板规划的人大概率都遇到过这个问题手里只有一颗颗牙齿的STL白模或者一整副带牙龈的牙颌模型得把每颗牙单独抠出来才能做后续的排列、设计和测量。直接按曲率切牙颈线过渡平缓曲率信号弱切出来的边界像狗啃用深度学习方法手头没有标注好的几千例数据训练根本推不动。所以在这个场景里投影算法也就是曲面的栅格化算法反而是最务实的一条路把三维网格先“拍扁”到二维栅格上在二维域里做分割再把分割结果回投到三维网格。这个过程稳定、可解释、参数可调而且对STL这种纯三角面片格式天然友好。这篇笔记就拆一份基于这个思路的牙齿分割实现覆盖从网格预处理、投影栅格化到牙龈外轮廓计算的完整链路并标出我在实际牙模数据上踩过的坑。适合刚接触网格分割、想绕过几何直接法和深度学习门槛的工程师或研究生。2. 网格预处理分割结果差八成是输入STL本身不干净2.1 原始牙模网格的病封闭性、法向和噪声三角面片牙齿STL模型来源通常是口内扫描仪或模型扫描仪输出三角网格质量参差不齐。最典型的问题有三个一是网格有孔洞牙龈底面或牙颈部扫描盲区会留洞区域增长算法碰到孔洞边缘就停住二是法向朝向不统一扫描仪重建时部分面片法向翻转导致投影高度图出现厚度翻倍或空洞三是牙面有离散噪声面片扫描时灰尘、唾液反光会造成飘在主体外的孤立碎面。对分割算法来说真正致命的不是网格粗糙而是“水密性”不达标。后续投影栅格化需要从网格往平面投射线如果网格不是封闭的射线在孔洞处直接穿透栅格上那一片就会变成无值或者错误值。所以预处理的第一件事不是平滑而是补洞和统一法向。2.2 预处理步骤补洞、去噪、法向统一的标准流水线实际工程里我一般用VTK或Trimesh库来做这套预处理。下面给出的是用Python Trimesh实现的流程适合在拿到STL后先做一遍清洗import trimesh import numpy as np # 读取原始STL mesh trimesh.load(tooth_arch.stl) # 1. 合并重复顶点去除退化三角面片 mesh.merge_vertices() mesh.update_faces(mesh.area 1e-10) # 2. 补洞Trimesh提供简单的孔洞修补能力 mesh.fill_holes() # 默认只填补边界较小的孔洞 # 3. 去除离群连通分量只保留最大连通块 ccs mesh.split(only_watertightFalse) ccs.sort(keylambda m: m.volume, reverseTrue) mesh ccs[0] # 4. 法向统一修正Trimesh的fix_normals会把法向调整为一致朝外 mesh.fix_normals() # 5. 可选拉普拉斯平滑但不要过度否则会消掉牙尖 smoothed trimesh.smoothing.filter_laplacian(mesh, iterations5, lamb0.5) # 导出清洗后的STL mesh.export(tooth_arch_clean.stl)这段代码的逻辑按顺序处理了四个问题。merge_vertices把扫描仪多次拼接产生的重复顶点合并减少后续邻接查询的歧义fill_holes把牙颈部的可见孔洞补齐但若扫描质量太差出现大面积缺失这一步补不动需要在上游补充扫描或交互式补洞split volume排序删除了扫描时产生的悬浮碎面——飘在外面的小连通块体积通常比主体小几个数量级最后fix_normals统一法向这一步对投影生成高度图来说至关重要。参数上注意两点。fill_holes对超过一定边界长度的孔洞默认不补工业级做法是在Meshlab里用“Close Holes”配合最大边长阈值来补一般阈值设置为网格平均边长的5-10倍。filter_laplacian的迭代次数建议控制在10次以内过重的平滑会把牙尖和牙窝这些分割关键特征抹平后面投影分割时边界反而更难定位。2.3 清洗后的判定标准输出什么才说明可以往下走清洗完别急着跑分割先检查三个指标。第一网格的is_wet水密性属性是否变为True或者补洞后边界边boundary edges数量是否归零第二把法向可视化渲染一遍确认所有法向锥指向外第三检查最大连通分量体积与整网格体积的比值应接近1.0。如果你的输入是只包含牙齿不含底座和牙龈的单颗牙STL跳过连通分量筛选和补洞也可以但法向统一这一步永远别省。投影栅格化假设所有面片法向与投影方向夹角关系是确定的法向一旦颠倒投影出的高度图会把牙窝变成牙尖分割直接翻车。3. 投影栅格化把三维曲面摊平成二维栅格分离牙冠与牙龈3.1 为什么选投影而不选曲面直分割栅格域里的邻域定义是可靠的直接在三网格上做区域增长和聚类邻域关系需要建立半边结构或拓扑处理非流形边时特别容易出错。投影法把三维问题转换为二维图像处理问题以后邻域系统退化为上下左右四邻域或八邻域连通性分析、形态学操作、边界追踪全都可以直接复用成熟的图像处理库。投影的物理意义是这样把牙颌模型沿一个合适的视角方向通常是咬合平面法向做正交投影每个三角面片被“拍”到二维平面上对应像素记录该位置处的深度、法向夹角或密度值。这样形成的栅格图像里牙冠区域因为凸起深度值比牙龈区域低或高取决于投影方向定义牙缝则表现为高低过渡带。这个做法的另外一个好处是栅格化天然是一个降采样过程。一个10万面片的STL投影到512x512栅格后每个像素可能覆盖多个三角面片噪声被平均掉后续分割算法的鲁棒性反而比直接处理顶点要高。3.2 投影方向与栅格分辨率的选取咬合面法向与体素尺寸的经验法则投影方向选不好牙齿在栅格上会互相遮挡分割结果就会出现牙间粘连。常见做法是用PCA对牙颌模型的全部顶点做主成分分析取最小特征值对应的特征向量作为牙颌的“扁平法向”再强制该法向朝上。栅格分辨率直接决定分割精度和计算量的平衡。经验法则是让单个牙齿在栅格上的最短边长不小于80像素否则牙缝过渡带只占两三个像素分割结果在回投三维时会出现严重的锯齿。还有一个更直接的标准栅格像素尺寸取网格平均边长的0.5-1.0倍。平均边长0.3mm的模型栅格分辨率就在0.15-0.3mm/像素对应的栅格尺寸约为包围盒投影宽度除以该像素尺寸。import numpy as np from scipy.spatial import ConvexHull def get_projection_params(mesh): # mesh: trimesh对象已经有vertices和faces verts mesh.vertices # 1. PCA主成分分析最小特征值向量即法向估计 cov np.cov(verts.T) eigvals, eigvecs np.linalg.eigh(cov) normal eigvecs[:, 0] # 最小特征值对应法向 if normal[2] 0: normal -normal # 强制法向朝上 # 2. 建立以法向为z轴的局部坐标系 z_axis normal x_axis np.array([1.0, 0.0, 0.0]) x_axis - (x_axis z_axis) * z_axis x_axis / np.linalg.norm(x_axis) y_axis np.cross(z_axis, x_axis) # 3. 顶点投影到该平面 proj_verts np.column_stack([ verts x_axis, verts y_axis ]) # 4. 计算包围盒和推荐栅格尺寸 min_corner proj_verts.min(axis0) max_corner proj_verts.max(axis0) bbox_size max_corner - min_corner # 5. 按平均边长估算栅格分辨率 edges mesh.edges_unique avg_edge_len np.mean(np.linalg.norm( verts[edges[:, 0]] - verts[edges[:, 1]], axis1)) pixel_size avg_edge_len * 0.8 grid_shape np.ceil(bbox_size / pixel_size).astype(int) return { normal: normal, x_axis: x_axis, y_axis: y_axis, min_corner: min_corner, pixel_size: pixel_size, grid_shape: grid_shape }PCA选法向这个做法的关键假设是牙颌模型在咬合平面方向上的分布最扁平。对于包含牙冠和牙龈的完整牙颌模型这个假设基本成立但如果模型只包含一排牙齿没有牙龈底座主轴方向会偏转这时应该手动指定一个近似咬合平面的法向宁可用粗略的眼观方向也不要盲目信PCA。pixel_size选平均边长0.8倍是基于一个直观考量面片投影到栅格上至少覆盖0.8个像素保证栅格值不是稀疏的。如果你需要分割结果更精细下调到0.5会导致栅格尺寸翻倍计算内存和耗时都涨但牙缝分辨能力提升有限。建议先在0.8倍下跑通全流程确认无误后再单独对感兴趣区域做高分辨率重投影。3.3 栅格值填充与牙冠-牙龈分离深度图、高度差图和种子区域增长栅格化完成后下一步是从栅格数据里把牙齿和牙龈分开。这里我采用的是深度图做种子区域增长法整体思路是牙冠从牙龈表面向上凸起在深度图上表现为局部极值区域从每个极值点出发沿下降方向生长直到遇到深度梯度突变的地方就停止那一条停止边界就是牙冠和牙龈的过渡带。具体实现上把每个栅格像素的值定义为其到投影平面的距离即深度。扫描仪重建出的牙颌模型牙龈部分相对平坦深度值变化平缓牙冠部分鼓起一个包深度变化剧烈。对深度图做形态学闭运算消除牙缝的小凹坑后用局部极大值检测找到候选牙尖点再从牙尖点向周围做梯度受限的区域增长。import cv2 import numpy as np def segment_teeth_by_projection(depth_map, min_grad_thresh0.8): # depth_map: 深度图float32单位mm牙齿区域深度小于牙龈区域 depth_map depth_map.astype(np.float32) # 1. 形态学闭运算填掉牙缝和细小的深度噪声 kernel cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5, 5)) closed cv2.morphologyEx(depth_map, cv2.MORPH_CLOSE, kernel) # 2. 高度差图原始深度与闭运算结果的差 # 闭运算会抹平突起的牙冠差图里牙冠区域显著 diff cv2.absdiff(depth_map, closed) # 3. 二值化高度差超过阈值的像素视为“牙冠种子候选” thresh_val np.percentile(diff, 75) seed_mask (diff thresh_val).astype(np.uint8) * 255 # 4. 形态学开运算去除噪点再找连通域作为种子区域 seed_mask cv2.morphologyEx(seed_mask, cv2.MORPH_OPEN, kernel) num_labels, labels, stats, _ cv2.connectedComponentsWithStats(seed_mask, 8) # 5. 对每个种子区域记录其质心和平均深度作为增长起点 seeds [] for i in range(1, num_labels): if stats[i, cv2.CC_STAT_AREA] 30: continue ys, xs np.where(labels i) seed_depth np.mean(depth_map[ys, xs]) seeds.append({ centroid: (np.mean(xs), np.mean(ys)), seed_depth: seed_depth, label: i }) # 6. 基于深度梯度约束的区域增长 # 从种子区域出发逐像素扩张深度值变化超过阈值的像素并入 # 梯度突变牙颈线位置处停止 gy, gx np.gradient(depth_map) grad_mag np.sqrt(gx**2 gy**2) seg_mask np.zeros_like(depth_map, dtypenp.uint8) for i, seed in enumerate(seeds): seed_y, seed_x np.where(labels seed[label]) stack list(zip(seed_y, seed_x)) visited set(zip(seed_y, seed_x)) while stack: y, x stack.pop() seg_mask[y, x] i 1 for dy, dx in [(-1,0),(1,0),(0,-1),(0,1)]: ny, nx y dy, x dx if (ny, nx) in visited: continue if ny 0 or ny depth_map.shape[0] or nx 0 or nx depth_map.shape[1]: continue if grad_mag[ny, nx] min_grad_thresh: continue if abs(depth_map[ny, nx] - seed[seed_depth]) 4.0: continue visited.add((ny, nx)) stack.append((ny, nx)) return seg_mask, seeds这段代码的核心是种子选取和增长条件的配合。闭运算生成的差图理论上在牙冠凸起处差异最大所以thresh_val取75分位数能筛掉一大片平坦牙龈区梯度阈值min_grad_thresh控制边界敏感度牙颈线处深度变化急剧梯度幅值会明显跳变数值调到0.8能压住内部凹陷的假边界又能在真边界处停住。参数上最值得注意的是4.0mm的深度窗口限制。这是给牙冠高度设的上限正常恒牙牙冠高度约7-12mm但投影后深度方向上牙尖到牙颈的距离一般不超4mm如果种子深度和待生长像素深度差超过4mm基本就是长到了邻牙或牙龈上。对于乳牙模型这个值要下调到2.5mm左右。次区域增长会在牙缝处粘连的主要原因是闭运算核尺寸不够大牙缝凹槽没有被完全填平改进做法是把核尺寸从5x5调到7x7或9x9代价是牙冠边缘会被稍微腐蚀掉一圈。4. 牙龈外轮廓计算从分割掩膜到三维等值线的完整链路4.1 为什么牙冠分割后还要单独算牙龈外轮廓有些场景比如做牙龈袖口形态分析、确定修复体颈缘线、或者打印手术导板时约束就位方向需要的不只是“哪颗牙是哪颗”还要一条明确的牙龈外轮廓也就是牙龈和牙颈的过渡边界线在三维空间中的走向。这条线在二维栅格上表现为分割掩膜的边界但直接把这个边界反投影回网格会得到锯齿严重、穿过三角面片内部的折线没法直接用于CAD建模或3D打印路径规划。更麻烦的一个细节是牙冠和牙龈在STL网格上是连续曲面牙龈外轮廓并不是网格上真实存在的一条边而是“牙冠区域掩膜边界”在网格表面上的插值轨迹。所以计算牙龈外轮廓本质上是三步二维边界提取、边界点三维化、在网格表面做平滑和投影优化。4.2 二维掩膜边界提取与像素坐标回三维承接上一节的分割掩膜seg_mask假设我们需要第k颗牙的牙龈轮廓。先把它单独抽出来做二值化再做一次形态学闭运算把牙冠区域的细小缺口补平用OpenCV的findContours提取外轮廓。然后对轮廓上的每个像素点用栅格化坐标系反变换回三维空间。这一步要注意的是反变换不能只做单点映射。像素中心对应的三维点是该像素覆盖的若干三角面片的“平均位置”直接取深度图的数值会丢掉面片的朝向信息导致后续网格表面投影时出现偏移。我通常在像素坐标转三维后沿投影法向方向做一次射线与网格求交取交点作为该轮廓点的最终位置方法和射线拾取一致。def contour_to_3d(contour_pixels, depth_map, projection_info, mesh): # contour_pixels: (N, 1, 2) int, opencv contour格式 # depth_map: 栅格深度图 # projection_info: get_projection_params的返回结果 # mesh: trimesh对象用于射线求交 pts_3d [] norm projection_info[normal] origin projection_info[min_corner] pixel_size projection_info[pixel_size] x_axis projection_info[x_axis] y_axis projection_info[y_axis] # 仿射变换矩阵栅格像素 - 三维空间 for px, py in contour_pixels.reshape(-1, 2): # 像素坐标转平面坐标 u origin[0] px * pixel_size v origin[1] py * pixel_size depth depth_map[py, px] # 平面坐标 深度 - 初始三维点 pt (u * x_axis v * y_axis depth * norm) # 沿法向做射线求交取网格表面交点 # 方向取法向负方向因为牙冠在牙龈上方 origin_ray pt direction -norm locations, index_ray, index_tri mesh.ray.intersects_location( [origin_ray], [direction], return_localFalse ) if len(locations) 0: # 取离起点最近的交点 dist np.linalg.norm(locations - origin_ray, axis1) pts_3d.append(locations[np.argmin(dist)]) else: # 求交失败则退回深度点 pts_3d.append(pt) return np.array(pts_3d)这段代码里射线求交是关键步骤。mesh.ray.intersects_location返回射线与网格的全部交点取最近的那个是因为模型表面可能有局部凹陷射线穿入穿出会产生多个交点最近交点就是可见表面。若求交失败比如轮廓点在牙缝空隙处射线穿出了网格外保守做法是退回到深度映射点后续平滑步骤会把它拉回来。4.3 轮廓平滑与重采样拉普拉斯回拉、样条拟合与均匀重采样三维轮廓点从二维边界来天然带锯齿。直接拿去CAD里扫掠生成曲面质量一定不过关。实际做法是先用拉普拉斯平滑做轻度去锯齿再用三次样条做参数化拟合最后在样条曲线上均匀重采样。拉普拉斯平滑注意不要迭代太多次否则轮廓会向内侧收缩牙颈线的偏颌特征被抹掉。我一般迭代5次权重0.3。样条拟合这一步因为轮廓是闭合曲线用样条拟合时需要在首尾处做周期延拓处理才能保证闭合处导数连续。from scipy.interpolate import CubicSpline def smooth_contour(pts_3d, num_samples200, smooth_iters5): # 1. 拉普拉斯平滑 for _ in range(smooth_iters): pts_new pts_3d.copy() for i in range(len(pts_3d)): prev pts_3d[(i - 1) % len(pts_3d)] nxt pts_3d[(i 1) % len(pts_3d)] pts_new[i] pts_3d[i] 0.3 * (prev nxt - 2 * pts_3d[i]) pts_3d pts_new # 2. 闭合曲线均匀参数化 diff np.linalg.norm(np.diff(pts_3d, axis0), axis1) t np.concatenate([[0], np.cumsum(diff)]) t_circ np.concatenate([t, [t[-1] 1.0]]) pts_circ np.vstack([pts_3d, pts_3d[0:1]]) # 3. 样条拟合周期边界 cs_x CubicSpline(t_circ, pts_circ[:, 0], bc_typeperiodic) cs_y CubicSpline(t_circ, pts_circ[:, 1], bc_typeperiodic) cs_z CubicSpline(t_circ, pts_circ[:, 2], bc_typeperiodic) # 4. 均匀重采样 t_new np.linspace(0, t[-1], num_samples, endpointFalse) smooth_pts np.column_stack([cs_x(t_new), cs_y(t_new), cs_z(t_new)]) return smooth_ptsbc_typeperiodic保证样条在闭合点处的切向量连续牙颈线不会出现折角。num_samples选200在10-20mm的牙颈周长上相当于每个采样点间距0.05-0.1mm足够支撑后续CAD扫掠和有限元网格生成。曲线均匀重采样之后牙龈外轮廓就被表示为一串有序三维点既可用于可视化也可导出为IGES或DXF曲线做进一步处理。5. 参数避坑投影方向、栅格分辨率和形态学核的四个翻车点5.1 翻车点一投影方向与真实咬合面偏差超过15度牙列全部重叠成一片现象投影图上所有牙齿在深度方向互相遮挡分割出来的每个牙冠区域都是“串糖葫芦”状无法分离单颗牙。原因PCA最小特征值法向在牙颌模型不对称比如单侧缺牙、牙龈萎缩时不稳定或者模型包含底座时主方向被底座偏置。解决不要盲信PCA用目视检查咬合平面的方法覆盖投影方向。做法是把模型渲染后用交互式工具选三个不在同一直线的点两个后牙颊尖、一个前牙切端计算三点平面的法向作为投影方向。若想全自动可以把PCA结果和三个解剖标志点法向之间做一个角度校验偏差超过12度就自动切换为标志点法向。5.2 翻车点二栅格像素尺寸小于平均边长0.3倍内存爆翻两倍现象grid_shape算出来4000x3000生成深度图的内存占用直接让16GB内存的机器开始swap分割进程卡死。原因像素尺寸设得太小栅格面积增大。分辨率提高对分割精度的贡献在平均边长0.5倍以下基本收敛再小只会增加计算量和内存消耗。解决把pixel_size限制在平均边长的0.5-1.0倍区间。如果确实需要更高精度的局部分割不要整体重投影而是对单颗牙的包围盒区域做局部高分辨率栅格化再把结果拼接回整体坐标系。5.3 翻车点三闭运算核尺寸固定9x9石膏模型牙缝照旧粘连现象形态学闭运算后牙缝仍有深度残留区域增长把邻牙连成一块分割掩膜里两牙共享一个连通域。原因牙缝宽度因牙弓位置而异——前牙牙缝窄后牙尤其是第一磨牙远中牙缝宽。固定尺寸的闭运算核很难同时适配所有牙缝。核太小填不满后牙宽缝太大了会把前牙牙冠边缘也吞掉。解决把闭运算核尺寸和投影分辨率挂钩。核的物理直径取牙缝最大宽度估计值的1.5倍换算到像素。临床测量牙缝宽度一般是0.5-1.5mm按0.8mm/像素的分辨率核直径约1-2个像素用7x7到11x11的核可覆盖。处理更稳妥的做法是分两步先用小核闭运算得到一个过分割掩膜再对每个连通域边缘做腐蚀和区域增长来精细定位牙缝边界。5.4 翻车点四区域增长梯度阈值0.8在低分辨率模型上牙龈轮廓缩成一团现象分割掩膜中牙冠区域偏小牙龈外轮廓位置明显低于真实牙颈线看起来像整颗牙被“砍了一截”。原因模型分辨率过低平均边长超过0.5mm时牙颈线处的深度梯度被栅格化平滑掉了梯度阈值0.8导致的停止位置比实际边界偏向牙冠导致轮廓内缩。解决对低分辨率模型把梯度阈值下调到0.4-0.5让区域增长更激进地向牙龈方向推进然后在后处理阶段用牙冠高度约束深度差不超过4mm来刹车。另一个补救措施是把分割结果回投到网格后沿网格法向做一次边界外扩扩张距离约2个平均边长再重新计算轮廓线。从那以后我在每个新模型上跑分割前都会先统计平均边长再决定梯度阈值的刻度。6. 结果验证与半自动修复把Dice系数和孔洞修补变成流程内的一步分割结果不能只看颜色渲染好看。我实际验证过肉眼觉得“分得还可以”的结果和真实牙冠边界对齐后Dice系数最低能掉到0.82——这在外科导板设计里是不可接受的。所以跑完分割后验证环节不能省至少要覆盖两方面数值指标和拓扑合法性。数值指标上我常用的两个验证手段是Dice系数和表面距离。Dice系数需要和标准分割结果对比如果没有真实标注可以退而求其次把分割出来的牙冠区域回投到三维网格上检查每个牙冠连通域的体积是否落在牙齿体积的合理范围内单颗恒牙牙冠体积约0.5-2.0cm³如果有明显偏离的连通域优先怀疑是牙缝粘连或牙龈残留。表面距离验证用MeshLab或VTK计算分割边界到牙颈线手工标注点集的距离95百分位小于0.3mm算通过。拓扑合法性检查看三点每个分割连通域是否连续边界是否闭合孔洞内是否还有未分配的孤立三角面片。这里我习惯在流程里加一段自动修复代码——把分割掩膜回映射到网格面片后对网格做一次孔洞检测发现小于30个三角形的孔洞就直接补掉大于30个的标记出来人工检查。这个阈值不是拍脑袋定的牙颈线附近的孔洞通常小于20个三角形而扫描盲区造成的孔洞动辄上百三角形两者在这个阈值下稳定区分。验证完之后遇到个别牙齿分割边界偏了的情况我通常不重新跑全流程而是在掩膜图上做局部笔刷修复把分割掩膜导成图像格式用标注工具在错分区域涂抹直接把掩膜改对再回投三维。这样做的好处是只影响被修正的那颗牙不会破坏其他牙齿已经稳定的分割边界。这套流程跑通的标志是从原始STL进、到牙龈外轮廓曲线出整个过程不需要人工干预且单颗牙分割的Dice系数稳定在0.90以上。我现在的习惯是每个新扫描设备的数据进来第一周先拿5例历史模型跑一遍全流程把投影方向、梯度阈值、闭运算核这三组参数记录成该设备的预设档之后每例模型直接套档运行遇到分割异常再回到参数调优。这个习惯帮我减少了大半重复调参时间也希望能帮你在自己的数据上少走这些弯路。本文还有配套的精品资源点击获取