简介：面向三维建模、计算机图形学与医疗影像处理场景，这份代码资源围绕牙齿STL网格模型分割问题，实现了基于投影算法（曲面栅格化）的牙齿与牙龈分离，并支持计算牙龈外轮廓。资源以MATLAB脚本为主要实现，覆盖坐标系统对齐、三角面片栅格化、边缘检测、阈值分割/区域生长以及外轮廓拟合等关键步骤，适合算法研究者和牙科数字化应用开发者参考。压缩包为zip格式，共含54个文件，整体约18.37MB，除40个.m脚本外，还附带少量JAVA源码、.mat数据文件和项目配置文件，能够支撑复现与二次开发。已有595人学习浏览。通过阅读源码中的主流程和辅助函数，可以完整理解从STL文件读取、模型预处理到分割结果可视化的技术链路；配合内置数据与示例脚本，有助于快速验证算法效果，节省自建实验环境的时间。

1. 牙齿STL网格分割，为什么先从投影算法说起

一次全牙弓口扫导出的STL网格，通常有七八万到二十万个三角形。要做单颗牙分割、牙龈外轮廓提取，最笨的办法是在三维空间里直接做区域生长或图割，但结果往往被邻面接触区、倒凹和牙龈缝隙带偏，调参到怀疑人生。标题里的“投影算法（曲面的栅格化算法）”换了个思路：把三维网格映射到一个规则的二维栅格上，生成深度图和掩膜图，在二维上完成牙齿分割和牙龈轮廓计算，再把结果反投影回三维网格。这套方案速度比纯三维分割快一个数量级，而且每一步都能用图像处理工具检查中间结果，适合口扫数据、模型扫描数据和大多数需要批量处理牙位分割的齿科CAD/CAM场景。注意，这个STL是网格文件格式，不是C++标准模板库那个STL，我第一次接手时也在命名上绕了一下。

2. 栅格化之前的网格卫生：清洗与预平滑的关键参数

2.1 先修法向、非流形边和退化面片，否则投影图全是噪声

投影栅格化是在对三角形逐个采样，任何一个三角形的法向翻转、面积为零或边被两个以上面片共用，都会直接在深度图上留下空洞或毛刺。所以拿到STL后我不会直接做投影，而是先跑一遍几何检查。用trimesh做这一步非常快：

import trimesh import numpy as np from collections import Counter mesh = trimesh.load("patient_arch.stl", process=False) # 法向统一朝外，typo会直接影响投影深度方向 mesh.fix_normals() # 非流形边检查：一条边被超过2个三角面共用就是不合法 edge_face_count = Counter() for face in mesh.faces: for edge in [(face[0], face[1]), (face[1], face[2]), (face[2], face[0])]: edge_face_count[tuple(sorted(edge))] += 1 non_manifold_edges = [e for e, c in edge_face_count.items() if c > 2] print("非流形边数量:", len(non_manifold_edges)) # 删除面积为0的退化面片并清理孤立顶点 mesh.update_faces(mesh.nondegenerate_faces()) mesh.remove_unreferenced_vertices() print("watertight:", mesh.is_watertight) print("winding consistent:", mesh.is_winding_consistent)

fix_normals()在trimesh里的实现是按连通分量统一法向方向，保证所有三角形法向朝外或朝内一致。如果模型本身有多个互相穿插的连通体，这一步只能保证每个连通体内部一致，无法保证整体方向，所以我后面还会再校验一次投影方向。nondegenerate_faces()会删掉面积接近零的面片，这类退化面片在牙齿扫描数据里常见于牙尖顶部和牙龈沟底部，扫描仪重建时会产生大量狭长三角形，如果不删，栅格化时它们的重心投影不稳定，会在深度图上产生孤立噪点。

2.2 Laplacian平滑的迭代次数与邻域半径：磨皮但不能磨掉牙尖

