简介针对数字化牙科与三维模型处理场景本项目以MATLAB语言实现了一套牙齿STL网格模型分割算法利用投影算法曲面栅格化完成牙齿与牙龈分离及牙龈外轮廓计算覆盖数据读取、栅格化投影、边缘检测、区域分割与轮廓拟合等完整流程。压缩包共54个文件约18.37MB以40个.m算法脚本/函数为核心辅以4个.mat测试数据集、4个Java辅助文件以及README说明文档和项目配置文件可直接用于算法复现与二次开发。目前已有594人浏览学习适合具备MATLAB和图形学基础的研究者、算法工程师及口腔医学工程相关学生参考。代码中包含STL文件读取、内外轮廓过滤、孔洞填充、自动边缘检测、牙龈轮廓凸包拟合等多个功能模块并附带demo.zip演示与多组mat/dat测试数据可帮助深入理解投影分割原理快速构建牙齿数字化分析工具。1. 牙齿STL分割为什么绕道二维投影栅格化解决了邻牙分离的哪个死结拿到一份牙齿STL网格模型最常见的诉求就是把每颗牙单独切出来顺手把牙龈外轮廓也算出来。直接对着三角网格做分割最大的痛点是邻牙之间几乎没有可测量的几何分界曲率接近、凹槽浅区域生长一跑就把两颗牙连成一块。这个标题给的路线恰恰相反——用投影算法把曲面往咬合平面上做栅格化得到一张深度图和一张面片索引图在二维图像上完成分割再把像素标签映射回三维面片。这就是投影算法的本质曲面的栅格化。它把三维空间里难做的邻牙分离转成二维图像上成熟的阈值、形态学和分水岭操作。适合正在做正畸分析、种植导板定位、牙冠设计的软件团队实现成本低比直接上三维深度网络稳得多。注意这里的STL是三角网格文件格式不是C STL容器那个STL。整个落地链路是摆正模型 → 网格预处理 → 栅格化 → 图像分割 → 映射回三维 → 牙龈外轮廓计算。下面按这条线拆开讲参数和坑都放在对应章节。2. 投影前先摆正模型咬合平面估计、姿态旋转与上下颌分离投影方向是整个栅格化的地基。地基歪了后面每一步都在补救而且补不回来。这一章解决的是投影方向怎么定、上下颌怎么拆两个前置问题。2.1 投影方向的坑为什么不能直接拿Z轴投影口扫设备给的STL坐标系是扫描仪自己定的Z轴大概率不指向咬合平面法向。直接拿原始Z轴投影牙冠颊侧和舌侧会大面积自遮挡深度图上出现成片空洞牙缝也被吞掉。投影算法的前提是先把模型旋转到咬合平面法向与Z轴重合投影方向才可靠。判断摆正是否成功的标准很简单旋转后所有顶点的Z坐标方差明显变小且牙尖对应的Z值集中在整个模型的高端。这里有人会问能不能用AABB的包围盒主轴代替咬合平面对单颗牙偶尔可以对全牙弓不行——牙弓是弯曲的包围盒主轴跟咬合平面法向偏差能到十几度。所以正规做法还是先拟合咬合平面。2.2 用PCA求咬合平面最小二乘拟合与旋转参数我一般直接对全部顶点做PCA最小特征值对应的特征向量就是咬合平面法向。这比用RANSAC拟合平面更省事因为牙颌模型顶点量大PCA一次矩阵分解就出结果而且对牙尖这种局部凸起不敏感。import numpy as np import trimesh mesh trimesh.load(jaw.stl, processFalse) verts mesh.vertices.astype(np.float64) centroid np.mean(verts, axis0) centered verts - centroid cov np.cov(centered.T) eigenvalues, eigenvectors np.linalg.eigh(cov) # eigh返回升序第0列是最小特征值对应的特征向量即咬合平面法向 normal eigenvectors[:, 0] if normal[2] 0: normal -normal # 统一朝牙尖方向避免模型被翻过来 target np.array([0.0, 0.0, 1.0]) axis np.cross(normal, target) if np.linalg.norm(axis) 1e-8: axis axis / np.linalg.norm(axis) angle np.arccos(np.clip(np.dot(normal, target), -1.0, 1.0)) rot trimesh.transformations.rotation_matrix(angle, axis, pointcentroid) mesh.apply_transform(rot)逻辑说明eigh得到的特征向量是正交的第0列对应最小特征值也就是顶点分布最集中的方向理论上就是咬合平面的法向。得到法向后用叉积求出旋转轴再用轴角公式构造旋转矩阵让法向转到Z轴正方向。processFalse是为了让顶点顺序保持原始读数做几何运算时不会触发trimesh的自动预处理。参数说明angle不需要手动设置方向Rodrigues公式会自动处理。关键判断是normal[2] 0时翻转法向否则旋转完下颌会倒扣。摆正后验证一下打印mesh.vertices[:, 2].std()牙弓模型一般从原始状态的几十毫米降到几毫米内。2.3 上下颌用一个STL文件时要先拆开很多口扫文件是上下颌一起导出的直接栅格化必然自遮挡。拆分用连通域分析把三角网格的面邻接关系建成图找独立分量。from trimesh.graph import connected_components # face_adjacency是面与面的共享边关系connected_components返回分组后的面索引 groups connected_components(mesh.face_adjacency, min_len64) groups.sort(keylambda g: len(g), reverseTrue) for i, g in enumerate(groups): sub mesh.submesh([g], maintain_vertsTrue) sub.export(fcomponent_{i}.stl)逻辑说明face_adjacency是(N,2)的边数组表示每对相邻面片。connected_components在面邻接图上做并查集把网格拆成互不相连的块。min_len64过滤掉游离面片和小碎块避免把扫描噪声当成独立颌体。参数说明min_len按面片数过滤口扫杂质通常只有几个面片设64比较稳。拆完判断上下颌用每个分量的质心Z值即可——摆正之后Z值大的是上颌小的是下颌。判断结果别省存到配置里后面栅格化和映射都要用它。3. 网格预处理别偷懒去噪、补洞与紧凑面片索引的构建STL网格直接拿来栅格化十有八九要翻车。口扫数据自带三类问题顶点重复存储造成面片缝隙、扫描噪点产生深度图椒盐、边缘破损导致深度图黑洞。这一章把栅格化之前必须做的处理讲透。3.1 拉普拉斯平滑与双边滤波的取舍深度图对Z值很敏感网格上一个0.1mm的噪声尖峰落到图像上就是一颗孤立亮像素后面分水岭会把它当成一颗牙的种子。平滑是必须的但方向很讲究。# 拉普拉斯平滑迭代次数和松弛因子是主要参数 mesh mesh.smoothed(iteration8, lam0.6) # 保留牙尖和窝沟特征的场景用Taubin平滑更稳 # trimesh没有内置Taubin时可以用open3d的filter_smooth_taubin # mesh_open3d.filter_smooth_taubin(number_of_iterations20, lambda_filter0.5, mu-0.53)逻辑说明拉普拉斯平滑是把每个顶点往邻域重心方向拉lam控制每次拉动幅度iteration控制次数。它对噪声敏感但也会让牙尖变钝。Taubin平滑是拉普拉斯的改进版一次正向拉一次反向拉能保住锐利特征齿科模型上用得多。如果只是做分割和轮廓计算拉普拉斯平滑凑合能用如果还要拿分割结果去设计导板建议上Taubin。参数说明lam0.6超过0.8会收缩模型iteration在5到10之间足够。深度图上如果还有残留噪点后面栅格化完再对深度图做一次3x3中值滤波比在网格上猛刷平滑安全。3.2 哪些洞必须补、哪些洞千万别补补洞是个双刃剑。牙冠边缘被扫描漏掉的洞不补的话深度图上就是一圈黑洞牙龈外轮廓会从这里断掉但牙缝区域如果出现洞千万别补——牙缝在深度图上是暗沟是分割邻牙的关键线索补了等于把边界抹掉。import trimesh # 补洞接口max_hole_area按模型尺度调 trimesh.repair.fill_holes(mesh, max_hole_area1.5) # 如果版本不支持max_hole_area先看每个洞的边界长度 for loop in mesh.boundaries: loop_verts mesh.vertices[list(loop)] edge_sum np.linalg.norm(np.diff(loop_verts, axis0), axis1).sum() # 边界总长小于阈值才补逻辑说明fill_holes会把网格上的孔洞用新增面片封起来max_hole_area限制只补小洞。牙缝区域的凹陷不是孔洞一般情况下不会触发补洞逻辑但有些口扫数据牙缝底部会形成细长孔洞一旦被判定为洞就补上了。所以补洞前先可视化所有mesh.boundaries确认哪些需要补。参数说明max_hole_area用相对值单位是模型尺寸平方。我执行的是全牙弓扫描时设1.5到2.0单颗牙模型可以放大到5。补洞后重新检查mesh.boundaries如果牙缝位置还残留边界说明没被误补可以继续。3.3 顶点哈希与半边结构为栅格化准备数据STL文件里顶点是重复存的同一个空间点可能出现在多个面片上。如果不做顶点去重栅格化时同一个位置会反复计算面片索引图也会出现错位。先用哈希把重复顶点合并成紧凑索引。vert_map {} new_verts [] new_faces [] for f in mesh.faces: face [] for vid in f: key tuple(np.round(mesh.vertices[vid], 6)) if key not in vert_map: vert_map[key] len(new_verts) new_verts.append(mesh.vertices[vid]) face.append(vert_map[key]) new_faces.append(face) mesh trimesh.Trimesh(verticesnp.array(new_verts), facesnp.array(new_faces), processFalse)逻辑说明把每个顶点坐标四舍五入到6位小数作为字典键保证同一个空间点只保留一份。new_faces里存的是合并后的新顶点索引。处理完后顶点数组无重复面片索引完全连续——这是后面逐三角形投影填充的基础。参数说明6位小数对应微米级精度STL文件坐标一般有6位小数再低会碰断邻接关系。做完后顺手构建半边邻接表即每个面片三个邻居面的索引后面第5章映射回三维时背面补齐要用。4. 把曲面栅格化成图深度图与面片索引图的生成参数这是整个标题的技术核心曲面的栅格化。所谓栅格化就是把连续的三角网格投影到XY平面按像素网格采样每个像素存两个值——该位置的Z坐标深度图以及覆盖该像素的三角形面片编号面片索引图。这两张图是后面所有操作的原料。4.1 栅格分辨率怎么定一个像素对应多少毫米分辨率直接决定牙缝能不能分开。牙缝宽度一般在0.2到0.5mm分水岭要想稳定工作牙缝在图像上至少占2到3个像素。分辨率设0.2mm/px时0.3mm的牙缝只有1.5像素分割基本靠猜。分辨率(mm/px)60x40mm模型图像尺寸0.3mm牙缝对应像素适用场景0.051200x8006px单颗牙、边缘精修0.1600x4003px全牙弓分割默认值0.2300x2001.5px快速预览不适合分割我默认全牙弓用0.1mm/px。图像分辨率再往上提计算量是平方增长的0.05mm/px时600x400的图变成1200x800遍历面片的耗时翻四倍收益只是边缘更平滑这个钱不值得花。单颗牙模型或要做备牙边缘线时再上0.05。4.2 逐三角形投影填充深度缓冲与面片索引同步写入栅格化的标准做法是遍历每个三角形计算它在XY平面的包围盒只对包围盒里的像素做点是否在三角形内的判断。命中像素写入深度Z和面片编号多个三角形覆盖同一个像素时取Z最小的那个——这等价于从投影方向看过去的遮挡关系。import numpy as np import trimesh def point_in_triangle(pt, tri): v0, v1, v2 tri def sign(a, b, c): return (c[0]-a[0])*(b[1]-a[1]) - (b[0]-a[0])*(c[1]-a[1]) d1 sign(pt, v0, v1) d2 sign(pt, v1, v2) d3 sign(pt, v2, v0) has_neg (d1 0) or (d2 0) or (d3 0) has_pos (d1 0) or (d2 0) or (d3 0) return not (has_neg and has_pos) def interpolate_z(pt, tri): A, B, C tri[0], tri[1], tri[2] v0, v1, v2 B - A, C - A, pt - A d00 np.dot(v0, v0); d01 np.dot(v0, v1) d11 np.dot(v1, v1); d20 np.dot(v2, v0) d21 np.dot(v2, v1) denom d00 * d11 - d01 * d01 v (d11 * d20 - d01 * d21) / denom w (d00 * d21 - d01 * d20) / denom u 1 - v - w return u * A[2] v * B[2] w * C[2] def rasterize_teeth(mesh, resolution0.1): verts mesh.vertices.astype(np.float64) faces mesh.faces x_min, x_max verts[:, 0].min(), verts[:, 0].max() y_min, y_max verts[:, 1].min(), verts[:, 1].max() width int((x_max - x_min) / resolution) 1 height int((y_max - y_min) / resolution) 1 depth np.full((height, width), np.inf, dtypenp.float32) face_id np.full((height, width), -1, dtypenp.int32) for fi, f in enumerate(faces): tri verts[f] px0 int((tri[:, 0].min() - x_min) / resolution) px1 int((tri[:, 0].max() - x_min) / resolution) 1 py0 int((tri[:, 1].min() - y_min) / resolution) py1 int((tri[:, 1].max() - y_min) / resolution) 1 px0, py0 max(px0, 0), max(py0, 0) px1, py1 min(px1, width), min(py1, height) for py in range(py0, py1): for px in range(px0, px1): # 取像素中心点做判断避免边缘像素抖动 p np.array([x_min (px 0.5) * resolution, y_min (py 0.5) * resolution]) if point_in_triangle(p, tri): z interpolate_z(p, tri) if z depth[py, px]: depth[py, px] z face_id[py, px] fi depth[np.isinf(depth)] np.nan return depth, face_id逻辑说明外层循环遍历全部面片对每个面片计算投影包围盒只处理包围盒内的像素避免全图扫描。point_in_triangle用叉积符号判断三条边同侧则点在三角形内。interpolate_z用重心坐标插值像素中心点的Z值。深度和面片索引同步写入保证每个像素既能拿到高度信息又能追溯到三维面片。参数说明resolution0.1是全牙弓默认值单颗牙可以调到0.05。depth初始化为无穷大最后把未覆盖区域改成NaN避免后续图像处理把空值当有效深度。这个Python原型对几十万面片比较慢生产环境把这层逻辑换成OpenGL深度缓冲或CUDA实现数学完全等价速度和正确性都能对齐。4.3 全牙弓的窗口裁剪别让后牙投影重叠全牙弓直接栅格化会遇到一个实际问题后牙向舌侧倾斜第二磨牙和第一磨牙的侧面投影会重叠深度图看不出牙缝。常见做法是沿牙弓分段栅格化。先用所有牙尖顶点拟合一条二次曲线作为牙弓中心线沿中心线以窗口滑动窗口宽度取单颗牙宽度加半个邻牙牙缝大约8到10mm。每个窗口独立栅格化、独立分割再把标签拼回全牙弓坐标。这一步不做后牙区分割结果大概率是错的。5. 从二维标签回到三维像素映射、牙龈外轮廓与避坑清单栅格化完成之后分割和映射就是标准的图像处理流程了。这一章把像素级标签怎么变成面片级标签、牙龈外轮廓怎么算、以及这条链路上最常见的几个坑一次说完。5.1 深度图上的牙齿分割阈值、形态学与分水岭的组合深度图上的牙齿分割单靠阈值必然不行——前牙和后牙的牙尖高度差很大一个固定阈值要么漏了尖牙要么把后牙的龈乳头一起切进来。我惯用的组合是Otsu给初值形态学去噪声距离变换加分水岭做邻牙分离。import cv2 import numpy as np from scipy import ndimage # NaN无效像素填到最底层让它们自然成为背景 z_min np.nanmin(depth) depth_ok np.where(np.isfinite(depth), depth, z_min - 1.0) depth_s ndimage.median_filter(depth_ok, size3) norm (depth_s - np.nanmin(depth_s)) / (np.nanmax(depth_s) - np.nanmin(depth_s)) * 255 depth_u8 np.clip(norm, 0, 255).astype(np.uint8) _, otsu cv2.threshold(depth_u8, 0, 255, cv2.THRESH_BINARY cv2.THRESH_OTSU) tooth_mask (otsu 0).astype(np.uint8) dist cv2.distanceTransform(tooth_mask, cv2.DIST_L2, 3) kernel np.ones((7, 7), np.uint8) local_max cv2.dilate(dist, kernel) dist local_max dist 8 # 种子点距离至少8像素压掉过分割 markers np.zeros(depth_u8.shape, dtypenp.int32) markers[tooth_mask 0] 1 # 背景固定为1 seed_id 1 for y, x in np.argwhere(local_max): seed_id 1 markers[y, x] seed_id labels cv2.watershed(cv2.merge([depth_u8] * 3), markers) labels[labels -1] 0 # 分水岭边界像素先置零后面用邻域标签补齐逻辑说明Otsu把图像分成高的牙冠区和低的牙龈区得到初始tooth_mask。距离变换计算每个前景像素到背景的距离局部极大值就是每颗牙的中心候选。local_max dist 8把距离太小的噪点种子删掉这是控制过分割最关键的一行。分水岭在markers基础上生长背景是1各牙种子从2开始编号边界像素输出为-1。参数说明size3的中值滤波去深度图椒盐不影响大结构。kernel7x7的膨胀比较距离局部极大值时用核越大种子越少全牙弓建议7x7单颗牙可以降到5x5。dist 8的阈值和分辨率挂钩0.1mm/px下8像素对应0.8mm低于这个距离的候选种子基本是噪声。5.2 像素标签映射回面片索引图的关键用法图像分割完每颗牙的像素区域有了唯一ID。要回到三维就要用第4章生成的面片索引图每个像素记录了它属于哪个三角形面片把像素标签按面片分组投票得到面片级标签。from collections import Counter face_labels np.full(len(mesh.faces), -1, dtypenp.int32) # 每个面片收集所有覆盖像素的标签投票决定归属 for fi in range(len(mesh.faces)): ys, xs np.where(face_id fi) if len(ys) 0: continue votes Counter(labels[ys, xs]) if votes: face_labels[fi] votes.most_common(1)[0][0] # 背面不可见面片沿邻接面片投票补齐 for fi in range(len(face_labels)): if face_labels[fi] -1: nbr adjacency[fi] cnt Counter(face_labels[nbr]) cnt.pop(-1, None) if cnt: face_labels[fi] cnt.most_common(1)[0][0]逻辑说明第一轮投票是像素到面片的映射。同一个面片可能被多个像素覆盖这些像素的标签一般一致但在牙齿边界处会跨标签投票能消除边界锯齿。第二轮补齐处理背面面片——投影时它们被挡住没有像素覆盖只能靠邻接面片的标签投票。两轮下来三维网格的每一块都有牙齿ID。这里的坑是第二轮投票会把邻牙标签串进来。如果模型摆正到位牙缝两侧的背面面片正好相邻投票很容易把上颌侧切牙的舌侧补成尖牙。我处理的办法是补齐时只取该面片法向与投影方向夹角大于150度的对面候选相当于只让同轴方向的邻居投票。5.3 牙龈外轮廓计算的完整链路牙龈外轮廓实际上就是牙冠掩膜和牙龈背景的公共边界。分水岭标签里背景是1牙齿是2到N把牙齿区域掩膜提取出来做一次开运算断开牙缝位置的细连接再提取掩膜边界映射回三维点最后样条拟合。teeth (labels 1).astype(np.uint8) teeth cv2.morphologyEx(teeth, cv2.MORPH_OPEN, np.ones((3, 3), np.uint8)) border teeth - cv2.erode(teeth, np.ones((3, 3), np.uint8)) ys, xs np.nonzero(border) margin_pts [] for y, x in zip(ys, xs): if np.isfinite(depth_ok[y, x]): margin_pts.append([ x_min (x 0.5) * resolution, y_min (y 0.5) * resolution, depth_ok[y, x] ]) margin_pts np.array(margin_pts) from scipy.interpolate import splprep, splev tck, u splprep([margin_pts[:, 0], margin_pts[:, 1], margin_pts[:, 2]], s1.0) smooth_pts np.array(splev(np.linspace(0, 1, 200), tck)).T逻辑说明开运算用3x3结构元素目的是断开牙缝方向可能残留的细桥避免外轮廓线钻进牙缝。边界像素通过深度图的坐标反算回三维点比取面片顶点更平滑。样条拟合把离散边界点连成闭合曲线s控制平滑程度。参数说明s1.0在0.1mm/px分辨率下基本保住轮廓细节曲线不会切进牙龈。如果轮廓抖动厉害把s提到3到5代价是切牙乳头位置会被磨平。输出smooth_pts后按需重采样200个点是导板设计和3D打印都够用的密度。5.4 避坑清单四个高频翻车现场现象一深度图边缘全是黑洞牙冠轮廓残缺。原因是模型没摆正侧面自遮挡严重。解决回到第2章检查旋转后Z坐标标准差超过3mm就重新估计咬合平面必要时把牙尖点聚类后再做PCA。现象二分水岭把一颗牙切成两半。原因是距离变换的局部极大值太多牙尖上的小隆起被当成独立种子。解决把dist 8提高到15或把7x7核改成9x9。先压种子再调平滑别反过来。现象三牙缝位置补洞后邻牙分割连成一片。原因是牙缝底部的细长孔洞被fill_holes判定为洞补掉了。解决补洞前用mesh.boundaries把孔洞可视化牙缝区域的边界环单独排除或者把max_hole_area降到0.5只补扫描破损的毛边。现象四映射回三维后舌侧面片标签错误。原因是投影方向看不到舌侧邻接投票把隔壁牙的标签拉进来了。解决先检查face_id里可见面片的覆盖率正常应该超过60%投票补齐时加法向夹角约束只接受来自同一颗牙方向邻居的标签。6. 验证与进阶Dice、边界误差以及深度学习分割怎么接进来分割做完了先说怎么验证。牙齿分割不像自然图像有现成GT我习惯让医生手动标一次或者用商用软件导出的结果当参考再算两个指标Dice系数和平均对称表面距离。Dice看区域重合度边界距离看轮廓贴不贴。def dice(pred, gt): p pred 0 g gt 0 inter np.logical_and(p, g).sum() total p.sum() g.sum() return 2.0 * inter / total if total 0 else 1.0 def assd(pred_pts, gt_pts): # 双向最近邻距离取平均 from scipy.spatial import cKDTree tree1 cKDTree(pred_pts) tree2 cKDTree(gt_pts) d1 tree1.query(gt_pts)[0] d2 tree2.query(pred_pts)[0] return (d1.mean() d2.mean()) / 2.0Dice在0.85以上基本可用ASSD小于0.3mm算优秀。低于这个线先去看是哪颗牙拖后腿通常出在第二磨牙的远中面或者尖牙舌侧。进一步这套投影栅格化管线天然适合接深度学习。深度图、面片索引图、法向图可以堆成三通道输入用U-Net这类语义分割算法直接在图像上做分割替代分水岭。训练标注不用从头画把分水岭结果让医生做一遍修正就是一份质量不错的GT。深度学习版跑通后0.1mm/px的深度图推理只需几十毫秒分水岭那套就可以退役了。我自己现在跑流程习惯把深度图、索引图、标签图、面片标签全部落盘每次参数改动都能对比出是哪一步出了问题。这个习惯比任何调参技巧都省钱。希望帮到你。本文还有配套的精品资源点击获取
牙齿STL分割绕道二维：投影栅格化解决邻牙分离
简介针对数字化牙科与三维模型处理场景本项目以MATLAB语言实现了一套牙齿STL网格模型分割算法利用投影算法曲面栅格化完成牙齿与牙龈分离及牙龈外轮廓计算覆盖数据读取、栅格化投影、边缘检测、区域分割与轮廓拟合等完整流程。压缩包共54个文件约18.37MB以40个.m算法脚本/函数为核心辅以4个.mat测试数据集、4个Java辅助文件以及README说明文档和项目配置文件可直接用于算法复现与二次开发。目前已有594人浏览学习适合具备MATLAB和图形学基础的研究者、算法工程师及口腔医学工程相关学生参考。代码中包含STL文件读取、内外轮廓过滤、孔洞填充、自动边缘检测、牙龈轮廓凸包拟合等多个功能模块并附带demo.zip演示与多组mat/dat测试数据可帮助深入理解投影分割原理快速构建牙齿数字化分析工具。1. 牙齿STL分割为什么绕道二维投影栅格化解决了邻牙分离的哪个死结拿到一份牙齿STL网格模型最常见的诉求就是把每颗牙单独切出来顺手把牙龈外轮廓也算出来。直接对着三角网格做分割最大的痛点是邻牙之间几乎没有可测量的几何分界曲率接近、凹槽浅区域生长一跑就把两颗牙连成一块。这个标题给的路线恰恰相反——用投影算法把曲面往咬合平面上做栅格化得到一张深度图和一张面片索引图在二维图像上完成分割再把像素标签映射回三维面片。这就是投影算法的本质曲面的栅格化。它把三维空间里难做的邻牙分离转成二维图像上成熟的阈值、形态学和分水岭操作。适合正在做正畸分析、种植导板定位、牙冠设计的软件团队实现成本低比直接上三维深度网络稳得多。注意这里的STL是三角网格文件格式不是C STL容器那个STL。整个落地链路是摆正模型 → 网格预处理 → 栅格化 → 图像分割 → 映射回三维 → 牙龈外轮廓计算。下面按这条线拆开讲参数和坑都放在对应章节。2. 投影前先摆正模型咬合平面估计、姿态旋转与上下颌分离投影方向是整个栅格化的地基。地基歪了后面每一步都在补救而且补不回来。这一章解决的是投影方向怎么定、上下颌怎么拆两个前置问题。2.1 投影方向的坑为什么不能直接拿Z轴投影口扫设备给的STL坐标系是扫描仪自己定的Z轴大概率不指向咬合平面法向。直接拿原始Z轴投影牙冠颊侧和舌侧会大面积自遮挡深度图上出现成片空洞牙缝也被吞掉。投影算法的前提是先把模型旋转到咬合平面法向与Z轴重合投影方向才可靠。判断摆正是否成功的标准很简单旋转后所有顶点的Z坐标方差明显变小且牙尖对应的Z值集中在整个模型的高端。这里有人会问能不能用AABB的包围盒主轴代替咬合平面对单颗牙偶尔可以对全牙弓不行——牙弓是弯曲的包围盒主轴跟咬合平面法向偏差能到十几度。所以正规做法还是先拟合咬合平面。2.2 用PCA求咬合平面最小二乘拟合与旋转参数我一般直接对全部顶点做PCA最小特征值对应的特征向量就是咬合平面法向。这比用RANSAC拟合平面更省事因为牙颌模型顶点量大PCA一次矩阵分解就出结果而且对牙尖这种局部凸起不敏感。import numpy as np import trimesh mesh trimesh.load(jaw.stl, processFalse) verts mesh.vertices.astype(np.float64) centroid np.mean(verts, axis0) centered verts - centroid cov np.cov(centered.T) eigenvalues, eigenvectors np.linalg.eigh(cov) # eigh返回升序第0列是最小特征值对应的特征向量即咬合平面法向 normal eigenvectors[:, 0] if normal[2] 0: normal -normal # 统一朝牙尖方向避免模型被翻过来 target np.array([0.0, 0.0, 1.0]) axis np.cross(normal, target) if np.linalg.norm(axis) 1e-8: axis axis / np.linalg.norm(axis) angle np.arccos(np.clip(np.dot(normal, target), -1.0, 1.0)) rot trimesh.transformations.rotation_matrix(angle, axis, pointcentroid) mesh.apply_transform(rot)逻辑说明eigh得到的特征向量是正交的第0列对应最小特征值也就是顶点分布最集中的方向理论上就是咬合平面的法向。得到法向后用叉积求出旋转轴再用轴角公式构造旋转矩阵让法向转到Z轴正方向。processFalse是为了让顶点顺序保持原始读数做几何运算时不会触发trimesh的自动预处理。参数说明angle不需要手动设置方向Rodrigues公式会自动处理。关键判断是normal[2] 0时翻转法向否则旋转完下颌会倒扣。摆正后验证一下打印mesh.vertices[:, 2].std()牙弓模型一般从原始状态的几十毫米降到几毫米内。2.3 上下颌用一个STL文件时要先拆开很多口扫文件是上下颌一起导出的直接栅格化必然自遮挡。拆分用连通域分析把三角网格的面邻接关系建成图找独立分量。from trimesh.graph import connected_components # face_adjacency是面与面的共享边关系connected_components返回分组后的面索引 groups connected_components(mesh.face_adjacency, min_len64) groups.sort(keylambda g: len(g), reverseTrue) for i, g in enumerate(groups): sub mesh.submesh([g], maintain_vertsTrue) sub.export(fcomponent_{i}.stl)逻辑说明face_adjacency是(N,2)的边数组表示每对相邻面片。connected_components在面邻接图上做并查集把网格拆成互不相连的块。min_len64过滤掉游离面片和小碎块避免把扫描噪声当成独立颌体。参数说明min_len按面片数过滤口扫杂质通常只有几个面片设64比较稳。拆完判断上下颌用每个分量的质心Z值即可——摆正之后Z值大的是上颌小的是下颌。判断结果别省存到配置里后面栅格化和映射都要用它。3. 网格预处理别偷懒去噪、补洞与紧凑面片索引的构建STL网格直接拿来栅格化十有八九要翻车。口扫数据自带三类问题顶点重复存储造成面片缝隙、扫描噪点产生深度图椒盐、边缘破损导致深度图黑洞。这一章把栅格化之前必须做的处理讲透。3.1 拉普拉斯平滑与双边滤波的取舍深度图对Z值很敏感网格上一个0.1mm的噪声尖峰落到图像上就是一颗孤立亮像素后面分水岭会把它当成一颗牙的种子。平滑是必须的但方向很讲究。# 拉普拉斯平滑迭代次数和松弛因子是主要参数 mesh mesh.smoothed(iteration8, lam0.6) # 保留牙尖和窝沟特征的场景用Taubin平滑更稳 # trimesh没有内置Taubin时可以用open3d的filter_smooth_taubin # mesh_open3d.filter_smooth_taubin(number_of_iterations20, lambda_filter0.5, mu-0.53)逻辑说明拉普拉斯平滑是把每个顶点往邻域重心方向拉lam控制每次拉动幅度iteration控制次数。它对噪声敏感但也会让牙尖变钝。Taubin平滑是拉普拉斯的改进版一次正向拉一次反向拉能保住锐利特征齿科模型上用得多。如果只是做分割和轮廓计算拉普拉斯平滑凑合能用如果还要拿分割结果去设计导板建议上Taubin。参数说明lam0.6超过0.8会收缩模型iteration在5到10之间足够。深度图上如果还有残留噪点后面栅格化完再对深度图做一次3x3中值滤波比在网格上猛刷平滑安全。3.2 哪些洞必须补、哪些洞千万别补补洞是个双刃剑。牙冠边缘被扫描漏掉的洞不补的话深度图上就是一圈黑洞牙龈外轮廓会从这里断掉但牙缝区域如果出现洞千万别补——牙缝在深度图上是暗沟是分割邻牙的关键线索补了等于把边界抹掉。import trimesh # 补洞接口max_hole_area按模型尺度调 trimesh.repair.fill_holes(mesh, max_hole_area1.5) # 如果版本不支持max_hole_area先看每个洞的边界长度 for loop in mesh.boundaries: loop_verts mesh.vertices[list(loop)] edge_sum np.linalg.norm(np.diff(loop_verts, axis0), axis1).sum() # 边界总长小于阈值才补逻辑说明fill_holes会把网格上的孔洞用新增面片封起来max_hole_area限制只补小洞。牙缝区域的凹陷不是孔洞一般情况下不会触发补洞逻辑但有些口扫数据牙缝底部会形成细长孔洞一旦被判定为洞就补上了。所以补洞前先可视化所有mesh.boundaries确认哪些需要补。参数说明max_hole_area用相对值单位是模型尺寸平方。我执行的是全牙弓扫描时设1.5到2.0单颗牙模型可以放大到5。补洞后重新检查mesh.boundaries如果牙缝位置还残留边界说明没被误补可以继续。3.3 顶点哈希与半边结构为栅格化准备数据STL文件里顶点是重复存的同一个空间点可能出现在多个面片上。如果不做顶点去重栅格化时同一个位置会反复计算面片索引图也会出现错位。先用哈希把重复顶点合并成紧凑索引。vert_map {} new_verts [] new_faces [] for f in mesh.faces: face [] for vid in f: key tuple(np.round(mesh.vertices[vid], 6)) if key not in vert_map: vert_map[key] len(new_verts) new_verts.append(mesh.vertices[vid]) face.append(vert_map[key]) new_faces.append(face) mesh trimesh.Trimesh(verticesnp.array(new_verts), facesnp.array(new_faces), processFalse)逻辑说明把每个顶点坐标四舍五入到6位小数作为字典键保证同一个空间点只保留一份。new_faces里存的是合并后的新顶点索引。处理完后顶点数组无重复面片索引完全连续——这是后面逐三角形投影填充的基础。参数说明6位小数对应微米级精度STL文件坐标一般有6位小数再低会碰断邻接关系。做完后顺手构建半边邻接表即每个面片三个邻居面的索引后面第5章映射回三维时背面补齐要用。4. 把曲面栅格化成图深度图与面片索引图的生成参数这是整个标题的技术核心曲面的栅格化。所谓栅格化就是把连续的三角网格投影到XY平面按像素网格采样每个像素存两个值——该位置的Z坐标深度图以及覆盖该像素的三角形面片编号面片索引图。这两张图是后面所有操作的原料。4.1 栅格分辨率怎么定一个像素对应多少毫米分辨率直接决定牙缝能不能分开。牙缝宽度一般在0.2到0.5mm分水岭要想稳定工作牙缝在图像上至少占2到3个像素。分辨率设0.2mm/px时0.3mm的牙缝只有1.5像素分割基本靠猜。分辨率(mm/px)60x40mm模型图像尺寸0.3mm牙缝对应像素适用场景0.051200x8006px单颗牙、边缘精修0.1600x4003px全牙弓分割默认值0.2300x2001.5px快速预览不适合分割我默认全牙弓用0.1mm/px。图像分辨率再往上提计算量是平方增长的0.05mm/px时600x400的图变成1200x800遍历面片的耗时翻四倍收益只是边缘更平滑这个钱不值得花。单颗牙模型或要做备牙边缘线时再上0.05。4.2 逐三角形投影填充深度缓冲与面片索引同步写入栅格化的标准做法是遍历每个三角形计算它在XY平面的包围盒只对包围盒里的像素做点是否在三角形内的判断。命中像素写入深度Z和面片编号多个三角形覆盖同一个像素时取Z最小的那个——这等价于从投影方向看过去的遮挡关系。import numpy as np import trimesh def point_in_triangle(pt, tri): v0, v1, v2 tri def sign(a, b, c): return (c[0]-a[0])*(b[1]-a[1]) - (b[0]-a[0])*(c[1]-a[1]) d1 sign(pt, v0, v1) d2 sign(pt, v1, v2) d3 sign(pt, v2, v0) has_neg (d1 0) or (d2 0) or (d3 0) has_pos (d1 0) or (d2 0) or (d3 0) return not (has_neg and has_pos) def interpolate_z(pt, tri): A, B, C tri[0], tri[1], tri[2] v0, v1, v2 B - A, C - A, pt - A d00 np.dot(v0, v0); d01 np.dot(v0, v1) d11 np.dot(v1, v1); d20 np.dot(v2, v0) d21 np.dot(v2, v1) denom d00 * d11 - d01 * d01 v (d11 * d20 - d01 * d21) / denom w (d00 * d21 - d01 * d20) / denom u 1 - v - w return u * A[2] v * B[2] w * C[2] def rasterize_teeth(mesh, resolution0.1): verts mesh.vertices.astype(np.float64) faces mesh.faces x_min, x_max verts[:, 0].min(), verts[:, 0].max() y_min, y_max verts[:, 1].min(), verts[:, 1].max() width int((x_max - x_min) / resolution) 1 height int((y_max - y_min) / resolution) 1 depth np.full((height, width), np.inf, dtypenp.float32) face_id np.full((height, width), -1, dtypenp.int32) for fi, f in enumerate(faces): tri verts[f] px0 int((tri[:, 0].min() - x_min) / resolution) px1 int((tri[:, 0].max() - x_min) / resolution) 1 py0 int((tri[:, 1].min() - y_min) / resolution) py1 int((tri[:, 1].max() - y_min) / resolution) 1 px0, py0 max(px0, 0), max(py0, 0) px1, py1 min(px1, width), min(py1, height) for py in range(py0, py1): for px in range(px0, px1): # 取像素中心点做判断避免边缘像素抖动 p np.array([x_min (px 0.5) * resolution, y_min (py 0.5) * resolution]) if point_in_triangle(p, tri): z interpolate_z(p, tri) if z depth[py, px]: depth[py, px] z face_id[py, px] fi depth[np.isinf(depth)] np.nan return depth, face_id逻辑说明外层循环遍历全部面片对每个面片计算投影包围盒只处理包围盒内的像素避免全图扫描。point_in_triangle用叉积符号判断三条边同侧则点在三角形内。interpolate_z用重心坐标插值像素中心点的Z值。深度和面片索引同步写入保证每个像素既能拿到高度信息又能追溯到三维面片。参数说明resolution0.1是全牙弓默认值单颗牙可以调到0.05。depth初始化为无穷大最后把未覆盖区域改成NaN避免后续图像处理把空值当有效深度。这个Python原型对几十万面片比较慢生产环境把这层逻辑换成OpenGL深度缓冲或CUDA实现数学完全等价速度和正确性都能对齐。4.3 全牙弓的窗口裁剪别让后牙投影重叠全牙弓直接栅格化会遇到一个实际问题后牙向舌侧倾斜第二磨牙和第一磨牙的侧面投影会重叠深度图看不出牙缝。常见做法是沿牙弓分段栅格化。先用所有牙尖顶点拟合一条二次曲线作为牙弓中心线沿中心线以窗口滑动窗口宽度取单颗牙宽度加半个邻牙牙缝大约8到10mm。每个窗口独立栅格化、独立分割再把标签拼回全牙弓坐标。这一步不做后牙区分割结果大概率是错的。5. 从二维标签回到三维像素映射、牙龈外轮廓与避坑清单栅格化完成之后分割和映射就是标准的图像处理流程了。这一章把像素级标签怎么变成面片级标签、牙龈外轮廓怎么算、以及这条链路上最常见的几个坑一次说完。5.1 深度图上的牙齿分割阈值、形态学与分水岭的组合深度图上的牙齿分割单靠阈值必然不行——前牙和后牙的牙尖高度差很大一个固定阈值要么漏了尖牙要么把后牙的龈乳头一起切进来。我惯用的组合是Otsu给初值形态学去噪声距离变换加分水岭做邻牙分离。import cv2 import numpy as np from scipy import ndimage # NaN无效像素填到最底层让它们自然成为背景 z_min np.nanmin(depth) depth_ok np.where(np.isfinite(depth), depth, z_min - 1.0) depth_s ndimage.median_filter(depth_ok, size3) norm (depth_s - np.nanmin(depth_s)) / (np.nanmax(depth_s) - np.nanmin(depth_s)) * 255 depth_u8 np.clip(norm, 0, 255).astype(np.uint8) _, otsu cv2.threshold(depth_u8, 0, 255, cv2.THRESH_BINARY cv2.THRESH_OTSU) tooth_mask (otsu 0).astype(np.uint8) dist cv2.distanceTransform(tooth_mask, cv2.DIST_L2, 3) kernel np.ones((7, 7), np.uint8) local_max cv2.dilate(dist, kernel) dist local_max dist 8 # 种子点距离至少8像素压掉过分割 markers np.zeros(depth_u8.shape, dtypenp.int32) markers[tooth_mask 0] 1 # 背景固定为1 seed_id 1 for y, x in np.argwhere(local_max): seed_id 1 markers[y, x] seed_id labels cv2.watershed(cv2.merge([depth_u8] * 3), markers) labels[labels -1] 0 # 分水岭边界像素先置零后面用邻域标签补齐逻辑说明Otsu把图像分成高的牙冠区和低的牙龈区得到初始tooth_mask。距离变换计算每个前景像素到背景的距离局部极大值就是每颗牙的中心候选。local_max dist 8把距离太小的噪点种子删掉这是控制过分割最关键的一行。分水岭在markers基础上生长背景是1各牙种子从2开始编号边界像素输出为-1。参数说明size3的中值滤波去深度图椒盐不影响大结构。kernel7x7的膨胀比较距离局部极大值时用核越大种子越少全牙弓建议7x7单颗牙可以降到5x5。dist 8的阈值和分辨率挂钩0.1mm/px下8像素对应0.8mm低于这个距离的候选种子基本是噪声。5.2 像素标签映射回面片索引图的关键用法图像分割完每颗牙的像素区域有了唯一ID。要回到三维就要用第4章生成的面片索引图每个像素记录了它属于哪个三角形面片把像素标签按面片分组投票得到面片级标签。from collections import Counter face_labels np.full(len(mesh.faces), -1, dtypenp.int32) # 每个面片收集所有覆盖像素的标签投票决定归属 for fi in range(len(mesh.faces)): ys, xs np.where(face_id fi) if len(ys) 0: continue votes Counter(labels[ys, xs]) if votes: face_labels[fi] votes.most_common(1)[0][0] # 背面不可见面片沿邻接面片投票补齐 for fi in range(len(face_labels)): if face_labels[fi] -1: nbr adjacency[fi] cnt Counter(face_labels[nbr]) cnt.pop(-1, None) if cnt: face_labels[fi] cnt.most_common(1)[0][0]逻辑说明第一轮投票是像素到面片的映射。同一个面片可能被多个像素覆盖这些像素的标签一般一致但在牙齿边界处会跨标签投票能消除边界锯齿。第二轮补齐处理背面面片——投影时它们被挡住没有像素覆盖只能靠邻接面片的标签投票。两轮下来三维网格的每一块都有牙齿ID。这里的坑是第二轮投票会把邻牙标签串进来。如果模型摆正到位牙缝两侧的背面面片正好相邻投票很容易把上颌侧切牙的舌侧补成尖牙。我处理的办法是补齐时只取该面片法向与投影方向夹角大于150度的对面候选相当于只让同轴方向的邻居投票。5.3 牙龈外轮廓计算的完整链路牙龈外轮廓实际上就是牙冠掩膜和牙龈背景的公共边界。分水岭标签里背景是1牙齿是2到N把牙齿区域掩膜提取出来做一次开运算断开牙缝位置的细连接再提取掩膜边界映射回三维点最后样条拟合。teeth (labels 1).astype(np.uint8) teeth cv2.morphologyEx(teeth, cv2.MORPH_OPEN, np.ones((3, 3), np.uint8)) border teeth - cv2.erode(teeth, np.ones((3, 3), np.uint8)) ys, xs np.nonzero(border) margin_pts [] for y, x in zip(ys, xs): if np.isfinite(depth_ok[y, x]): margin_pts.append([ x_min (x 0.5) * resolution, y_min (y 0.5) * resolution, depth_ok[y, x] ]) margin_pts np.array(margin_pts) from scipy.interpolate import splprep, splev tck, u splprep([margin_pts[:, 0], margin_pts[:, 1], margin_pts[:, 2]], s1.0) smooth_pts np.array(splev(np.linspace(0, 1, 200), tck)).T逻辑说明开运算用3x3结构元素目的是断开牙缝方向可能残留的细桥避免外轮廓线钻进牙缝。边界像素通过深度图的坐标反算回三维点比取面片顶点更平滑。样条拟合把离散边界点连成闭合曲线s控制平滑程度。参数说明s1.0在0.1mm/px分辨率下基本保住轮廓细节曲线不会切进牙龈。如果轮廓抖动厉害把s提到3到5代价是切牙乳头位置会被磨平。输出smooth_pts后按需重采样200个点是导板设计和3D打印都够用的密度。5.4 避坑清单四个高频翻车现场现象一深度图边缘全是黑洞牙冠轮廓残缺。原因是模型没摆正侧面自遮挡严重。解决回到第2章检查旋转后Z坐标标准差超过3mm就重新估计咬合平面必要时把牙尖点聚类后再做PCA。现象二分水岭把一颗牙切成两半。原因是距离变换的局部极大值太多牙尖上的小隆起被当成独立种子。解决把dist 8提高到15或把7x7核改成9x9。先压种子再调平滑别反过来。现象三牙缝位置补洞后邻牙分割连成一片。原因是牙缝底部的细长孔洞被fill_holes判定为洞补掉了。解决补洞前用mesh.boundaries把孔洞可视化牙缝区域的边界环单独排除或者把max_hole_area降到0.5只补扫描破损的毛边。现象四映射回三维后舌侧面片标签错误。原因是投影方向看不到舌侧邻接投票把隔壁牙的标签拉进来了。解决先检查face_id里可见面片的覆盖率正常应该超过60%投票补齐时加法向夹角约束只接受来自同一颗牙方向邻居的标签。6. 验证与进阶Dice、边界误差以及深度学习分割怎么接进来分割做完了先说怎么验证。牙齿分割不像自然图像有现成GT我习惯让医生手动标一次或者用商用软件导出的结果当参考再算两个指标Dice系数和平均对称表面距离。Dice看区域重合度边界距离看轮廓贴不贴。def dice(pred, gt): p pred 0 g gt 0 inter np.logical_and(p, g).sum() total p.sum() g.sum() return 2.0 * inter / total if total 0 else 1.0 def assd(pred_pts, gt_pts): # 双向最近邻距离取平均 from scipy.spatial import cKDTree tree1 cKDTree(pred_pts) tree2 cKDTree(gt_pts) d1 tree1.query(gt_pts)[0] d2 tree2.query(pred_pts)[0] return (d1.mean() d2.mean()) / 2.0Dice在0.85以上基本可用ASSD小于0.3mm算优秀。低于这个线先去看是哪颗牙拖后腿通常出在第二磨牙的远中面或者尖牙舌侧。进一步这套投影栅格化管线天然适合接深度学习。深度图、面片索引图、法向图可以堆成三通道输入用U-Net这类语义分割算法直接在图像上做分割替代分水岭。训练标注不用从头画把分水岭结果让医生做一遍修正就是一份质量不错的GT。深度学习版跑通后0.1mm/px的深度图推理只需几十毫秒分水岭那套就可以退役了。我自己现在跑流程习惯把深度图、索引图、标签图、面片标签全部落盘每次参数改动都能对比出是哪一步出了问题。这个习惯比任何调参技巧都省钱。希望帮到你。本文还有配套的精品资源点击获取
