基于形態(tài)學(xué)方法的牙頜模型自動(dòng)分割算法

陳 虎 馬珂楠 胡畢富

1(北京大學(xué)口腔醫(yī)院口腔醫(yī)學(xué)數(shù)字化研究中心 北京 100081)

2(北京航空航天大學(xué)機(jī)械工程及自動(dòng)化學(xué)院飛行器制造工程系 北京 100191)

0 引 言

在口腔醫(yī)學(xué)領(lǐng)域中,CAD/CAM技術(shù)的引入正逐漸改變傳統(tǒng)口腔治療方式,而牙科CAD系統(tǒng)作為口腔數(shù)字化醫(yī)療中的重要環(huán)節(jié),對(duì)其設(shè)計(jì)策略和算法的研究都具有重要的意義[1]?？谇会t(yī)學(xué)包括口腔修復(fù)、正畸、種植以及頜面修復(fù)等主要方面,每種不同治療方式需要構(gòu)建的CAD系統(tǒng)差異性很大,但很多情況下都需要從牙頜的數(shù)字化模型分離開(kāi)牙齒和牙齦部分,實(shí)現(xiàn)對(duì)牙頜模型的分割,例如口腔正畸中的虛擬排牙等。受不同牙頜的形狀差異性,以及牙齒畸形、牙列缺損等因素影響,目前的牙科CAD軟件大多是以繁瑣、低效的手動(dòng)交互方式完成分割的,對(duì)牙頜模型分割的自動(dòng)化程度較低。