扫描网格表面通常有高频噪点，尤其是牙龈边缘和牙颈部。直接在噪点上做投影，深度图的梯度会乱跳，后续找牙龈轮廓时会找到一堆假边界。我一般会做一轮拉普拉斯平滑，但不会用大数据包自带的默认参数，而是自己控制迭代次数和松弛因子：

import numpy as np from scipy import sparse from scipy.sparse.linalg import spsolve def laplacian_smooth(vertices, faces, iterations=10, lam=0.6): v = vertices.copy() n = len(v) rows, cols = [], [] for f in faces: rows.extend([f[0], f[1], f[1], f[2], f[2], f[0]]) cols.extend([f[1], f[0], f[2], f[1], f[0], f[2]]) data = np.ones(len(rows)) A = sparse.coo_matrix((data, (rows, cols)), shape=(n, n)).tocsr() deg = np.asarray(A.sum(axis=1)).flatten() deg[deg == 0] = 1 D_inv = sparse.diags(1.0 / deg) P = D_inv @ A for _ in range(iterations): v = v + lam * (P @ v - v) return v

lam取0.6、迭代10次是我在口腔扫描数据上的默认值。lam大于1时迭代不收敛，网格会震荡发散；小于0.3时收敛太慢，10次迭代基本看不出平滑效果。这里要提醒一下，纯Laplacian平滑会让网格整体收缩，牙尖高度会稍微变矮，对分割影响不大，但如果你后续要做修复体边缘贴合，收缩量就不能忽略。想避免收缩，就用Taubin平滑，形式上是一样的，只是迭代时交替用正的lam和负的mu，一正一负相互抵消收缩趋势。

2.3 投影前先把坐标系摆正，后面的方向二义性少一半

这一步容易被跳过，但很关键。投影算法需要指定一个投影方向，也就是把网格往哪个平面上压。齿科扫描数据导出时坐标系往往来自扫描仪本身，不同设备差异很大。我一般先做一次粗配准：用所有顶点的PCA分析，把最小特征值对应的特征向量作为近似牙体长轴方向，旋转到Z轴上。这样后续做咬合平面拟合时初始值就准了。

points = mesh.vertices centroid = points.mean(axis=0) cov = np.cov((points - centroid).T) eigvals, eigvecs = np.linalg.eigh(cov) # eigvecs[:, 0] 对应最小特征值，近似牙体长轴方向 axis = eigvecs[:, 0] if axis[2] < 0: axis = -axis

这里有一个常见翻车点：PCA的eigh返回的向量符号是任意的，你无法确定axis指向牙尖还是牙根方向。我加了一个启发式规则：齿科扫描数据导出后Z轴向上是大概率事件，所以如果axis[2] < 0就反转。这个规则对多数扫描仪靠谱，但对本身坐标系就奇怪的模型不适用，所以我在下一章做了更严格的投影方向校验，而不是完全信任这一部。

3. 曲面的栅格化算法：投影方向、像素尺寸与深度场生成

3.1 投影方向的两条路线：PCA拟合与咬合平面约束

投影方向是整个栅格化算法的灵魂。方向选错，相邻牙齿在投影图上会重叠，牙龈沟和牙间乳头全部糊在一起，后面再怎么调分水岭参数都救不回来。

常见做法有两种。第一种是直接用PCA得到的最小特征向量作为投影方向，优点是全自动、不需要任何解剖先验，缺点是遇到牙弓形态特殊（比如重度磨损）时，PCA主轴会被整体形状带偏。第二种是拟合咬合平面，用咬合平面上远离牙龈区域的那些顶点做最小二乘平面拟合，取平面法向作为投影方向。我一般在预处理阶段已经把模型粗旋转到牙体长轴与Z轴平行，所以直接取Z轴附近方向就行，但如果你想处理更一般的数据，投影方向要用咬合平面约束来得稳。

这里有个工程直觉：投影方向决定了相邻牙接触区在图像上的表现形态。如果从咬合面垂直往下投影，邻面接触区是一个窄条，在栅格分辨率不够时容易连成一片；如果从偏颊侧方向投影，邻面会拉开，更容易分割，但牙龈轮廓会变形。我的做法是分开处理：做单颗牙分离时用偏颊侧5到10度的投影方向，做牙龈外轮廓提取时回到咬合面垂直方向。这样两个任务各自用最有利的视角。

