基于三角网格演化的CBCT牙齿图像分割方法

【作 者】刘枭雄，石峰，张继武

1 上海交通大学生物医学工程系，上海，200240

2 Carestream Health 全球研发中心，上海，201206

在现代牙科领域，锥形束计算机断层扫描(CBCT，Cone Beam Computed Tomography)越来越广泛采用，使三维测量牙齿的形态及大小成为可能。应用三维图像分割方法，获得患者三维个性化牙冠及牙根的解剖信息，可以帮助临床医生作出更准确的判断，确定个性化的治疗方案。将牙齿的几何信息与计算机辅助制造相结合时，可以提高牙齿修补材料制造的精度和效率，因此从CBCT中准确地分割出牙齿是非常关键的一步。然而，由于CBCT牙齿图像存在以下特点，使牙齿分割成为一项极具挑战的工作：(1)相邻牙齿间隙小，造成图像中相邻牙齿接触处轮廓线丢失；(2)牙齿的密度从牙冠到牙根各不相同，引起CBCT图像中单个牙齿的灰度不均匀；(3)由于牙齿根部埋在牙槽骨里，二者密度相近，故边缘不清晰；(4)牙齿的轮廓在牙根和牙冠部分存在拓扑结构的变化。

迄今为止，主要的牙齿分割算法有：(1)手动分割法，例如Rahimi等人[2]由牙齿的二维组织图像和显微CT图像构建出三维模型，其缺点是需要多人参与，分割时间长；(2)自适应阈值分割法，由于需要解决牙齿灰度值不均匀以及相邻牙齿边界模糊的问题，该方法常采用大阈值，因此导致出现过分割现象[3]；(3) 参数曲线演化分割法，该方法的缺点是无法解决牙齿拓扑结构改变以及相邻牙齿间边界模糊的问题[4]；(4) 水平集分割法，该方法对初始化非常敏感，并且缺少对演化表面的局部控制，导致在牙齿的不同部分进行自适应演化很难实现。

本文提示了一种利用三角网格演化分割CBCT牙齿图像的新方法。Smith等人将这种方法用在脑组织的MRI图像分割上并取得了较好的效果。本文在借鉴文献[6]的基础上提出了适合分割CBCT牙齿图像的方法。该方法通过控制三角网格模型上每个顶点的运动，使模型表面不断地迭代变化，最终达到演化后的模型能够与牙齿表面相吻合。在这个演化过程中，主要受到两个方面条件的制约：一方面，模型表面需要保持良好的平滑度，并且能够与牙齿表面的平滑度相匹配；另一方面，模型能够准确地演化到牙齿表面。该方法的优点是实现过程简单，人工参与少，并且鲁棒性好。

1 算法与实现

可变形表面方法是可变形模型方法在三维上的应用。其基本思想是将可变模型表示为参数曲线或曲面V(s)，模型的形状由满足以下能量泛函的极小化条件所决定：

其中Eint为模型的内能，代表了对模型形状的约束，使模型保持一定的光滑性和连续性；Eext为模型的外能，来自目标图形特征，引导可变模型轮廓演化到目标特征，当模型达到图像特征处时该能量达到最小。因此，利用可变形表面方法分割图像的一个关键问题，是如何定义控制模型演化的内力和外力。Smith等人[6]提出利用三角网格建立初始模型，通过与模型节点相关的一组向量控制模型演化，并利用模型不断迭代得到最终结果。

在定义内力和外力前，需先引入两个基本向量。首先，定义三角网格初始模型的每个顶点vi的法向量7]：

其中：N(vi)是顶点vi的相邻三角形集合；#N(vi)为集合N(vi)中的元素个数；是第k个三角形的法向量。

其次，得到每个顶点临近所有顶点的平均位置，然后将由当前顶点指向平均位置的向量，分解为相互垂直的两个向量和，分别平行和垂直于向量（如图1），定义：

图1 二维情况下将向量 分解为 和Fig.1 For the 2D case decomposing the vector into and

1.1 内力的定义

为保持模型表面三角网格大小均匀，且具有较好的平滑度和连续性，定义内力由两部分组成：(1)向量=0.5，作用是保持模型表面三角网格均匀分布；(2)向量=f，使当前节点向与邻近节点共面的2方向移动，达到增加模型表面平滑度的目的。为使模型在曲率较大的位置具有较高的内力，在当前顶点处定义：

图2 局部曲率半径r，距离l和Sn向量的关系。Fig.2 The relation between local curvature r, vertex spacing l and vector

由于牙冠和牙根部位曲率极值差异较大，因此需分别估计这两部分的参数rmin和rmax。为简化算法，以牙冠和牙根在z轴方向的位置区间为依据，确定参数rmin和rmax的取值。这样不仅有利于提高模型演化结果的准确性，而且保持了算法迭代过程的稳定性。

1.2 外力的定义

通常，外力由图像灰度的梯度向量决定，将形变模型拉向灰度梯度值高的位置，即目标物体表面。但这样定义存在缺陷，即模型难以演化到目标物体的弱边界，且易受图像噪声影响。为此，引入局部阈值tl，通过沿当前顶点其局部法向量方向所搜索到的图像灰度值信息确定tl值。形变模型演化的外力定义为：

f3定义为：

其中，tl为局部阈值，通过沿当前顶点其局部法向量方向所搜索到的图像灰度值信息确定，其定义为：

Imin和Imax为从当前顶点开始，沿法向量向模型内部延伸距离d所搜索到的最小灰度值和最大灰度值，分别定义为:

为限制图像中过“亮”或过“暗”的体素点所产生的影响，使用灰度值t2，t和tm加以限定。其中t2为位于灰度累计直方图中2%处的灰度值；t为通过Otsu算法[8]得到能够大致区分目标物体和背景的阈值；tm为在初始模型内部目标图像灰度值的中位数。

参数tm是控制模型演变的重要参数。但由于牙根部位与邻近组织的灰度值偏差较小，导致牙根部位的图像灰度值的平均值高于图像中的其他部分。所以，为得到有效的外力，需要单独计算该范围的灰度值中位数tm。同样，为简化算法，以牙根在z轴方向的位置区间为依据，确定参数tm的取值。

1.3 模型初始化

根据牙齿的形状特点，我们采用椭球形的初始模型，其表面由三角网格构成。首先，由一个20面体通过迭代将每个三角形细分为四个更小三角形，同时调整新产生的顶点到中心的距离，使所有顶点均处于该二十面体对应的球面上。然后，根据图像信息大致确定初始模型在x，y，z三个方向的半径、模型的中心点和旋转角度，从而对该球形模型进行空间变化，建立椭球形模型。参数设定的原则是使初始模型尽量位于牙齿的内部，如图1所示。

图3 初始椭球模型。Fig. 3 The initial ellipsoid model

2 实验结果与分析

为了验证本文算法的有效性，测试在CPU为Intel Core2 Duo T9600 2.8 GHz，内存为4GB的PC机上运行，算法开发采用Visual C++ 2008实现。算法验证所使用的图像数据来源于Carestream Kodak K9000型CBCT。操作者首先在每颗牙齿的近似中心点上放置初始模型，然后模型的演化自动开始直至收敛。分割结果如图4所示。

图5演示了算法在牙齿根部的分割有效性。图5（a）和5（c）为牙根部CBCT图像的横截面，可以看出牙齿的密度和槽骨的密度已经非常接近，边缘也非常模糊。图5(b)、5(d)分别为5(a)、5(c)的分割结果，可以看出牙齿的边缘被完整地提取出来了。

图4 牙齿分割结果（横截面）Fig. 4 Teeth segmentation result (Axial view)

图5 牙齿分割结果（横截面）Fig. 5 Teeth segmentation result (Axial view)

图6(a)为一颗牙齿的矢状面截图，在牙冠部位两颗邻近的牙齿几乎相连，并且牙根和牙槽骨的密度十分接近。图6(b)为本算法的分割结果，在牙冠处能正确地区分两颗牙齿，在牙根处牙齿与牙槽骨也得到了很好地分离。

图6 牙齿分割结（矢状面）Fig. 6 Teeth segmentation result (Sagittal view)

为进一步证明本方法的有效性和准确性，将本文算法的分割结果和专家的分割结果进行了量化对比。真阳性区域TP表示为手工分割和算法分割的公共区域；假阳性区域FP表示在算法分割区域内但在手工分割区域范围外；假阴性区域FN表示包含在手工分割区域内，但在算法分割区域外。由以上三个区域定义的敏感度Cs和精确度Ca为[9]：

其中Am= TR+FN，表示手工分割所包含的区域。

本文使用了7组数据来进行分割对比，以验证算法的精度和效率。结果的平均敏感度Cs为89.3%， 平均精确度Ca为86.5%。本文算法分割单个牙齿的平均时间为37 s。

3 结论

本文提出了一种利用三角网格模型演化进行牙齿CBCT图像分割的方法。方法首先由操作者在牙齿中心放置一个椭球形初始模型，模型在自身内力约束和图像信息外力约束下不断演化到牙齿边缘，直至收敛。其优点是模型在演化的过程始终保持光滑和连续，避免了目标轮廓上的噪声和不规则，同时也有效地保持了牙齿结构拓扑不变，防止分割泄露到邻近的牙齿。与专家手工分割结果的对比，说明了方法的准确性和有效性，能够为医生的临床诊断和治疗提供辅助手段。

[1] Hui Gao, Oksam Chae. Individual tooth segmentation from CT images using level set method with shape and intensity prior [J].Pattem Recognition, 2010, 43(7):2406-2417.

[2] A.Rahimi,L.Keilig, et al. 3D Reconstruction of dental specimens from 2D histological images and μCT-Scans[J]. Computer Methods in Biomechanics and Biomedical Engineering, 2005,8(3):167-176.

[3] H.Heo, O.S.Chae, et al. Segmentation of tooth in CT images for the 3D reconstruction of teeth[C], Proceedings of SPIE-IS&T Electronic Imaging, 2004, 455-466.

[4] X.L.Wu, H.Gao, H.Heo, et al. Improved B-spline contour fitting using genetic algorithm for the segmentation of dental computerized tomography image sequences[J], Imaging Science and Technology, 2007, 51(4): 328-336.

[5] Sh. Keyhaninejad, R.A. Zoroo fi, et al. Automated segmentation of teeth in multi-slice CT images[C], IET International Conference on Visual Information Engineering, 2006, 339-344.

[6] Stephen M.Smith. Fast Robust Automated Brain Extraction[J],Human Brain Mapping, 2002, 17(3): 143-155.

[7] 罗良峰. 离散三角网格上的法向量和曲率估计[D], 大连理工大学, 2007

[8] Nobuyuki Otsu. A threshold selection method from gray-level histogram [J], IEEE Transactions on System, Man and Cybemetics,1979, 9(1): 62-66.

[9] 杨振森, 李传富, 周康源, 等. 基于改进测地线模型的前列腺超声图像分割 [J], 中国医疗器械杂志, 2008: 32(5): 316-318.