牙頜模型的分割屬于數(shù)字幾何處理領(lǐng)域中的三維網(wǎng)格分割問(wèn)題,近年來(lái)三維網(wǎng)格分割已成為一個(gè)研究熱點(diǎn)[2],但已有的網(wǎng)格模型分割方法不能實(shí)現(xiàn)對(duì)任意形狀的模型的自動(dòng)分割。目前針對(duì)具有復(fù)雜形狀的牙頜模型,已提出了較多的分割方法。文獻(xiàn)[3-4]從全景圖像中構(gòu)建了兩幅距離圖像實(shí)現(xiàn)了牙齒的自動(dòng)分割,避免了直接處理復(fù)雜的三維網(wǎng)格數(shù)據(jù)問(wèn)題,但也使得這種分割方式不適合于處理三維牙頜模型的牙科CAD系統(tǒng)。目前比較常見(jiàn)的方法是通過(guò)分析牙頜模型的離散曲率實(shí)現(xiàn)對(duì)牙齒和牙齦的分離。Kronfeld等[5]利用牙齒和牙齦邊界附近頂點(diǎn)的曲率值較大的特點(diǎn),構(gòu)建了活動(dòng)輪廓線,之后通過(guò)能量法局部?jī)?yōu)化邊界線實(shí)現(xiàn)牙頜分割。但該方法不易處理活動(dòng)輪廓線出現(xiàn)不連續(xù)的牙頜模型。Kumar等[6]通過(guò)設(shè)置曲率閾值選取出高曲率值的頂點(diǎn)獲得齦緣線輪廓,在分離開(kāi)牙齦和牙齒后,將齦緣線上的尖點(diǎn)作為相鄰牙的分割邊界點(diǎn)分離開(kāi)相鄰牙齒。該方法也難以處理牙齒畸形等情形,比較依賴于標(biāo)準(zhǔn)理想的牙頜模型。為了提高分割的質(zhì)量,在曲率獲取齦緣線的基礎(chǔ)上,利用形態(tài)學(xué)操作獲取更精細(xì)的牙齒邊界線。Zhao等[7]在加入較多交互的方式下,利用形態(tài)學(xué)操作實(shí)現(xiàn)了牙頜模型的分割。Yuan等[8]通過(guò)形態(tài)學(xué)操作分離開(kāi)牙齒和牙齦曲面,并實(shí)現(xiàn)了對(duì)相鄰牙齒的分離以及對(duì)牙齒曲面的修補(bǔ),該方法同樣需要較多的手動(dòng)交互操作。Wu等[9]在形態(tài)學(xué)操作分割牙頜模型的基礎(chǔ)上,實(shí)現(xiàn)了少量交互下對(duì)牙齒牙齦以及相鄰牙齒的分割。除了利用牙頜網(wǎng)格模型的曲率特性實(shí)現(xiàn)分割外,還有其他反映網(wǎng)格曲面特征的信息可用于牙頜模型的分割。Zou等[10]通過(guò)計(jì)算網(wǎng)格上的調(diào)和場(chǎng)來(lái)反映網(wǎng)格的幾何細(xì)節(jié),再手動(dòng)選取具有約束作用的特征點(diǎn),實(shí)現(xiàn)牙頜模型的分割,在之后的改進(jìn)方法中[11],通過(guò)對(duì)牙頜模型的多層切分實(shí)現(xiàn)了自動(dòng)選擇特征點(diǎn),從而完成對(duì)牙頜模型的自動(dòng)分割,但該方法在處理相鄰牙齒貼靠緊密的情況時(shí)難以將牙齒分割開(kāi)。Hwang等[12]也基于調(diào)和場(chǎng)進(jìn)行牙頜模型的分割,并采用狄利克雷邊界條件,凸分割等后處理方法,在具有噪聲和網(wǎng)格碎片的臨床實(shí)際的牙頜模型上取得了一定的效果。Fan等[13]采用離散Morse理論分析曲面結(jié)構(gòu),把分割問(wèn)題轉(zhuǎn)化成從曲面上冗余特征線族中提取分割邊界線,再通過(guò)優(yōu)化方法從底層分割中聚合產(chǎn)生子部分,獲得了較好的分割效果。Kim等[14]通過(guò)對(duì)牙頜模型進(jìn)行截面投影,隨后應(yīng)用區(qū)域增長(zhǎng)算法實(shí)現(xiàn)牙齒的自動(dòng)分割。He等[15]則基于曲率信息和局部分布密度信息初步劃分牙齒-牙齦邊界和牙齒間邊界,然后通過(guò)區(qū)域增長(zhǎng)算法實(shí)現(xiàn)牙齒分割。Pupykina等[16]利用可視化工具包進(jìn)行交互,基于用戶選擇的少量關(guān)鍵特征點(diǎn),通過(guò)模糊綜合評(píng)判的方法確定牙齒的邊界。近年來(lái)也出現(xiàn)了許多基于機(jī)器學(xué)習(xí)的方法。Xu等[17]建立了牙齒分割的深度卷積神經(jīng)網(wǎng)絡(luò),將牙頜網(wǎng)格面片的特征作為輸入,經(jīng)過(guò)兩層神經(jīng)網(wǎng)絡(luò),第一層用于區(qū)分牙齒和牙齦區(qū)域,第二層用于區(qū)分不同牙齒區(qū)域,最終通過(guò)模糊聚類對(duì)結(jié)果進(jìn)行優(yōu)化,取得了較好的分割效果。Lian等[18]建立了MeshSegNet,實(shí)現(xiàn)了牙頜模型多尺度的特征提取和學(xué)習(xí),采用圖約束(Graph-Constrained)學(xué)習(xí)模塊進(jìn)行局部特征提取,并通過(guò)密集融合(dense fusion)策略進(jìn)行更高層次的特征學(xué)習(xí),最后采用圖割算法對(duì)分類結(jié)果進(jìn)行優(yōu)化調(diào)整。Sun等[19]基于圖卷積神經(jīng)網(wǎng)絡(luò)實(shí)現(xiàn)了牙齒和牙齦的同時(shí)分割。Zhang等[20]將牙頜三維模型映射到二維的調(diào)和參數(shù)空間,轉(zhuǎn)換為圖像信息,從而利用卷積神經(jīng)網(wǎng)絡(luò)進(jìn)行牙齒分割,隨后將二維圖像上的分割結(jié)果映射回三維實(shí)體模型中,并通過(guò)改進(jìn)的模糊聚類和分割方法完成邊界的調(diào)整和優(yōu)化。這些基于機(jī)器學(xué)習(xí)的方法雖然在一些數(shù)據(jù)集上取得了較高的準(zhǔn)確率,但是受訓(xùn)練數(shù)據(jù)規(guī)模、質(zhì)量等影響較大,在處理數(shù)據(jù)質(zhì)量較差的臨床實(shí)際模型時(shí)可能無(wú)法達(dá)到理想的效果。另一方面,由于訓(xùn)練數(shù)據(jù)通常需要由專業(yè)人員進(jìn)行大規(guī)模和高質(zhì)量的標(biāo)注,成本高昂,訓(xùn)練和測(cè)試時(shí)對(duì)機(jī)器運(yùn)算能力也有非常高的要求。

現(xiàn)有的牙科CAD系統(tǒng)在處理牙頜模型的分割時(shí),都是以手動(dòng)交互的方式分出牙齒再進(jìn)行其他處理,頻繁地交互操作難以提高CAD系統(tǒng)的使用效果。但受到牙頜的形狀差異性的影響,以及牙齒排列畸形、牙列缺損等口腔病癥的出現(xiàn),使得牙頜模型的自動(dòng)分割一直沒(méi)有好的解決方案,還沒(méi)有通用的方法,已有方法都或多或少存在交互和準(zhǔn)確性問(wèn)題。

