基于網格抽取的牙齒模型快速分割算法

馬 勇 柯永振 楊 帥

(天津工業大學計算機科學與軟件學院 天津 300387)

0 引 言

隨著社會的快速發展，人們的生活水平得到了很大的提高。與此同時，很多人開始重視自身的外表，而牙齒作為外表的一個重要組成部分，直接影響著個人的儀容和生活質量。由于一般的金屬牙套存在形狀突出、難清洗等缺點，所以近年來興起的隱形牙套受到了人們的廣泛歡迎。制作隱形牙套主要有牙齒數據獲取以及三維重建、牙齒分割、繪制牙弓線、虛擬排牙以及隱形牙套生成等過程，而牙齒分割是非常重要的一個過程。

近十年，在牙齒模型分割領域中，專家和學者提出了許多牙齒分割的算法。文獻[1]針對牙齦上的輪廓線，一種基于能量方程的牙齒分割算法被提出，通過求解能量方程，在條件收斂之前，該算法對牙齦上的輪廓線不斷地進行迭代并進行移動，這種算法雖然沒有進行人工交互，但是對于很多排列不整齊的牙齒并不能進行精確的分割。文獻[2]采用洪泛法提取牙齒上的信息，通過這些信息識別牙齒之間的凹區域，即牙齒間的邊界，該算法對于排列比較整齊的牙齒有較精確分割結果，但對于排列不規則的牙齒沒有良好的分割效果，而且分割效率非常低。文獻[3]提出了交互標記控制算法，該算法交互選取種子面片，然后根據相鄰兩個面的曲率進行區域生長，該算法分割速度較慢，而且存在過度分割現象。文獻[4]提出一種基于改進后的交互標記控制算法，該方法在分割過程中采用人為的設定閾值作為控制參數來控制隊列的操作，這種算法使牙齒的分割速度得到了提高，但對于一些不規則牙齒并不能進行精確的分割。文獻[5]提出一種基于調和場的牙齒分割交互算法，該算法在分割過程中采用調和場的屬性提取等值線，進而提取牙齒邊界，然后再交互進行牙齒分割，這種方法提高了分割的精度，但對于一些不規則牙齒并不能得到良好的分割效果。文獻[6]通過計算牙齒三角網格的曲率以及進行一系列的形態學的開閉操作，提出了一種對于三角網格進行形態學操作的分割算法，這種算法對牙齒可以進行比較精確的分割，但是對于牙齒之間的融合區域，并沒有提出有效的解決方法，而且分割的速度較慢。文獻[7]針對牙齒之間的融合區域，找到了一種解決的方法，該方法首先剔除牙齒之間的融合區域然后再對沒有融合區域的牙齒重新建模，那樣做破壞了原有的牙齒網格模型上的拓撲結構。文獻[8]對提取的牙齒特征區域進行膨脹和腐蝕操作，然后利用形態學開閉操作對牙齒的特征區域進行細化，將牙齒分離，該方法對排列規則的牙齒有良好的分割結果，而對排列不規則的牙齒，分割結果較差，而且分割效率較低。文獻[9]提出了一種基于調和場的自動分割算法，該算法能有效的一次性獲取牙齒網格模型上的所有牙齒區域，但該算法并不能進行單個牙齒分離。

本文在交互標記控制算法的基礎上，提出了一種基于網格抽取的牙齒模型快速分割算法。該算法主要改進了三個方面：第一，進行網格抽取，可以適當減少網格的數量，加快網格的分割速度；第二，進行網格平滑，消除網格中的噪聲，提高網格分割的精度；第三，改變數據結構，取消排序算法，減少算法的運行時間。

1 交互標記控制算法

交互標記控制算法[4]是在分水嶺算法[10]的基礎上提出來的。該算法遵循Hoffman等[11]提出的極小值法則，并按照Fast Marching Watershed[12]理論進行區域擴展從而確定分割區域的位置和范圍[13-14]。

該算法的思想描述如下：

交互標記控制算法首先在每個牙齒上人工拾取一個三角網格，并對其進行標記。接著利用堆棧作為存放三角網格的數據結構，對堆棧進行初始化操作，即把每個牙齒上拾取并被標記的三角網格的相鄰的三角網格全部壓入堆棧中。然后把堆棧里的三角網格按照相對彎曲程度進行從大到小的排序，每次從堆棧中取出位于棧頂的三角網格，即相對彎曲程度最大的三角網格，將該網格進行標記，并且將與該面片相鄰的且沒有被標記的三角面片按照相對彎曲程序的大小順序進入堆棧中。此過程循環進行，直到堆棧為空。至此，牙齒上所有的三角網格都被標記了一次，通過標記，可以確定每個網格屬于某個區域，即屬于某個牙齒上。

