简介这是一份面向口腔医学图像处理与计算机图形学学习者的C课程设计资源基于VS 2017、Qt5.13.2与Vtk8.2.0开发实现三维牙齿模型自动化预处理系统涵盖牙齿分割、计数、编号、轴向标定及缺失识别等功能并将结果渲染至屏幕辅助牙医开展诊断。资源包共21个文件约8.42MB除5个STL牙齿扫描模型和8张处理效果示意图外还包含3个头文件、2个C源文件及README说明文档其中源码覆盖连通分量提取、去除牙龈、建立OBB包围盒、牙齿轴向计算等关键流程模型示例涵盖上下颌多组数据便于对照验证算法。目前已有265人浏览学习适合具备一定C与VTK基础、希望深入理解三维模型处理管线或完成类似课程设计的本科高年级学生及研究人员。通过阅读源码和效果图可快速掌握从原始扫描数据到牙齿分割、编号与缺失判断的完整处理思路并直接复用或扩展相关模块到自己的项目中。1. 三维牙齿预处理的自动化价值在哪里在口腔正畸和修复的数字化流程里口内扫描得到的整口牙 STL 模型是一个连续三角网格牙龈把所有牙齿连成一个整体后续排牙、模拟移动、设计牙套都需要先把单颗牙齿分离出来。手动在网格上画分割线一位医生处理一口牙往往要花十几分钟而且重复性操作容易在牙颈部边界出现主观差异。这个基于 C 的自动化预处理系统能够在读取完 STL 文件后自动完成牙齿分割、计数、编号、轴向标定和缺失识别并用不同颜色渲染到 Qt 界面中。开发环境是 VS2017 Qt 5.13.2 VTK 8.2.0核心处理全部在 VTK 管线中完成适合医学图像处理、VTK 二次开发和三维网格分析方向的工程师参考。2. 数据读取与牙弓方向预对齐2.1 STL 网格的读取与质量校验STL 文件是口内扫描最常见的网格格式VTK 提供了开箱即用的读取器。读取后首先要检查点数、单元数以及网格是否闭合因为后续的曲率计算和连通性分析都依赖高质量输入。#include vtkSTLReader.h #include vtkPolyData.h #include vtkCleanPolyData.h vtkNewvtkSTLReader reader; reader-SetFileName(jaw.stl); reader-Update(); vtkPolyData* mesh reader-GetOutput(); std::cout Points: mesh-GetNumberOfPoints() Triangles: mesh-GetNumberOfPolys() std::endl; vtkNewvtkCleanPolyData cleaner; cleaner-SetInputData(mesh); cleaner-SetTolerance(1e-6); cleaner-Update(); mesh cleaner-GetOutput();这段代码通过SetFileName指定输入文件Update触发读取管线。vtkCleanPolyData用来去除重复点和退化三角形SetTolerance(1e-6)控制合并点的距离阈值——单位是模型坐标单位过大会把正常顶点合并掉过小则无法清除重复点。网格规模通常会从几十万到上百万三角形打印日志能第一时间发现文件是否读错或坐标单位是否异常。2.2 用 PCA 确定咬合平面法向整口牙模型中所有牙冠的咬合面大致落在同一平面内而每颗牙齿的长轴方向与该平面近似垂直。对整个模型点云做主成分分析最小特征值对应的特征向量就是咬合面的法向。这个方向是所有后续操作的基础分割后的牙齿需要按该方向判断上下缺失识别也需要在咬合平面投影中计算间距。#include vtkMath.h double center[3]; vtkMath::ComputeCenter(mesh, center); int n mesh-GetNumberOfPoints(); double cov[3][3] {{0.0}}; for (vtkIdType i 0; i n; i) { double p[3]; mesh-GetPoint(i, p); double x p[0] - center[0]; double y p[1] - center[1]; double z p[2] - center[2]; cov[0][0] x * x; cov[0][1] x * y; cov[0][2] x * z; cov[1][1] y * y; cov[1][2] y * z; cov[2][2] z * z; } cov[1][0] cov[0][1]; cov[2][0] cov[0][2]; cov[2][1] cov[1][2]; for (int i 0; i 3; i) for (int j 0; j 3; j) cov[i][j] / n; double eigenval[3], eigenvec[3][3]; vtkMath::Jacobi(cov, eigenval, eigenvec); // eigenvec[2] 是最小特征值方向应为咬合面法向这里直接累加协方差矩阵然后调用vtkMath::Jacobi做特征分解。eigenval数组按特征值从大到小排列eigenvec[2]对应最小特征值方向。注意模型坐标一般是毫米协方差需要除以点数获得统计意义上的方差如果模型没有平移特征分解容易受数值误差影响所以要先用ComputeCenter去中心化。2.3 法线一致化与网格顶点的坐标归一化扫描得到的 STL 法线有时方向不一致三角片之间会出现法线翻转这会影响曲率计算符号。VTK 的vtkPolyDataNormals可以基于相邻三角片做一致化处理。同时把模型质心平移到原点会让后续 PCA 和连通分量提取更稳定。#include vtkPolyDataNormals.h vtkNewvtkPolyDataNormals normals; normals-SetInputData(mesh); normals-SplittingOff(); // 不按尖锐角切分 normals-AutoOrientNormalsOn(); // 自动翻转不一致的法线 normals-Update(); mesh normals-GetOutput();SplittingOff是关键很多扫描模型在尖锐边缘有较多折痕如果不关闭 Splitting法线计算会在锐边处产生分裂的拓扑导致后续曲率不连续。归一化时对每个点减去质心坐标即可注意原坐标需要保留一份用于最终渲染避免反复变换累积误差。3. 基于曲率阈值与区域生长的牙齿分割3.1 为什么不能直接做连通分量分析这个问题在课程答疑中经常出现。整口牙 STL 模型在拓扑上是一个连通分量因为牙龈把每颗牙齿的颈部连成一整片。直接调用vtkPolyDataConnectivityFilter只会得到一个区域根本分不开牙齿。必须先构造一个“牙齿区域”的标记图让牙龈沟处的高曲率网格成为天然屏障然后再对标记区域做连通性分析。3.2 顶点曲率估算与阈值选取VTK 给出了多种曲率估计方式这里使用平均曲率Mean Curvature。牙龈沟是表面凹陷最剧烈的地方平均曲率绝对值明显高于牙冠其它区域合理的阈值可以把牙冠上的较平坦区域保留下来。#include vtkCurvatures.h vtkNewvtkCurvatures curv; curv-SetInputData(mesh); curv-SetCurvatureTypeToMean(); curv-Update(); vtkPolyData* withCurv curv-GetOutput(); withCurv-GetPointData()-SetActiveScalars(Curvatures);计算完成后曲率值保存在点数据数组Curvatures中。阈值选取需要结合网格分辨率口腔扫描模型平均边长约 0.1~0.2 mm牙冠平滑区平均曲率通常在 0.3~0.8 之间牙龈沟平均曲率常超过 1.5。如果模型是从 CBCT 重建的体素化造成的噪声会让曲率整体偏大此时阈值需要上调 50% 以上。建议先在软件中绘制曲率直方图取双峰的谷底作为阈值。3.3 用阈值和连通成分提取得到每颗牙齿的候选区域曲率过滤的核心思想保留平均曲率低于阈值的三角片这些三角片大部分属于牙冠牙龈沟区域三角片被过滤掉之后牙齿之间在拓扑上不再相连。然后对剩余区域提取连通分量。#include vtkThreshold.h #include vtkGeometryFilter.h #include vtkPolyDataConnectivityFilter.h vtkNewvtkThreshold threshold; threshold-SetInputData(withCurv); // 只保留曲率低于 1.0 的区域 threshold-SetInputArrayToProcess( 0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, Curvatures); threshold-ThresholdByLower(1.0); vtkNewvtkGeometryFilter geom; geom-SetInputConnection(threshold-GetOutputPort()); vtkNewvtkPolyDataConnectivityFilter conn; conn-SetInputConnection(geom-GetOutputPort()); conn-SetExtractionModeToAllRegions(); conn-ColorRegionsOn(); conn-Update(); int regionCount conn-GetNumberOfExtractedRegions();ThresholdByLower(1.0)会保留曲率值小于等于 1.0 的单元格vtkThreshold输出非结构化网格因此需要用vtkGeometryFilter转回vtkPolyData。SetExtractionModeToAllRegions让每个连通区域都在输出的RegionId中分配自己的颜色通过该 ID 可以分离出每一颗牙的网格。但要注意如果相邻牙齿紧贴、牙龈沟狭窄阈值过滤后牙齿之间可能仍有小桥连接这时需要用形态学开运算或增大阈值。一个补救办法是在区域生长时提高阈值同时用区域面积过滤掉小于 500 mm² 的碎块因为牙冠面积通常大于该值。下面是区域生长替代方案的核心思路以每个区域质心为种子在原始网格上做 BFS 扩展相邻点的曲率必须低于阈值才会加入当前牙齿。这样即使存在细小桥梁生长过程也会被高曲率顶点阻断。不过 BFS 需要自己维护三角网格的邻接关系开销较大优先推荐阈值连通分量方案。4. 单颗牙齿的 OBB 轴向标定与牙龈去除4.1 用协方差矩阵构建设定长轴分割后每个连通区域可以看作一颗独立的牙齿候选。要计算每颗牙齿的轴向和归一化尺寸需要对区域内点云再次做 PCA最大特征值对应的特征向量就是牙齿长轴方向。注意前臼齿和磨牙的长轴方向并不严格平行于咬合面法向而是有一定倾斜因此不能用整口牙的全局方向替代。#include vector #include vtkPolyData.h #include vtkMath.h double ComputeLongAxis(vtkPolyData* tooth, double axis[3]) { double center[3]; vtkMath::ComputeCenter(tooth, center); double cov[3][3] {{0.0}}; int n tooth-GetNumberOfPoints(); for (vtkIdType i 0; i n; i) { double p[3]; tooth-GetPoint(i, p); double x p[0] - center[0]; double y p[1] - center[1]; double z p[2] - center[2]; cov[0][0] x * x; cov[0][1] x * y; cov[0][2] x * z; cov[1][1] y * y; cov[1][2] y * z; cov[2][2] z * z; } for (int i 0; i 3; i) for (int j 0; j 3; j) cov[i][j] / n; double eigenval[3], eigenvec[3][3]; vtkMath::Jacobi(cov, eigenval, eigenvec); for (int i 0; i 3; i) axis[i] eigenvec[0][i]; return eigenval[0]; }返回的最大特征值代表了牙齿沿长轴方向点分布的方差可以用它判断牙齿的“细长程度”。普通门牙的该值约为侧向方差的 5~8 倍磨牙则相对接近 2~3 倍。如果该比值小于 1.5说明分割出来的可能是一块牙龈而不是完整牙齿可以将其过滤掉。4.2 用 OBB 方向做轴向矫正和裁剪得到长轴后将每颗牙齿的局部坐标旋转到以长轴为 Z 轴的坐标系然后计算 OBB。OBB 中心取点云中心三个半边长分别从旋转后坐标的最大最小值得出。#include vtkMatrix4x4.h #include vtkTransform.h // axis 为长轴方向将长轴对齐到全局 Z 轴 vtkNewvtkMatrix4x4 mat; vtkMath::Perpendiculars(axis, mat-GetElement(0), mat-GetElement(1), 0); mat-SetElement(0, 3, 0); mat-SetElement(1, 3, 0); mat-SetElement(2, 3, 0); mat-SetElement(2, 0, 0); mat-SetElement(2, 1, 0); mat-SetElement(2, 2, 1); vtkNewvtkTransform transform; transform-SetMatrix(mat); transform-TransformPoint(center, center);vtkMath::Perpendiculars会根据给定长轴生成两个正交轴组成旋转矩阵的前两列。之后每个顶点投影到新坐标系记录 X/Y/Z 三个方向的最小最大值就得到了 OBB。这个 OBB 不仅用于渲染框体更重要的是后续裁剪牙龈。4.3 自动去除残余牙龈组织的比例裁剪曲率分割后部分牙龈仍然附着在牙冠颈部尤其在后牙的舌侧和颊侧。比较可靠的自动去除方法是在牙齿长轴方向OBB 局部 Z 轴上牙龈总是分布在靠近牙根一侧而牙冠分布在相反一侧。把每颗牙旋转到长轴竖直后统计所有顶点 Z 值的分布保留从最大 Z 往下 80%~85% 范围内的三角片底部 15%~20% 作为牙龈残余剔除。#include vtkClipPolyData.h double zRange[2]; tooth-GetBounds(); // 全局坐标系下取 z 范围 double lowerPlane zMin (zMax - zMin) * 0.15; vtkNewvtkPlane clipPlane; clipPlane-SetOrigin(0, 0, lowerPlane); clipPlane-SetNormal(0, 0, 1); vtkNewvtkClipPolyData clipper; clipper-SetInputData(tooth); clipper-SetClipFunction(clipPlane); clipper-Update();注意裁剪前必须先应用上一步的旋转让长轴对齐 Z 轴否则SetNormal(0,0,1)无法作用在正确的方向上。裁剪比例 15% 是经验默认值对于磨牙可以提高到 20%因为磨牙颈部牙龈包绕更明显。裁剪后还需要检查剩余三角片数量如果少于原始数量的 30%说明可能整颗牙齿被误删此时需要回退阈值。5. 牙齿计数、编号与缺失识别5.1 牙弓投影与空间排序分割和裁剪完成后每颗牙齿得到一组顶点和一个质心。要完成编号首先需要确定牙弓方向。把质心投影到第 2 章求得的咬合平面得到二维坐标。由于牙弓呈 U 形直接用 X 坐标排序在尖牙转弯处会产生错误应使用最近邻路径从最左端的一个质心出发每次选择距离当前点最近且未访问过的点形成顺序链。std::vectorPoint2D pts; // 从质心投影 std::vectorint order(pts.size(), -1); int current 0; // 可选根据 X 坐标最小值确定起点 for (int i 0; i pts.size(); i) { order[i] current; double minDist 1e10; int next -1; for (int j 0; j pts.size(); j) { if (order[j] -1 j ! current) { double d (pts[j] - pts[current]).norm(); if (d minDist) { minDist d; next j; } } } current next; }这段贪心路径适用于正常牙弓形态。如果患者有牙齿拥挤或扭转最近邻可能出现跨行跳变。因此建议在排序前先拟合一条二次曲线把点到曲线的投影距离作为排序依据。5.2 牙位编号映射牙位编号采用临床常用的 FDI 两位数字系统。实际项目中最直接的做法是从患者视角把下颌牙从患者左侧到右侧按 1 到 16 编号然后映射到 FDI 对应象限。但是开发时往往不知道哪侧是患者左因此先要区分上颌和下颌通过整口模型的重心和咬合面法向判断 Z 方向的正负。// 根据质心坐标与咬合面法向夹角判断上下颌 int quadrantBase (centerToPlaneSign 0) ? 10 : 20; // 再根据横坐标排序结果赋予 1~8 或 9~16 int fdiCode quadrantBase (sideIndex 8 ? 8 - sideIndex : sideIndex 1);这段是简化示意。真正稳妥的方案是制作一个 28 个标准牙位的模板模型用点云配准如 ICP把输入牙齿质心对齐到模板位置然后根据最近邻关系得到每个牙位编号。配准方案对缺失牙齿不敏感因为缺失位置的模板点没有对应输入点不会参与匹配。5.3 缺失识别数量与相邻间距离群有了编号顺序后缺失识别就有两个依据。第一是数量完整恒牙列不含智齿有 28 颗牙少于 28 直接说明有缺失。第二是位置即使数量恰好为 28也可能出现多生牙取代正常牙位的情况因此要通过相邻牙齿中心距离判断。std::vectordouble spacings; for (int i 0; i sortedCenters.size() - 1; i) { double d (sortedCenters[i1] - sortedCenters[i]).norm(); spacings.push_back(d); } double mean accumulate(spacings.rbegin(), spacings.rend(), 0.0) / spacings.size(); double var ...; // 计算方差 double threshold mean 2.0 * sqrt(var); for (int i 0; i spacings.size(); i) { if (spacings[i] threshold) { // 在 i 与 i1 之间标记缺失一颗牙 } }这里采用“均值2 倍标准差”作为离群阈值。正常牙弓中相邻牙齿中心距大约在 5~8 mm标准差约 0.8 mm因此阈值大约为 8~9 mm。若间距超过该值就在这两个牙位之间插入一个缺失标记。处理阶段默认参数说明曲率阈值1.0低于 1.0 认定为平滑牙冠区面积过滤500 mm²过滤碎小连通体OBB 裁剪比例15%~20%去除长轴末端牙龈残留缺失间距阈值均值 2σ大于该间距则判定缺牙6. Qt 与 VTK 渲染分割结果的几个关键技巧6.1 在 QVTKOpenGLWidget 中搭建渲染管线Qt 5.13.2 搭配 VTK 8.2.0 时建议使用QVTKOpenGLWidget而非旧的QVTKWidget。前者支持 OpenGL 3.2 上下文渲染性能更好。为了让 VTK 事件循环与 Qt 事件循环兼容必须在main函数中先设置默认 OpenGL 格式#include QApplication #include QSurfaceFormat #include QVTKOpenGLWidget.h int main(int argc, char* argv[]) { QSurfaceFormat format; format.setRenderableType(QSurfaceFormat::OpenGL); format.setVersion(3, 2); QSurfaceFormat::setDefaultFormat(format); QApplication app(argc, argv); // 创建主窗口将 QVTKOpenGLWidget 放入 centralWidget return app.exec(); }6.2 并行渲染每颗牙齿并分配颜色每颗牙齿输出一个vtkActor通过vtkLookupTable按编号分配颜色。实测发现当牙齿数量超过 20 颗时逐个创建 Actor 的渲染帧率仍可达到 60 FPS。需要注意设置vtkRenderer的ResetCameraClippingRange否则裁剪面可能吃掉部分牙齿。vtkNewvtkLookupTable lut; lut-SetNumberOfTableValues(32); lut-Build(); for (int i 0; i toothActors.size(); i) { double rgb[3]; lut-GetColor(i % 28 1, rgb); toothActors[i]-GetProperty()-SetColor(rgb); renderer-AddActor(toothActors[i]); } renderer-GetRenderWindow()-Render();6.3 拾取牙齿并联动显示编号在界面中点击某一颗牙齿时使用vtkPointPicker返回被选中网格所在的 Actor再通过 Actor 与牙齿编号的映射关系显示信息。这里有一处容易踩坑如果多个 Actor 数据共享一个 Mapper拾取后无法区分是哪个 Actor。务必让每颗牙齿独立创建一个 Mapper。vtkNewvtkPointPicker picker; picker-Pick(x, y, 0, renderer); vtkActor* pickedActor picker-GetActor(); if (pickedActor) { int toothId actorToId[pickedActor]; statusLabel-setText(QString(牙位编号: %1).arg(toothId)); }验证分割正确率时可以加一个键盘快捷键切换整个模型的半透明显示观察每颗牙齿的边界是否贴合牙龈沟。如果相邻牙齿颜色出现穿越说明曲率阈值偏低如果牙齿边缘过度收缩说明 OBB 裁剪比例过大需要调小 5%。本文还有配套的精品资源点击获取