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

发布时间：2026/10/26 23:59:01
投影栅格化实现牙齿STL网格分割与牙龈轮廓提取

简介：面向三维建模、计算机图形学与口腔医疗影像处理领域的牙齿STL网格模型分割算法资源，核心通过投影（曲面栅格化）实现牙齿与牙龈分离，并计算牙龈外轮廓，适合算法研究者、MATLAB开发者及数字化牙科技术爱好者参考实践。压缩包共54个文件，以40个m格式的MATLAB源码为主，涵盖STL模型读取、坐标对齐、栅格化、边缘检测、区域填充、轮廓拟合等模块，并附带4个java辅助文件、4个mat数据文件以及md说明文档，整体大小18.37MB。已有595人学习下载，可见其实用价值。使用者可获得完整的牙齿分割算法实现思路与可运行脚本，包括投影栅格化、Canny/Sobel边缘检测、阈值/区域生长分割以及外轮廓计算等关键步骤的代码示范，并借助内置数据与演示脚本快速验证效果，便于进一步研究或迁移到相关三维网格分割任务。代码结构清晰，模块划分明确，方便按需复用与二次开发。

1. 投影式牙齿分割：为什么栅格化比纯几何切割省一个量级

拿到一颗带牙龈的牙齿 STL 网格模型，第一反应往往是直接做三维几何分割——用曲率、法向聚类、区域生长去切牙冠和牙龈。真跑起来就发现不对：牙齿之间的邻接面是凹曲面，法向变化平缓，曲率噪声大，纯几何方法在牙缝处几乎必然过分割或欠分割。这个项目换了个思路：把三维网格投影到二维平面做栅格化，把分割问题降维成图像处理问题，最后再反投影回三维取牙龈外轮廓。投影分割算法比直接三维几何切割在工程上省一个量级，代码好写、参数好调、边界稳定，尤其适合 Qt5.15.2 + C++ 这类桌面工具链。下面按我的拆解顺序，把 STL 读取、曲面栅格化、区域生长分割、轮廓反投影和踩坑记录全部展开。

2. STL 读取与邻接重建：只靠三角面片也能算牙龈轮廓

2.1 STL 格式的隐藏缺陷

STL 文件本质是三角形面片的无序集合，每个面片包含三个顶点坐标和法向量，顶点按面片重复存储，没有拓扑关系。做牙龈外轮廓计算时，你必须知道哪些边是边界边、哪些顶点是相邻顶点，这个信息 STL 文件里不直接给，得自己重建。常见做法是读入后先做顶点去重，再构建半边结构或邻接表。

我在 Qt 环境里一般直接用QFile读二进制 STL，因为二进制比 ASCII 省一半空间，解析也简单：

struct STLHeader { char info[80]; quint32 triCount; }; struct STLVector { float x, y, z; }; struct STLTriangle { STLVector normal; STLVector v[3]; quint16 attr; }; bool readBinarySTL(const QString& path, std::vector<STLTriangle>& tris) { QFile f(path); if (!f.open(QIODevice::ReadOnly)) return false; QDataStream in(&f); in.setByteOrder(QDataStream::LittleEndian); STLHeader hdr; if (in.readRawData(hdr.info, 80) != 80) return false; in >> hdr.triCount; tris.resize(hdr.triCount); for (quint32 i = 0; i < hdr.triCount; ++i) { in.readRawData((char*)&tris[i], sizeof(STLTriangle)); } return true; }

这段逻辑很直白：跳过 80 字节文件头，读三角形数量，然后一次性把三角形数组读入。需要注意的坑是 STL 二进制格式在 Windows 和 Linux 下字节序都是小端，但QDataStream默认在大端机器上可能出错，所以显式setByteOrder是必须的。而且STLTriangle结构体可能有 padding 对齐问题，稳妥做法是用QDataStream逐个字段读，或加#pragma pack(push,1)。我从 Qt 里读取 STL 后通常会顺手检查法向量长度，离 1 偏差超过 0.01 的面片直接标记，否则后续栅格化时深度值会偏。

2.2 顶点去重与邻接表构建

