趙夫群,耿國華
(1.咸陽師范學(xué)院 教育科學(xué)學(xué)院,陜西 咸陽 712000;2.西北大學(xué) 信息科學(xué)與技術(shù)學(xué)院,陜西 西安 710127)
顱面復(fù)原是一項(xiàng)對(duì)人類顱骨的面部容貌進(jìn)行復(fù)原的技術(shù),它以顱骨的形狀特征和顱面復(fù)原技術(shù)為基礎(chǔ),以人類的面部軟組織統(tǒng)計(jì)厚度為依據(jù),采用一定的算法將軟組織添加到顱骨上,從而實(shí)現(xiàn)顱骨面貌復(fù)原。常用的顱骨軟組織添加方法有兩種[1]:一是顱骨變形法,即采用硬組織填充軟組織;二是三維體積變形法,即從與待復(fù)原顱骨相似的顱骨中獲得軟組織厚度。這兩種方法都涉及未知顱骨與參考顱骨的配準(zhǔn)問題,因此,顱骨配準(zhǔn)是顱面復(fù)原的一個(gè)重要步驟,其正確性對(duì)顱骨面貌的復(fù)原起著關(guān)鍵性的作用。
基于孔洞輪廓線的顱骨配準(zhǔn)的基本思路為:對(duì)于一個(gè)待復(fù)原顱骨U(也叫未知顱骨U),采用一定的配準(zhǔn)算法從顱骨數(shù)據(jù)庫中找出與U最為相似的一個(gè)或多個(gè)參考顱骨S,那么顱骨S的面貌即可作為U的參考面貌,從而為未知顱骨U的顱面復(fù)原提供依據(jù)。目前,顱骨配準(zhǔn)已經(jīng)在考古、醫(yī)學(xué)研究以及刑事案件偵破等領(lǐng)域[2-4]得到了一定的應(yīng)用。
由于顱骨的三維數(shù)據(jù)模型復(fù)雜,含噪聲和外點(diǎn)較多,因此對(duì)其配準(zhǔn)精度要求較高。目前,三維顱骨配準(zhǔn)大多采用基于特征的配準(zhǔn)方法,即全局特征配準(zhǔn)方法和局部特征配準(zhǔn)方法[5-7]。全局特征描述了整個(gè)顱骨模型,而局部特征只描述顱骨的關(guān)鍵特征點(diǎn)或點(diǎn)的鄰域特征。由于三維顱骨模型的點(diǎn)或線特征較為明顯,因此局部特征比全局特征更適用于部分覆蓋的三維顱骨模型的配準(zhǔn)。在基于局部特征的顱骨配準(zhǔn)方法中,特征點(diǎn)標(biāo)定法[8-9]是使用較多的方法,但配準(zhǔn)結(jié)果并不十分理想。法向、曲率、凹或凸的特征區(qū)域等也是顱骨局部特征描述的重要方法,能夠在一定程度上實(shí)現(xiàn)三維顱骨模型的配準(zhǔn)。此外,迭代最近點(diǎn)(Iterative closest point, ICP)[10]及其改進(jìn)算法[11-14]也被用在了三維顱骨模型的配準(zhǔn)中。但是由于ICP算法對(duì)待配準(zhǔn)模型的初始相對(duì)位置要求較高,因此一般要先進(jìn)行顱骨的粗配準(zhǔn),然后再采用ICP或其改進(jìn)的算法來實(shí)現(xiàn)顱骨細(xì)配準(zhǔn)。
針對(duì)三維顱骨模型數(shù)據(jù)量大、分辨率差異大等問題,本文提出一種基于孔洞輪廓線的三維顱骨模型配準(zhǔn)方法。首先,提取顱骨的眼眶、鼻框、顳骨、上頜骨以及下頜骨等孔洞輪廓線,并通過對(duì)輪廓線的匹配來實(shí)現(xiàn)顱骨粗配準(zhǔn);然后,再采用PICP算法將顱骨進(jìn)行細(xì)配準(zhǔn),從而實(shí)現(xiàn)顱骨的最終精確配準(zhǔn)。
這里提取的顱骨孔洞輪廓線(簡稱輪廓線)主要包括眼眶、鼻框、顳骨、上頜骨以及下頜骨輪廓線等5種類型,如圖1所示。
圖1 顱骨輪廓線的類型Fig.1 Types of skull contour lines
對(duì)于顱骨的三角網(wǎng)格數(shù)據(jù)模型,若一條邊只被一個(gè)三角形使用,則稱該邊為邊界邊,邊界邊上的點(diǎn)稱作邊界點(diǎn)。多條邊界邊首尾相連則構(gòu)成一條輪廓線。定義一條輪廓線的長度為該輪廓線包含的邊界點(diǎn)的數(shù)目,而兩條輪廓線l1i和l2j之間的最短距離為min{distance(pm,qn)|pm∈l1i,qn∈l2j},distance(pm,qn)為輪廓線上邊界點(diǎn)pm和qn的歐氏距離。
根據(jù)文獻(xiàn)[15]統(tǒng)計(jì)的輪廓線的長度以及輪廓線間的最短距離的均值和標(biāo)準(zhǔn)差,即可對(duì)三維顱骨模型的輪廓線類型進(jìn)行自動(dòng)識(shí)別。具體步驟如下:
1)由于眼眶、鼻框和顳骨的輪廓線長度相近,上頜骨和下頜骨的輪廓線長度與這3種輪廓線長度差異很大,因此可直接從中區(qū)分出上頜骨和下頜骨輪廓線。
2)如果顱骨包含上頜骨輪廓線,則先任意確定其左右方位。
3)對(duì)于剩下的待識(shí)別輪廓線,根據(jù)其長度的均值和標(biāo)準(zhǔn)差便可確定其可能的輪廓類型。通過該步驟判斷出的每條輪廓線可能有多種類型,而且某些類型還要進(jìn)一步細(xì)化為左右兩種。比如,眼眶輪廓線又有左眼眶和右眼眶輪廓線之分。
4)對(duì)于每一種輪廓線類型(左眼眶、右眼眶、鼻框、左顳骨、右顳骨),執(zhí)行步驟3)便可獲得其對(duì)應(yīng)的輪廓線集合。
5)采用兩步法從所有輪廓線組合中篩選出正確的輪廓線組合。首先,剔除元素有重復(fù)的輪廓線組合,再在剩下的每個(gè)組合中添加上/下頜骨輪廓線;然后,計(jì)算所有輪廓線間的最短距離,并判斷每個(gè)距離值是否滿足輪廓線間最短距離的條件,若滿足則記1分,否則記0分,最終,得分最高的組合即為輪廓線的分類結(jié)果。
一條輪廓線上的邊界點(diǎn)構(gòu)成一條空間離散曲線,本文采用四次B樣條曲線將其擬合成光滑的空間曲線。B樣條曲線的定義如下[16]:
給定m=n+k+1個(gè)頂點(diǎn),則可以定義n+1段k次參數(shù)曲線,第i段B樣條曲線函數(shù)可以表示為
(1)
其中,s=0,1,…,k,i=0,1,…,n;pi+s為第i段曲線特征多邊形的k+1個(gè)頂點(diǎn);fs,k(t)為B樣條基底函數(shù),可表示為
(2)
對(duì)于四次B樣條曲線,k=4,即s=0,1,2,3,4,其基底函數(shù)為
(3)
(4)
(5)
(6)
(7)
那么,第i段B樣條曲線的矩陣表達(dá)式可寫為
(8)
對(duì)于第1部分提取的輪廓線上每一個(gè)邊界點(diǎn)pi,選取與其相鄰的前后各兩個(gè)點(diǎn),即pi-2,pi-1和pi+1,pi+2,對(duì)pi及其相鄰點(diǎn)共5個(gè)點(diǎn)采用四次B樣條曲線對(duì)其擬合,假設(shè)得到擬合輪廓線l(t)=(x(t),y(t),z(t))。
輪廓線采用邊界點(diǎn)的曲率和撓率表示,即將輪廓線上邊界點(diǎn)的曲率和撓率組成特征串,通過計(jì)算兩條輪廓線的相似度進(jìn)行輪廓線的匹配。
設(shè)輪廓線l(t)=(x(t),y(t),z(t))的一階和二階導(dǎo)數(shù)分別為l′(t)=(x′(t),y′(t),z′(t))和l″(t)=(x″(t),y″(t),z″(t)),于是輪廓線l(t)的曲率k(t)和撓率τ(t)分別為
(9)
(10)
那么,基于輪廓線特征串的顱骨配準(zhǔn)方法的具體步驟描述如下:
1)設(shè)置初值i=1,j=1,i≤m,j≤n。
2)取未知顱骨U的第i條輪廓線l1i。
3)取參考顱骨S中的第j條輪廓線l2j,用式(11)計(jì)算l1i和l2j的相似度ξ1。若ξ1大于給定閾值,則用四元數(shù)法計(jì)算l1i和l2j的旋轉(zhuǎn)矩陣R和平移矩陣t,將l1i和l2j對(duì)齊;否則,轉(zhuǎn)到步驟4)。
(11)
4)執(zhí)行j=j+1,若j≤n,轉(zhuǎn)到步驟3),否則,轉(zhuǎn)到步驟5)。
5)執(zhí)行i=i+1操作,若i≤m,轉(zhuǎn)到步驟2),否則,轉(zhuǎn)到步驟6)。
6)判斷,若顱骨U中的輪廓線都能在S中找到匹配的輪廓線,且顱骨S中的輪廓線都能在U中找到匹配的輪廓線,并且一一對(duì)應(yīng),那么顱骨U和S配準(zhǔn)成功,否則配準(zhǔn)失敗。
通過顱骨輪廓線的匹配,兩個(gè)顱骨已經(jīng)基本對(duì)齊。接下來再采用PICP 算法[14]將兩個(gè)顱骨進(jìn)行細(xì)配準(zhǔn)。PICP 算法是在ICP 算法的基礎(chǔ)上,通過添加高斯概率模型實(shí)現(xiàn)的,該算法具有較強(qiáng)的抗噪性,適用于顱骨的細(xì)配準(zhǔn)。
實(shí)驗(yàn)數(shù)據(jù)采用西北大學(xué)可視化技術(shù)研究所提供的CT掃描的顱骨三角網(wǎng)格數(shù)據(jù)模型。對(duì)于一個(gè)未知顱骨U,在顱骨數(shù)據(jù)庫中配出一個(gè)或者幾個(gè)相似顱骨,作為顱面復(fù)原的參考顱骨。采用本文配準(zhǔn)方法,首先提取待配準(zhǔn)顱骨的眼眶、鼻框、顳骨、上頜骨以及下頜骨等孔洞輪廓線,然后通過特征串的匹配實(shí)現(xiàn)顱骨輪廓線的匹配,由此實(shí)現(xiàn)顱骨的配準(zhǔn),最后采用PICP算法將顱骨進(jìn)一步細(xì)配準(zhǔn)。
通過將未知顱骨U(如圖2(a)所示)與顱骨庫中265個(gè)顱骨S1~ S265進(jìn)行配準(zhǔn),找到了顱骨U的一個(gè)最佳匹配參考顱骨S1(如圖2(b)所示),其配準(zhǔn)結(jié)果如圖3所示。U與剩下的264個(gè)顱骨S2~ S265均配準(zhǔn)失敗,部分配準(zhǔn)失敗的結(jié)果如圖4所示。
圖2 待配準(zhǔn)顱骨Fig.2 Skulls need to be registered
圖3 U和S1配準(zhǔn)結(jié)果的正、側(cè)面Fig.3 The front and side registration result of U and S1
從圖3可見,采用本文配準(zhǔn)方法,未知顱骨U和參考顱骨S1能夠得到良好的配準(zhǔn)結(jié)果,顱骨S1可以作為顱骨U的參考顱骨,可以為未知顱骨U提供可能的復(fù)原面貌參考。
圖4 U和S2~S4的配準(zhǔn)結(jié)果Fig.4 Registration results of U and S2~S4
為了進(jìn)一步說明本文顱骨配準(zhǔn)方法的性能,對(duì)未知顱骨U和參考顱骨S1再分別單獨(dú)采用PICP算法和可信區(qū)域配準(zhǔn)算法[17]進(jìn)行配準(zhǔn),PICP算法的配準(zhǔn)結(jié)果如圖5所示,可信區(qū)域配準(zhǔn)算法的配準(zhǔn)結(jié)果如圖6所示。
圖5 PICP算法配準(zhǔn)結(jié)果的正、側(cè)面Fig.5 Registration result of PICP algorithm
圖6 可信區(qū)域配準(zhǔn)算法配準(zhǔn)結(jié)果的正、側(cè)面Fig.6 Registration result of trust region registration algorithm
從圖5和圖6的配準(zhǔn)結(jié)果來看,雖然PICP算法具有較強(qiáng)的抗噪性,但是單獨(dú)采用PICP算法并不能實(shí)現(xiàn)兩個(gè)顱骨的正確配準(zhǔn)。這是由于PICP算法對(duì)兩個(gè)待配準(zhǔn)顱骨的初始位置要求較高,一般需要先進(jìn)行粗配準(zhǔn),再采用PICP算法進(jìn)行細(xì)配準(zhǔn),這樣才能達(dá)到精確的配準(zhǔn)結(jié)果。而可信區(qū)域配準(zhǔn)能夠?qū)蓚€(gè)顱骨進(jìn)行配準(zhǔn),但是跟本文顱骨配準(zhǔn)算法相比,可信區(qū)域配準(zhǔn)的配準(zhǔn)結(jié)果的精度有所降低。
對(duì)于未知顱骨U和參考顱骨S1~S4,分別采用PICP算法、可信區(qū)域配準(zhǔn)以及本文配準(zhǔn)算法進(jìn)行配準(zhǔn),其配準(zhǔn)結(jié)果如表1所示。
從表1的配準(zhǔn)結(jié)果來看,本文配準(zhǔn)算法的配準(zhǔn)精度最高、耗時(shí)最短。跟PCIP算法相比,本文算法的配準(zhǔn)精度和耗時(shí)分別提高了約60%和50%,跟可信區(qū)域配準(zhǔn)相比,本文算法的配準(zhǔn)精度和耗時(shí)分別提高了約30%和30%。因此,本文提出的基于孔洞輪廓線的配準(zhǔn)方法是一種速度更快、精度更高的三維顱骨模型配準(zhǔn)方法。
表1 算法的配準(zhǔn)結(jié)果Tab.1 Registration result of the algorithms
顱骨配準(zhǔn)是計(jì)算機(jī)輔助顱面復(fù)原的重要研究內(nèi)容之一。由于顱骨數(shù)據(jù)量大、含噪聲多、結(jié)構(gòu)復(fù)雜,因此對(duì)配準(zhǔn)精度要求較高。鑒于此,本文提出了基于孔洞輪廓線的三維顱骨模型配準(zhǔn)方法。首先,提取顱骨的孔洞輪廓線;然后,采用輪廓線的特征串實(shí)現(xiàn)顱骨輪廓線的匹配;最后,采用PICP算法將顱骨進(jìn)行進(jìn)一步的細(xì)配準(zhǔn),從而實(shí)現(xiàn)顱骨精確配準(zhǔn)的目的。雖然該方法達(dá)到了較高的顱骨配準(zhǔn)精度,但是不能正確實(shí)現(xiàn)孔洞缺失的顱骨的配準(zhǔn),因?yàn)樵摲椒ㄒ髢蓚€(gè)待配準(zhǔn)的顱骨孔洞輪廓線之間存在一一對(duì)應(yīng)的關(guān)系。在今后的研究中,要進(jìn)一步解決缺失顱骨的配準(zhǔn)問題,并將顱骨配準(zhǔn)結(jié)果應(yīng)用到顱面復(fù)原的研究中,提高顱面復(fù)原技術(shù)在人類考古、刑事案件偵破等領(lǐng)域的應(yīng)用價(jià)值。