石沁祎，閆 方，楊 陽，陳玥甫，林曉浪，王遠軍

上海理工大學 醫(yī)療器械與食品學院，上海 200093

引 言

錐形束計算機斷層掃描（Cone Beam Computed Tomography，CBCT）即錐形束CT，因具有掃描時間短、輻射劑量小、設備輕便等優(yōu)點而被廣泛應用于口腔診斷領域.通過對CBCT所得二維圖像進行重建，實現(xiàn)牙齒及牙槽骨等骨性結構的三維顯示，這將大大減少醫(yī)生的診斷難度，同時為治療方案的確定和虛擬手術的實現(xiàn)提供便利.但由于 CBCT具有輻射性，磁共振成像（Magnetic Resonance Imaging，MRI）逐漸成為了更安全的口腔檢測成像手段.因為成像原理的不同，牙齒及牙槽骨在CBCT圖像中為高密度影像，而在MRI中為低密度影像，且由于MRI主要適用于軟組織成像，骨性組織分辨率較差，因此雖能基本完整地展現(xiàn)牙齒與牙槽骨的形態(tài)與結構，但依然存在相鄰牙齒間黏連、牙根處拓撲結構改變、牙齒與牙槽骨皮質連接處等復雜結構無法清晰分辨的問題.利用MRI檢測下頜骨病，并將牙齒從圖像分割出來后，可以對牙槽骨和下頜骨病變作出更準確地診斷.

把牙齒和牙槽骨分別從二維圖像中分割出來，是實現(xiàn)骨性結構三維重建的基本前提；并且分割的精確度將直接影響重建的效果.因此，如何實現(xiàn)二維圖像中牙齒及牙槽骨的精確分割就成為一個關鍵問題.水平集方法由Osher和Sethian[1]在1988年提出，因其在描述復雜拓撲結構上的獨特優(yōu)勢，水平集方法近年來被廣泛應用于圖像識別、檢測等領域.國內外研究人員對于改進水平集算法有過許多研究.Gao等[2]首次將基于邊界、基于區(qū)域及基于先驗形狀約束能量項加入到水平集中，提出了一種混合的水平集模型.Zhang等[3]基于反應擴散正則化提出了一種新的水平集演化方法，佐以兩步裂項法，可大大降低計算的復雜性.由于單一的能量模型很難準確分割目標，Wang等[4]提出了一種由多尺度局部似然圖像擬合（Local Likelihood Image Fitting，LLIF）能量項、自適應先驗形狀約束能量項和反應擴散（Reaction Diffusion，RD）正則化能量項組成的混合水平集模型，能夠有效提高分割精度.

Li等[5]于2010年提出一種避免初始化水平集函數(shù)的距離正則化水平集演化（Distance Regularized Level Set Evolution，DRLSE）模型，其通過加入距離正則項（Distance Regularized Term）使零水平集函數(shù)在向目標物體邊緣移動的過程中，保持了其符號距離特性.但通過觀察其距離正則項所對應的勢阱函數(shù)，發(fā)現(xiàn)該勢阱函數(shù)在水平集函數(shù)梯度?φ為 0時會出現(xiàn)演化速度變至無窮的現(xiàn)象，這將直接導致水平集函數(shù)越過目標物體的某些較弱的邊緣，出現(xiàn)曲線侵入.而牙齒的牙冠部分通常情況下存在緊貼或擠壓，這種情況在CBCT與磁共振圖像上即表現(xiàn)為弱邊緣，故無法將該種勢阱函數(shù)所得到的距離正則項直接應用于對CBCT與磁共振圖像中牙齒的分割.基于此，本文在該勢阱函數(shù)的基礎上進行了改進，提出了一種新的勢阱函數(shù).該勢阱函數(shù)所產生的正則項克服了原勢阱函數(shù)對弱邊緣不敏感的缺點、并加快水平集函數(shù)演化速度，更適合用于對牙齒牙槽骨圖像的分割.

1 水平集模型