STL 文件里同一个顶点会在多个三角形中重复存储，比如一个典型的牙齿模型里，顶点重复次数平均是 6 次。如果不去重，牙缝两侧的面片会被当成完全独立的网格，牙龈轮廓根本连不上。我们去重的方式是用一个哈希表记录顶点坐标到索引的映射：

struct Vec3Hash { size_t operator()(const STLVector& v) const { size_t h1 = std::hash<float>()(v.x); size_t h2 = std::hash<float>()(v.y); size_t h3 = std::hash<float>()(v.z); return h1 ^ (h2 << 1) ^ (h3 << 2); } bool operator()(const STLVector& a, const STLVector& b) const { return std::abs(a.x - b.x) < 1e-6 && std::abs(a.y - b.y) < 1e-6 && std::abs(a.z - b.z) < 1e-6; } }; std::unordered_map<STLVector, int, Vec3Hash, Vec3Hash> vertexMap; std::vector<STLVector> vertices; std::vector<std::array<int,3>> triVerts; for (auto& t : tris) { std::array<int,3> idx; for (int k = 0; k < 3; ++k) { auto it = vertexMap.find(t.v[k]); if (it == vertexMap.end()) { int newIdx = (int)vertices.size(); vertices.push_back(t.v[k]); vertexMap[t.v[k]] = newIdx; idx[k] = newIdx; } else { idx[k] = it->second; } } triVerts.push_back(idx); }

去重之后再做邻接表就轻松了：对每条边(a, b)建立一个键，映射到它所属的三角形列表，如果一个边只被一个三角形使用，它就是网格边界边。牙龈外轮廓本质上就是边界边的集合，这个计算依赖顶点去重和边邻接。这里有个参数值得注意：哈希比较器里1e-6的容差，扫描仪生成的 STL 顶点坐标一般能精确到 0.001mm，容差设 1e-6 是安全的，但如果输入模型经过网格简化，顶点坐标会有微小漂移，容差可以放宽到 1e-5，再大就可能把相邻牙冠的顶点错误合并，后面分割就会翻车。

2.3 法向量统一：投影分割的前置条件

STL 文件里的法向量经常是乱的，有的面片朝外、有的朝内，特别是从口腔扫描软件导出的模型，法线不一致是常事。栅格化投影虽然不关心法向量方向，但后面区域生长分割要用法向量夹角做判据，法线乱成一团就没法算了。

解决办法是遍历所有三角形，检查它的法向量和三个顶点的平均法向是否一致，不一致就翻转。常见做法是从任意一个三角形开始 BFS，沿公共边传播，每到一个新三角形，检查它与前一个三角形在公共边上的顶点顺序是否一致。如果两个三角形共享边(a,b)，但一个按a->b方向、另一个按b->a方向，说明邻接三角形法线翻面了，把它反转。这一步做完后，后续分割和轮廓提取才会稳定。这也是为什么我每次拿到 STL 都比比例子花多十分钟先做法线矫正——血泪经验，省不掉。

3. 曲面栅格化与深度投影：把三维网格变成二维问题

3.1 投影平面的选择与建立

栅格化的核心是把三维曲面降维到二维平面。对牙齿模型，最自然的投影方向是咬合面法向。你可以用全部顶点的协方差矩阵求主成分方向，这个方向通常就是牙弓的竖直轴。也可以让用户在界面上点三个牙尖点，拟合一个平面作为投影面，后者更可控，因为每颗牙的生理角度有差异，全自动拟合偶尔会被歪斜的磨牙带偏。

我采用的方案是先算全体顶点协方差矩阵的最小特征值对应的特征向量作为投影方向，再取模型质心作为平面原点。投影平面的二维基向量就取另外两个特征向量。这样得到的是一个正交投影映射P : R3 -> R2。栅格分辨率的确定要结合牙冠实际尺寸：常见做法是把投影平面划分成W x H的网格，每个格子对应的物理尺寸取 0.15mm 到 0.25mm。分辨率设太粗，牙缝细节丢失；太细，后文深度图的噪声会变得很大，分割算法容易产生碎块。

3.2 三角形栅格化与深度缓冲

把三角形投影到二维平面后，需要遍历它覆盖的像素格子，写入深度值和面片索引。这一步是算法最核心的计算密集部分。我用扫描线法填充三角形，对每个像素用重心坐标判断是否在三角形内，并且维护一个深度缓冲：

