简介面向牙科数字化模型处理这套资源提供了基于投影与曲面栅格化的STL网格分割实现核心是牙龈区域外轮廓的自动提取。方法先将三维牙颌网格曲面投影到二维平面经栅格化离散、局部曲率计算与区域生长聚类识别牙齿与牙龈交界最终把轮廓映射回三维空间完成分割。资源共65个文件、约27.99MB以40个MATLAB脚本为主覆盖曲率计算、轮廓拟合、孔洞修复、区域生长等模块另有4个Java源码用于STL读取与模型构建4个mat数据文件、README说明及调试备份的zbak文件便于对照实现思路与参数调整。已有46人学习下载。整套代码流程完整、参数可控适合口腔数字化与三维网格处理方向学习者复现研究借助投影精度与栅格分辨率设置即可适配不同扫描质量的牙颌数据为修复体设计或病理分析提供可靠基础。1. 牙齿STL网格模型分割为什么投影和栅格化能解决分牙问题在正畸、种植和隐形牙套方案设计里口腔扫描件导出的通常是完整牙弓的STL网格模型牙齿和牙龈粘在同一个三角网格上没有语义标签。要做单颗牙的移动模拟、牙冠备牙量分析或龈缘线设计第一步就是把混合网格拆成“每颗牙一个独立网格”再把牙龈区域的外轮廓算出来。这个需求听着简单直接上手却发现通用网格分割算法根本不好用牙齿之间的邻接面平缓、凹痕浅基于曲率的聚类会把两颗邻牙并成一个区域而基于区域生长的做法又容易从牙颈部漏到牙龈上去。我拆过好几套这类资源最后稳定跑通的方案就是标题里这条路线——先用投影把三维网格摊到二维平面上做规则栅格化在二维栅格里做连通域分割再映射回三维网格最后在牙龈带上做外轮廓计算。这个思路适合谁主要是做口腔数字化项目、医用网格后处理、以及想用Python/ C 处理STL但不想碰重型网格库的朋友原理不复杂代码量也不大能落地。2. 网格预处理与坐标系矫正分割前必须做对的三件事2.1 顶点焊接与法向统一扫描得到的STL有两种常见问题重复顶点和法向不一致。重复顶点会让邻接关系计算出现孤立面片——看上去两个三角形挨着实际它们的顶点不是同一个内存对象区域生长根本走不过去。法向不一致则会让基于法向差的边界检测失效因为相邻面片一个朝外一个朝内法向夹角会算出一个假边界。我一般会先用trimesh把网格载入然后做顶点焊接和法向统一代码大概是这样的import trimesh import numpy as np # 载入STLmerge_verticesTrue会把距离小于tol的顶点合并 mesh trimesh.load(tooth_arch.stl, forcemesh, merge_verticesTrue) # 法向统一到朝向网格外部 mesh.fix_normals() # 强制计算顶点邻接关系后续区域生长会频繁用到 mesh.vertex_adjacency_graph # 触发邻接图构建这里的merge_vertices默认容差是tol1e-8如果扫描件本身顶点有微小的坐标抖动我会放宽到1e-6不然焊不干净。fix_normals()做的是让相邻面片法向夹角尽量小再配合有向包围盒把整体翻转修正到朝外。法向统一之后建议顺手把退化三角形面积接近零的面片剔除否则后续投影时会出现除零和异常栅格点。2.2 质心对齐与牙弓方向估计STL模型从扫描设备出来后坐标系是随机的可能旋转了任意角度。直接对原始坐标做投影基本上是灾难——牙齿会斜着摊开邻牙在投影平面上重重叠叠栅格化之后连成一整片。所以预处理里最关键的一步是把牙弓摆正先计算质心并把模型平移到原点再用PCA估计牙弓的主方向把牙弓长轴旋转到X轴方向。# 计算顶点质心并平移 vertices mesh.vertices centroid vertices.mean(axis0) vertices_centered vertices - centroid # PCA估计主方向 cov np.cov(vertices_centered.T) eig_val, eig_vec np.linalg.eigh(cov) # 最大特征值对应的特征向量是牙弓长轴方向旋转到X轴 long_axis eig_vec[:, -1] rotation trimesh.geometry.align_vectors(long_axis, [1, 0, 0]) mesh.apply_transform(rotation)这里用eigh而不是eig是因为协方差矩阵是实对称的eigh数值上更稳。旋转后牙弓长轴对应X方向Z轴指向咬合方向Y方向则是颊舌向。这个坐标系的建立直接决定了后续投影平面怎么选——我会把投影平面定义为XY平面也就是把牙弓从咬合方向往下看这样每颗牙的冠面会摊开成一个不重叠的椭圆区域。如果原始模型是单颗牙而不是整副牙弓PCA的主方向可能是牙长轴方向需要根据实际需求改成把牙长轴对齐到Z轴我一般做个开关参数控制。2.3 预处理失败的观察窗口预处理做没做对最直接的观察办法是把网格渲染出来看坐标系轴或者输出质心、长轴方向向量打印出来检查。常见错误是PCA方向选反了——牙弓长轴指向负X方向旋转矩阵跟着反投影结果整体镜像。这个问题我遇到过两次现在的做法是固定一个规则令旋转后模型的X轴方向与牙弓末端的磨牙侧一致简单点说取牙弓两侧最远点的连线方向做X轴再通过叉乘保证Y轴朝向舌侧。坐标方向统一之后后面栅格化的每一个坐标判断才不会出幺蛾子。3. 投影与曲面栅格化分割从三维网格到二维掩膜的映射计算3.1 投影平面的选择与牙弓展开投影不是简单地把顶点压到XY平面上就完事。牙弓是弯曲的前端切牙区弧线陡、后段磨牙区弧线平如果直接竖直投影切牙区的邻牙在平面上会挤在一起栅格化后连通域边界模糊。因此要先把牙弓“拉直”相当于把三维曲面展开到二维参数域。常见做法是沿牙弓中心线做弧长参数化。先拾取牙弓上的一个参考点序列通常是手工标几个点或从顶点密度大的区域自动拟合一条三次样条曲线然后对每个网格顶点计算它在中心线上的投影点和弧长位置。这样展开后X坐标是牙弓弧长方向Y坐标是到中心线的垂直距离Z坐标仍然保留原始咬合高度信息。展开之后的网格再往XY平面投影邻牙之间就有清晰的沟壑。# 假设centerline是弧长参数化的牙弓中心线numpy数组形状(N,3) # vertices_centered是预处理后的顶点 s_coords np.zeros(len(vertices_centered)) perp_dists np.zeros(len(vertices_centered)) for i, v in enumerate(vertices_centered): # 找中心线上最近的两个点做插值得到弧长和垂距 diff centerline - v dist np.linalg.norm(diff, axis1) idx np.argmin(dist) # 近似直接用最近点的弧长作为该顶点的弧长坐标 s_coords[i] arc_length[idx] # 垂直距离顶点到中心线最近点的切线方向的叉积投影 perp_dists[i] np.dot(v - centerline[idx], perp_dir)这段代码是近似实现精度够用。arc_length是预先对中心线采样点累计求和得到的弧长数组perp_dir是中心线在最近点处的法向。对牙弓这种形态这种逐点最近邻近似就能得到稳定的展开效果。展开后观察点云切牙、尖牙、磨牙应该各自分离成孤立的簇如果还有重叠说明中心线拟合偏了要回头检查参考点。3.2 栅格化的分辨率设定与连通域标记展开后的点云没有拓扑关系栅格化就是把点云变成规则的二维网格。我一般把牙弓展开区域划分成网格每个格子统计是否有顶点投影进去得到一个二值掩膜。关键参数是栅格分辨率——太粗会把邻牙之间的沟壑直接抹平太细则掩膜上出现大量空洞连通域分析会把同一颗牙切成好几块。# 栅格化把展开后的XY坐标映射到栅格 resolution 0.2 # 单位mm按平均边长2-3倍取 x_min, y_min s_coords.min(), perp_dists.min() x_max, y_max s_coords.max(), perp_dists.max() width int((x_max - x_min) / resolution) 1 height int((y_max - y_min) / resolution) 1 mask np.zeros((height, width), dtypenp.uint8) for s, p in zip(s_coords, perp_dists): col int((s - x_min) / resolution) row int((p - y_min) / resolution) if 0 row height and 0 col width: mask[row, col] 1 # 形态学闭运算填补小孔 from scipy import ndimage mask ndimage.binary_closing(mask, iterations2).astype(np.uint8) # 连通域标记 labeled, num_features ndimage.label(mask)分辨率我通常取网格平均边长的2~3倍。牙齿STL的扫描精度普遍在0.05~0.1mm平均边长在0.2mm左右所以resolution0.5会太粗0.2左右比较稳妥。binary_closing迭代两次能把扫描噪声造成的空洞补上如果迭代太多会把邻牙之间的缝隙也填掉这个参数要针对具体数据微调。ndimage.label默认用4连通我习惯改成8连通因为牙齿投影边缘是斜的4连通容易产生锯齿断裂。3.3 从二维掩膜映射回三维网格连通域标记完成后每个二维栅格单元有了一个标签号。接下来要把标签映射回三维网格顶点然后通过顶点标签聚合出每颗牙的子网格。做法是建立一个从栅格坐标到顶点索引的映射关系遍历每个顶点找到它落在哪个栅格单元里读取对应标签。# 建立栅格坐标到顶点索引的映射 cell_to_verts {} for i, (s, p) in enumerate(zip(s_coords, perp_dists)): col int((s - x_min) / resolution) row int((p - y_min) / resolution) key (row, col) label labeled[row, col] if label 0: continue cell_to_verts.setdefault(key, []).append((i, label)) # 给每个顶点赋标签 vert_labels np.zeros(len(vertices_centered), dtypeint) for (row, col), items in cell_to_verts.items(): for vert_idx, label in items: vert_labels[vert_idx] label # 按标签抽取子网格 from collections import defaultdict label_to_faces defaultdict(list) for face_idx, face in enumerate(mesh.faces): labels vert_labels[face] if labels[0] labels[1] labels[2]: label_to_faces[labels[0]].append(face_idx)这里有个关键细节只有三个顶点标签一致的面片才归入对应牙齿因为一个三角形如果跨越两个标签它大概率位于牙齿边界上直接丢弃能避免相邻牙齿面片互相污染。丢弃边界三角面片会在牙齿边缘留下一圈小缺口但这不影响后续外轮廓计算轮廓提取本来就基于边界边而不是完整表面。如果缝隙明显可以再用形态学膨胀把边界顶点吸附到相邻牙上——不过这一步我通常不做边界缺口对轮廓测量的影响可以忽略。4. 牙龈外轮廓计算边界环提取与曲线闭合4.1 牙龈区域的曲率刻画与阈值选取分割完成后模型里除了牙齿剩下的就是牙龈区域。牙龈外轮廓要算的是牙龈与牙槽黏膜交界的地方临床上叫膜龈联合线但在网格上并没有明确的解剖标记只能用几何特征近似。牙龈区域靠近牙颈部的曲率变化明显向外过度到平坦的牙槽区域这个过渡带就是轮廓所在。计算网格曲率我用的是二环邻域法——对每个顶点取它的一环邻域和两环邻域拟合二次曲面然后从曲面系数导出主曲率。这个做法比直接调库稳因为我需要知道每个顶点的曲率张量方向而不只是最大值。from scipy.linalg import lstsq def fit_quadric(verts): # 拟合 z a*x^2 b*y^2 c*xy d*x e*y f A np.c_[verts[:, 0]**2, verts[:, 1]**2, verts[:, 0]*verts[:, 1], verts[:, 0], verts[:, 1], np.ones(len(verts))] coef, _, _, _ lstsq(A, verts[:, 2]) return coef def vertex_curvature(mesh, v_idx): neighbors mesh.vertex_neighbors[v_idx] # 一环不够扩展到二环 for n in list(neighbors): neighbors list(set(neighbors) | set(mesh.vertex_neighbors[n])) pts mesh.vertices[neighbors] local pts - mesh.vertices[v_idx] # 转成局部坐标系法向为z轴 normal mesh.vertex_normals[v_idx] z_axis normal x_axis np.cross(z_axis, [1, 0, 0]) if np.linalg.norm(x_axis) 1e-6: x_axis np.cross(z_axis, [0, 1, 0]) x_axis x_axis / np.linalg.norm(x_axis) y_axis np.cross(z_axis, x_axis) local_coords np.c_[local x_axis, local y_axis, local z_axis] coef fit_quadric(local_coords[:, :2] local_coords[:, 2:]) # 补齐 # 平均曲率近似 H coef[0] coef[1] return H这段代码里mesh.vertex_neighbors是顶点邻接表需要预处理阶段构建。二环邻域能覆盖更广的曲面趋势比单纯一环对噪声鲁棒。平均曲率H的阈值我一般取所有牙龈顶点曲率分布的前5%分位数作为边界起始点再往下生长到曲率降到中位数附近停止。不同扫描仪出来的网格平滑度不一样阈值最好做成参数暴露在外层配置里。4.2 外轮廓边界环的抽取与修补有了曲率标记的牙龈网格外轮廓就是牙龈网格的边界边。在网格拓扑里边界边是只属于一个三角面片的边这个性质跟曲率无关直接遍历面片统计边出现次数即可。from collections import Counter edge_count Counter() for face in mesh.faces: for i in range(3): e tuple(sorted((face[i], face[(i1)%3]))) edge_count[e] 1 boundary_edges [e for e, cnt in edge_count.items() if cnt 1] # 把边界边连成环 def build_loops(boundary_edges): adj defaultdict(list) for u, v in boundary_edges: adj[u].append(v) adj[v].append(u) loops [] visited set() for start in adj: if start in visited: continue loop [] curr start prev None while curr not in visited: visited.add(curr) loop.append(curr) nbrs [n for n in adj[curr] if n ! prev] if not nbrs: break prev, curr curr, nbrs[0] loops.append(loop) return loops边界环抽出来之后往往不是一条干净曲线会有小锯齿、分叉、断点。原因有两个分割阶段丢弃了跨标签面片导致边界上出现凹陷部分区域牙龈和牙齿靠得太近曲率阈值没切干净。我的处理顺序是先去掉过短的环——长度小于总边界长度1%的环直接扔掉然后把主环上的相邻点做滑动平均平滑最后对环上的缺口做线性插值补点。4.3 轮廓点序列的输出与下游对接外轮廓计算的最终产物是一组有序点序列每个点有原始三维坐标。输出格式上我做两种一种是按原始网格顶点索引编号方便下游直接引用网格顶点属性另一种是导出成独立的多段线点云文件可以直接加载进CAD或者MeshLab做对比。对于牙龈外轮廓实际使用时还会要求把轮廓点按逆时针方向排序并且每隔一个固定间距采样。我对比过不同算法最后确定的方案是先把边界环点投影到最佳拟合平面在二维平面里按角度排序再映射回三维坐标。这样就算原始边界环的顶点顺序是乱的也能重新整理成有序轮廓。输出时记录每个点在原网格上的顶点索引这样如果需要算轮廓长度或曲率可以直接读取。5. 牙齿分割避坑指南五条高频翻车记录5.1 投影重叠导致相邻牙粘连现象栅格化掩膜上两颗相邻牙连成一个连通域标签数比实际牙数少一颗。原因牙弓展开不彻底前牙区弧度大的位置展开后仍然存在局部重叠投影点落进了同一个栅格单元。解决不要用全局PCA拟合牙弓方向改成沿牙弓中心线做逐段弧长展开如果展开后仍然粘连把栅格分辨率从0.2mm降到0.15mm同时减少形态学闭运算的迭代次数。5.2 退化三角形让面积计算直接报错现象程序跑到一半报ZeroDivisionError或者输出nan定位后发现是某些三角形的面积为零。原因扫描件里存在共线或重复顶点的退化三角形计算法向时归一化除零。解决预处理阶段用np.linalg.norm(np.cross(v1-v0, v2-v0))检查每个三角形面积小于1e-12的面片直接剔除如果剔除后面片数量损失超过2%说明网格质量问题严重应回到原始扫描件重新做网格修复。5.3 栅格分辨率选错导致断裂与误合并现象同一颗牙的牙冠在掩膜上碎成三四个连通域或者相邻牙之间明明有沟壑掩膜上却连成一片。原因分辨率对网格密度不匹配。分辨率太粗沟壑被栅格单元直接吞没分辨率太细扫描噪声在掩膜上形成大量离散点把牙齿分成碎片。解决先统计网格边长的中位数median_edge_len把分辨率设为它的2倍作为初始值然后用一个交互式滑块快速预览掩膜效果调到所有牙齿刚好分离且无碎片时锁定参数。这个参数建议写入配置文件不要硬编码。5.4 曲率阈值一刀切导致轮廓断裂现象牙龈外轮廓环在尖牙和磨牙区域出现断裂人工数了数缺口有四五处。原因牙龈形态本身曲率差异大切牙区过渡平缓、磨牙颊侧过渡陡峭固定阈值只能切出一部分边界。解决改用相对阈值——先计算所有牙龈顶点的曲率分布取75%分位数作为高阈值、中位数作为低阈值在高低阈值之间做区域生长形态陡峭区域从高阈值出发向下生长平缓区域从低阈值出发向上生长双向逼近后衔接边界环。5.5 标记号映射回网格时索引错位现象分割出来的牙齿网格面片数量异常有的牙只剩一半有的牙混进了邻牙面片。原因投影阶段顶点坐标被修改过但顶点索引没有同步或者cell_to_verts字典里同一个顶点被多次赋予不同标签后写入的标签覆盖了先前的。解决在预处理阶段就固定顶点索引不变投影展开计算出的坐标单独存数组不写回网格标签映射时用if vert_labels[vert_idx] 0: vert_labels[vert_idx] label做首次赋值保护避免覆盖。添加一行断言检查每个顶点至少有一个标签且标签为0的顶点数量不超过总顶点数的0.1%。6. 分割质量与轮廓的验证方法用几何指标代替肉眼检查6.1 边界环闭合性检查分割和轮廓算法跑完之后第一件事不是肉眼盯着看渲染效果而是做量化检查。闭合性是最基础的验证把每条轮廓环的起点和终点连起来计算几何距离如果大于平均边长的2倍说明环没闭合需要回到边界环修补环节重新跑。我习惯把检查逻辑写成一个独立脚本输入是分割结果目录输出检查报告。# 闭合性检查 def check_loop_closed(loop, mesh, tol0.5): start mesh.vertices[loop[0]] end mesh.vertices[loop[-1]] dist np.linalg.norm(start - end) return dist tol # tol单位mm这个检查每次分割后强制跑一遍。尤其当边界环经过平滑处理后首尾点会被拉离原位置检查能提前发现参数是否把轮廓改得变形了。6.2 曲率与法向一致性校验分割后的单颗牙网格法向应该全部朝外且过渡连续。我的验证做法是统计相邻面片法向夹角超过30度的面片对数占总面片数的比例。比例高于5%说明网格表面存在明显褶皱多半是分割边界切得不齐导致的。同时计算每颗牙的平均曲率和牙龈区域的曲率均方差如果两颗相邻牙的曲率方差差距太大说明其中一颗可能混入了牙龈面片。6.3 结果导出的标准化做法分割完成后我统一导出三个文件原始网格的标签顶点色每个顶点一个整数标签、分离好的单颗牙STL列表、牙龈外轮廓点云。标签顶点色用mesh.visual.vertex_colors写入这样在MeshLab里可以直接按颜色快速检查分割结果。单颗牙STL命名用牙位编号方便对接下游排牙算法。轮廓点云导出成CSV加一行文件头说明坐标系和单位这能省掉后续和别的模块对接时的沟通成本。这套流程跑通的第一个版本也是各种翻车后面把预处理、投影、栅格化、轮廓修补每个环节的检查点都固化成了脚本才真正做到“跑完就能交付”。从那以后我每次处理牙齿STL分割都强制走一遍闭合性检查和法向一致性校验分不清是算法问题还是参数问题的情况少了大半。希望帮到你。本文还有配套的精品资源点击获取