一種對錐束CT圖像進行牙列分割的方法與流程
2023-04-26 03:19:11 2
本發明涉及計算機視覺和圖像處理技術領域，尤其涉及一種對錐束ct(cbct)圖像進行牙列分割的方法。
背景技術：
錐束ct(cbct)圖像常用於頜面與口腔正畸外科輔助手術預測以及牙列對齊計劃的制定。cbct圖像可以提供病人特定的解剖結構信息，其中包含完整的牙列信息。而傳統的基於光學的方式，例如三維雷射掃描以及立體視覺設備，僅僅能夠獲取牙冠的幾何信息，其中由於牙根被埋藏在牙齦以及頜骨中無法獲取其幾何形態。在臨床口腔外科中，cbct圖像由於其低放射劑量以及較低的採集成本，相對於傳統的ct成像技術具有極大的優勢。但是，較低的放射劑量以及信噪比一般會造成cbct圖像質量較差，常常出現圖像退化現象，例如口腔治療以及牙齒種植在牙頜結構中放置的金屬物體造成的線束硬化(beamhardening)問題，有限視域的截斷，以及由於牙齒與周圍牙槽骨相似的灰度而造成的模糊的輪廓。同時在圖像採集過程中，病人可能有微小的頭部運動，也會造成牙齒輪廓的模糊。特別是為了診斷牙齒畸形，cbct在採集過程中通常上下牙列處於咬合狀態，這使得上下牙列在咬合部分的分割變得更加困難。
近年來有大量技術被應用於從ct以及cbct圖像中進行牙列分割，其中包括使用積分灰度投影、基於閾值和區域增長的方法、圖割法、基於水平集活動輪廓的方法等。積分灰度投影方法僅僅能獲取牙列的粗略分割，例如不同牙齒之間的分割平面。閾值與區域增長方法對於cbct中可能出現的牙列與牙槽骨以及其它組織之間模糊的輪廓難以有效處理，其分割與真實的牙齒輪廓的一致性較差。基於交互的分割方法例如圖割法非常依賴初始前景區域的定義，對於不同的初始前景定義所得到的分割常常存在差異，難以重複分割結果。為了獲取可靠的分割，活動輪廓方法諸如水平集技術集成形狀以及灰度先驗、從切片圖像分割獲取的三維形狀先驗、以及隱參數表達模型，逐層對cbct圖像進行分割。在基於水平集的方法中，對於約束項需要進行精心設計以避免分割過程中可能出現的收縮(shrinkage)以及洩露(leakage)問題。此外，在基於活動輪廓方法進行逐層分割中還可能有誤差累計的問題。此外，三維牙列模型被引入分割系統，包括三維形狀圖集以及三維形狀統計模型。圖集以及統計形狀模型一般來自大量的ct圖像，這勢必會增加對數據的需求以及數據處理的負擔。
技術實現要素：
為了克服上述現有技術的不足，本發明提供一種對錐束ct(cbct)圖像進行牙列分割的方法，基於可變形模型的隨機遊走方法，從cbct圖像中分割得到完整牙列。
本發明的原理是：將體圖像中的分割問題定義為圖結構中的標籤擴散問題，其中，預先對少量切片圖像進行標註，並使用擴散機制對感興趣體圖像(volumeofinterest，voi)中其它像素進行自動標註。首先，利用原始的隨機遊走方法對cbct圖像中牙列進行初始分割，基於原始的隨機遊走的初始分割通常會包含分割錯誤，例如由於cbct圖像退化出現的線束硬化、圖像中模糊的牙列輪廓所造成的標記擴散誤差。因此，引入三維可變形模型用於擬合由標記擴散得到的體圖像分割，利用組合logistic函數，並基於經驗牙本質厚度定義柔化約束，將柔化約束引入隨機遊走算法中更新體素的類別標籤的估計。此外，為了改善從單步柔化約束下的標記擴散所得到的圖像分割精度，引入迭代修正算法，針對系統的兩個方面，即柔化約束下標記擴散所對應的大規模稀疏線性系統和三維模型到圖像分割所得到的牙列表面體素的擬合，進行迭代求解。在每步迭代中，三維可變形模型的配準可以看作是對體素標記的正則化，該過程可以有效地消除分割誤差。
本發明提供的技術方案是：
一種對錐束ct圖像進行牙列分割的方法，針對錐束ct圖像中的感興趣體圖像進行牙列分割，得到牙列表面模型；包括如下步驟：
1)針對錐束ct圖像中的感興趣體圖像區域定義圖結構；所述圖結構中，結點集合對應感興趣體圖像區域中所有的體素；邊連接集合為感興趣體圖像區域中相鄰體素之間的邊連接；
2)根據牙列分布情況定義類別標籤集合，針對步驟1)所述感興趣體圖像中所有的體素，採用隨機遊走方法進行牙列的初始分割，得到所述所有體素的類別，對所述所有體素進行類別標記；
3)根據步驟2)得到的所有體素的類別，利用其中屬於牙列的表面體素定義牙列的表面形態，生成初始分割得到的表面體素集；進一步利用三維可變形模型，通過非剛性變形配準擬合所述牙列表面體素集合，得到三維可變形模型的非剛性變換參數，用於擬合所述初始分割得到的表面體素集；
4)基於三維可變形模型定義柔化約束，並基於柔化約束下的隨機遊走方法進行牙列分割，得到感興趣體圖像區域中所有體素的類別標籤；
5)通過迭代修正過程對牙列分割進行迭代修正，結合步驟4)所述柔化約束下的類別標籤擴散和步驟3)所述三維可變形模型的非剛性配準過程，改進單次柔化約束下隨機遊走方法所得到的體圖像分割結果。
針對上述對錐束ct圖像進行牙列分割的方法，進一步地，步驟1)所述圖結構中，採用基於上下文的體素特徵描述子，通過體素的表觀計算相鄰體素之間的相似度，用於表示體素的上下文特徵差異；所述體素特徵描述子向量中的元素對應當前體素的包圍圖像塊與模式中的採樣體素的包圍圖像塊之間的像素灰度差異累計；圖結構中邊連接上的權重通過體素上下文特徵描述差異和灰度差異的加權組合得到。
上述上下文的體素描述子在當前體素的一個半徑為的包圍球體中通過隨機採樣生成一個模式p，所述模式p由採樣得到的體素集合表示；所述體素特徵描述子向量f(vi)中的元素對應當前體素的包圍圖像塊與模式p中的採樣體素的包圍圖像塊之間的像素灰度差異累計，表示為式11：
f(vi)＝{dp(c(vi),c(vi+γ))|vi+γ∈p}(式11)
式11中，f(vi)表示當前體素vi的特徵描述子向量；c(vi)表示體素vi的包圍圖像塊，dp為圖像塊的灰度差異；γ為當前體素與模式p中採樣體素之間的位移；c(vi+γ)表示模式中的採樣體素的包圍圖像塊；其中，圖像塊的灰度差異dp定義為式12：
式12中，dp為圖像塊的灰度差異；c(vi)與c(vk)分別表示體素vi與vk的包圍圖像塊；δv為包圍體素塊中體素的位移；vi+δv屬於圖像塊c(vi)；i(vi+δv)為體素vi+δv的灰度；i(vk+δv)為體素vk+δv的灰度；
所述圖像塊的灰度差異可通過圖像的線性卷積算子計算得到：對差異圖像v-v(γ)進行卷積，通過式13得到：
dp＝(v-v(γ))2*g(式13)
式13中，dp為圖像塊的灰度差異；γ為體素vi與模式中一個採樣體素之間的位移向量；對體圖像進行平移操作得到v(γ)；v-v(γ)為差異圖像；g為卷積核；
圖結構中邊連接上的權重通過上下文特徵描述差異與灰度差異的加權組合得到，如式14：
式14中，wij為體素i與相鄰體素j之間邊連接的權重；α為常數係數用於調節上下文特徵以及體素灰度對於權重計算的影響；ii與ij對應第i個與第j個體素的灰度，fi與fj對應第i個與第j個體素的上下文描述；ψ為圖結構的邊連接集合。
針對上述對錐束ct圖像進行牙列分割的方法，進一步地，步驟2)採用隨機遊走方法對感興趣體圖像區域中所有的體素的類別進行標記，得到牙列的初始分割結果；具體通過最小化式21的能量函數得到感興趣體圖像中體素類別標記：
式21中，erw(x)表示基於隨機遊走算法的圖像分割能量；|s|為系統類別標籤的個數；表示第i個結點被分為第s類的概率；系統中每個結點對均對應一個|s|維的概率向量，其中第i個結點所屬的類別定義為表示預先由用戶交互定義的結點概率，當第i個結點被用戶交互定義為第s類，則取值為1，否則為0；nl是用戶預先交互標註的體素個數；nv為感興趣體圖像中所有體素的個數；μ1與μ2為常數；l為拉普拉斯矩陣；h為對角線指示矩陣，當感興趣體圖像中的第i個結點被預先標註時，h的元素hii＝1；x與表示感興趣體圖像中所有體素對應的類別概率以及所組成的矩陣；。
其中，將所述式21的目標函數轉化為矩陣形式，通過將所述式21的一階導數設為0，將目標函數求解轉化為求解式22所示的大規模稀疏的線性系統：
式22中，表示單位矩陣；l為拉普拉斯矩陣；h為對角線指示矩陣，當感興趣體圖像中的第i個結點被預先標註時，h的元素hii＝1；x與表示感興趣體圖像中所有體素對應的類別概率以及所組成的矩陣；
通過求解式22的線性系統，得到感興趣體圖像中所有體素的類別，作為牙列的初始分割結果。
針對上述對錐束ct圖像進行牙列分割的方法，進一步地，步驟3)所述三維可變形模型用於處理錐束ct圖像自身的局部退化造成的邊界混淆問題；所述通過非剛性變形配準擬合所述牙列表面體素集合，具體包括如下步驟：
31)利用procrustes分析計算三維模型的全局變換使其最大程度地與表面體素集yc一致；
32)隨後利用非剛性配準獲取特定的三維模型，使得所述特定的三維模型與初始分割的表面體素集一致；
33)為降低非剛性變換的參數空間，利用嵌入變形方法，將所述特定的三維模型表面的非剛性變形定義為一個稀疏基網格頂點變換的線性組合，三維模型非剛性配準的能量函數定義如式31：
式31中，為變形前後表面頂點的位移；變形後的三維模型ye的頂點為其中，權重βij對應頂點與其在基網格b上的近鄰頂點ηj，並基於兩點之間的歐式距離定義為其中參數tj為第j個基網格頂點上的三維變換矩陣；能量函數ere的第一項用於最小化變形後的三維模型表面與體圖像分割獲取的表面體素集之間的距離；dc表示變換後的頂點與標籤擴散得到的表面體素集yc之間的距離,其中為表面體素；第二項用於保持三維模型表面平滑，該項最小化頂點與其近鄰的空間變換之間的差異；ω為常數係數，ne為三維模型表面頂點個數。
34)通過最小化能量函數ere，得到三維模型的非剛性變換參數，用於配準擬合所述牙列表面體素集合。
針對上述對錐束ct圖像進行牙列分割的方法，進一步地，步驟4)所述基於柔化約束下的隨機遊走方法進行牙列分割，包括如下步驟：
41)通過設定與體圖像分割結果擬合後的三維模型表面內部體素更有可能屬於牙列，根據體素結點概率定義柔化約束；
42)設定柔化約束下的感興趣體圖像的標籤擴散能量為式43：
erwsc(x)＝λ1erw(x)+λ2esc(x)(式43)
其中erw與esc式21式42所定義的隨機遊走以及柔化約束的能量。λ1與λ2為常數係數。
將式43中的能量函數的一階導數設為0，得到關於體素結點概率向量的線性系統如式44，
其中l為拉普拉斯矩陣；h為對角線指示矩陣，當voi中的第i個結點被預先標註時，h的元素hii＝1；x與表示voi中所有體素對應的類別概率以及用戶定義的先驗概率所組成的矩陣。對應感興趣區域中所有體素的柔化約束所構成的矩陣；λ1與λ2為常數係數,實驗中分別設置為1、0.1；μ1與μ2為常數，本實驗中分別設置為1與0.001；表示單位矩陣。
43)利用共軛梯度方法求解線性系統，該線性系統需要求解|s|-1次，其中|s|為系統中類別標籤的個數；通過求解線性系統得到感興趣體圖像中所有體素的類別標籤。
進一步地，步驟41)採用以下步驟定義所述柔化約束：
410)根據牙齒表面內部的灰度分布，對牙本質賦以更大的屬於牙齒的概率，定義基於三維可變形模型屬於前景牙齒的概率為式41：
式41中，pd(r,θ)為基於三維可變形模型屬於前景牙齒的概率；r1(θ)與r2(θ)表示牙齒的外輪廓以及牙髓腔輪廓，其中r2(θ)＝r1(θ)-a，a為經驗的牙本質厚度；參數η用於控制如上概率函數經過牙齒輪廓時函數的形狀；
411)針對沒有牙髓腔的分層圖像，基於三維可變形模型屬於前景牙齒的概率定義為式411：
式411中，
412)針對牙齒可能出現的多牙根以及多牙尖結構情況，牙齒在切片上呈多輪廓狀態時基於三維可變形模型屬於前景牙齒的概率表示為式412：
式412中，nc為當前牙齒在切片圖像中輪廓的個數；pd為單外輪廓的基於三維可變形模型屬於前景牙齒的概率(式41與式411)；r,θ為體素在極坐標系中的坐標；
413)利用卡方核的一對多支持向量機生成分類器，所述分類器的輸出被歸一化到0-1之間，表示基於體素表觀屬於前景(牙列)的概率pc，柔化約束定義為式413：
式413中，pc表示由體素vi表觀估計的體素屬於前景(牙列)的概率；表示在體素vi基於三維可變形模型屬於前景牙齒的概率；r,θ為體素在極坐標系中的坐標；是對體素vi相對於類別s所定義的柔化約束。
針對上述對錐束ct圖像進行牙列分割的方法，進一步地，步驟5)所述通過迭代修正過程對牙列分割進行迭代修正，具體包括如下步驟：
51)將柔化約束下的標籤擴散和三維模型的非剛性配準能量結合，定義能量如式51：
e(x,t)＝λ1erw+λ2esc+λ3ere(式51)
式51中，e(x,t)是基於三維可變形模型隨機遊走算法進行牙列分割的總能量函數；erw與esc為式21與式42所定義的隨機遊走以及柔化約束的能量函數；ere為式31所定義三維模型非剛性配準的能量；λ1，λ2與λ3為常數，實驗中分別設置為1、0.1、1；；
52)在第一個階段，設定基於標籤擴散得到的體圖像分割獲取牙列表面體素集，由於表面體素集由第i-1次迭代中的cbct圖像分割得到，所以將與柔化約束下的隨機遊走相關的項erw和esc從能量函數式51中去掉，通過最小化能量e(x(i-1),t)對三維模型進行變形，用於擬合牙列的表面體素集；此階段對於三維模型變形參數的求解和步驟3)相同；
53)在第二個階段，去掉能量函數中與三維模型配準相關的項ere，最小化能量e(x,t(i-1))，得到感興趣圖像中所有體素的標籤；此階段通過求解一個大規模稀疏的線性系統求解所有體素的標籤，和步驟4)所述柔化約束下的隨機遊走方法相同；
54)迭代優化進行步驟52)和步驟53)，當相鄰兩步迭代中三維模型更新||ye(i)-ye(i-1)||2小於預先給定的閾值時終止迭代，得到感興趣圖像中所有體素的標籤，作為迭代修正後的牙列分割結果。其中ye(i)與ye(i-1)分別為第i步與第i-1步三維表面模型。
針對上述對錐束ct圖像進行牙列分割的方法，進一步地，所述得到的牙列表面模型與手工分割得到的表面模型的平均距離誤差在0.3mm以下。
與現有技術相比，本發明的有益效果是：
本發明提供一種對錐束ct(cbct)圖像進行牙列分割的方法，基於可變形模型的隨機遊走方法，從cbct圖像中分割得到完整牙列。本發明方法結合半監督的隨機遊走算法和三維可變形模型配準定義的柔化約束，獲取牙列的可靠分割。其中，利用三維可變形模型引入體圖像的柔化約束，可以看作為正則化項用於處理基於半監督標籤擴散的體圖像分割中的噪聲。本發明還採用迭代修正方法，對柔化約束下標記擴散以及三維模型對表面體素集的擬合問題，進行迭代求解，可有效地消除分割誤差，改進單次標籤擴散所獲取的牙列分割，提高分割結果的精度。
利用本發明提供的方法，可以對醫學臨床cbct圖像進行自動的牙列分割，分割得到的三維牙列模型可用於頜面與口腔正畸外科輔助手術預測以及牙列對齊計劃的制定，該牙列分割可滿足口腔臨床對於精度的要求。度量自動分割所得到的牙列表面模型與手工分割對應的表面模型的平均距離，該距離誤差在0.3mm以下。
附圖說明
圖1是本發明提供方法的流程框圖。
圖2是本發明中採用的基於上下文的體素描述子的結構示意圖；
其中，vi為當前體素；為當前體素的一個包圍球體的半徑；γ為當前體素與模式p中採樣體素之間的位移；c(vi)為當前體素vi的包圍圖像塊；c(vi+γ)為模式中的採樣體素的包圍圖像塊。
具體實施方式
下面結合附圖，通過實施例進一步描述本發明，但不以任何方式限制本發明的範圍。
cbct圖像是三維體圖像，由一系列的二維圖像組成，本說明將其中二維圖像稱作切片圖像(分層圖像)。本發明實施例針對醫學臨床cbct圖像，基於cbct圖像中感興趣體圖像區域中定義的圖結構和少量用戶交互標註的分層圖像(切片圖像)，對cbct圖像中的牙列進行分割。其中，採用三維可變形模型定義柔化約束，通過基於柔化約束下的隨機遊走算法更新體圖像中的牙列分割，再以迭代修正方法獲取可靠的分割結果。
圖1是本發明提供方法的流程框圖。具體地，本發明提供方法包括如下步驟：
步驟1：設定cbct圖像中的感興趣體圖像(voi)，定義voi中的圖結構；
本方法以一種半監督的方式對cbct圖像進行牙列分割，其中，將用戶所指定的少量切片圖像中的標記擴散到感興趣體圖像的整個區域。該擴散過程在voi中定義的圖結構中進行。針對voi定義的圖結構中，結點集合對應voi中所有的體素，邊連接集合為相鄰體素之間的邊連接。
一般地，圖結構的邊連接集合包含兩類，一類對應切片圖像內部的邊連接，該連接對應切片圖像中規則的像素網格；另一類對應切片圖像之間的邊連接，依據由稠密光流算法獲取的切片圖像之間的對應體素建立邊連接。值得注意的是，由稠密光流算法建立的對應並非一一對應，即第i層切片圖像中的一個結點可能與第(i+1)層切片圖像中不止一個結點對應。該情況常常發生在牙齒區域分叉的切片圖像，例如磨牙的多牙根區域，或者牙冠中的多牙尖區域。
系統基於體素的表觀計算相鄰體素之間的相似度(為上下文特徵描述差異與灰度差異的加權組合)。使用體素的灰度i計算體素的相似度是一種較直觀的方式，但是僅僅利用體素灰度定義體素相似度對於低輻射劑量採集的cbct圖像並不適合。這是由於cbct圖像中可能存在退化區域，其中僅由灰度相似並不能說明對應體素相似。因此系統引入基於上下文的體素描述子。
圖2是基於上下文的體素描述子；其中，當前體素的一個半徑為的包圍球體中，γ為當前體素與模式p中採樣體素之間的位移；通過隨機採樣生成模式p；特徵描述子向量f(vi)中的元素對應當前體素vi的包圍圖像塊c(vi)與模式中的採樣體素的包圍圖像塊c(vi+γ)之間的像素灰度差異累計。
為了降低計算的複雜度，該上下文的體素描述子並未對當前體素的所有近鄰體素進行計算，而是在當前體素的一個半徑為的包圍球體中通過隨機採樣生成一個模式p。該模式由採樣得到的體素集合表示。特徵描述子向量f(vi)中的元素對應當前體素的包圍圖像塊與模式中的採樣體素的包圍圖像塊之間的像素灰度差異累計，表示為式11：
f(vi)＝{dp(c(vi),c(vi+γ))|vi+γ∈p}(式11)
式11中，f(vi)表示當前體素vi的特徵描述子向量；c(vi)表示體素vi的包圍圖像塊，dp為圖像塊的灰度差異；γ為當前體素與模式p中採樣體素之間的位移；c(vi+γ)表示模式中的採樣體素的包圍圖像塊。其中，圖像塊的灰度差異dp定義為式12：
式12中，dp為圖像塊的灰度差異；c(vi)與c(vk)分別表示體素vi與vk的包圍圖像塊；δv為包圍體素塊中體素的位移；vi+δv屬於圖像塊c(vi)；i(vi+δv)為體素vi+δv的灰度；i(vk+δv)為體素vk+δv的灰度；
圖像塊的灰度差異可以看作圖像塊之間對應體素灰度差異的求和，該計算可以通過圖像的線性卷積算子實現。如果體素vi與模式中一個採樣體素之間的位移向量為γ，則可預先對體圖像進行平移操作得到v(γ)。圖像塊的灰度差異dp則可以對差異圖像v-v(γ)進行卷積得到，即通過式13得到：
dp＝(v-v(γ))2*g(式13)
式13中，dp為圖像塊的灰度差異；g為卷積核,v(γ)為對原始體圖像v進行平移γ後得到的體圖像。
結合體素的上下文特徵描述，圖結構中邊連接上的權重定義為上下文特徵描述差異與灰度差異的加權組合，如式14：
式14中，wij為體素i與相鄰體素j之間邊連接的權重；α為常數係數用於調節上下文特徵以及體素灰度對於權重計算的影響；ii與ij對應第i個與第j個體素的灰度，fi與fj對應第i個與第j個體素的上下文描述。ψ為圖結構的邊連接集合；本實施例中，系統預先將cbct圖像的灰度範圍變換到[0,255]區間，其中α被設置為0.1。通過式14計算得到體素的相似度矩陣。
步驟2：採用隨機遊走算法進行牙列的初始分割；
系統根據牙列分布情況定義標籤集合s，其中對不同的牙齒分配不同的標籤，正常的牙列包含32顆牙齒，其中第三磨牙，即智齒通常未發育完全，系統將其排除在外，還餘下28顆牙齒，因而加上背景類別，系統的標籤集合中將定義29類標籤。但在實際cbct數據中，由於牙齒缺失等問題，實際cbct圖像中包含的類別可能少於該類別個數。在初始牙列分割中，系統採用原始的隨機遊走算法對感興趣體圖像區域進行標記。voi中體素標記可以通過最小化如下的能量函數(式21)得到。
式21中，erw(x)表示基於隨機遊走算法的圖像分割能量；|s|為系統類別標籤的個數；表示第i個結點被分為第s類的概率；系統中每個結點均對應一個|s|維的概率向量，其中第i個結點所屬的類別定義為表示預先由用戶交互定義的結點概率，當第i個結點被用戶交互定義為第s類，則取值為1，否則為0；nl是用戶預先交互標註的體素個數；nv為voi中所有體素的個數；μ1與μ2為常數，本實驗中分別設置為1與0.001。上述目標函數(式21)可以轉化為矩陣形式，其中l為拉普拉斯矩陣，在給定步驟1)中通過公式14得到的相似度矩陣w後，拉普拉斯矩陣l＝d-w；其中，對角線矩陣d中元素定義為h為對角線指示矩陣，當voi中的第i個結點被預先標註時，h的元素hii＝1；x與表示voi中所有體素對應的類別概率以及所組成的矩陣。通過將上面的能量函數的一階導數設為0，該能量函數可以轉化為求解如下大規模稀疏的線性系統，表示為式22：
式22中，表示單位矩陣；l為拉普拉斯矩陣；h為對角線指示矩陣，當voi中的第i個結點被預先標註時，h的元素hii＝1；μ1與μ2為常數，本實驗中分別設置為1與0.001。給定用戶定義的形狀先驗可通過求解上面的線性系統(式22)得到voi中所有體素的類別，即得到牙列的初始分割結果。
步驟3：利用三維可變形模型，通過非剛性變形配準擬合步驟2進行圖像分割所得到的牙列表面體素集合，得到三維可變形模型的非剛性變換參數；
步驟2得到voi中所有體素的類別，利用其中屬於牙列表面的體素，可以定義牙列的表面形態；考慮到cbct圖像中可能存在的退化現象，以及處於咬合狀態上下頜接觸的牙齒，基於原始的隨機遊走技術獲取的初始分割與真實牙齒輪廓之間會存在較大的差異。系統引入三維可變形模型處理cbct圖像自身的局部退化問題所造成的邊界混淆問題。三維可變形模型來自對實體牙齒模型的三維光學掃描，其具有良好的拓撲結構定義。從cbct圖像分割可以定義牙列表面的體素集合，即屬於牙列的體素如果其近鄰體素屬於背景，則該體素被放入表面體素集yc。對三維可變形模型施加非剛性變換，以擬合初始分割得到的表面體素集。與體圖像分割結果擬合後的三維模型表面內部體素更有可能屬於牙列，並基於該假設定義柔化約束。
三維模型非剛性擬合分為兩個階段，首先利用procrustes分析計算三維模型的全局變換使其最大程度地與表面體素集yc一致；隨後利用非剛性配準獲取特定病人的三維模型，使其與初始分割的表面體素集一致。為了降低非剛性變換的參數空間，系統利用嵌入變形技術，將三維模型表面的非剛性變形定義為一個稀疏基網格頂點變換的線性組合，其中變形後的三維模型頂點為其中，權重βij對應頂點與其在基網格b上的近鄰頂點ηj，並基於兩點之間的歐式距離定義為其中參數tj為第j個基網格頂點上的三維變換矩陣，三維模型非剛性配準的能量函數定義如式31：
式31中，即變形前後表面頂點的位移；能量函數ere的第一項用於最小化變形後的三維模型表面與體圖像分割獲取的表面體素集之間的距離；dc表示變換後的頂點與標籤擴散得到的表面體素集yc之間的距離,其中為表面體素；第二項用於保持三維模型表面平滑，該項最小化頂點與其近鄰的空間變換之間的差異。ω為常數係數，ne為三維模型表面頂點個數。通過最小化能量函數ere可以得到三維模型的非剛性變換參數。
步驟4：基於柔化約束下的隨機遊走算法進行牙列分割，得到voi中所有體素的類別標籤，即哪些體素屬於牙列(前景)，哪些體素屬於背景；
可以直觀看到，與體圖像分割得到的表面體素集配準後的三維模型表面內部的體素更有可能屬於牙列。但是，在體圖像中牙齒表面內部體素並不具有均勻的灰度，其中牙本質具有較亮的灰度，而牙髓腔灰度較低，考慮到牙齒表面內部的灰度分布，系統對牙本質賦以更大的屬於牙齒的概率。由於三維模型來自對牙齒模型的光學掃描，該模型不包括內部的牙髓腔的表面，利用經驗的牙本質厚度定義a，以及組合logistic函數依據步驟3得到的變形後的三維模型表面定義柔化約束。極坐標系的中心定義為變形後三維模型與切片圖像相交輪廓的中心，令r1(θ)與r2(θ)表示牙齒的外輪廓以及牙髓腔輪廓，其中r2(θ)＝r1(θ)-a。基於三維可變形模型屬於前景牙齒的概率pd定義為式41：
式41中，pd(r,θ)為基於三維可變形模型屬於前景牙齒的概率；r,θ為體素在極坐標系中的坐標；r1(θ)與r2(θ)表示牙齒的外輪廓以及牙髓腔輪廓，其中r2(θ)＝r1(θ)-a，a為經驗的牙本質厚度；參數η用於控制如上概率函數經過牙齒輪廓時函數的形狀。在實驗中，參數η的值設為0.99。
如果考慮沒有牙髓腔的分層圖像，基於三維可變形模型屬於前景牙齒的概率定義為式411：
其中，pd(r,θ)為基於三維可變形模型屬於前景牙齒的概率；r,θ為體素在極坐標系中的坐標；r1(θ)表示牙齒的外輪廓；其中參數η用於控制如上概率函數經過牙齒輪廓時函數的形狀。在實驗中，參數η的值設為0.99。
上述的基於三維可變形模型屬於前景牙齒的概率針對切片圖像中的單輪廓定義，考慮到牙齒可能出現的多牙根以及多牙尖結構，該定義可以直接擴展到牙齒在切片上的多輪廓狀態，表示為式412：
其中，nc為當前牙齒在切片圖像中輪廓的個數；pd(r,θ)為式41與式411定義的在坐標(r,θ)處基於三維可變形模型屬於前景牙齒的概率；為擴展到牙齒在切片上的多輪廓狀態時基於三維可變形模型屬於前景牙齒的概率。
由於三維模型變形後與牙齒真實表面之間可能存在空隙，系統還利用體素的表觀定義體素屬於牙列的分類器。從交互標註的切片圖像中學習該分類器，利用卡方(chi-square)核的一對多(one-vs-rest)支持向量機(svm)生成分類器，該分類器的輸出被歸一化到0-1之間，表示體素屬於前景(牙列)的概率pc。柔化約束定義為式413：
式413中，pc表示由體素表觀估計的體素vi屬於前景(牙列)的概率；表示基於三維可變形模型vi屬於前景牙齒的概率；r,θ為體素vi在極坐標系中的坐標。是對體素vi相對於類別s所定義的柔化約束。
柔化約束的能量項定義為式42：
其中，表示第i個結點被分為第s類的概率；x表示voi中所有體素對應的類別概率所組成的矩陣；|s|表示系統中類別個數；nv為感興趣區域中體素的個數；為式413所定義的柔化約束；對應感興趣區域中所有體素的柔化約束所構成的矩陣；給定柔化約束，組合式21與式42，將感興趣體圖像的標籤擴散能量定義為式43：
erwsc(x)＝λ1erw+λ2esc(式43)
其中erw與esc式21式42所定義的隨機遊走以及柔化約束的能量。λ1與λ2為常數係數。將該能量函數(式43)的一階導數設為0，可得到關於體素結點概率向量的線性系統，利用共軛梯度方法求解如下線性系統(式44)。
其中l為拉普拉斯矩陣；h為對角線指示矩陣，當voi中的第i個結點被預先標註時，h的元素hii＝1；x與表示voi中所有體素對應的類別概率以及用戶定義的先驗概率所組成的矩陣。對應感興趣區域中所有體素的柔化約束所構成的矩陣；λ1與λ2為常數係數,實驗中分別設置為1、0.1；μ1與μ2為常數，本實驗中分別設置為1與0.001；表示單位矩陣。由於xi是|s|維的向量，其中|s|為系統中類別標籤的個數，該線性系統需要求解|s|-1次，其中求解線性系統可以得到voi中所有體素的類別標籤。
步驟5：對牙列分割迭代修正；
為了改進單次柔化約束下隨機遊走算法所得到的體圖像分割，系統引入迭代修正過程，對柔化約束下標記擴散，以及三維模型對表面體素集的擬合問題，進行迭代求解，以改進單次標籤擴散所獲取的牙列分割。
將柔化約束下的標籤擴散以及三維模型的非剛性配準能量結合，定義能量如下：
e(x,t)＝λ1erw+λ2esc+λ3ere(式51)
式51中，e(x,t)是基於三維可變形模型隨機遊走算法進行牙列分割的總能量函數；erw與esc為式21式42所定義的隨機遊走以及柔化約束的能量函數；ere為式31所定義三維模型非剛性配準的能量；λ1，λ2與λ3為常數，實驗中分別設置為1、0.1、1。
該能量函數的優化可以分解為兩個子問題。在第一個階段，給定基於標籤擴散得到的體圖像分割定義牙列表面體素集，對三維模型進行變形，用於擬合更新後的表面體素集。該變形通過最小化能量e(x(i-1),t)得到，其中表面體素集由第i-1次迭代中的cbct圖像分割得到，因而與柔化約束下的隨機遊走相關的項，即erw與esc可以從能量函數中去掉。對於三維模型變形參數的求解與步驟3)的描述一致。
在第二個階段，最小化能量e(x,t(i-1))以得到voi中所有體素的標籤。由於在標籤擴散過程中使用的柔化約束來自第i-1步迭代中的三維模型配準，因而能量函數中與三維模型配準相關的項，即ere可被去掉，該階段對應步驟4)中柔化約束下的隨機遊走算法，通過求解一個大規模稀疏的線性系統求解所有體素的標籤。
隨著兩個子問題的優化迭代進行，變形後的三維物體表面與基於標籤擴散所得到的表面體素集會越來越接近。當相鄰兩步迭代中三維模型更新||ye(i)-ye(i-1)||2(其中ye(i)與ye(i-1)分別為第i步與第i-1步三維表面模型)小於預先給定的閾值時，終止迭代，得到voi中所有體素的標籤，作為迭代修正後的牙列分割結果。
為了驗證基於三維可變形模型隨機遊走算法對牙列進行分割的精度，利用dice測度計算自動分割與手工分割的相似度，該相似度大於0.95。同時，度量自動分割所得到的牙列表面模型與手工分割對應的表面模型的平均距離，該距離誤差在0.3mm以下，可滿足口腔臨床對於精度的要求。
利用本發明的方法，可以對臨床採集cbct圖像進行牙列分割。該方法結合半監督的隨機遊走算法，以及由三維可變形模型配準定義的柔化約束獲取牙列的可靠分割。其中利用三維可變形模型引入柔化約束，也可以看作為正則化項用於處理基於半監督標籤擴散的體圖像分割中的噪聲。該發明引入迭代修正方法，其中對系統的兩個子問題，即柔化約束下標記擴散，以及三維模型對表面體素集的擬合問題，進行迭代求解以改進單次標籤擴散所獲取的牙列分割。
需要注意的是，公布實施例的目的在於幫助進一步理解本發明，但是本領域的技術人員可以理解：在不脫離本發明及所附權利要求的精神和範圍內，各種替換和修改都是可能的。因此，本發明不應局限於實施例所公開的內容，本發明要求保護的範圍以權利要求書界定的範圍為準。