水平集方法能夠簡明地表示復雜拓撲結構的基礎思想在于其將n維曲面的演化問題轉化為n+1維空間中隱式方程解集合，即水平集的演化.因此二維曲線運動將轉化為三維曲面運動，即該曲線每一個時刻的變化均可利用三維曲面的零水平面來表達.圖1為水平集方法示意圖.

但利用水平集進行曲線演化的缺點是要進行周期性初始化，為避免復雜的初始化過程，Li等[5]提出了歸一化符號距離函數(shù).通過在能量泛函中構造勢阱函數(shù)決定的正則項來保證每一次演化后符號距離函數(shù)的值滿足一定條件，從而避免了對水平集函數(shù)即該符號距離函數(shù)的重新初始化.水平集函數(shù)定義為：

圖1 水平集方法示意圖Fig. 1 Diagram of level set method

其中，d為正常數(shù)，Ω表示初始輪廓區(qū)域.初始輪廓內部設定為負值，初始輪廓外部設定為正值.為驅動該水平集函數(shù)表示的曲線向著感興趣區(qū)域的邊界移動，定義水平集能量泛函E(φ)為：

(2)式第一項中，Rp(φ)為距離正則項，由(3)式表示：

其中P(·)為勢阱函數(shù)，?表示梯度算子，?φ為水平集函數(shù)梯度，根據勢阱函數(shù)最小點處?φ的取值，分為單勢阱函數(shù)和雙勢阱函數(shù)，在 DRLSE模型中，選用了雙勢阱函數(shù).由于構造勢函數(shù)的最終目的是保證符號距離的值，因此采用這種勢函數(shù)的正則項被稱為距離正則項.由于對勢阱函數(shù)取積分后所得的距離正則項可驅動曲線演化，所以針對不同的圖像情況，選擇恰當?shù)膭葳搴瘮?shù)是得到準確結果的關鍵.

(2)式第二項中，Lg(φ)為長度項，由(4)式表示：

其中δ(·)是Dirac函數(shù)，此處采用(5)式 δε(·)逼近δ(·)，ε常取1.5[5]：

邊緣指示函數(shù)g用(6)式表示，其中I0為待分割圖像，Gσ為標準偏差為σ的高斯核，Gσ*I0為高斯核與待分割圖像作卷積平滑后的圖像，?為微分算子符號：

在圖像真實的邊界處，梯度指向牙齒內部，然而相鄰牙齒、邊界的梯度和牙齒內部邊界（牙釉質和牙本質的邊界）的梯度都是指向牙齒外部.水平集符號距離函數(shù)的梯度方向總是指向輪廓的內部，即與真實邊界的梯度方向夾角成銳角.因此，只有當輪廓邊界的梯度方向和水平集符號距離函數(shù)的梯度方向指向同側時，得到的輪廓為目標邊界[6].根據文獻[7]的分析，g的范圍在0~1之間，當接近目標區(qū)域邊緣時取值較小，接近于0；當處于相同背景的區(qū)域時，取值較大，接近于1.

(2)式第三項中，Ag(φ)為面積項，由(7)式表示：

其中g為邊緣指示函數(shù)，H(·)為Heaviside函數(shù)，此處采用(8)式Hε(·)逼近H(·)，ε常取1.5，Hε(·)與 δε(·)的關系為

(2)式中，μ為正則權重系數(shù)；λ為長度權重系數(shù)；α為面積權重系數(shù)，當α為正時，輪廓向內收縮，α為負時，輪廓向外擴張，從而驅動零水平集曲線向邊界演化[8].在曲線演化的過程中，若曲線尚未到達待分割圖像目標區(qū)域的輪廓邊界，能量泛函E(φ)將驅動曲線向邊界靠攏；當水平集函數(shù)位于目標位置時，能量泛函E(φ)取得最小值，曲線停止演化，得到分割結果.

2 基于混合水平集模型的牙齒牙槽骨分割算法

2.1 DRLSE勢阱函數(shù)