void rasterizeTriangle(const Tri2D& t2d, int imgW, int imgH, float depthA, float depthB, float depthC, int triIdx, float* depthBuf, int* faceBuf) { // 包围盒裁剪 int minX = clamp((int)floor(std::min({t2d.x0, t2d.x1, t2d.x2})), 0, imgW-1); int maxX = clamp((int)ceil (std::max({t2d.x0, t2d.x1, t2d.x2})), 0, imgW-1); int minY = clamp((int)floor(std::min({t2d.y0, t2d.y1, t2d.y2})), 0, imgH-1); int maxY = clamp((int)ceil (std::max({t2d.y0, t2d.y1, t2d.y2})), 0, imgH-1); float invDenom = 1.0f / ((t2d.x1 - t2d.x0) * (t2d.y2 - t2d.y0) - (t2d.y1 - t2d.y0) * (t2d.x2 - t2d.x0)); for (int y = minY; y <= maxY; ++y) { for (int x = minX; x <= maxX; ++x) { float px = x + 0.5f, py = y + 0.5f; float w0 = ((t2d.x1 - px) * (t2d.y2 - py) - (t2d.y1 - py) * (t2d.x2 - px)) * invDenom; float w1 = ((t2d.x2 - px) * (t2d.y0 - py) - (t2d.y2 - py) * (t2d.x0 - px)) * invDenom; float w2 = 1.0f - w0 - w1; if (w0 < 0 || w1 < 0 || w2 < 0) continue; float depth = w0 * depthA + w1 * depthB + w2 * depthC; int idx = y * imgW + x; // 只保留离投影面最近的片元 if (depth < depthBuf[idx]) { depthBuf[idx] = depth; faceBuf[idx] = triIdx; } } } }

这份代码里depth用的是投影方向上的有向距离，而不是到投影面的绝对距离。因为牙齿的牙尖高、牙龈低，深度越小代表越靠近观察者，种子点也就是区域生长的出发点，一般是深度局部极小值。重心坐标的 w0/w1/w2 会算成负数，说明像素不在三角形内，直接跳过。

这里有两个关键参数：一是采样点取像素中心x+0.5, y+0.5，避免像素角点判断造成锯齿；二是深度比较用<，这意味着如果两个三角形在像素上有重叠，只保留最靠外的那个，天然得到一个深度缓冲，这就是曲面栅格化算法的精髓。计算结果就是一张单通道深度图和一张面片索引图，索引图是后面区域生长的输入。如果想在 QOpenGLWidget 里实时查看效果，可以把深度图上传成纹理，灰度映射后用着色器平滑显示。栅格化这一步大约处理 20 万面片的模型，单线程在 0.5 秒内能完成，如果模型到百万面片级别，可以加 OpenMP 并行化，但注意深度缓冲写入有竞争，要按像素行分区或使用原子操作。

3.3 空洞填充与形态学预处理

扫描得到的 STL 模型在牙缝深处、牙龈舌侧经常有空缺，栅格化之后深度图会带孔洞。这些孔洞如果不处理，区域生长分割时会把一个牙冠分裂成多个区域，或者牙龈轮廓线断成几截。

常见做法是对深度图做形态学闭运算。先对深度图做膨胀，把边缘像素扩展到孔洞里，再用腐蚀恢复。闭运算的核大小建议用 3x3，因为 5x5 会把牙缝里仅有的一个像素宽的深沟填平，导致相邻牙冠在分割时连在一起。做闭运算之前，最好先做一次中值滤波去掉孤立噪声点。中值滤波窗口设为 5x5 比较稳妥，能滤掉单像素的椒盐噪声，同时保留牙尖和牙龈的过渡边缘。

形态学预处理完成后，深度图从「原始浮点数据」变成「干净的二维地形图」。后续处理都在这张图上操作，处理结果直接决定牙龈外轮廓的准确性。这一步也是为什么我说投影算法比纯几何切割好调参——图像处理那套滤波、形态学方法全都复用了，OpenCV 的cv::morphologyEx可以直接上，但要注意 STL 网格的深度值尺度是毫米级，做形态学前把深度图归一化到 0-255 反而会损失精度，我都是直接对浮点图操作。

4. 区域生长分割与牙龈外轮廓：从种子点到边界线

4.1 种子点选取：牙尖自动定位

分割的第一步是找种子点。在深度图上，牙冠的牙尖是局部极小值点（深度最小，离投影面最近）。直接检测局部极小值会得到一大堆伪种子——牙缝底、牙龈凹陷处也符合这个条件。我用一个三重条件过滤：

  1. 深度值要小于该牙弓的平均深度减去两倍标准差，防止在过高位置选种子。
  2. 局部窗口尺度取 7x7 像素，窗口中心必须严格小于窗口内 90% 的像素，排除平坦区。
  3. 两个候选种子点之间的距离小于牙冠最小宽度（比如 5mm 或 25 像素）时，保留深度值更小的那个，抑制同一个牙冠上出多个种子。

种子点定位后，区域生长就在深度图上展开。代码思路：

struct Seed { int x, y; float depth; }; std::vector<Seed> findSeeds(const float* depthImg, int w, int h, float meanDepth, float stdDepth) { std::vector<Seed> seeds; float threshold = meanDepth - 2.0f * stdDepth; for (int y = 3; y < h - 3; ++y) { for (int x = 3; x < w - 3; ++x) { float d = depthImg[y * w + x]; if (d > threshold) continue; int count = 0; for (int dy = -3; dy <= 3; ++dy) for (int dx = -3; dx <= 3; ++dx) { float nd = depthImg[(y+dy)*w + (x+dx)]; if (nd <= d) ++count; } if (count < 36) continue; // 需要窗口内90%像素深度更低 // 检查已选种子距离 bool dup = false; for (auto& s : seeds) { int dist2 = (s.x - x)*(s.x - x) + (s.y - y)*(s.y - y); if (dist2 < 25 * 25) { // 约5mm间距 if (s.depth <= d) { dup = true; break; } else s = {x, y, d}; // 替换旧种子 } } if (!dup) seeds.push_back({x, y, d}); } } return seeds; }

窗口 7x7、阈值 90% 这两个参数要结合栅格分辨率调整。若分辨率是 0.2mm/像素，7x7 对应 1.4mm 范围，对正常磨牙的牙尖来说足够。如果是前牙区比较窄的牙冠，窗口可以缩到 5x5。种子点之间的去重半径 25 像素对应 5mm，比最小牙冠宽度略小——为的是避免磨牙的远中颊尖和近中舌尖两个牙尖被合并成一个，如果去重半径设太大，会丢失真正的牙尖种子。种子点找得准不准直接决定分割的好坏，我会在调试时把种子点叠加到原始 STL 上可视化，用 QOpenGLWidget 加载 STL 渲染半透明网格，种子点画成红色小球。这一步很值得做，因为深度图上的极小值未必是三维空间里的牙尖，投影方向偏了会出现投影重叠，看着像牙尖的像素其实不是。

4.2 二维区域生长的约束条件

种子点到位后，区域生长的判据比纯深度阈值要复杂。我只用单一深度阈值时，牙龈区域和牙冠底部经常连成一片，原因是牙冠到牙龈的过渡在深度图上是一个平缓斜坡，没有清晰的断崖。解决方案是加一个法向量约束——从面片索引图里取出当前像素对应的三角形法向量，如果两个相邻像素的法向量夹角超过 30 度，判定为边界，不生长。

区域生长用广度优先遍历：

struct PixelNode { int x, y; }; void regionGrow(const float* depthImg, const int* faceIdx, const std::vector<STLVector>& triNormals, const Seed& seed, int w, int h, float depthTol, float angleTolDeg, std::vector<int>& labelImg, int label) { std::queue<PixelNode> q; q.push({seed.x, seed.y}); labelImg[seed.y * w + seed.x] = label; const float cosTol = std::cos(angleTolDeg * M_PI / 180.0f); while (!q.empty()) { auto [cx, cy] = q.front(); q.pop(); int cidx = cy * w + cx; int cf = faceIdx[cidx]; STLVector cn = triNormals[cf]; float cd = depthImg[cidx]; for (int dy = -1; dy <= 1; ++dy) { for (int dx = -1; dx <= 1; ++dx) { if (dx == 0 && dy == 0) continue; int nx = cx + dx, ny = cy + dy; if (nx < 0 || ny < 0 || nx >= w || ny >= h) continue; int nidx = ny * w + nx; if (labelImg[nidx] != 0) continue; float nd = depthImg[nidx]; if (std::abs(nd - cd) > depthTol) continue; int nf = faceIdx[nidx]; if (nf < 0 || nf >= (int)triNormals.size()) continue; STLVector nn = triNormals[nf]; float dot = cn.x*nn.x + cn.y*nn.y + cn.z*nn.z; if (dot < cosTol) continue; labelImg[nidx] = label; q.push({nx, ny}); } } } }

两个约束条件配合效果是最好的：深度差限制在 0.5mm 内，法向夹角限制在 30 度内。牙冠牙尖处法向平滑，可以继续生长；到了牙冠和牙龈交界处，虽然深度变化不大，但法向突然转了一个大角度，生长就停了。牙龈区域本身法向变化平缓，但深度值比牙冠大不少，深度差约束会把牙龈排除在外。这两个阈值都要处以栅格分辨率换算成像素深度差，如果某颗牙有严重磨损或牙冠倾斜角度大，夹角阈值可以放宽到 45 度，但放宽了牙龈和牙冠的边界会模糊，后面提取出来的外轮廓会偏到牙龈上去。

4.3 轮廓提取与三维反投影

二维分割完成后，每个像素都有了一个标签，牙龈外轮廓就在分割区域和非分割区域的交界处。用 Marching Squares 从标号图提取等值线，或者更直接——用边界跟踪算法沿标签边界走一圈，得到像素级多边形。得到二维轮廓之后，把它反投影到三维平面：

STLVector unprojectPoint(float u, float v, const STLVector& origin, const STLVector& axisU, const STLVector& axisV) { STLVector p; p.x = origin.x + u * axisU.x + v * axisV.x; p.y = origin.y + u * axisU.y + v * axisV.y; p.z = origin.z + u * axisU.z + v * axisV.z; // 沿投影方向修正回网格表面 // 用深度图上存储的深度值，沿着投影方向移动 return p; }

反投影后，每个像素点的三维坐标都在投影平面上，但牙龈外轮廓实际上是三维曲面上的曲线，投影平面上的点会偏离真实网格表面。所以最常见的做法是从每个二维轮廓点出发，沿投影方向反向发射射线，与 STL 网格求交，交点才是真正的牙龈外轮廓点。这个求交可以用光栅化时记录的 faceIdx 加速：先在面片索引图里找到该点所在的面片，然后在该面片上做点-三角形最近距离投影。拉开牙冠之后，轮廓线就是一条绕牙颈部的封闭曲线，这正是需要导出的牙龈外轮廓。轮廓线平滑可以用 Chaikin 曲线磨角算法，迭代两次就够了，迭代太多次会往牙冠方向缩，量起来偏小。

轮廓提取出来之后，如果做正畸排牙或牙冠设计，这条封闭曲线就是下一步的裁剪边界。可以把轮廓线画回 QOpenGLWidget 里做验证，半透明显示 STL 网格，叠加一条高亮曲线，确认它确实绕在牙龈边缘，而不是跑到了牙冠中央或牙龈底端。

5. 常见问题排查：锯齿、空洞与过分割的五个踩坑记录

5.1 牙缝处出现过分割

现象：相邻两颗牙被分割成三四块，多出来的是牙缝里的小碎片区域。原因是深度差阈值设置太小，牙缝底部本身的深度比两侧牙冠低，单个像素的深度噪声就足以触发边界判定，生长被提前终止。解决方法是把深度差阈值从 0.3mm 上调到 0.5mm，同时在生长完成后加一步最小区域过滤——面积小于 2mm 的碎块直接并入相邻最大标签区域。不要用形态学腐蚀再膨胀去修，会让牙尖被磨平。

5.2 牙龈轮廓线锯齿严重

现象：轮廓线在牙颈部附近呈明显的阶梯状，一个像素一个台阶。原因有两层：栅格化分辨率不够，以及 Marching Squares 未做插值直接按像素格走线。解决方法是把栅格分辨率从 0.25mm 提升到 0.15mm，同时轮廓提取时用双线性插值计算等值线，对每个格子内的跨边界点插值到亚像素精度。我一般还会加一步 Savitzky-Golay 滤波，对轮廓坐标做平滑，窗口取 11 点，既保形又去锯齿。

5.3 牙龈舌侧轮廓出现断口

现象：投影后的深度图在舌侧有孔洞，分割到孔洞边缘停止，轮廓线不闭合。原因是扫描模型本身舌侧牙龈有缺损，或者投影方向选择导致舌侧凹面自遮挡。解决方法是栅格化之后不立刻使用深度图，先做 3x3 闭运算填孔，如果孔洞直径超过 20 像素，则用口腔内典型的牙龈厚度 1.5-2mm 做插值修补——取孔洞边缘的深度均值补进去，再做高斯平滑衔接。如果 STL 模型本身缺损严重，二维修补只能救急，最终轮廓精度还是得靠原始模型质量。

5.4 区域生长把牙龈和牙冠连在一起

现象：分割结果显示牙冠和牙龈是一个连通区域，轮廓线跑到牙龈根部去了。原因通常是法向量夹角阈值设得过大，牙颈部的法向偏转被放过了。解决方法分两步：先检查投影方向是否与牙齿长轴平行，如果投影方向偏斜超过 10 度，深度差约束会失效；然后把法向夹角阈值从 45 度收回到 25 度。如果还不能分离，再额外加一个深度绝对阈值——只允许生长在距离投影面 6mm 以内的区域，牙龈底部深度超过这个值就切掉。

5.5 同一颗磨牙出现两个分割区域

现象：一颗磨牙的颊尖和舌尖被分割成两个独立区域，多出来的轮廓线横穿牙冠中央。原因是我的种子点去重半径设小了，颊尖和舌尖各自生成了种子，而深度图上两个牙尖之间有一条浅沟，生长过程在那里汇合后又被深度约束卡断。解决方法不是调深度阈值——那是治标——而是把种子点去重半径加大到 8mm，让形态学上同属一颗牙的种子合并。更稳的做法是分割结束后做一次拓扑检查：如果两个区域的公共边界长度占各自边界长度的 30% 以上，强制合并。从那以后我每次做多牙弓分割都强制走一遍这个合并检查，磨牙再也不劈两半了。

6. 参数标定与验证技巧：投影分辨率与阈值组的一键标定

参数太多容易乱，我的习惯是把所有参数集中到一个结构体里，并且按物理尺寸定义，不用像素值，这样换模型时不会因为分辨率变了就失控。先定投影分辨率，方法是量 STL 模型的包围盒，取最长边对应的物理长度除以期望的像素数。期望像素数怎么定？看牙冠的最小特征——牙颈部的曲面转折宽度，这个宽度在真实模型里大约 2mm，投影后至少要有 8 个像素才能被法向约束识别出来，所以分辨率下限是 0.25mm。上限取 0.15mm，因为扫描模型本身的网格密度有限，投影超过网格密度就是纯插值出来的假细节。

深度差阈值和法向夹角阈值我建了一个小规模的快速回归集。用 5 颗不同形态的牙齿模型作为标准测试集，跑完全部流程后人工检查分割结果和轮廓线是否合理。调整的时候先锁定分辨率，再调深度差，最后微调法向夹角。每次调整都记录分割区域的数目和轮廓线的闭合性，区域数目偏离人工标注值超过 10% 就回滚。

验证轮廓精度时，我习惯把提取出的牙龈外轮廓投影回三维网格，计算轮廓上每个采样点到最近三角面片的距离，统计平均距离和最大距离。一个合格的轮廓提取，平均距离应该在 0.1mm 以内，最大距离不超过 0.3mm。如果最大距离超过 0.5mm，基本可以断定轮廓线穿过了牙龈表面到了软组织区域，需要调整投影方向或者加强法向约束。这套验证流程我现在每次跑分割都会执行一遍，已经变成肌肉记忆了。写代码的最后一步，是导出一条 CSV 记录分割后的牙冠编号、体积和轮廓线弧长，方便后续做正畸分析时直接对照。希望这篇拆解能帮你在做牙齿 STL 网格分割时少走几步弯路。

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