3.2 栅格分辨率：先统计网格边长分布，再定像素尺寸

像素尺寸是栅格化里最需要认真对待的参数。分辨率定太高，单个三角形覆盖不到像素，掩膜上全是孔洞；定太低，牙尖窝沟细节丢失，邻面边界模糊。我不拍脑袋设固定值，而是先统计STL网格的平均边长：

edges = mesh.edges_unique edge_lengths = np.linalg.norm(mesh.vertices[edges[:, 0]] - mesh.vertices[edges[:, 1]], axis=1) mean_edge = edge_lengths.mean() res = 1.5 * mean_edge # 经验系数，后面按需微调

像素尺寸与网格边长最好保持同一个量级。齿科口扫数据平均边长通常在0.05到0.15mm之间，所以res落在0.08到0.2mm比较合理。给你一组我在全牙弓模型上测试过的对照：

注意，STL网格边长反映的是重建时的采样密度，不代表扫描仪的物理精度。有些扫描仪导出网格很细（0.02mm），但真实表面精度只有0.1mm，这时把像素设成0.02mm纯属浪费算力，还会让噪声更突出。我通常以最终需要的轮廓精度为准来定分辨率，牙龈外轮廓提取用0.1mm左右，单颗牙分离用0.15mm左右就够。

注意：分辨率不是越高越好，像素比三角形还小的时候，栅格化就退化成对顶点的随机采样，掩膜会出现大量单像素孔洞。闭运算可以补孔，但也会把细小的牙缝吃掉，这是分辨率与形态学操作之间最常见的矛盾。

3.3 三角形栅格化的最小实现：深度场、命中掩膜与像素-面片索引

理解了方向和分辨率，下面就是栅格化本身。我在这里不只存深度图，还额外维护一个“像素到面片”的索引列表，这是后面从二维结果反投影回三维网格的数据桥梁，很多人做投影分割时把这个信息丢了，最后反投影只能靠最近邻搜索，慢且不准。

import numpy as np def edge_function(p, a, b): return (b[0] - a[0]) * (p[1] - a[1]) - (b[1] - a[1]) * (p[0] - a[0]) def rasterize_triangle(pix_pts, fid, res, depth_img, mask_img, origin_pix, face_hits): xmin = max(0, int(np.floor(pix_pts[:, 0].min()))) xmax = min(depth_img.shape[1] - 1, int(np.ceil(pix_pts[:, 0].max()))) ymin = max(0, int(np.floor(pix_pts[:, 1].min()))) ymax = min(depth_img.shape[0] - 1, int(np.ceil(pix_pts[:, 1].max()))) area = edge_function(pix_pts[0], pix_pts[1], pix_pts[2]) for y in range(ymin, ymax + 1): for x in range(xmin, xmax + 1): p = np.array([x + 0.5, y + 0.5]) w0 = edge_function(p, pix_pts[1], pix_pts[2]) w1 = edge_function(p, pix_pts[2], pix_pts[0]) w2 = edge_function(p, pix_pts[0], pix_pts[1]) if area > 0: inside = (w0 >= 0 and w1 >= 0 and w2 >= 0) else: inside = (w0 <= 0 and w1 <= 0 and w2 <= 0) if not inside: continue # 重心坐标，用于插值深度 w0 /= area; w1 /= area; w2 /= area depth_val = w0 * pix_pts[0, 2] + w1 * pix_pts[1, 2] + w2 * pix_pts[2, 2] idx = y * depth_img.shape[1] + x face_hits[idx].append(fid) if depth_val > depth_img[y, x]: depth_img[y, x] = depth_val mask_img[y, x] = 1

