陳 虎 馬珂楠 胡畢富

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

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

0 引 言

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

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

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

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

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

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

對于三角網(wǎng)格模型M通常可表示為M={V,T},其中V是網(wǎng)格頂點(diǎn)的集合,T是網(wǎng)格上三角面片面的集合。

(1)

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

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

(2)

在牙頜模型的特征線附近都是凹下去的“谷線”區(qū)域,這些區(qū)域頂點(diǎn)的平均曲率值都是偏小的負(fù)值,通常直接設(shè)置曲率閾值ε來獲得初始的特征點(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ū)域上,其中使用的幾種算子主要包括：膨脹、腐蝕、開閉操作、骨架化以及裁剪操作等。該方法主要在于利用頂點(diǎn)的鄰接關(guān)系來表示形態(tài)學(xué)中的結(jié)構(gòu)元,為此定義了任意頂點(diǎn)vi的鄰接關(guān)系nhd,對于頂點(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ò)展,比原來的特征區(qū)域連續(xù)性更好。對特征區(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á)到簡化區(qū)域,使得特征結(jié)構(gòu)更清晰的目的。對特征區(qū)域F以nhdn為結(jié)構(gòu)元的腐蝕可表達(dá)為：

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

(7)

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

1.2.2骨架化和裁剪操作

骨架化操作可以把網(wǎng)格上的特征區(qū)域細(xì)化為單個(gè)頂點(diǎn)寬度的特征線,骨架化操作的實(shí)質(zhì)是在判斷應(yīng)去除特征區(qū)域中的哪些頂點(diǎn)。特征區(qū)域中的頂點(diǎn)根據(jù)其復(fù)雜度的不同可以分為三類[13]。對于頂點(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)的過程,可用數(shù)學(xué)式表示為：

(9)

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

(10)

2 牙頜模型自動分割算法

2.1 算法概述

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

圖1 算法概述

Step1初始化特征點(diǎn)集。計(jì)算每個(gè)頂點(diǎn)v的平均曲率值κH(v),對于給定的曲閾值ε,若滿足κ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é)處理。對特征點(diǎn)集依次做膨脹F=dilate(F),腐蝕F=erode(F),骨架化F=skeletonize(F),修剪F=prune(F)四種操作,得到單個(gè)頂點(diǎn)寬度的特征線F。

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

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

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

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

2.2 自動閉合特征線

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

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

(11)

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

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

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

圖2 本文算法自動分割效果

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

圖3 本文算法在5種常見模型上的測試

4 結(jié) 語

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

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