該算法的具體流程為：

(1) 堆棧初始化，即把所有被人工交互標記的三角網格的相鄰網格按照曲率的大小壓入堆棧中。

(2) 區域生長，即每次從堆棧中取出位于棧頂的三角網格，記為，則將網格face′標記為flag′，并將與該網格face′相鄰的并且沒有被標記的網格gi、flag′和相應的曲率H(face′,gi)構成的結構體按曲率大小排序，然后順序壓入堆棧中。

2 基于網格抽取的牙齒模型快速分割算法

由于之前的交互標記控制算法分割時間長且存在欠分割或過分割的現象，為了提高算法的分割速度和精度，本文提出了一種運行速度較快，分割精度較高的牙齒模型快速分割算法。

基于網格抽取的牙齒模型快速分割算法思想描述如下：

(1) 網格抽取 對于輸入的牙齒三維網格模型，進行網格抽取，適當縮減網格的數量。

(2) 網格平滑 利用拉普拉斯平滑算法對三角網格進行平滑處理，以消除網格的噪聲。

(3) vector容器初始化 交互操作，使已經標記的三角網格相鄰的網格進入vector容器中。

(4) 區域生長 訪問牙齒網格所有面片，根據極小值法則[11]，對網格面片進行分組，進而實現牙齒分離。

基于網格抽取的牙齒模型快速分割算法的流程圖如圖1所示。

圖1 基于網格抽取的牙齒模型快速分割算法流程圖

2.1 網格抽取

一般經過激光掃描得到的STL格式的牙齒模型的網格數量非常多，為了能夠加快算法的運行速度，本文采用了文獻[15]提出的網格縮減算法進行牙齒模型網格抽取，即適當減少牙齒模型網格的數量。多次實驗表明，網格抽取10%，即網格數量減少10%時，使用本文提出的算法，既不會影響牙齒網格分割的精度，又提高了網格的分割速度。該算法首先將網格G=中的每個頂點Vi分類，并插入到優先隊列S1中。優先級隊列是基于誤差來刪除頂點并重新三角化，無法刪除或三角化的頂點將被跳過。接著處理優先級隊列S中每個頂點Vi(即使用邊折疊來刪除頂點，然后進行三角剖分)，繼續迭代，直到優先級隊列S1為空。最后處理所有剩余的頂點，沿著銳利的邊緣將網格劃分成單獨的片段，并且將其插入到優先級隊列S2中。再次處理優先級隊列S2，直到隊列S2為空。算法的偽代碼如下：

網格抽取算法偽代碼：

Input: G

//輸入網格的幾何和拓撲信息

classify Viand Push S1

//識別網格頂點并放入隊列S1

While(S1!=empty) do

If(delete Viand re-triangulated )

//刪除不符合條件的頂點并且重新三角化

Vitriangulated

End if

Else ViPush S2

//不能被刪除的頂點進入隊列S2

End while

While(S2!=empty)

Vitriangulated

//剩余的頂點三角化

End while

Return G1

2.2 網格平滑

經過激光掃描得到的STL文件，往往存在著噪聲，為了消除噪聲，提高分割的精度，本文把拉普拉斯平滑算法用在了對牙齒三角網格的平滑處理上。該算法首先找到網格G=的每個頂點Vi的1鄰域頂點集合，其1鄰域定義為與點Vi所有共邊的點的合集：

L1(Vi)={Vi|?edge(Vi,Vj)}

(1)

然后為每個頂點Vi構建1鄰域數組Array[n]，n是每個頂點Vi的1鄰域頂點的個數。最后迭代所有頂點，對于每個頂點Vi，根據1鄰域數組Array[n]里所有頂點坐標的平均值，即采用傘狀算子：

(2)

來修改Vi的坐標。重復多次迭代(根據多次實驗的結果，迭代次數為20次時，效果最好)。算法的偽代碼如下：

網格平滑算法偽代碼：

Input: G

//輸入網格的幾何和拓撲信息

While(count < 20) do

//迭代20次平滑

For i < vertex number of G do

n = L1(Vi)

//每個頂點執行方程(2)

new array[n]

modify coordinate of Viby U(Vi)

//根據方程(3)修改每個頂點的坐標

End for

End while

Return G1

2.3 替換堆棧結構