这段代码有三个设计点需要说明。第一，像素中心取x+0.5, y+0.5而不是x, y，否则大量像素中心恰好落在三角形边缘上，会产生细密锯齿。第二，area的正负反映了三角形顶点在投影平面的绕向，法向翻转的三角形area为负，判断条件要相应反转，这也是为什么前面要先做法向一致化。第三，深度更新用的是max策略，也就是保存该像素上所有覆盖三角形中“最高”的那个深度值，这相当于一个Z-Buffer，含义是沿着投影方向看过去，表面最外层的高程；如果任务是检测牙龈沟等凹陷结构，就改用min策略存最低深度。

face_hits这个列表数组是反投影的命脉。每个像素记录了哪些三角形覆盖过它，后面做分割得到像素标签后，只需按面片投票就能把二维标签还原到三维网格上，一次遍历搞定，不需要任何空间索引或最近邻查询。

4. 二维深度图上的分割算法：牙齿分离与牙龈外轮廓提取

4.1 牙体与牙龈的掩膜分割：深度梯度与曲率辅助

投影栅格化之后得到的是一张二维深度图，接下来第一步不是急着分牙，而是先把牙体区域和牙龈区域分开。在深度图上，牙体与牙龈的交界处有明显的深度梯度变化：牙颈部以下深度突然变浅（从咬合面视角看，牙龈表面比牙冠低），而且这个过渡带在颊侧和舌侧都有。

我一般会先对深度图做高斯模糊，然后计算梯度模长，再用Otsu自动取阈值：

from scipy import ndimage as ndi from skimage.filters import threshold_otsu, gaussian blur_depth = gaussian(depth_img, sigma=1.5) grad = np.hypot(*np.gradient(blur_depth)) thresh = threshold_otsu(grad[grad > 0]) body_mask = grad < thresh body_mask = ndi.binary_opening(body_mask, iterations=2) body_mask = ndi.binary_closing(body_mask, iterations=3)

sigma=1.5是经验值，单位是像素。如果像素尺寸0.1mm，这个模糊半径大约是0.15mm，刚好能压住表面噪点而不抹掉牙龈沟。梯度阈值用Otsu自动求，比手工给百分比经验值可靠，因为不同病例的磨损程度、牙龈形态差异很大。形态学开运算是去掉牙尖和牙龈边缘的孤立噪点，闭运算是把邻面接触区的细缝和牙颈部过渡带的小断裂接起来。

这个掩膜就是“牙体区域”。后面分水岭只在body_mask里做，牙龈区域直接归为背景，这样可以避免牙龈缝隙把牙齿分割结果带偏。

4.2 单颗牙分离：距离变换分水岭与最小牙尖间距

得到牙体掩膜后，要做的是把每一颗牙从牙弓上切下来。这里我用距离变换加分水岭，思路是：先把牙体掩膜做距离变换，每个牙冠中心是距离值局部极大值，把这些极大值点作为标记，然后对距离场的负值做分水岭分割。

