楊 玲 李 宇 王中科 饒妮妮
1(電子科技大學(xué)生命科學(xué)與技術(shù)學(xué)院,成都 6 10054)2(成都信息工程學(xué)院電子工程學(xué)院,成都 6 10225)3(成都信息工程學(xué)院網(wǎng)絡(luò)工程學(xué)院,成都 6 10225)
牙齒的三維建模在牙種植領(lǐng)域中有著十分重要的意義。目前,基于特征點(diǎn)的三維重建算法由于其數(shù)據(jù)量少而精和可快速重建目標(biāo)物的結(jié)構(gòu)特征,具有較強(qiáng)的實(shí)用性[1]。目標(biāo)物的特征點(diǎn)從幾何方面可分為點(diǎn)、線、面三種,其中點(diǎn)特征是最常采用的一種,如灰度的局部最大值點(diǎn)、邊緣點(diǎn)、角點(diǎn)和線交叉點(diǎn)等[2-8]。本研究中提及的牙齒特征點(diǎn)擬指在牙齒CT各斷層圖像中表示牙齒形狀輪廓的邊緣特征點(diǎn)。曲線或曲面演化是解決圖像分割和目標(biāo)檢測(cè)的一種有效方法。曲線演化是曲線在自身內(nèi)力與外力(代表圖像信息)的共同作用下,變化收斂到圖像的邊緣或其它特征處的過程。因此,曲面演化常在醫(yī)學(xué)上用于特征點(diǎn)的提取。曲線演化的方式有多種,如Sethian等提出的Level Set算法[9],Kass等提出的Snake方法[10],Chan等提出的基于Mumford-Shah理論的CV模型[11],以及近期提出的窄帶法和無(wú)重復(fù)初始化速度函數(shù)等方法[12-13]。上述方法為處理不同拓?fù)浣Y(jié)構(gòu)物體的變形提供了方便,且抗噪性能較強(qiáng),然而這些曲線演化方法仍然存在許多的不足之處。例如,基于水平集的曲線演化依賴于圖像梯度的局部信息,在梯度不明顯處或?qū)Ρ榷炔淮蟮膮^(qū)域容易產(chǎn)生曲線泄漏[14];而CV模型沒有利用圖像局部的梯度信息,當(dāng)初始水平集位于平滑區(qū)域或者極度凹陷區(qū)域進(jìn)行分割時(shí),水平集曲線的收斂速度會(huì)非常緩慢[15]。盡管在近期提出的算法中,對(duì)速度函數(shù)、參加運(yùn)算的數(shù)據(jù)量進(jìn)行了改善,但運(yùn)算效率仍是該類算法需要重點(diǎn)優(yōu)化的問題。
本研究結(jié)合曲線演化算法抗噪能力強(qiáng)的優(yōu)勢(shì)和牙齒CT序列圖像數(shù)據(jù)的特點(diǎn),改進(jìn)離散曲線演化(discrete curve evolution,DCE)模型,并基于改進(jìn)的DCE提出一種提取牙齒邊緣特征點(diǎn)的新方法。新方法通過控制不同層圖像的特征點(diǎn)數(shù),能快速、準(zhǔn)確地提取牙齒的三維特征點(diǎn),既保證了特征點(diǎn)提取的信息量足以重建出整顆牙齒,又能降低數(shù)據(jù)的冗余量,簡(jiǎn)化了演化過程。最后,采用真實(shí)的磨牙CT圖像數(shù)據(jù)驗(yàn)證新方法的可行性和計(jì)算效率。
DCE算法是Latecki等為提取目標(biāo)物骨架而提出的一種快速獲得復(fù)雜多邊形拓?fù)浣Y(jié)構(gòu)的算法[16]。由于圖像中物體的輪廓常受噪聲影響而形變,DCE實(shí)際上是輪廓預(yù)處理過程,包括兩方面的含義:一是除去邊界噪聲;二是保留視覺感官上的重要部分。它對(duì)復(fù)雜多邊形進(jìn)行離散曲線演化,既能降低噪聲的影響,又能排除不重要的形狀特征,從而簡(jiǎn)化形狀。文獻(xiàn)[17,18]證明了這個(gè)結(jié)論。
DCE算法在獲取邊緣信息時(shí),首先采用形態(tài)學(xué)算法對(duì)圖像進(jìn)行去噪、二值化,提取出感興趣的連通域,然后采用梯度算子和輪廓跟蹤法獲取連通域的封閉輪廓線,以此作為DCE演化的初始特征點(diǎn),然后從最大凸(或凹)弧線的起點(diǎn)開始演化,獲取圖像邊緣上的突變點(diǎn),并在每一步演化中去除相關(guān)性最小的頂點(diǎn),最后保留該圖形的拓?fù)浣Y(jié)構(gòu)。其中初始輪廓點(diǎn)的選取對(duì)曲線演化結(jié)果有著重要的影響,算法中采用了多種去噪方法去除圖像的偽邊緣,降低了初始輪廓點(diǎn)對(duì)噪聲的敏感性,以便獲取單一的閉合輪廓曲線,確保曲線演化結(jié)果的可靠性。
由于任意數(shù)字曲線可以看作一個(gè)含有大量頂點(diǎn)的多邊形,所以能將多邊形的演化用數(shù)字曲線來代替。設(shè)C為實(shí)平面上的一閉合多邊形,S1…,Sm為構(gòu)成C的連續(xù)線段,P1…Pm,為C的頂點(diǎn)。其中,為Si和Si+1的公共頂點(diǎn)。為得到閉合曲線,可約定Pm+i=Pi。在曲線演化的每一步中,找出相關(guān)度最小的頂點(diǎn)Pi,用一段連接Pi+1和Pi-1的線段代替原來的兩段相鄰線段Si和Si+1。這樣每一步將去掉多邊形的一個(gè)頂點(diǎn)Pi。
由此可見,離散曲線演化的基本思想是:在演化的每一步,對(duì)連續(xù)線段S1、S2用單一線段S1∪S2連接端點(diǎn)替代,如圖1所示。圖中是在曲線演化過程中將相關(guān)性最小的Pi點(diǎn)去掉。
圖1 DCE演化示意圖Fig.1 The evolution scheme of DCE
每個(gè)頂點(diǎn)的相關(guān)度參數(shù)K的計(jì)算式為
式中,β(S1,S2)是線段 S1和 S2的拐角,以逆時(shí)針為正;根據(jù)曲線的凹凸特性,β(S1,S2)的計(jì)算方法有所差異。
圖2所示為曲線凹凸對(duì)拐角的影響,其中b點(diǎn)為當(dāng)前頂點(diǎn),a、c分別為b點(diǎn)前后→相鄰的兩個(gè)頂點(diǎn),S1表示b和a兩點(diǎn)構(gòu)成的線段—ba、S2表示 b 和 c 兩點(diǎn)構(gòu)成的線段→—bc,β表示從線段 S 逆時(shí)針旋轉(zhuǎn)至線1段S2的夾角。圖(a)是凸曲線的拐角,小于90°,則β(S1,S2)的計(jì)算公式為
式中,l為線段的長(zhǎng)度函數(shù),計(jì)算公式為
式中,(ΔxS1,ΔyS2),(ΔxS2,ΔyS1)分別表示線段 S1和S2兩端點(diǎn)x和y的坐標(biāo)差值。
圖(b)是凹曲線的拐角,大于 9 0°,則 β(S1,S2)的計(jì)算公式為
根據(jù)式(1)計(jì)算的K值越高,則表明弧S1∪S2對(duì)形狀的貢獻(xiàn)越大。DCE算法將K作為曲線演化的控制參數(shù)。如對(duì)采樣后的閉合多邊形C,若P1,…,Pm為C的頂點(diǎn),可定義M(p)為其頂點(diǎn)個(gè)數(shù),也即特征點(diǎn)個(gè)數(shù),設(shè)定其閾值為T,使DCE算法演化時(shí),通過K值的計(jì)算,依次去掉K較小的點(diǎn),使M(p)=T,從而結(jié)束曲線的演化,產(chǎn)生相應(yīng)圖像的骨架。
DCE算法可以保留圖形的拓?fù)浣Y(jié)構(gòu)??芍?,采用DCE提取的圖像特征點(diǎn)個(gè)數(shù)越多,表達(dá)的圖像信息就越準(zhǔn)確。然而,由于DCE算法停止演化是由設(shè)定的控制點(diǎn)數(shù)來確定,因此,對(duì)序列圖像而言,這會(huì)造成不同結(jié)構(gòu)具有相同的停止條件,也會(huì)造成數(shù)據(jù)點(diǎn)冗余或不足,最終導(dǎo)致三維重建達(dá)不到滿意的效果。
圖3(a)為一磨牙的第330層CT圖像(512像素×512像素)。用DCE算法對(duì)其右上方連通域進(jìn)行特征提取,結(jié)果如圖3中(b)~(e)所示,其中(b)和(d)將拐角定義為凸點(diǎn)提取,分別提取了10和20個(gè)特征點(diǎn);(c)和(e)將拐角定義為凹點(diǎn)提取,分別提取了10和20個(gè)特征點(diǎn)。圖3結(jié)果表明,DCE算法會(huì)因拐角設(shè)置不同而導(dǎo)致提取的特征點(diǎn)存在差異,同時(shí)也表明特征點(diǎn)個(gè)數(shù)增加到某種程度時(shí),拐角定義為凹或凸對(duì)特征點(diǎn)的提取影響不大。
為了克服DCE算法的特征點(diǎn)表達(dá)信息不夠的問題,現(xiàn)有兩種解決方案:一是定義兩種拐角,對(duì)一幅圖像先進(jìn)行凹(凸)點(diǎn)提取,再進(jìn)行凸(凹)點(diǎn)提?。欢菍⒚糠鶊D像特征點(diǎn)的數(shù)目提高到一定程度。然而,上述兩種方案會(huì)使程序的執(zhí)行效率降低,甚至產(chǎn)生特征點(diǎn)的冗余問題。
本研究提出改變DCE的停止準(zhǔn)則M(p)=T,使之成為可控點(diǎn)的操作方式;并提出曲線特征因子量,使DCE特征點(diǎn)數(shù)量的選取具有自適應(yīng)性和合理性。
設(shè)閉合曲線Y=F(X)在某一區(qū)間I上連續(xù),曲線上共有N個(gè)拐點(diǎn),第j個(gè)拐點(diǎn)為(x(j),F(xiàn)(x(j)))(j=1,2,…,N),則將 N相對(duì)于構(gòu)成該曲線所有點(diǎn)的個(gè)數(shù)即曲線長(zhǎng)度進(jìn)行歸一化處理,用來表示曲線的復(fù)雜程度,并稱之為曲線特征因子量(characteristic factor of curve,CFOC)。即
式中,L為曲線上數(shù)據(jù)點(diǎn)總數(shù)。
如圖4所示是一牙齒CT圖像(分辨率為512像素×512像素)的13至483斷層的曲線特征因子量。橫坐標(biāo)表示不同斷層圖像,縱坐標(biāo)表示特征因子量。圖4中(a)是13至483斷層的曲線特征因子量,每一層圖像對(duì)應(yīng)的多個(gè)離散點(diǎn)分別表示多個(gè)連通域的特征因子量。取其各層曲線特征因子量中最大值對(duì)應(yīng)的點(diǎn),得到13至483斷層的最大曲線特征因子量,如圖4(b)所示。為自適應(yīng)地確定不同斷層圖像的特征點(diǎn)數(shù)Z,可將特征因子量CFOC分成若干等級(jí),不同等級(jí)提取不同數(shù)目的特征點(diǎn)。
在保留其拓?fù)淝疤嵯?,改進(jìn)算法對(duì)不同等級(jí)的特征因子量設(shè)定特征點(diǎn)數(shù)Z的原則是:若各斷層曲線相差較大,即CFOC曲線隨CT層數(shù)變化陡峭,則各層的Z值可設(shè)定為跨度較大的值;反之,當(dāng)各斷層曲線相差較小,即CFOC曲線隨CT層數(shù)變化比較平緩,各層的Z值可取跨度較小的值。
分別將CFOC分成了5個(gè)和8個(gè)等級(jí),并按式(7)設(shè)定了各斷層的特征點(diǎn)數(shù)Z與最大曲線特征因子量等級(jí)n的關(guān)系:
從圖 4(b)可以看出,當(dāng)?。?.0、0.2],[0.2、0.4],[0.4、0.6],[0.6、0.8],[0.8、1.0]5 個(gè)等級(jí)區(qū)間時(shí),曲線變化相對(duì)平坦的區(qū)域歸為同一區(qū)間,而不同區(qū)間變換相對(duì)陡峭,此時(shí)不同區(qū)間可采用跨距較大的特征點(diǎn)數(shù) Z,即為 3、5、9、17和 33個(gè)。而當(dāng)分成8個(gè)均勻等級(jí)區(qū)間時(shí),就會(huì)降低各區(qū)間之間的陡峭性,使Z值跨距減少,冗余點(diǎn)增多。
圖4 牙齒CT序列切片的斷層曲線特征因子量。(a)各斷層曲線特征因子量;(b)各斷層最大曲線特征因子量Fig.4 CFOC of Dental CT slices.(a)CFOC of each slice;(b)The max CFOC of each slice
該方案既能保證牙齒圖像特征點(diǎn)提取的有效性,又能減少數(shù)據(jù)的冗余度。
圖5為一磨牙的第29、119、229和329層 CT圖像。采用改進(jìn)的DCE算法對(duì)圖5中4層CT圖像進(jìn)行特征點(diǎn)提取,結(jié)果如圖6所示。圖6中(a)~(c)分別是從圖5的第29層CT圖像的3個(gè)連通域中提取的特征點(diǎn),且從每個(gè)連通域中提取的特征點(diǎn)數(shù)均為5;圖6(d)是從圖5的第119層CT圖中提取的特征點(diǎn),數(shù)目為17個(gè);圖6(e)是從圖5的第229層CT圖中提取的特征點(diǎn),數(shù)目為17個(gè);圖6中(f)~(h)分別是從圖5的第329層圖像的3個(gè)連通域中提取的特征點(diǎn),每個(gè)連通域提取了9個(gè)特征點(diǎn)。從圖6的結(jié)果可以看出,改進(jìn)的DCE算法能夠從磨牙CT中準(zhǔn)確有效的提取出特征點(diǎn)。
接著可利用提取的特征點(diǎn),重建牙齒的三維曲面,具體方法是:對(duì)每一層牙齒的特征點(diǎn)采用不同比率的樣條插值得到各斷層曲線數(shù)據(jù),再將每層曲線數(shù)據(jù)擴(kuò)展到三維空間,用樣條曲線插值生成對(duì)應(yīng)的曲面數(shù)據(jù)。圖7是在VTK(Visualization ToolKit)開發(fā)平臺(tái)下,利用圖6的特征點(diǎn),采用樣條曲線重建出的牙齒表面圖??梢钥闯?,牙齒整個(gè)表面重建完整,牙冠和牙根明顯,與原牙齒曲面結(jié)構(gòu)吻合,證實(shí)了所提出算法的可行性和有效性。
圖5 一磨牙的4層CT圖像(512像素×512像素)。(a)第29層;(b)第119層;(c)第229層;(d)第329層Fig.5 Four CT slices of a molar(512 pixel×512 pixel).(a)The 29thslice;(b)The 119thslice;(c)The 229thslice;(d)The 329thslice
表1給出了采用DCE和改進(jìn)的DCE算法提取圖5磨牙CT圖像特征點(diǎn)的計(jì)算效率情況。DCE算法沒有經(jīng)過特征點(diǎn)個(gè)數(shù)判定,直接在每一層牙齒CT圖的各個(gè)連通域提取固定的12個(gè)特征點(diǎn),所需總時(shí)間為7 min 35 s,每一層平均所需時(shí)間為0.97 s。改進(jìn)的DCE算法能夠自適應(yīng)地確定每一層牙齒CT圖的各個(gè)連通域需要提取的特征點(diǎn)數(shù)目,降低了數(shù)據(jù)存儲(chǔ)的冗余量,同時(shí)提高對(duì)目標(biāo)特征表達(dá)的精確性。其中特征點(diǎn)數(shù)目確定方法所需的總時(shí)間為2 min 56 s和3 min 11 s,每一層平均所需時(shí)間分別為0.37 s和0.41 s??梢姼倪M(jìn)算法在每層上提取特征點(diǎn)所花時(shí)間約為原算法的50%。此外從提取的總特征量而言,改進(jìn)算法提取的特征量約為原算法的80%左右,因此改進(jìn)DCE算法的計(jì)算效率明顯優(yōu)于DCE算法。
圖6 對(duì)圖5中的4層CT圖像進(jìn)行特征點(diǎn)提取的結(jié)果。(a)對(duì)第29層圖像的連通域1提取5點(diǎn);(b)對(duì)第29層圖像的連通域2提取5點(diǎn);(c)對(duì)第29層圖像的連通域3提取5點(diǎn);(d)對(duì)第119層圖像的連通域提取17點(diǎn);(e)對(duì)第229層圖像的連通域提取17點(diǎn);(f)對(duì)第329層圖像的連通域1提取9點(diǎn);(g)對(duì)第329層圖像的連通域2提取9點(diǎn);(h)對(duì)第329層圖像的連通域3提取9點(diǎn)Fig.6 The extracted feature points of Fig.5’s four CTslices.(a)5pointsofthefirst connected area of the 29thslice;(b)5 points of the second connected area of the 29thslice;(c)5 points of the third connected area of the 29th slice;(d)17 points of the connected area of the 119thslice;(e)17 points of the connected area of the 229thslice;(f)9 points of the first connected area of the 329thslice;(g)9 points of the second connected area of the 329thslice;(h)9 points of the third connected area of the 329th slice
圖7 利用改進(jìn)DCE算法提取的特征點(diǎn)重建的牙齒表面圖Fig.7 The reconstructed tooth surface with the extracted feature points by improved DCE algorithm
表1 DCE和改進(jìn)DCE算法提取牙齒13~483層CT圖特征點(diǎn)效率的對(duì)比Tab.1 The contrast of efficiency off eature points extracted by DCE and improved DCE for 13rd~483rd CT slices
從表1中也可看出,當(dāng)選擇不同的特征因子量等級(jí)劃分方法和特征點(diǎn)數(shù),改進(jìn)算法效率會(huì)有所差異,故為發(fā)揮改進(jìn)算法的最優(yōu)效率,還需進(jìn)一步優(yōu)化特征因子量的等級(jí)劃分和特征點(diǎn)數(shù)的選擇方法。
本研究利用改進(jìn)后的離散曲線演化算法提取牙齒CT圖像的邊緣特征點(diǎn),并在VTK平臺(tái)上用樣條插值方法對(duì)牙齒特征點(diǎn)進(jìn)行曲面重建,驗(yàn)證了算法的有效性。所提出的特征點(diǎn)提取算法的優(yōu)勢(shì)主要體現(xiàn)在:一是由于采用基于曲線演化模型,所以算法的抗噪能力強(qiáng),特征點(diǎn)的定位準(zhǔn)確;二是算法在提取特征點(diǎn)之前自適應(yīng)地計(jì)算所要提取的特征點(diǎn)數(shù)量。由于不同信息量的連通域所需的特征點(diǎn)數(shù)不同,因此,本算法既保證了獲取的特征點(diǎn)信息能足以表達(dá)出整幅牙齒信息,也降低了特征點(diǎn)的冗余量。三是本算法相比傳統(tǒng)曲線演化模型,效率較高,對(duì)多幅圖像的處理具有優(yōu)勢(shì)。實(shí)驗(yàn)結(jié)果表明,該算法能夠準(zhǔn)確地提取出牙齒圖像邊界的特征點(diǎn),且重建出的牙齒三維曲面結(jié)構(gòu)完整清晰,其性能優(yōu)于相關(guān)傳統(tǒng)算法。
所提出算法還存在計(jì)算速度較慢的問題,今后的工作將致力于算法執(zhí)行效率的提高。初步思路是進(jìn)一步優(yōu)化各斷層特征點(diǎn)數(shù)與特征因子量等級(jí)的函數(shù)關(guān)系,并對(duì)同一層包含多個(gè)連通域的情況,可分別根據(jù)每個(gè)連通域的自身的曲線特征因子量選取不同的特征點(diǎn)數(shù),以克服目前的統(tǒng)一處理模式,從而進(jìn)一步降低特征點(diǎn)的冗余度,提高算法的執(zhí)行效率。
[1]王付新,黃毓瑜,孟偲.三維重建中特征點(diǎn)提取算法的研究與實(shí)現(xiàn)[J].工程圖學(xué)學(xué)報(bào),2007,3:91-96.
[2]謝萍,馬小勇,張憲民.一種快速的復(fù)雜多邊形匹配算法[J].計(jì)算機(jī)工程,2003,29(16):176-181.
[3]彭文,童若鋒,錢歸平.使用特征點(diǎn)與曲線配準(zhǔn)醫(yī)學(xué)圖像[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2007,19(9):1126-1131.
[4]Baboulaz L,Dragotti PL.Exact feature extraction using finite rate of innovation principles with an application to image superresolution[J].IEEE Transactions on Image Processing,2009,18(2):281-298.
[5]邵巍,朱圣英,陳靈芝.一種快速多尺度特征點(diǎn)匹配算法[J].中國(guó)圖像圖形學(xué)報(bào),2009,14(12):2572-2576.
[6]雷琳,王壯,粟毅.基于多尺度Gabor濾波器組的不變特征點(diǎn)提取新方法[J].電子學(xué)報(bào),2009,37(10):2314-2319.
[7]李竹林,王文東,趙宗濤.改進(jìn)的邊緣特征點(diǎn)提取算法[J].計(jì)算機(jī)工程與應(yīng)用,2009,45(2):63-65.
[8]王建文,李青.基于點(diǎn)和邊緣相結(jié)合特征提取的圖像配準(zhǔn)算法[J].計(jì)算機(jī)工程與設(shè)計(jì),2009,30(4):928-930.
[9]Osher S,Sethian J.Fronts propagation with curvature dependent speed:Algorithms based on Hamilton-Jacobi formulations[J].Journal of Comp.Phys.1988,79:12-49.
[10]Kass M,Witkin A,Terzopoulos D.Snakes:Active contour models[J].International Journal of Computer Vision,1988,1(4):321-331.
[11]Chan TF,Vese LA.Active Contours Without Edges[J].IEEE Transactions on Image Processing,2001,10(2):266-276.
[12]Wang Hongyuan,Zhou Zeming.Study on the Narrow Band Level Set Method[J].Journal of Communication and Computer,2006,3(9):37-40.
[13]Li Chunming,Xu Chenyang,Gui Changfeng,et al.Level Set Evolution Without Re-initialization:A New Variational Formulation[C]//Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition.San Diego:IEEE,2005,1:430-436.
[14]溫鐵祥,楊豐.一種新的曲線演化混合模型圖像分割算法[J].電路與系統(tǒng)學(xué)報(bào),2007,12(4):48-52.
[15]Vese LA,Chan TF.A multiphase level set framework for image segmentation using the Mumford and Shah model[J].International Journal of Computer Vision,2002,50(3):271-293.
[16]Latecki LJ,Lakamper R.Shape similarity measure based on correspondence of visual parts[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2000,22(10):1185-1190.
[17]Bai Xiang,Latecki LJ,Liu Wenyu.Skeleton Pruning by Contour Partitioning with Discrete Curve Evolution[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2007,29(3):449-462.
[18]Krinidis S,Chatzis V.A skeleton family generator via physicsbased deformable models[J].IEEE Transactions on Image Processing,2009,18:1-11.