正則項的作用不僅是平滑水平集函數(shù)，更重要的一點是令至少在零水平集附近區(qū)域內滿足|?φ|=1，這樣才可以保證曲線演化過程中相關計算的準確性.在勢阱方程中用s表示?φ，可以通過構造在s=1處取得最小值的勢函數(shù)p(s)來實現(xiàn)以上目的.根據勢阱函數(shù)最小點處s的取值，分為單勢阱函數(shù)和雙勢阱函數(shù).單勢阱函數(shù)表達式簡單，但在某些環(huán)境下，為得到最小化能量而衍生出的水平集演化會產生很嚴重的副作用，故多采用雙勢阱函數(shù)構造正則項.

Li等[5]在 DRLSE模型中采用了一種雙勢阱函數(shù)來保持水平集函數(shù)的符號距離特性即|?φ|=1，在該模型中，勢阱函數(shù)p(s)以及擴散速率函數(shù)dp(s)表達如下：

其中，s表示水平集梯度?φ.不難發(fā)現(xiàn)，該雙勢阱函數(shù)在零勢阱附近的演化速度接近于1，達到演化速度的最大值，使得水平集急速向前演化，最終使弱邊緣處的s減小至0，導致零水平集侵入被分割目標內部，無法實現(xiàn)對弱邊緣的準確檢測；這種不足使得將其運用在圖像分割的算法中時條件受限.

2.2 改進的雙勢阱函數(shù)

針對Li提出的雙勢阱的不足，孫超等[9]提出了V勢阱函數(shù)，該函數(shù)統(tǒng)一了演化速度的數(shù)學表達形式，在避免零勢阱修正的同時，降低了零勢阱附近的演化速度.本文在雙勢阱的基礎上提出了一種新的勢阱函數(shù)—單勾勢阱函數(shù).勢阱函數(shù)p(s)以及擴散速率函數(shù)dp(s)表達如下：

如圖2(a)所示，本文的勢阱函數(shù)同屬于雙勢阱函數(shù)，由零勢阱和壹勢阱構成.擴散速率函數(shù)為水平集梯度?φ的勢阱函數(shù)的導數(shù)除以?φ，它能表示?φ的變化方向和速度，為方便表示，用s代替水平集梯度?φ.如圖2(b)所示，當0≤s≤0.5時，勢阱函數(shù)的演化速度大于零，水平集向前演化，s向 0靠近，保持水平集為常數(shù)；當0.5＜s＜1，勢阱函數(shù)的演化速度小于零，水平集向后演化，s增大，向 1靠近；s≥1時，勢阱函數(shù)的演化速度大于零，水平集向前演化，s減小，同樣使s向1靠近.這使得水平集函數(shù)保持符號距離特性.

圖2 單勾勢阱函數(shù)圖.(a)勢阱函數(shù)；(b)擴散速率函數(shù)Fig. 2 Diagram of single hook well function. (a) Potential well function; (b) Diffusion rate function

與Li等提出的DRLSE模型中的雙勢阱函數(shù)相比，本文提出的單勾勢阱函數(shù)有如下優(yōu)勢：（1）擴散速率函數(shù)在s=0處時，不存在分數(shù)計算，因此不需要對分母取0時的函數(shù)進行修正，避免了算法的下降；（2）如圖3所示，擴散速率函數(shù)在s=0處時的值較小，即勢阱函數(shù)在零勢阱附近的演化速度偏小，不會導致s急速減小，避免了曲線前進速度太快導致水平集侵入被分割目標內部，因此可以增強對弱邊緣的檢測準確性.

在以往的水平集模型中，為了抑制由噪聲引起的非均勻性，常常使用大尺度高斯濾波對圖像進行預處理[10].而圖像中的紋理變化常常是非均勻性的，恒定方差的高斯濾波不能平滑圖像區(qū)域的非均勻性.若方差較大，對應的高斯濾波器模板尺度也較大，圖像邊緣損失較多，容易導致目標區(qū)域出現(xiàn)過分割；反之，若方差較小，高斯濾波器模板尺度較小，圖像非均勻性程度過大，容易導致目標區(qū)域出現(xiàn)欠分割[10].另外，DRLSE模型是基于邊緣的分割方法，此類方法的主要缺點有：（1）初始輪廓會影響分割的精確度；（2）對圖像噪聲很敏感；（3）在弱邊緣情況下演化結果不好[11].而邊緣指示函數(shù)g依賴光滑圖像的梯度，所以在分割含噪圖像時，演化曲線受噪聲干擾容易陷入虛假邊緣[12].因此，當目標之間的距離很近時，此時若使用大高斯核濾波，會使邊界處的灰度值被其鄰域的加權灰度值取代，而其鄰域中像素多為目標像素，故原先低灰度邊界灰度上升，到一定程度后，邊界消失，目標邊界可能相互連接，導致分割失敗[13].