本文在現(xiàn)有形態(tài)學(xué)方法處理牙頜模型分割的基礎(chǔ)上[21],提出了一種自動(dòng)閉合牙頜模型輪廓特征線的方法。牙頜三角網(wǎng)格模型在經(jīng)過(guò)形態(tài)學(xué)處理得到單個(gè)頂點(diǎn)寬度的特征線后,有時(shí)存在不連續(xù)的情形。首先通過(guò)計(jì)算特征頂點(diǎn)的復(fù)雜度[8,21]判斷出特征線的端點(diǎn),之后通過(guò)距離約束需要配對(duì)連接的端點(diǎn)以及無(wú)效的端點(diǎn),再逐層分別搜索端點(diǎn)周圍的非特征頂點(diǎn)并標(biāo)記,直到兩配對(duì)的端點(diǎn)都搜索到同一非特征頂點(diǎn),以該頂點(diǎn)作為兩個(gè)配對(duì)端點(diǎn)的中點(diǎn),根據(jù)最短路徑原則分別與兩端點(diǎn)相連,從而使得特征線封閉。該方法解決了原來(lái)需要手動(dòng)處理特征線不封閉的問(wèn)題,實(shí)現(xiàn)了對(duì)牙頜模型的自動(dòng)分割。

1 網(wǎng)格曲面上的形態(tài)學(xué)操作

1.1 三角網(wǎng)格的描述

對(duì)于三角網(wǎng)格模型M通?？杀硎緸镸={V,T},其中V是網(wǎng)格頂點(diǎn)的集合,T是網(wǎng)格上三角面片面的集合。

(1)

三角網(wǎng)格拓?fù)湫畔⒌慕⒅饕紤]網(wǎng)格M的頂點(diǎn)V、邊E和三角面片T的信息以及相互之間的遍歷和訪問(wèn)。

牙頜模型曲面的特征提取主要可歸結(jié)為三角網(wǎng)格分割問(wèn)題,特征線一般是分割的邊界線。實(shí)現(xiàn)對(duì)牙頜模型的分割關(guān)鍵在于齦緣線的提取,在牙頜模型中,齦緣線及相鄰牙齒之間的特征線附近的頂點(diǎn)平均曲率值比其他區(qū)域偏大,可以通過(guò)平均曲率值初步區(qū)分出齦緣線附近的特征區(qū)域。為了計(jì)算網(wǎng)格頂點(diǎn)上的離散曲率,對(duì)于每個(gè)頂點(diǎn)v∈V的主曲率可分別表示為κ1(v)和κ2(v),從而平均曲率表示為：

(2)

在牙頜模型的特征線附近都是凹下去的“谷線”區(qū)域,這些區(qū)域頂點(diǎn)的平均曲率值都是偏小的負(fù)值,通常直接設(shè)置曲率閾值ε來(lái)獲得初始的特征點(diǎn)集F：

F={v∈V|κH(v)≤ε}

(3)

根據(jù)曲率閾值初步篩選出的特征點(diǎn)集有較多的噪聲點(diǎn),而且能夠表示齦緣線和相鄰牙齒分割線的特征點(diǎn)也存在缺陷,需要做后續(xù)處理,才能獲得適用于牙頜模型分割的只有一個(gè)頂點(diǎn)寬度的特征線。

1.2 形態(tài)學(xué)操作

形態(tài)學(xué)方法可以處理初始特征點(diǎn)集并獲得只有一個(gè)頂點(diǎn)寬度可用于模型分割的特征線。R?ssl等[21]最早將形態(tài)學(xué)操作應(yīng)用于三角網(wǎng)格曲面的特征線提取,并將用于圖像處理的形態(tài)學(xué)操作轉(zhuǎn)化到處理三角網(wǎng)格曲面特征區(qū)域上,其中使用的幾種算子主要包括：膨脹、腐蝕、開(kāi)閉操作、骨架化以及裁剪操作等。該方法主要在于利用頂點(diǎn)的鄰接關(guān)系來(lái)表示形態(tài)學(xué)中的結(jié)構(gòu)元,為此定義了任意頂點(diǎn)vi的鄰接關(guān)系nhd,對(duì)于頂點(diǎn)vi的1-ring鄰接關(guān)系定義為：

nhd{i}={i}∪{j|?edge(vi,vj)}

(4)

鄰接關(guān)系可以逐次擴(kuò)展為n-ring鄰接關(guān)系nhdn：