傳統的交互標記控制算法，利用堆棧存儲牙齒模型的面片信息，并根據直接插入排序算法對堆棧里的三角面片進行排序。鑒于堆棧元素后進先出的特點，本文采用vector容器代替堆棧存儲牙齒模型的面片信息。由于直接插入排序算法比較耗時，本文不再使用排序算法進行排序，充分利用vector容器自身的特點，可以直接在容器內找到曲率最大的元素，這樣節省了大量時間。

2.4 vector容器初始化

在每個牙齒上人工拾取一個三角網格并對其進行標記f，找出每個拾取的三角網格f相鄰的三個網格fi，把每個網格f的ID，標記m，以及曲率c(f，fi)存入對應的結構體數組中，并放入vector容器中；而曲率c(f，fi)的計算方法如下：

設兩個三角網格f1、f2，它們的單位法向量分別為n1和n2，這里設f1、f2的共享邊界為AC，AB為其中一個三角面片的非公共邊(如圖2所示)，則標記的兩個三角網格f1,f2，定義式(3)作為它們的相對彎曲程度的函數。

H(f1,f2)=H(f2,f1)=

(3)

圖2 相鄰兩個網格之間相對彎曲程度的計算

2.5 區域生長算法

如果vector容器不為空，每次從vector中取出曲率最大的三角面片fmax，判斷fmax的訪問標記Visited是否為真，若Visited為真，將fmax直接從vector容器中刪除，繼續迭代；若Visited為假，設置fmax的訪問標記Visited為TRUE，將fmax插入到對應的數組Array[m]中(每個牙齒對應著一個數組，存放訪問標記Visited為真的三角面片信息)，把fmax相鄰的且訪問標記Visited為假的面片fi放入vector容器中，繼續迭代。其偽代碼如下所示：

區域生長算法的偽代碼：

Input:G

//輸入網格的幾何和拓撲信息

While(vector!=empty)

Select fmax

//選擇曲率最大的三角面片

If(Visited==true)

Delete fmaxfrom vector

//刪除已經被標記的三角面片

End if

ElseVisited=true

//設置訪問標記為真

fmaxPop Array[m]

Neighboring face of fmaxand Visited is false Pop vector

//把相鄰的且訪問標記為假的面片放入vector中

End while

Return Array[m]

3 算法實驗及分析

為驗證本文提出的基于網格抽取的牙齒模型快速分割算法的精確性和快速性，利用激光掃描得到的STL文件格式的牙齒模型進行實驗，所有實驗均在處理器為Intel(R) Core(TM) i5-3470 CPU @ 3.20 GHz，內存4 GB，64位Windows7系統的電腦上運行的。開發軟件為Visual Studio 2012和Qt5.5.1，編程語言為C++，所用的工具包為VTK可視化工具包。

圖3是兩個激光掃描得到的STL文件格式的牙齒網格模型(這兩個網格模型均來自天津某醫院)。模型一的三角網格數量為365 092，頂點數量為182 544，大小為18 290 KB，體積為49 464.945 m3，面積為11 487.523 m2；模型二的三角網格數量為506 702，頂點數量為253 353，大小為24 742 KB，體積為78 645.768 m3，面積為20 356.143 m2。圖4是兩個STL文件格式的網格牙齒模型經過交互標記控制算法分割后，不同角度的分割結果圖。圖5是兩個STL文件格式的牙齒網格模型經過本文的算法分割后，不同角度的分割結果圖。從圖中可以看出，交互標記控制算法分割牙齒存在明顯的欠分割(單顆牙齒沒有被完全分割出來)現象，特別是在雙尖牙處(如圖4圓圈圍住的牙齒)。而本文的算法分割牙齒，基本上解決了欠分割現象的存在(如圖5圓圈圍住的牙齒)。其主要原因是交互標記控制算法的區域生長是根據直接插入排序算法從堆棧中每次取曲率最大的三角面片，而雙尖牙上的凹區域比較明顯，直接插入排序算法則會把雙尖牙上凹區域的三角面片排到堆棧的底部，這樣區域生長時無法取出該區域的三角面片，所以極容易產生欠分割的現象。而本文的算法利用vector代替堆棧結構，并且沒有利用直接插入排序算法，而是利用了vector容器自帶的取最大值的方法，這樣區域生長時可以非常準確地取出雙尖牙凹區域的三角面片，所以本文的算法解決了交互標記控制算法存在的欠分割的現象。

圖3 兩個牙頜模型

圖4 交互標記控制算法分割結果

圖5 本文算法分割結果

