☰
牙齿STL网格分割:投影与栅格化的高效分牙方案
2026/10/3 4:53:53 网站建设 项目流程

简介:面向牙科数字化模型处理,这套资源提供了基于投影与曲面栅格化的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 # 载入STL,merge_vertices=True会把距离小于tol的顶点合并 mesh = trimesh.load("tooth_arch.stl", force="mesh", merge_vertices=True) # 法向统一到朝向网格外部 mesh.fix_normals() # 强制计算顶点邻接关系,后续区域生长会频繁用到 mesh.vertex_adjacency_graph # 触发邻接图构建

这里的merge_vertices默认容差是tol=1e-8,如果扫描件本身顶点有微小的坐标抖动,我会放宽到1e-6,不然焊不干净。fix_normals()做的是让相邻面片法向夹角尽量小,再配合有向包围盒把整体翻转修正到朝外。法向统一之后,建议顺手把退化三角形(面积接近零的面片)剔除,否则后续投影时会出现除零和异常栅格点。

2.2 质心对齐与牙弓方向估计

STL模型从扫描设备出来后,坐标系是随机的,可能旋转了任意角度。直接对原始坐标做投影,基本上是灾难——牙齿会斜着摊开,邻牙在投影平面上重重叠叠,栅格化之后连成一整片。所以预处理里最关键的一步是把牙弓摆正:先计算质心并把模型平移到原点,再用PCA估计牙弓的主方向,把牙弓长轴旋转到X轴方向。

# 计算顶点质心并平移 vertices = mesh.vertices centroid = vertices.mean(axis=0) 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, axis=1) 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), dtype=np.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, iterations=2).astype(np.uint8) # 连通域标记 labeled, num_features = ndimage.label(mask)

分辨率我通常取网格平均边长的2~3倍。牙齿STL的扫描精度普遍在0.05~0.1mm,平均边长在0.2mm左右,所以resolution=0.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), dtype=int) 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[(i+1)%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, tol=0.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分割,都强制走一遍闭合性检查和法向一致性校验,分不清是算法问题还是参数问题的情况少了大半。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询