from scipy import ndimage as ndi from skimage.feature import peak_local_max from skimage.segmentation import watershed dist = ndi.distance_transform_edt(body_mask) dist_smooth = gaussian(dist, sigma=1.0) coords = peak_local_max( dist_smooth, min_distance=int(round(2.5 / res)), # 2.5mm约等于前磨牙颊舌径的一半 exclude_border=False, ) markers = np.zeros_like(dist_smooth, dtype=np.int32) for i, (r, c) in enumerate(coords): markers[r, c] = i + 1 labels = watershed(-dist_smooth, markers, mask=body_mask)

min_distance是我在牙弓数据上最敏感的参数，它的含义是两个标记点之间的最小像素距离，换算成物理尺寸约等于两颗牙中心的最小间距。口扫数据里相邻牙中心距一般在4到6mm，我取2.5mm作为下限，太大会漏掉前牙区的相邻牙，太小会把同一颗牙的牙尖和牙窝误判成两颗牙。如果你在前牙区遇到粘连分割不开，优先调大min_distance而不是去改分水岭本身。

分水岭切不开邻面接触区是常见问题。原因是真实扫描里相邻牙的邻面接触区完全贴合，深度图上没有凹陷，分水岭没有“岭”可走。我这边有两条后路：一是对body_mask做一次轻微腐蚀，让邻面区域先断开，再分水岭，代价是会吃掉牙间乳头的细节；二是先检测邻面接触区的位置（曲率突变的狭长区域），把那里在距离场上强制拉低，人为制造一条分水岭，这样不损失几何细节。后者效果更好，但实现量略大。

4.3 牙龈外轮廓提取：掩膜边界追踪与牙间乳头的闭运算

牙齿分离完成后，牙龈外轮廓计算就变成图像处理里的边界追踪问题。牙体掩膜的边界就是牙体与牙龈的交界线，也就是临床说的牙龈边缘线。直接find_contours出来会有两个问题：一是边界锯齿严重，二是牙间乳头区域掩膜可能断裂，导致轮廓不闭合。我的处理流程是先闭运算保证掩膜连通，再取0.5等值线作为轮廓：

from skimage.measure import find_contours closed = ndi.binary_closing(body_mask, iterations=5) contours = find_contours(closed.astype(float), 0.5) # contours 是像素坐标折线，按长度取最外层那条作为牙弓外轮廓 outer = max(contours, key=len)

取0.5等值线而不是1.0，是为了让轮廓落在掩膜边界像素的中间位置，这样反投影回三维网格时误差是对称的，后续光顺不会系统性内缩。binary_closing迭代次数我控制在3到5次，太多会把牙间乳头整个填平，太少又补不上真实的断裂。这里有个判断标准：闭运算后轮廓数量应该等于1或2（牙弓内外两条），如果还是多条，说明牙间乳头断得太厉害，这时候要回去调栅格分辨率而不是继续暴力闭运算。

5. 投影分割的五个翻车现场：现象、原因与解决方案

5.1 法向与投影方向导致的两种“伪深度场”

现象一：深度图上出现黑色空洞和孤立小岛，掩膜像撒了一把芝麻。

原因不是算法bug，而是模型里存在法向不统一的三角形。fix_normals()只保证连通分量内部一致，但牙冠和牙龈如果因为扫描拼接误差被分成了多个独立连通体，部分连通体的法向可能整体朝下。栅格化时我用max策略更新深度，朝下的三角形更新的是“最低深度”而不是“最高深度”，结果在牙体表面形成了凹陷洞。

解决办法是栅格化前先做一次“投影方向一致性校验”：把所有三角形的法向投影到投影方向上，如果超过2%的三角形的法向与主方向夹角大于90度，就把整个模型的法向翻转。这个校验比任何修复函数都直接。

现象二：整个深度图左右颠倒，牙冠看起来像牙槽骨。

原因是PCA求主方向时的符号二义性。eigh返回的特征向量方向是任意的，如果最小特征值对应的向量指向了牙根方向而不是牙尖方向，栅格化后深度场就反了，后续梯度阈值和轮廓提取全部错位。

解决办法是加一个语义校验小步骤：深度图生成后，统计牙体掩膜中心的深度均值与边界区域的深度均值，正常情况下牙体中心比边缘高（深），如果反了就把投影方向取反重新栅格化。这个校验成本极低，但能避免后面所有步骤白做。

5.2 分辨率、邻面接触区与反投影断裂

现象三：像素尺寸设置过小，掩膜出现大量单像素孔洞，闭运算后邻面又被填上了。

原因我在3.2里提过：像素比三角形还小的时候，栅格化等价于对顶点随机采样，覆盖率天然达不到100%。闭运算确实能补孔，但补完的代价是牙缝消失，分水岭更切不开。

解决办法是回到平均边长的统计上，让res不小于1.0倍平均边长，同时微调闭运算的迭代次数。如果必须用小像素捕捉细节，那就反过来用大的res生成粗掩膜、用小的res生成细节图，两者叠加，而不是指望单一分辨率一步到位。

现象四：相邻牙之间没有分水岭，整片牙弓连在一起。

原因是邻面接触区在投影方向上完全重合，深度图上没有凹陷。这跟5.2现象三不一样，是几何形态问题而不是分辨率问题。强行增大min_distance会把前牙也漏掉，强行腐蚀掩膜又会吃掉牙尖。

我的做法是检测邻面接触区的脊线。先计算深度图的局部主曲率，在牙体掩膜内找曲率负值最强（即凸脊）且方向垂直于牙弓走向的像素带，把这些像素在距离场上赋零值，强制分水岭从这里走。这个方案对拥挤牙列也有效，但曲率计算对噪声敏感，要先做高斯滤波。

现象五：分水岭结果反投影回三维网格后，边界是开环的，单颗牙的边界不闭合。

原因不在二维分割，而在反投影投票逻辑。边界像素上命中的三角形，有些是牙冠的，有些是牙龈的，有些是两者过渡区共用的。按像素标签直接投票时，过渡区面片票数往往接近五五开，argmax后归属随机，导致网格边界呈锯齿状且不闭合。

解决办法是面片投票前先做一个过渡面片标记：设某面片的三个顶点投影到像素后，既有牙冠标签又有牙龈标签，就把它标记为“边界过渡面片”，不参与投票。最后沿着过渡面片带，用Dijkstra求一条连接牙冠区域边界起终点的最短路径，这样闭合边界一定存在且连通。

6. 把二维轮廓映回三维空间：坐标复用与牙龈线样条光顺

6.1 坐标复用：一张变换表把像素点变成三维点

栅格化时保存的投影平面原点、两个正交基向量和像素尺寸，在这一步派上大用场。从轮廓像素坐标恢复三维点时不需要任何最近邻查询，直接做线性变换：

def pixel_to_3d(r, c, uv_origin, u_vec, v_vec, n_vec, res): u = (c * res) + res / 2.0 v = (r * res) + res / 2.0 return uv_origin + u * u_vec + v * v_vec

这里的uv_origin是投影平面原点，u_vec和v_vec是平面内两个正交单位向量。轮廓点从二维像素坐标映射回三维空间后，得到的是一条三维折线，但它位于投影平面上，严格说还不在网格表面上，需要再沿投影方向投回去。常见做法是把折线上的每个点沿n_vec方向投影，和网格求交。这个过程的精度取决于前面栅格化的分辨率，通常能控制在0.1mm量级，对牙龈外轮廓足够。

6.2 等弧长参数化与B样条光顺：平滑度的取舍

三维折线直接输出会有锯齿，必须光顺。我习惯用等弧长参数化加三次B样条拟合，而不是简单的高斯平滑，因为B样条可以控制平滑度与贴合度的权重：

from scipy.interpolate import splprep, splev points = contour_3d # 三维折线点集 tck, u = splprep(points.T, s=len(points) * 0.01) smooth_line = splev(np.linspace(0, 1, 200), tck)

s是平滑因子，s越大曲线越光滑但离原始点越远。我一般取len(points) * 0.01到len(points) * 0.05，具体看牙龈沟的形态。牙龈沟深、牙间乳头尖锐的病例，s调小，保留曲率；牙龈平坦的病例，s调大，得到更流畅的边界线。注意牙间乳头处不要全局光顺，那里曲率变化剧烈，全局平滑会把乳头尖直接削平，所以在牙间乳头区域我会降低权重，或者把轮廓按牙位拆段后逐段光顺。

这套流程走到这里就是一个可用的最小闭环：网格清洗、投影栅格化、二维分割、反投影、轮廓光顺。我最早做反投影时省掉了face_hits索引，改成每拿到一个轮廓点就去网格上搜最近三角形，全牙弓轮廓几百个点跑了将近一分钟；后来老老实实把像素到面片的索引留下来，整个反投影降到毫秒级。这个教训让我养成了习惯，任何投影类方案都先把“坐标映射关系”当成一等数据保存，后悔药没吃过第二次。希望帮到你。

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