nhdn+1=nhd(nhdn{i})n>1

(5)

1.2.1膨脹和腐蝕操作

膨脹操作可以將比較接近但不連續(xù)的特征區(qū)域合并在一起,并且特征區(qū)域向周圍擴(kuò)展,比原來(lái)的特征區(qū)域連續(xù)性更好。對(duì)特征區(qū)域F以nhdn為結(jié)構(gòu)元的膨脹可表達(dá)為：

dilaten(F)={j|?i∈F：j∈nhdn{i}}

(6)

膨脹操作是把特征區(qū)域F中頂點(diǎn)v的n-ring鄰接點(diǎn)全部都加入到特征集F中。

腐蝕操作的作用和膨脹操作相反,可以將特征區(qū)域細(xì)化,達(dá)到簡(jiǎn)化區(qū)域,使得特征結(jié)構(gòu)更清晰的目的。對(duì)特征區(qū)域F以nhdn為結(jié)構(gòu)元的腐蝕可表達(dá)為：

eroden(F)={j|nhdn{j}?F}

(7)

膨脹和腐蝕操作通常是組合來(lái)應(yīng)用的,網(wǎng)格上的特征區(qū)域需要經(jīng)過(guò)一系列的膨脹或腐蝕操作,兩者不同的組合會(huì)出現(xiàn)不同的操作效果,其中包括常用的開(kāi)操作、閉操作等,這些操作會(huì)使得特征區(qū)域在保持原有基本形狀,獲得較清晰的骨架特征。

1.2.2骨架化和裁剪操作

骨架化操作可以把網(wǎng)格上的特征區(qū)域細(xì)化為單個(gè)頂點(diǎn)寬度的特征線,骨架化操作的實(shí)質(zhì)是在判斷應(yīng)去除特征區(qū)域中的哪些頂點(diǎn)。特征區(qū)域中的頂點(diǎn)根據(jù)其復(fù)雜度的不同可以分為三類[13]。對(duì)于頂點(diǎn)vi的復(fù)雜度CP(i)定義為：

(8)

式中：Cnei(i)表示頂點(diǎn)vi依次按順時(shí)針的1-ring鄰接頂點(diǎn){v1,v2,…,vm},如果vk∈Cnei(i)且vk∈F,就令Cnei(i)k=1,否則Cnei(i)k=0。當(dāng)CP(i)≥4時(shí),頂點(diǎn)vi就是復(fù)雜點(diǎn)(Complex Vertex),構(gòu)成的點(diǎn)集記作C;當(dāng)CP(i)=0時(shí)是中心點(diǎn)(Center Vertex),構(gòu)成的點(diǎn)集記作⊙;其他的特征點(diǎn)為外側(cè)點(diǎn)(Disk Vertex),其點(diǎn)集記作○,通常CP(i)=2。

骨架化操作就是去除外側(cè)點(diǎn)并保留復(fù)雜點(diǎn)的過(guò)程,可用數(shù)學(xué)式表示為：

(9)

骨架化操作產(chǎn)生短的無(wú)關(guān)的分支可以通過(guò)裁剪操作消除。令S?F表示獲得的骨架的特征點(diǎn)集合,裁剪操作為：

(10)

2 牙頜模型自動(dòng)分割算法

2.1 算法概述

本文的牙頜模型分割算法中使用的牙頜模型是通過(guò)三維掃描儀根據(jù)牙頜的石膏模型獲取的點(diǎn)云數(shù)據(jù)生成的,牙頜數(shù)字化模型是一張三角網(wǎng)格曲面。本文的算法首先使用形態(tài)學(xué)方法生成齦緣線以及相鄰牙交線這兩種能用于牙齒分割的特征線,再通過(guò)端點(diǎn)的距離約束用最短路徑原則連接特征線斷開(kāi)的地方,得到封閉特征線用于牙頜模型的分割[22]。如圖1所示,算法的主要步驟有：

圖1 算法概述

Step1初始化特征點(diǎn)集。計(jì)算每個(gè)頂點(diǎn)v的平均曲率值κH(v),對(duì)于給定的曲閾值ε,若滿足κH(v)≤ε,則把該點(diǎn)v記為特征點(diǎn)加入到集合F。

Step2刪除噪點(diǎn)。將集合F中特征點(diǎn)按是否相連分類得到m個(gè)特征區(qū)域{F1,F2,…,Fm},記num(Fi)為特征區(qū)域Fi中頂點(diǎn)個(gè)數(shù),刪除num(Fi)/F<εd的特征區(qū)域,得到新的特征點(diǎn)集F,εd是設(shè)置的刪除比例。