圖3 單勾勢阱函數(shù)與Li等提出的雙勢阱函數(shù)對比圖.(a)勢阱函數(shù)對比；(b)擴散速率函數(shù)對比Fig. 3 Comparison of single hook well function and double well function proposed by Li et al. (a) Comparison of potential well function; (b) Comparison of diffusion rate function

為了得到較好的分割效果，將序貫濾波應用于本文的水平集模型，利用小方差高斯濾波的優(yōu)勢，對圖像進行多次小方差濾波的疊加，平滑過程可抽象為：

其中，I0為輸入圖像，I為高斯濾波后的圖像，g(σ1)、g(σ2)是方差分別為σ1、σ2的高斯核，根據高斯濾波函數(shù)的數(shù)學表達式，可將該濾波過程表達為：

(14)式表明，兩次高斯函數(shù)方差σ1、σ2的序貫濾波可等效為方差為σ3的一次高斯濾波.則可理解為，兩次較小方差的高斯濾波疊加產生的平滑效果可等效為一次較大方差的高斯濾波.基于此，將疊加較小方差高斯濾波的次數(shù)與水平集迭代次數(shù)相關聯(lián)，在設定好小方差的初始值σ0后，可由前面的公式推得第k次序貫濾波的方差σk為：

為避免水平集迭代次數(shù)過多導致目標區(qū)域過度平滑失去邊緣特征而出現(xiàn)曲線侵入，文獻[10]提出用區(qū)域置信度作為判定的標準之一，根據不同等效方差的高斯函數(shù)分割結果的區(qū)域置信度決定序貫濾波次數(shù)，其區(qū)域置信度Pr定義為：

其中，Ak-1為第(k-1)次小方差高斯濾波后，即第(k-1)次水平集演化的目標分割區(qū)域，S(Ak-1)為其對應區(qū)域的面積.同理，Ak為第k次迭代得出的目標區(qū)域，S(Ak)為其對應的面積.(16)式的分子為兩次分割區(qū)域的交集，在算法的實現(xiàn)上，判定區(qū)域的交集會稍微復雜一些.本文結合 DRLSE算法的優(yōu)勢，在對牙齒進行分割時，常常將初始輪廓設定在目標區(qū)域外部，故水平集多保持向內收斂，故上式中的S(Ak-1∩Ak)即為 m in{S(Ak-1),S(Ak)}，意為S(Ak)與S(Ak-1)中的最小值.這樣，不需要在引入判定交集增加算法的運行時間，只需簡單的選擇結構語句就可實現(xiàn)區(qū)域置信度的判定.故可得區(qū)域置信度為：

在以往的研究中，對于同一患者的不同牙齒多采用相同尺寸的牙盒作為初始輪廓，這樣雖然在選取時較為方便，但卻是以增加水平集迭代次數(shù)、增加數(shù)據提取時間作為代價的.本文將針對不同牙齒的特點，設置不同尺寸、不同迭代次數(shù)，以提高分割效率.為了進一步利用相鄰序列圖片中同一結構的相似性，本文將先驗信息[14]加入對下層圖像的分割中，以當前層對單一牙齒的分割結果作為下層圖像中同一牙齒的初始輪廓，以求在相對少的迭代次數(shù)下得到準確的結果.考慮到牙齒的部分圖像會出現(xiàn)當前層輪廓小于下層輪廓的現(xiàn)象，故在定義下層初始輪廓時先對當前層最終輪廓做形態(tài)學膨脹處理，以避免出現(xiàn)曲線侵入而使得結果不準確.

2.3 基于閾值法和孔洞填充的牙槽骨分割