表1列出了兩個牙齒模型在兩種算法中的分割時間，可以看出基于網格抽取的牙齒網格快速分割算法的分割時間有了明顯的提高，和原來的交互標記控制算法相比，分割時間的加速比提高了6到9倍左右，而且牙齒模型的三角面片數量越少，分割的時間越快。其主要原因是：第一，本文提出的算法包含了網格抽取算法，在牙齒網格分割前，進行了網格抽取，適當地減少了牙齒網格的數量，從而加快了網格的分割速度；第二，交互標記控制算法中包含了一個插入排序算法，而插入排序算法的時間復雜度是O(n2)，所以交互標記控制算法的時間復雜度是O(n2)，而本文提出的算法使用了vector容器自帶的尋找最大值的方法，沒有使用插入排序算法，因為vector容器自帶的尋找最大值的方法的時間復雜度為O(n)，所以本文提出的算法的時間復雜度是O(n)，因此本文提出的算法的分割時間明顯快于交互標記控制算法的分割時間。從實驗中也可以看出，網格模型的頂點數越多，體積和面積越大，算法執行的速度越慢，但分割的精度越高，對后續操作的影響也越小。

表1 兩種算法分割時間比較

4 結 語

本文提出了一種基于網格抽取的牙齒模型快速分割算法。多次實驗的結果表明，該算法主要包含三個優點：1) 分割速度快。和傳統的交互標記控制算法相比較，本文提出的算法牙齒分割速度有了明顯的提高。2) 分割精度高。傳統的交互標記算法存在欠分割，本文提出的算法一般不會出現欠分割。3) 適應性強。對于大部分的三維牙齒網格模型，本文提出的算法都可以進行精確的分割。同時，本文的算法除了適用于牙齒網格分割外，對于其他網格模型的分割，也可以用本文提出的算法。當然，本文提出的算法也有一定的局限性，對于一些牙齒網格模型，體積和面積太大，本文的算法不能進行精確的分割，而且對于網格數量超過80萬的三維牙齒模型，分割速度也有點慢。

參 考 文 獻

[1] Kronfeld T,Brunner D,Brunnett G.Snake-Based Segmentation of Teeth from Virtual Dental Casts[J].Computer-Aided Design and Applications,2010,7(2):221-233.

[2] Kumar Y,Janardan R,Larson B,et al.Improved Segmentation of Teeth in Dental Models[J].Computer-Aided Design and Applications,2013,8(2):211-224.

[3] 李成軍,張弛,汪國平.交互標記控制的快速網格分割[J].北京大學學報(自然科學版),2006,42(5):662-667.

[4] 吳松和.牙齒隱形矯正CAD軟件開發[D].石家莊:河北科技大學,2013.

[5] Li Z,Wang H.Interactive Tooth Separation from Dental Model Using Segmentation Field[J].Plos One,2016,11(8):e0161159.

[6] Wu K,Chen L,Li J,et al.Special Section on CAD/Graphics 2013:Tooth segmentation on dental meshes using morphologic skeleton[J].Computers & Graphics,2014,38(1):199-211.

[7] Yuan T,Liao W,Dai N,et al.Single-tooth modeling for 3D dental model[J].International Journal of Biomedical Imaging,2010,2010(9):1029-1034.

[8] Mouritsen D A.Automatic segmentation of teeth in digital dental models[D].Birmingham: the university of alabama at birmingham,2013.

[9] Liao S H,Liu S J,Zou B J,et al.Automatic Tooth Segmentation of Dental Mesh Based on Harmonic Fields[J].Biomed Research International,2015,2015(3):1-10.

[10] Mangan A P,Whitaker R T.Partitioning 3D surface meshes using watershed segmentation[J].IEEE Educational Activities Department,1999,5(4):308-321.

[11] Hoffman D D,Singh M.Salience of visual parts[J].Cognition,1997,63(1):29-78.

[12] Koschan A F.Perception-based 3D triangle mesh segmentation using fast marching watersheds[C]//Computer Vision and Pattern Recognition,2003.Proceedings.2003 IEEE Computer Society Conference on.IEEE,2003,2:27-32.

[13] Liu R,Zhang H.Segmentation of 3D Meshes through Spectral Clustering[C]//Computer Graphics and Applications,Pacific Conference.IEEE Computer Society,2004:298-305.

[14] Katz S,Tal A.Hierarchical mesh decomposition using fuzzy clustering and cuts[J].Acm Transactions on Graphics,2003,22(3):954-961.

[15] Schroeder W J.Decimation of triangle meshes[J].Acm Siggraph Computer Graphics,1992,26(2):65-70.