Step3形態(tài)學(xué)處理。對(duì)特征點(diǎn)集依次做膨脹F=dilate(F),腐蝕F=erode(F),骨架化F=skeletonize(F),修剪F=prune(F)四種操作,得到單個(gè)頂點(diǎn)寬度的特征線F。

Step4識(shí)別特征線上的端點(diǎn)。對(duì)于頂點(diǎn)vi∈F,若vi的復(fù)雜度CP(vi)≤2,則把vi放入連接端點(diǎn)的集合ConnectP。

Step5對(duì)集合ConnectP中的點(diǎn)配對(duì),計(jì)算連接中點(diǎn)MidP,根據(jù)最短路徑連接配對(duì)的端點(diǎn),得到新的特征點(diǎn)集F。

Step6根據(jù)封閉的特征線F分割牙頜模型。

下面具體介紹Step5中的方法。

2.2 自動(dòng)閉合特征線

特征線的端點(diǎn)集合ConnectP中的頂點(diǎn)既有用于連接間斷處的端點(diǎn),也有一些未被完全刪除的分支上的末端點(diǎn),會(huì)造成干擾,如圖1(h)所示。需要解決如何從ConnectP中找出需要相互連接的頂點(diǎn),可以通過(guò)設(shè)置一個(gè)距離范圍δmax來(lái)約束端點(diǎn)的連接,避免錯(cuò)誤連接的出現(xiàn)。

令δ(v)表示頂點(diǎn)u到v的最短距離：

(11)

式中：Pmin(u,v)表示頂點(diǎn)u到v的最短路徑。查找相互連接的端點(diǎn)的過(guò)程是不斷計(jì)算每個(gè)頂點(diǎn)的vi∈ConnectP與其n-ring鄰域內(nèi)頂點(diǎn)u在δmax范圍內(nèi)的最短距離δ(v),對(duì)n依次取n=1,2,…,n,每計(jì)算一次都要判斷頂點(diǎn)u是否在其他端點(diǎn)vj也被訪問(wèn),如果訪問(wèn)過(guò),則把頂點(diǎn)u存入連接中點(diǎn)的集合MidP中,并將頂點(diǎn)u與端點(diǎn)vi、vj根據(jù)最短路徑相連,這個(gè)過(guò)程如圖1(i)、(j)所示。對(duì)ConnectP中的端點(diǎn)同時(shí)計(jì)算,查找中間點(diǎn),如果是特征線上的分支端,通常在δmax范圍內(nèi)不會(huì)找到與之配對(duì)的端點(diǎn),在特殊情況下如果存在,可以在形成封閉輪廓后,對(duì)網(wǎng)格模型做分割時(shí)剔除這些小的區(qū)域,不影響分割效果。

3 實(shí)驗(yàn)結(jié)果

本研究中用于調(diào)參的牙頜模型如圖2所示,這些模型是由激光掃描儀對(duì)牙頜的石膏模型掃描獲取的曲面信息生成的三角網(wǎng)格曲面,圖2(a)、(b)、(c)中模型的頂點(diǎn)數(shù)量分別為1.53×105、1.32×105、1.56×105。程序中的平均曲率閾值取ε=0.66,在刪除噪聲點(diǎn)過(guò)程中的刪除比例取εd=0.01,在閉合特征線過(guò)程中的距離范圍取δmax=3.5 mm。牙齒與牙齦的分割以及相鄰牙齒之間的分割都比較準(zhǔn)確完整,表明自動(dòng)的牙頜模型分割算法具有良好的分割效果。

圖2 本文算法自動(dòng)分割效果

本研究收集了另外5種常見(jiàn)模型(圖3)共49顆牙齒進(jìn)行測(cè)試,本算法成功實(shí)現(xiàn)了46顆牙齒的自動(dòng)分割,成功率達(dá)到93.88%。

圖3 本文算法在5種常見(jiàn)模型上的測(cè)試

4 結(jié) 語(yǔ)

本文改進(jìn)了基于形態(tài)學(xué)方法的牙頜模型分割算法,解決了依靠交互方式選取牙頜模型特征線端點(diǎn)的問(wèn)題,實(shí)現(xiàn)了牙頜模型的自動(dòng)分割。實(shí)驗(yàn)表明,該算法的分割效果較好,適用于不同牙頜模型的自動(dòng)分割。

相關(guān)的后續(xù)工作包括：研究其他特殊的牙病患者模型對(duì)算法的影響,進(jìn)一步提高算法的適用性;建立對(duì)牙頜模型分割效果的評(píng)估方法,判斷分割后的模型能否用于患者的治療。