采用閾值法[15]對牙槽骨分割時，由于牙髓區(qū)域位于牙齒內部，并且灰度值較低，通常被視為背景區(qū)域，導致分割結果產生孔洞.Wang等[16]通過孔洞填充校正閾值分割得到的掩模，再利用自適應擴散流（Active Diffusion Flow，ADF）主動輪廓模型對單顆牙齒進行分割.本文提出一種基于邊緣信息檢測的填充算法，如圖4所示，其主要思想為：對閾值分割后已被識別為背景的像素點進行十字形擴散檢測，若以該像素點為中心向外擴散的同半徑十字上均存在非背景點，則認為此處為孔洞即中心像素點為非背景點，依次檢測非背景點后即可將牙槽骨分割過程中出現(xiàn)的孔洞全部填充，并避免出現(xiàn)邊界泄露的情況.

圖4 十字形擴散檢驗示意圖Fig. 4 Diagram of cross diffusion test

2.4 牙齒牙槽骨分割算法步驟

牙齒牙槽骨分割算法步驟如圖5所示，具體步驟如下：

Step1. 參數(shù)初始化：設定步長、正則項系數(shù)μ、長度項系數(shù)λ、面積項系數(shù)α、內外迭代次數(shù)

Step2. 讀取DICOM圖像，對其進行灰度均一化處理和對比度增強

Step3. 牙齒分割

a.水平集初始化，對第一層圖像手動取框作為初始水平集曲線

b.在區(qū)域置信度內進行小尺度高斯濾波，并計算圖像梯度，生成邊緣停止函數(shù)

c.計算單勾勢阱函數(shù)演化速度

d.計算水平集能量項最小化迭代公式，使其驅使輪廓線停止在牙齒邊緣處得到分割結果

e.讀取下一層圖像，對上一層的分割結果小尺度膨脹后作為下一層的初始水平集曲線

f.重復步驟b~e，直到分割完全部圖層，完成此顆牙齒的序列分割

g.重復步驟a~f，進行下一顆牙齒的分割，直到完成所有牙齒的序列分割

Step4. 牙槽骨分割

a.通過閾值法對牙槽骨進行分割，并對孔洞進行填充

b.讀取下一層圖像

c.重復步驟a~c，直到分割完全部圖層，完成牙槽骨分割

圖5 本文提出的牙齒牙槽骨圖像分割算法的流程圖Fig. 5 Algorithm flow chart for segmentation of tooth and alveolar bone proposed in this research

3 實驗結果

本文實驗數(shù)據為10組CBCT口腔圖像以及3組口腔磁共振圖像，采集于上海市第九人民醫(yī)院，其中CBCT掃描電壓為90 kV，曝光時間為10.8 s，獲得數(shù)據矩陣大小為595×595×358，像素分辨率為 0.25×0.25×0.25 mm3；MRI掃描序列為 T1W_TSE，TR為 507 ms，層厚為 4 mm，圖片尺寸為336×336.圖像的數(shù)據類型為16位無符號整型，保存的格式為DICOM文件.

本文的實驗環(huán)境為：CPU處理器為Inter(R) Core(TM) i5 CPU(2.30GHz)，內存8.0GB，MATLAB R2016b軟件下編寫程序完成實驗，實驗參數(shù)為μ=0.04、λ=5、α=1.5、內迭代8次、外迭代16次.

3.1 牙齒分割效果

將16位數(shù)據轉換為8位的灰度圖像，手動畫取合適大小的相同初始框，并選擇相同的實驗參數(shù)后，分別采用DRLSE模型和本文提出的水平集模型對CBCT圖像及磁共振圖像中的牙齒進行分割.圖6為CBCT圖像分割結果，可以發(fā)現(xiàn)Li等的方法可能會把相鄰牙齒的輪廓錯分為目標牙齒的輪廓，導致欠分割.而本文提出的算法模型在相鄰牙齒相互粘連、牙齒和牙槽骨灰度值相近的情況下，能夠驅使水平集函數(shù)表示的曲線到達目標牙齒的邊界處，而不會停留在相鄰牙齒和牙槽骨的分界處.由此可見，本文的算法對牙齒有精確的分割能力，可以很好地改善DRLSE模型的欠分割的情況.圖7為文本算法運用于磁共振圖像的實驗結果，亦可準確有效地分割出磁共振圖像中的單顆牙齒.

