针对三维计算机断层扫描(CT)体数据的牙齿分割问题, 本文提出了一种基于区域自适应形变模型的CT图像牙齿结构测量方法。本文方法结合了自动阈值分割、CV活动轮廓模型和图割方法, 利用自动阈值分割实现牙冠的分割与定位, 然后利用牙冠分割结果作为初始轮廓逐层分割牙齿。在分割难度最大的牙根上采用CV活动轮廓和图割互补的方法实现了牙根的准确分割。实验结果表明本文提出的牙齿结构测量方法能够准确地自动分割出牙齿的牙冠部分, 进而在牙冠分割基础上快速准确地分割出牙颈和牙根。本文提出的牙齿结构测量方法能够准确地从临床CT牙齿数据中分割提取牙齿结构, 鲁棒性强、精度高, 可以有效辅助医生的临床治疗。
引用本文： 王立新, 刘新新, 刘希云, 杨志海, 杨健, 艾丹妮, 王涌天. 基于区域自适应形变模型的CT图像牙齿结构测量方法研究. 生物医学工程学杂志, 2016, 33(2): 308-314. doi: 10.7507/1001-5515.20160052 复制
引言
计算机断层成像(computed tomography, CT)技术为牙科医生提供了丰富的牙齿图像信息，有效地拓展了医生的临床诊断视野，而三维(three dimension, 3D)图像处理技术可以有效地利用图像数据，辅助牙科医生高效准确地治疗口腔疾病[1]。在临床牙齿种植手术中，植入点位置和方向的选择是影响牙齿种植手术成功的关键因素之一，牙齿轮廓的正确分割是牙种植手术规划和导板设计的基础，牙齿分割结果的精准度会直接影响后续的应用[2]。
研究表明基于3D CT数据可以分割创建准确的牙齿三维模型，通过牙齿模型测量得到的牙齿数据是非常可靠的[3]。牙齿分割技术将感兴趣的牙齿轮廓提取出来，基于分割结果可以对牙齿进行特征提取和参数测量，进而为医生提供牙冠及牙根的完整解剖信息，从而制定高效准确的治疗方案[4]。因此，牙齿结构精准分割技术的研究具有非常重要的临床意义。然而，基于CT图像的牙齿分割目前仍然是一项非常具有挑战性的工作，原因如下：牙齿通常紧密相连(尤其是牙冠)，这造成了在CT图像中相邻牙齿轮廓合并在一起；牙齿的不同结构具有不同的牙齿密度，CT图像表现为牙齿灰度分布不均匀；牙根被颚骨完全包围，中间只有薄薄一层软组织，造成牙齿边界模糊且不连续[5]。
迄今为止，主要的牙齿分割算法有四种：①人工分割法，比如Rahimi等[6]利用CT图像人工分割构建牙齿三维模型。但是这种方法的缺点是主要依靠人工干预，效率低下，耗时长；②自动阈值分割法，由于牙齿灰度不均匀并且相邻牙齿边界模糊，难于分割，该方法常采用较大阈值进行分割，过分割现象严重；③活动轮廓模型，该方法在边界模糊或不连续时很容易被错误的边界吸引，造成分割结果的错误[7]；④水平集分割法，该方法对初始化轮廓非常敏感，针对不同CT图像或同一数据的不同层图像很难保持良好的准确率。
本文提出了一种基于区域自适应形变模型的CT图像牙齿结构测量方法，该方法结合了自动阈值分割、Chan-Vess活动轮廓分割方法(CV-AC)和图割方法，能够实现从CT体数据中自动分割出牙齿的完整轮廓。该方法首先利用自动阈值分割方法得到覆盖牙齿外层的牙釉质，牙釉质主要覆盖牙冠，具有较高的CT值，通过分析牙釉质的分割结果可以对牙齿进行初始定位；然后将自动阈值分割结果作为活动轮廓分割方法的初始轮廓，通过CV活动轮廓模型可以在牙釉质的基础上完整地分割出牙冠；最后按照从牙冠到牙根的顺序，对每个牙齿横断面分别进行基于CV-AC的牙齿分割，每层横断面的分割过程都用上一层分割结果作为CV-AC的初始轮廓。当牙齿与骨骼连接紧密，灰度没有明显区别，造成CV-AC模型发生错误时，本文利用图割分割方法替换CV-AC寻找一条最优的牙齿轮廓，最终实现了牙齿轮廓的准确提取，算法流程如图 1所示。
图1
牙齿分割算法流程
Figure1.
Flow chart of the tooth segmentation algorithm
1 算法与实现
1.1 图像预处理
本文测试数据来源于首都医科大学附属北京世纪坛医院口腔科，牙齿数据本身的选择没有经过任何筛选和处理，按照姓名首字母排序的第一个患者数据。该患者数据采集时间为2014年9月19日，图像大小为512×512×512，空间分辨率为0.16 mm×0.16 mm×0.16 mm。
牙齿结构分为牙冠、牙颈和牙根三部分，其中牙根与颌骨相连，在它们之间有一层薄薄的牙周膜。在临床CT图像中，颌骨和牙齿都具有较高的灰度值，而牙周膜属于软组织，CT值远远低于骨骼和牙齿。由于牙周膜厚度很小，在CT图像中牙周膜非常模糊，牙齿和颌骨的分界线清晰度低且不连贯，增加了分割的难度。为了提高牙齿分割的准确度，本文通过直方图拉伸的方法处理CT图像[8]，进而增强牙齿、颌骨和牙周膜的对比度。本文采用的是线性函数实现图像对比度的拉伸，函数如下：
| $ y=\left\{ \begin{array}{l} x + {b_1}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} x \le {n_1}\ a \times x + {b_2}{\kern 1pt} {\kern 1pt} {\kern 1pt} {n_1}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} < x < {n_2}\ x + {b_3}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} x \ge {n_2} \end{array} \right., $ |
其中x是原始CT灰度值，y是经过增强后的灰度值，a、b1、b2和b3为增强参数，n1和n2是骨骼灰度值分布范围的高低双阈值。
1.2 牙齿轮廓初始化
牙釉质是牙冠外层的白色半透明的钙化程度最高的坚硬组织，在CT图像中其CT值远高于颌骨[9]。利用牙釉质高CT值的特性，我们可以使用阈值分割从CT图像中提取牙釉质图像。由于牙釉质主要位于牙齿的牙冠部分，所以牙釉质的分割结果可以帮助我们确定牙冠的位置和方向。本文中牙釉质分割采用了基于灰度直方图的自动阈值分割方法，牙釉质高CT值的特性在灰度直方图中表现为在灰度大约1 800的位置有一个明显的波峰，位于图 2(a)矩形区域，图 2(b)是矩形区域的放大图，图 2(b)中波峰就是牙釉质峰，波峰位置即为牙釉质分割的最佳阈值。
图2
牙齿CT图像的直方图和牙釉质体积随横断面的变化
(a)灰度直方图；(b)局部放大直方图；(c)牙釉质体积随横断面的变化
Figure2. Histogram of the tooth CT image and volume difference of the tooth enamel between different slices(a) histogram; (b) local histogram; (c) volume difference of the tooth enamel between different slice
通过自动阈值分割我们得到牙釉质在牙齿上的分布，我们统计了牙釉质在所有横断面上的体积，统计结果如图 2所示。在图 2(c)中，横轴表示CT数据的512个横断面，纵轴表示每个横断面牙釉质的体积。在图中有两个波峰，左侧的波峰DT表示下颌骨牙齿的牙冠上的牙釉质，右侧的波峰UT表示上颌骨牙齿的牙冠上的牙釉质，波峰之间的波谷表示上排牙齿和下排牙齿的分界线。牙冠的分割是为后续分割提供精准的定位和准确的初始轮廓，所以我们选择范围在DT-UT之间进行牙冠的分割。
1.3 牙冠和牙颈分割
在牙冠部分，牙齿暴露在空气中，不与任何软组织和骨组织相连，牙齿边界非常清晰。我们利用牙釉质分割结果作为初始轮廓，利用CV-AC方法对牙齿的牙冠部分进行分割。Chan-Vese模型是Chan和Vese于2001年提出的基于区域的水平集活动轮廓模型，具有对噪声、杂波和对初始轮廓的设定不敏感等优点[10]。牙釉质(只覆盖在牙冠表面薄薄一层)作为初始轮廓，CV-AC模型通过较高次数的迭代演化后可以完美地实现牙冠的整体分割。
由于暴露在空气当中，牙冠相比牙颈和牙根是最容易分割的部分。在牙冠分割完成后，牙颈以下的部分因其直径减小，通常情况下是相互分离的，我们利用连通域运算获得独立的牙齿，然后分别完成对每个牙齿的整体分割。如果有相连的牙齿，我们利用CV-AC对牙齿继续向牙根分割，直到每个牙齿都相互独立。本文的CT数据在横轴方向上有非常高的分辨率(0.16 mm)，相邻横断面之间的牙齿轮廓不会有非常大的改变，所以我们可以把上一横断面的分割结果作为下一横断面分割的先验轮廓。这样每次分割的初始轮廓都非常接近于待分割的真实牙齿轮廓，我们只需要很少的迭代次数，活动轮廓就能找到牙齿的边界，进而可以极大地提高分割的精准度。但是除牙冠外，牙齿其他部位处于软组织和骨组织的包围中，尤其牙根部分与骨骼组织的边界非常狭窄且连贯性较差。为了更好地对牙齿进行分割，我们在CV-AC模型的基础上进行针对性改进，为CV-AC模型加入了灰度限制项和曲线演进先验限制项(PKL)，对于任意曲线变量C定义如下的“拟合能量函数”：
| $ \begin{array}{l} {F_1}\left(C \right) + {F_2}\left(C \right)=\int_{inside\left(C \right)} {{{\left| {I\left(x \right)-{c_1}} \right|}^2} + \omega } \ {\left| {I\left(x \right)-Threshold-I\left(x \right)} \right|^2}dx + \int_{outside\left(C \right)} {{{\left| {I\left(x \right)-{c_2}} \right|}^2} + } \ \omega \times {\left| {Threshold-I\left(x \right)} \right|^2}dx + PKL \end{array} $ |
式(2)中C是任意曲线变量，c1和c2是依赖于C的常数，它们分别是曲线C内部、外部区域的灰度平均值，Threshold是CT图像中骨骼和非骨骼组织的全局最佳骨骼阈值，该值也可以通过CT图像直方图获得，ω是骨骼灰度限制项的权重系数。其中，PKL是基于上一横断面先验知识的限制项，该项由上一横断面的牙齿轮廓经过距离变换生成。限制力大小与曲线移动距离成正比，当曲线移动距离超出一定范围时限制力变为无穷大，阻止曲线过大地偏离先验曲线。PKL的加入可以限制曲线的移动范围，在有限的范围内寻找最优的牙齿边界线。
1.4 牙根分割
由于牙根完全被骨骼包围，牙根相比牙冠和牙颈更加难以分割。基于1.3中图像增强和CV模型的分割方法在大部分情况下可以实现牙齿的准确分割，但是在一些情况下(犬齿和门牙)牙根会和颌骨紧密相连，边界处的灰度差别会非常小。当出现这种情况时，CV-AC模型在寻找牙齿边界过程中很容易在某一段边界发生错误。为了避免错误在下一横断面累积，我们采用图割的方法接替CV-AC继续分割牙齿。判断是否发生错误的方法是当活动轮廓分割错误时，活动曲线会在某一段曲线严重偏离先验曲线，重叠率(Dice系数)和质心位置会产生较大变化。图割是一种基于图论的组合优化方法，它将一幅图像映射成一个网络图，并建立关于标号的能量函数，运用最大流/最小割算法对网络图进行切割，得到网络图的最小割，即目标函数的最小值[11]。图割方法需要设定前景和背景种子点，种子点由上一横断面的牙齿分割结果作为先验知识生成，如图 3(a)所示；在牙齿边界内外分别选择两个像素加上边界线组成宽度为5个像素的窄带(灰色部分)，窄带内部区域(白色)设定为图割前景像素，窄带外部区域(黑色)设定为图割背景像素，灰色部分为待分割区域的划分结果，如图 3(b)所示。由此，可利用图割算法在初始轮廓区域内寻找一条最优的牙齿边界线，该方法在保证分割准确性的同时可避免程序陷入局部最优而无法找到真实边界的情况。
图3
待分割区域划分示意图
(a)牙齿先验轮廓；(b)待分割区域划分结果
Figure3. Diagram of the segmentation region determination(a) prior outline of the tooth; (b) result of the segmentation region
在牙冠分割完成后，图割方法就可以准确分割出大部分牙齿结构。但是，对于牙根分裂的情况，图割方法不能快速准确地识别出牙根的轮廓，这里主要是由于种子点的选取是基于上一横断面设定的，而种子点设定对图割的分割结果又有直接的影响。然而，CV-AC模型可以很灵活地进行闭合曲线的分裂和合并，能够准确处理牙根分裂的部位。所以，本文中只有当CV-AC模型出现错误才会使用图割方法对牙齿继续进行分割。
2 实验和结果
为了验证本文算法的有效性，本文在临床牙齿CT数据上进行了测试，测试在CPU为Intel Core I7-2600 3.4 GHz，内存为16 GB的PC机上运行，算法在Matlab平台上实现。
首先，我们对CT数据进行自动阈值分割得到了牙釉质在牙齿上的分布，如图 4所示。在图 4(a)中，牙釉质位于牙冠的外层形成一个环形，并且牙釉质主要位于牙齿的顶端。然后，我们把牙釉质阈值分割结果作为初始轮廓，利用CV-AC方法实现牙冠的完整分割。CV-AC方法有对初始轮廓不敏感和对噪声不敏感等优点，并且牙冠周围没有被覆其他身体组织，以牙釉质分割结果作为初始轮廓，经过较大次数的迭代(本文采用迭代600次)就可以获得牙冠完美的分割结果。图 4(b)中红色部分是牙釉质分割结果，白色部分是CV-AC经过600次迭代后得到的牙冠分割结果，图 4(c)是对DT和UT之间的所有横断面进行分割后牙冠分割结果。而如图 4(d)和图 4(e)所示则展示了在实际CT图像中CV-AC方法利用牙釉质分割结果作为初始轮廓完整分割出牙冠的过程。
图4
牙冠分割结果
(a)牙釉质阈值分割；(b)CV-AC分割；(c)牙冠分割三维渲染结果；(d)初始轮廓；(e)最终分割结果
Figure4. Results of the tooth crowns segmentation(a) result of the enamel threshold segmentation; (b) result of the CV-AC segmentation; (c) 3D rendering of the dental crown segmentation result; (d) initial contour; (e) final segmentation result
在完成牙冠分割后我们开始对牙颈和牙根进行分割，由于相邻横断面之间牙齿轮廓线只有很小的差别，这里我们利用上一横断面分割结果作为本横断面的先验轮廓。这时通过CV-AC方法对牙齿进行分割我们只需要很少的迭代次数(本文采用25次)就可以实现牙齿的精准分割。
相比牙颈被软组织(牙龈)包围，牙根完全被骨骼包围，牙齿和骨骼之间只有薄薄的一层牙周膜，并且在CT图像中非常模糊，所以分割难度很高。如果牙齿和颚骨的边界模糊不清，活动轮廓就会在骨骼和牙齿连接处被骨骼吸引，产生错误的结果。当CV-AC方法发生错误时，我们采用图割方法来寻找牙齿最优边界，如图 5所示。在图 5中，图 5(a)是牙齿原始CT图像，图 5(b)是作为先验轮廓的上一横断面牙齿分割结果，图 5(c)是基于先验轮廓设定图割的前景和背景(白色表示前景区域，黑色表示背景区域，灰色表示待分割区域)，图 5(d)为图割方法的分割结果。图 5显示无论牙齿边界比较清晰还是非常模糊，图割方法都能较好地分割出牙齿的轮廓。
图5
牙齿分割结果展示(3例)
(a)原始图像；(b)初始轮廓；(c)种子区域；(d)分割结果
Figure5. Tooth segmentation results based on the proposed method (three examples)(a) original image; (b) initial contour; (c) seed point region; (d) segmentation result
运用自动阈值分割、CV-AC模型和图割方法后，我们得到了牙齿的整体分割结果，如图 6所示。从牙齿分割结果的3D渲染图中我们可以看出本文提出的牙齿分割方法可以很好地分割出所有牙齿结构，对于分裂的牙根也保持了良好的准确度。图 6(b)是犬齿的单独3D渲染图，犬齿植入骨骼较深，分割难度相对较大，而采用本文提出的牙齿分割方法有效地分割出了整个牙齿结构。
图6
牙齿分割结果3D渲染图
(a)下排牙齿3D渲染图；(b)犬齿3D渲染图
Figure6. 3D rendering of the tooth segmentation results(a) 3D rendering of the segmented lower teeth; (b) 3D rendering of the canine tooth
3 结论
随着计算机图像处理技术的快速发展，医学图像处理技术在现代医学、生物医学研究和临床诊断等领域得到广泛的应用，并取得了非常好的效果。本文提出了一种基于区域自适应形变模型的CT图像牙齿结构测量方法。该方法首先利用自动阈值分割方法对牙齿进行分割，获得了牙釉质在牙冠上的分布情况，进而统计牙釉质在不同横断面上的分布，获得牙冠的位置和牙齿的方向。然后，我们利用牙釉质分割结果作为初始轮廓，使用CV-AC模型分割牙齿得到完整的牙冠轮廓。牙冠的分割结果作为先验轮廓用于后续牙齿轮廓的分割，这样可极大地减少迭代的次数，提高分割方法的运行效率。当活动轮廓发生错误时，可使用图割替代CV-AC方法进行牙齿分割。实验证明通过结合CV-AC和图割方法可以很好地实现牙齿整体的准确分割，解决了牙齿拓扑结构改变以及相邻牙齿间边界模糊的问题。对于临床数据的分割结果验证了方法的准确性和有效性，能够为医生的临床诊断和治疗提供辅助手段，在牙齿种植、牙齿正畸、医学信息学等多个方面有着重要的临床应用前景。
引言
计算机断层成像(computed tomography, CT)技术为牙科医生提供了丰富的牙齿图像信息，有效地拓展了医生的临床诊断视野，而三维(three dimension, 3D)图像处理技术可以有效地利用图像数据，辅助牙科医生高效准确地治疗口腔疾病[1]。在临床牙齿种植手术中，植入点位置和方向的选择是影响牙齿种植手术成功的关键因素之一，牙齿轮廓的正确分割是牙种植手术规划和导板设计的基础，牙齿分割结果的精准度会直接影响后续的应用[2]。
研究表明基于3D CT数据可以分割创建准确的牙齿三维模型，通过牙齿模型测量得到的牙齿数据是非常可靠的[3]。牙齿分割技术将感兴趣的牙齿轮廓提取出来，基于分割结果可以对牙齿进行特征提取和参数测量，进而为医生提供牙冠及牙根的完整解剖信息，从而制定高效准确的治疗方案[4]。因此，牙齿结构精准分割技术的研究具有非常重要的临床意义。然而，基于CT图像的牙齿分割目前仍然是一项非常具有挑战性的工作，原因如下：牙齿通常紧密相连(尤其是牙冠)，这造成了在CT图像中相邻牙齿轮廓合并在一起；牙齿的不同结构具有不同的牙齿密度，CT图像表现为牙齿灰度分布不均匀；牙根被颚骨完全包围，中间只有薄薄一层软组织，造成牙齿边界模糊且不连续[5]。
迄今为止，主要的牙齿分割算法有四种：①人工分割法，比如Rahimi等[6]利用CT图像人工分割构建牙齿三维模型。但是这种方法的缺点是主要依靠人工干预，效率低下，耗时长；②自动阈值分割法，由于牙齿灰度不均匀并且相邻牙齿边界模糊，难于分割，该方法常采用较大阈值进行分割，过分割现象严重；③活动轮廓模型，该方法在边界模糊或不连续时很容易被错误的边界吸引，造成分割结果的错误[7]；④水平集分割法，该方法对初始化轮廓非常敏感，针对不同CT图像或同一数据的不同层图像很难保持良好的准确率。
本文提出了一种基于区域自适应形变模型的CT图像牙齿结构测量方法，该方法结合了自动阈值分割、Chan-Vess活动轮廓分割方法(CV-AC)和图割方法，能够实现从CT体数据中自动分割出牙齿的完整轮廓。该方法首先利用自动阈值分割方法得到覆盖牙齿外层的牙釉质，牙釉质主要覆盖牙冠，具有较高的CT值，通过分析牙釉质的分割结果可以对牙齿进行初始定位；然后将自动阈值分割结果作为活动轮廓分割方法的初始轮廓，通过CV活动轮廓模型可以在牙釉质的基础上完整地分割出牙冠；最后按照从牙冠到牙根的顺序，对每个牙齿横断面分别进行基于CV-AC的牙齿分割，每层横断面的分割过程都用上一层分割结果作为CV-AC的初始轮廓。当牙齿与骨骼连接紧密，灰度没有明显区别，造成CV-AC模型发生错误时，本文利用图割分割方法替换CV-AC寻找一条最优的牙齿轮廓，最终实现了牙齿轮廓的准确提取，算法流程如图 1所示。
图1
牙齿分割算法流程
Figure1.
Flow chart of the tooth segmentation algorithm
1 算法与实现
1.1 图像预处理
本文测试数据来源于首都医科大学附属北京世纪坛医院口腔科，牙齿数据本身的选择没有经过任何筛选和处理，按照姓名首字母排序的第一个患者数据。该患者数据采集时间为2014年9月19日，图像大小为512×512×512，空间分辨率为0.16 mm×0.16 mm×0.16 mm。
牙齿结构分为牙冠、牙颈和牙根三部分，其中牙根与颌骨相连，在它们之间有一层薄薄的牙周膜。在临床CT图像中，颌骨和牙齿都具有较高的灰度值，而牙周膜属于软组织，CT值远远低于骨骼和牙齿。由于牙周膜厚度很小，在CT图像中牙周膜非常模糊，牙齿和颌骨的分界线清晰度低且不连贯，增加了分割的难度。为了提高牙齿分割的准确度，本文通过直方图拉伸的方法处理CT图像[8]，进而增强牙齿、颌骨和牙周膜的对比度。本文采用的是线性函数实现图像对比度的拉伸，函数如下：
| $ y=\left\{ \begin{array}{l} x + {b_1}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} x \le {n_1}\ a \times x + {b_2}{\kern 1pt} {\kern 1pt} {\kern 1pt} {n_1}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} < x < {n_2}\ x + {b_3}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} x \ge {n_2} \end{array} \right., $ |
其中x是原始CT灰度值，y是经过增强后的灰度值，a、b1、b2和b3为增强参数，n1和n2是骨骼灰度值分布范围的高低双阈值。
1.2 牙齿轮廓初始化
牙釉质是牙冠外层的白色半透明的钙化程度最高的坚硬组织，在CT图像中其CT值远高于颌骨[9]。利用牙釉质高CT值的特性，我们可以使用阈值分割从CT图像中提取牙釉质图像。由于牙釉质主要位于牙齿的牙冠部分，所以牙釉质的分割结果可以帮助我们确定牙冠的位置和方向。本文中牙釉质分割采用了基于灰度直方图的自动阈值分割方法，牙釉质高CT值的特性在灰度直方图中表现为在灰度大约1 800的位置有一个明显的波峰，位于图 2(a)矩形区域，图 2(b)是矩形区域的放大图，图 2(b)中波峰就是牙釉质峰，波峰位置即为牙釉质分割的最佳阈值。
图2
牙齿CT图像的直方图和牙釉质体积随横断面的变化
(a)灰度直方图；(b)局部放大直方图；(c)牙釉质体积随横断面的变化
Figure2. Histogram of the tooth CT image and volume difference of the tooth enamel between different slices(a) histogram; (b) local histogram; (c) volume difference of the tooth enamel between different slice
通过自动阈值分割我们得到牙釉质在牙齿上的分布，我们统计了牙釉质在所有横断面上的体积，统计结果如图 2所示。在图 2(c)中，横轴表示CT数据的512个横断面，纵轴表示每个横断面牙釉质的体积。在图中有两个波峰，左侧的波峰DT表示下颌骨牙齿的牙冠上的牙釉质，右侧的波峰UT表示上颌骨牙齿的牙冠上的牙釉质，波峰之间的波谷表示上排牙齿和下排牙齿的分界线。牙冠的分割是为后续分割提供精准的定位和准确的初始轮廓，所以我们选择范围在DT-UT之间进行牙冠的分割。
1.3 牙冠和牙颈分割
在牙冠部分，牙齿暴露在空气中，不与任何软组织和骨组织相连，牙齿边界非常清晰。我们利用牙釉质分割结果作为初始轮廓，利用CV-AC方法对牙齿的牙冠部分进行分割。Chan-Vese模型是Chan和Vese于2001年提出的基于区域的水平集活动轮廓模型，具有对噪声、杂波和对初始轮廓的设定不敏感等优点[10]。牙釉质(只覆盖在牙冠表面薄薄一层)作为初始轮廓，CV-AC模型通过较高次数的迭代演化后可以完美地实现牙冠的整体分割。
由于暴露在空气当中，牙冠相比牙颈和牙根是最容易分割的部分。在牙冠分割完成后，牙颈以下的部分因其直径减小，通常情况下是相互分离的，我们利用连通域运算获得独立的牙齿，然后分别完成对每个牙齿的整体分割。如果有相连的牙齿，我们利用CV-AC对牙齿继续向牙根分割，直到每个牙齿都相互独立。本文的CT数据在横轴方向上有非常高的分辨率(0.16 mm)，相邻横断面之间的牙齿轮廓不会有非常大的改变，所以我们可以把上一横断面的分割结果作为下一横断面分割的先验轮廓。这样每次分割的初始轮廓都非常接近于待分割的真实牙齿轮廓，我们只需要很少的迭代次数，活动轮廓就能找到牙齿的边界，进而可以极大地提高分割的精准度。但是除牙冠外，牙齿其他部位处于软组织和骨组织的包围中，尤其牙根部分与骨骼组织的边界非常狭窄且连贯性较差。为了更好地对牙齿进行分割，我们在CV-AC模型的基础上进行针对性改进，为CV-AC模型加入了灰度限制项和曲线演进先验限制项(PKL)，对于任意曲线变量C定义如下的“拟合能量函数”：
| $ \begin{array}{l} {F_1}\left(C \right) + {F_2}\left(C \right)=\int_{inside\left(C \right)} {{{\left| {I\left(x \right)-{c_1}} \right|}^2} + \omega } \ {\left| {I\left(x \right)-Threshold-I\left(x \right)} \right|^2}dx + \int_{outside\left(C \right)} {{{\left| {I\left(x \right)-{c_2}} \right|}^2} + } \ \omega \times {\left| {Threshold-I\left(x \right)} \right|^2}dx + PKL \end{array} $ |
式(2)中C是任意曲线变量，c1和c2是依赖于C的常数，它们分别是曲线C内部、外部区域的灰度平均值，Threshold是CT图像中骨骼和非骨骼组织的全局最佳骨骼阈值，该值也可以通过CT图像直方图获得，ω是骨骼灰度限制项的权重系数。其中，PKL是基于上一横断面先验知识的限制项，该项由上一横断面的牙齿轮廓经过距离变换生成。限制力大小与曲线移动距离成正比，当曲线移动距离超出一定范围时限制力变为无穷大，阻止曲线过大地偏离先验曲线。PKL的加入可以限制曲线的移动范围，在有限的范围内寻找最优的牙齿边界线。
1.4 牙根分割
由于牙根完全被骨骼包围，牙根相比牙冠和牙颈更加难以分割。基于1.3中图像增强和CV模型的分割方法在大部分情况下可以实现牙齿的准确分割，但是在一些情况下(犬齿和门牙)牙根会和颌骨紧密相连，边界处的灰度差别会非常小。当出现这种情况时，CV-AC模型在寻找牙齿边界过程中很容易在某一段边界发生错误。为了避免错误在下一横断面累积，我们采用图割的方法接替CV-AC继续分割牙齿。判断是否发生错误的方法是当活动轮廓分割错误时，活动曲线会在某一段曲线严重偏离先验曲线，重叠率(Dice系数)和质心位置会产生较大变化。图割是一种基于图论的组合优化方法，它将一幅图像映射成一个网络图，并建立关于标号的能量函数，运用最大流/最小割算法对网络图进行切割，得到网络图的最小割，即目标函数的最小值[11]。图割方法需要设定前景和背景种子点，种子点由上一横断面的牙齿分割结果作为先验知识生成，如图 3(a)所示；在牙齿边界内外分别选择两个像素加上边界线组成宽度为5个像素的窄带(灰色部分)，窄带内部区域(白色)设定为图割前景像素，窄带外部区域(黑色)设定为图割背景像素，灰色部分为待分割区域的划分结果，如图 3(b)所示。由此，可利用图割算法在初始轮廓区域内寻找一条最优的牙齿边界线，该方法在保证分割准确性的同时可避免程序陷入局部最优而无法找到真实边界的情况。
图3
待分割区域划分示意图
(a)牙齿先验轮廓；(b)待分割区域划分结果
Figure3. Diagram of the segmentation region determination(a) prior outline of the tooth; (b) result of the segmentation region
在牙冠分割完成后，图割方法就可以准确分割出大部分牙齿结构。但是，对于牙根分裂的情况，图割方法不能快速准确地识别出牙根的轮廓，这里主要是由于种子点的选取是基于上一横断面设定的，而种子点设定对图割的分割结果又有直接的影响。然而，CV-AC模型可以很灵活地进行闭合曲线的分裂和合并，能够准确处理牙根分裂的部位。所以，本文中只有当CV-AC模型出现错误才会使用图割方法对牙齿继续进行分割。
2 实验和结果
为了验证本文算法的有效性，本文在临床牙齿CT数据上进行了测试，测试在CPU为Intel Core I7-2600 3.4 GHz，内存为16 GB的PC机上运行，算法在Matlab平台上实现。
首先，我们对CT数据进行自动阈值分割得到了牙釉质在牙齿上的分布，如图 4所示。在图 4(a)中，牙釉质位于牙冠的外层形成一个环形，并且牙釉质主要位于牙齿的顶端。然后，我们把牙釉质阈值分割结果作为初始轮廓，利用CV-AC方法实现牙冠的完整分割。CV-AC方法有对初始轮廓不敏感和对噪声不敏感等优点，并且牙冠周围没有被覆其他身体组织，以牙釉质分割结果作为初始轮廓，经过较大次数的迭代(本文采用迭代600次)就可以获得牙冠完美的分割结果。图 4(b)中红色部分是牙釉质分割结果，白色部分是CV-AC经过600次迭代后得到的牙冠分割结果，图 4(c)是对DT和UT之间的所有横断面进行分割后牙冠分割结果。而如图 4(d)和图 4(e)所示则展示了在实际CT图像中CV-AC方法利用牙釉质分割结果作为初始轮廓完整分割出牙冠的过程。
图4
牙冠分割结果
(a)牙釉质阈值分割；(b)CV-AC分割；(c)牙冠分割三维渲染结果；(d)初始轮廓；(e)最终分割结果
Figure4. Results of the tooth crowns segmentation(a) result of the enamel threshold segmentation; (b) result of the CV-AC segmentation; (c) 3D rendering of the dental crown segmentation result; (d) initial contour; (e) final segmentation result
在完成牙冠分割后我们开始对牙颈和牙根进行分割，由于相邻横断面之间牙齿轮廓线只有很小的差别，这里我们利用上一横断面分割结果作为本横断面的先验轮廓。这时通过CV-AC方法对牙齿进行分割我们只需要很少的迭代次数(本文采用25次)就可以实现牙齿的精准分割。
相比牙颈被软组织(牙龈)包围，牙根完全被骨骼包围，牙齿和骨骼之间只有薄薄的一层牙周膜，并且在CT图像中非常模糊，所以分割难度很高。如果牙齿和颚骨的边界模糊不清，活动轮廓就会在骨骼和牙齿连接处被骨骼吸引，产生错误的结果。当CV-AC方法发生错误时，我们采用图割方法来寻找牙齿最优边界，如图 5所示。在图 5中，图 5(a)是牙齿原始CT图像，图 5(b)是作为先验轮廓的上一横断面牙齿分割结果，图 5(c)是基于先验轮廓设定图割的前景和背景(白色表示前景区域，黑色表示背景区域，灰色表示待分割区域)，图 5(d)为图割方法的分割结果。图 5显示无论牙齿边界比较清晰还是非常模糊，图割方法都能较好地分割出牙齿的轮廓。
图5
牙齿分割结果展示(3例)
(a)原始图像；(b)初始轮廓；(c)种子区域；(d)分割结果
Figure5. Tooth segmentation results based on the proposed method (three examples)(a) original image; (b) initial contour; (c) seed point region; (d) segmentation result
运用自动阈值分割、CV-AC模型和图割方法后，我们得到了牙齿的整体分割结果，如图 6所示。从牙齿分割结果的3D渲染图中我们可以看出本文提出的牙齿分割方法可以很好地分割出所有牙齿结构，对于分裂的牙根也保持了良好的准确度。图 6(b)是犬齿的单独3D渲染图，犬齿植入骨骼较深，分割难度相对较大，而采用本文提出的牙齿分割方法有效地分割出了整个牙齿结构。
图6
牙齿分割结果3D渲染图
(a)下排牙齿3D渲染图；(b)犬齿3D渲染图
Figure6. 3D rendering of the tooth segmentation results(a) 3D rendering of the segmented lower teeth; (b) 3D rendering of the canine tooth
3 结论
随着计算机图像处理技术的快速发展，医学图像处理技术在现代医学、生物医学研究和临床诊断等领域得到广泛的应用，并取得了非常好的效果。本文提出了一种基于区域自适应形变模型的CT图像牙齿结构测量方法。该方法首先利用自动阈值分割方法对牙齿进行分割，获得了牙釉质在牙冠上的分布情况，进而统计牙釉质在不同横断面上的分布，获得牙冠的位置和牙齿的方向。然后，我们利用牙釉质分割结果作为初始轮廓，使用CV-AC模型分割牙齿得到完整的牙冠轮廓。牙冠的分割结果作为先验轮廓用于后续牙齿轮廓的分割，这样可极大地减少迭代的次数，提高分割方法的运行效率。当活动轮廓发生错误时，可使用图割替代CV-AC方法进行牙齿分割。实验证明通过结合CV-AC和图割方法可以很好地实现牙齿整体的准确分割，解决了牙齿拓扑结构改变以及相邻牙齿间边界模糊的问题。对于临床数据的分割结果验证了方法的准确性和有效性，能够为医生的临床诊断和治疗提供辅助手段，在牙齿种植、牙齿正畸、医学信息学等多个方面有着重要的临床应用前景。
