針對三維計算機斷層掃描(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和圖割方法可以很好地實現牙齒整體的準確分割，解決了牙齒拓撲結構改變以及相鄰牙齒間邊界模糊的問題。對于臨床數據的分割結果驗證了方法的準確性和有效性，能夠為醫生的臨床診斷和治療提供輔助手段，在牙齒種植、牙齒正畸、醫學信息學等多個方面有著重要的臨床應用前景。