圖6 基于CBCT圖像的(a) DRLSE模型牙齒分割結果和(b)本文模型牙齒分割結果Fig. 6 Tooth segmentation with (a) DRLSE model, and (b) the proposed model in this research based on CBCT images

圖7 基于磁共振圖像的本文模型牙齒分割結果Fig. 7 Tooth segmentation with the proposed model in this research based on magnetic resonance image

本文從分割的精確性及運算時間兩個角度來評價算法的優(yōu)劣.采用積重疊誤差（Volumetric Overlap Error，VOE）[17]來進行基于區(qū)域的精確性描述，VOE定義如下：

其中A為算法自動分割的結果，B為手動勾畫結果并作為標準.通過計算水平集分割算法和手動分割結果的交集和并集的體積比得到的就是兩者的真實重疊度，VOE越小，說明分割的結果與標準的重疊度越高，即分割結果越準確.基于CBCT圖像，本文算法與DRLSE分割結果的定量比較如表1所示，可見本文算法在分割精確度上相比DRLSE有明顯的提升，算法效率上也有一些提高，但由于序貫濾波增加了算法的時間，因此該優(yōu)勢并不明顯.

表1 基于CBCT圖像，本文算法與DRLSE模型進行牙齒分割結果的定量比較Table 1 The quantitative comparison for tooth segmentation between the algorithm proposed in this research and DRLSE model based on CBCT images

利用本文提出的水平集算法對同一樣本中同一顆牙齒的CBCT多層序列圖像進行分割，結果如圖8所示.對第一層圖像進行手動畫框后，即可通過本文的水平集算法對目標牙齒進行準確的分割，將此層分割結果曲線進行小尺度膨脹后作為下一層的水平集初始輪廓，即可無中斷地進行序列圖像的牙齒分割，避免了重復初始化.實驗結果表明，利用本文的方法能過對CBCT序列圖像完成準確的牙齒分割，基本不會出現(xiàn)欠分割或是過分割的情況.

圖8 基于CBCT序列圖像的牙齒分割結果Fig. 8 Tooth segmentation based on CBCT sequence image

3.2 牙槽骨分割效果

采用本文提出的水平集算法可對一層圖像中的所有牙齒依次完成分割，得到的牙齒輪廓即為牙槽骨的內輪廓.再采用閾值分割法可得到牙槽骨的外輪廓，由于牙槽骨和牙齒內部存在灰度值較小的區(qū)域，閾值法得到的結果會存在孔洞，通過采用前文所述的孔洞填充的方法，選擇尺度為 15個像素點對其進行填充，即可分割出完整的牙槽骨，CBCT圖像牙槽骨分割結果如圖9所示，紅色曲線為各顆牙齒的輪廓，綠色曲線即為牙槽骨的外輪廓，紅、綠曲線之間的區(qū)域即為分割出的牙槽骨.對下牙槽骨進行序列分割后，將結果進行三維重建，結果如圖10所示.

圖9 (a) CBCT原圖；(b)使用本文算法的牙槽骨分割結果Fig. 9 (a) Original CBCT image; (b) Segmentation results of alveolar bone using the algorithm proposed in this research

圖10 下牙槽骨三維重建結果Fig. 10 3D reconstruction results of lower alveolar bone

4 結論

本文針對CBCT及磁共振圖像中牙齒的特點，對以往水平集模型中不適合應用于此類圖像的部分做出了改進，并將其應用在牙齒CBCT及磁共振圖像分割中，得到了更為準確的分割結果.其中創(chuàng)新的單鉤勢阱函數(shù)克服了現(xiàn)有勢阱函數(shù)使用限制多、弱邊緣難以檢測等缺點；將序貫濾波與水平集模型結合改善了原有水平集模型中單次使用大尺度高斯核而使圖像信息丟失的情況；根據斷層序列圖像特點，將當前層分割結果作為下層圖像初始框，實現(xiàn)先驗信息的充分利用.

利益沖突

無