結(jié)合分區(qū)和改進(jìn)ICP的三維顱骨自動配準(zhǔn)算法

史重陽,劉曉寧,羅星海,胡曉靜,耿國華

(西北大學(xué) 信息科學(xué)與技術(shù)學(xué)院,西安 710127) E-mail:xnliu@nwu.edu.cn

1 引 言

顱面復(fù)原是對未知顱骨根據(jù)統(tǒng)計(jì)先驗(yàn)知識恢復(fù)其生前面貌的一種技術(shù)[1].為此需要建立科學(xué)、有效的參考顱骨與所有樣本顱骨形狀關(guān)系的統(tǒng)計(jì)計(jì)算模型,即顱骨形態(tài)學(xué)統(tǒng)計(jì).為了對顱骨形態(tài)進(jìn)行統(tǒng)計(jì)分析,在建立顱骨數(shù)據(jù)庫時(shí),需要對所有樣本顱骨點(diǎn)云數(shù)據(jù)與參考顱骨進(jìn)行配準(zhǔn),將它們統(tǒng)一到一個(gè)坐標(biāo)系下,從而對顱骨形態(tài)進(jìn)一步分析.此外,面對一個(gè)未知顱骨,基于配準(zhǔn)的顱骨復(fù)原方法也需要從顱骨數(shù)據(jù)庫中檢索出最相似的顱骨.因此,顱骨配準(zhǔn)是顱面復(fù)原的一個(gè)關(guān)鍵步驟.

由于獲取的顱骨是三維的,因此三維顱骨配準(zhǔn)屬于三維模型配準(zhǔn)的具體應(yīng)用.目前,常用的三維配準(zhǔn)算法是由Besl等提出的迭代最近點(diǎn)(Iterated Closest Point,ICP)算法,該方法本質(zhì)是基于最小二乘法的最優(yōu)配準(zhǔn)方法,運(yùn)用Euclidean距離計(jì)算得到與特征點(diǎn)最近的匹配點(diǎn),估算最優(yōu)的變換矩陣參數(shù),直到滿足正確配準(zhǔn)的收斂精度要求[2].但是該算法對初始位置和匹配點(diǎn)對的要求比較高,容易形成局部迭代收斂,造成迭代收斂速度低.為此,國內(nèi)外學(xué)者對此算法進(jìn)行了改進(jìn).文獻(xiàn)[3]提出了一種基于點(diǎn)到平面的改進(jìn)ICP算法,改進(jìn)了ICP算法中點(diǎn)對點(diǎn)歐氏距離的評價(jià)函數(shù).文獻(xiàn)[4]提出了對點(diǎn)云對應(yīng)點(diǎn)賦予權(quán)重,剔除權(quán)重大于閾值的匹配點(diǎn)對,通過對目標(biāo)函數(shù)引入M-估計(jì)(M-estimation),剔除異常點(diǎn),有效的確定了初始位置和匹配點(diǎn)對.liu等提出了基于遺傳(Genetic Algorithm,GA)算法和ICP的配準(zhǔn)算法[5],GA是一種全局的搜索算法,很好的解決了ICP容易形成局部迭代的問題,但是GA算法時(shí)間復(fù)雜度較大.KO等提出基于曲面幾何特征性質(zhì)的配準(zhǔn)算法[6],與迭代配準(zhǔn)算法相比,該算法的配準(zhǔn)時(shí)耗明顯降低.文獻(xiàn)[7]提出采用幾何哈希方法找出使得最多數(shù)量的點(diǎn)對法矢一致的變換,運(yùn)用該變換對散亂點(diǎn)云作初始配準(zhǔn),再運(yùn)用ICP算法進(jìn)行精確配準(zhǔn).文獻(xiàn)[8]提出了一種基于曲率特征的配準(zhǔn)算法,該方法根據(jù)曲率相似度約束進(jìn)行配準(zhǔn).文獻(xiàn)[9,10]提出一種基于有界旋轉(zhuǎn)角的配準(zhǔn)算法,該方法在ICP算法的基礎(chǔ)上通過設(shè)置旋轉(zhuǎn)角度的邊界來改進(jìn)算法,從而提高了配準(zhǔn)的準(zhǔn)確度.

也有很多學(xué)者對三維顱骨配準(zhǔn)進(jìn)行了研究.Du等提出了一種基于ICP和薄板樣條函數(shù)(Thin-plate Spline,TPS)的三維顱骨配準(zhǔn)算法[11],該方法引入積分不變量的多尺度約束,選定最佳匹配點(diǎn),然后采用ICP實(shí)現(xiàn)初始配準(zhǔn),利用TPS實(shí)現(xiàn)精確配準(zhǔn).該方法提高了配準(zhǔn)的準(zhǔn)確度,但是ICP和TPS都是迭代算法,因此增加了配準(zhǔn)的時(shí)間復(fù)雜度.He等提出了以TPS算法實(shí)現(xiàn)顱骨的初始配準(zhǔn),然后根據(jù)歐氏距離與局部幾何特征的加權(quán)距離實(shí)現(xiàn)顱骨的精確配準(zhǔn)的自動配準(zhǔn)算法[12].M.Berar等提出了一種新的3D網(wǎng)格彈性配準(zhǔn)方法[13],該方法是基于多分辨率,對定義頂點(diǎn)的兩個(gè)網(wǎng)格實(shí)現(xiàn)對稱距離最小化.該方法有效的提高了配準(zhǔn)的準(zhǔn)確度,但是對于多頂點(diǎn)的網(wǎng)格模型,實(shí)現(xiàn)對稱距離最小化,使得計(jì)算量過大,從而增加了配準(zhǔn)的時(shí)耗.Yusuf Sahillioglu等提出了一種基于三維立體形狀體積的配準(zhǔn)算法[14],該方法根據(jù)顱骨四面體網(wǎng)格模型,獲取數(shù)據(jù)的體積性質(zhì)進(jìn)行配準(zhǔn).與網(wǎng)格模型相比,能更好的捕獲數(shù)據(jù)的體積性質(zhì),利用體積性質(zhì)進(jìn)行配準(zhǔn),對于畸形的顱骨具有一定的魯棒性.但是該方法在網(wǎng)格劃分的過程中,會丟失一些數(shù)據(jù),從而降低配準(zhǔn)的準(zhǔn)確度.

綜上所述,目前顱骨配準(zhǔn)具有配準(zhǔn)效率不高,配準(zhǔn)準(zhǔn)確度低,對缺損或者畸形顱骨沒有一定的魯棒性等問題.本文針對配準(zhǔn)效率和準(zhǔn)確度不高等問題,提出一種在區(qū)域劃分的基礎(chǔ)上,引入動態(tài)估計(jì)(Destimation)的ICP配準(zhǔn)算法.該算法避免的點(diǎn)云數(shù)據(jù)量過大而導(dǎo)致的計(jì)算量過大的問題,為精確配準(zhǔn)階段提供了良好的初始值,根據(jù)引入的Destimation有效剔除了誤匹配點(diǎn)對,降低了匹配點(diǎn)對的數(shù)量,有效提高了配準(zhǔn)的準(zhǔn)確度以及效率,使得迭代收斂性更快.

2 RANSAC算法

隨機(jī)抽樣一致(Random Sample Consensus,RANSAC)算法思想是從一組包含局外點(diǎn)的觀測數(shù)據(jù)集中,通過迭代的方法計(jì)算出數(shù)據(jù)的數(shù)學(xué)模型參數(shù),產(chǎn)生最大一致性數(shù)據(jù)集.該算法的基本假設(shè)是樣本中包含正確數(shù)據(jù)(inliers,可以被模型描述的數(shù)據(jù)),也包含異常數(shù)據(jù)(outliers,偏離正常范圍很遠(yuǎn),無法自適應(yīng)數(shù)學(xué)模型的數(shù)據(jù)),即數(shù)據(jù)集中含有噪聲.這些異常數(shù)據(jù)產(chǎn)生的原因是錯(cuò)誤的測量、假設(shè)或計(jì)算等.本文數(shù)據(jù)集中的異常數(shù)據(jù)是在初始配準(zhǔn)中尋找匹配點(diǎn)對時(shí)產(chǎn)生的誤匹配點(diǎn)對.該算法通過反復(fù)選擇樣本數(shù)據(jù)中的一組隨機(jī)子集來達(dá)成目標(biāo)[15].

RANSAC基本思想描述如下:

1)存在一個(gè)最小抽樣集,其個(gè)數(shù)為n(n為初始化模型參數(shù)所需的最小樣本數(shù))和一個(gè)樣本集N,其中|N|>n,從N中隨機(jī)抽取n個(gè)樣本,構(gòu)成N的子集S,用來初始化模型W.

2)余集SC=NS中與模型W的誤差小于某一設(shè)定閾值ε的樣本集合和集合S構(gòu)成集合S′.S′是內(nèi)點(diǎn)集,它們構(gòu)成S的一致集.

3)若內(nèi)點(diǎn)集的個(gè)數(shù)大于樣本集N的個(gè)數(shù),則得到正確的模型參數(shù),并根據(jù)內(nèi)點(diǎn)集S′,采用最小二乘等方法重新設(shè)計(jì)新的模型W′.重新隨機(jī)抽取新的樣本子集S,否則返回步驟1)進(jìn)行迭代.

4)在完成一定的抽樣次數(shù)后,若未找到一致集,則算法失敗,否則選取抽樣后得到的最大一致集判斷內(nèi)外點(diǎn),算法結(jié)束.

3 顱骨區(qū)域劃分

根據(jù)基于Euclidean距離的Voronoi圖對待匹配點(diǎn)集P和參考點(diǎn)集Q進(jìn)行區(qū)域劃分.首先對顱骨數(shù)據(jù)進(jìn)行隨機(jī)均勻采樣,得到采樣點(diǎn)集{v1,v2,v3,…,vn}和{u1,u2,u3,…,um}.然后以采樣點(diǎn)為初始聚類點(diǎn),得到每個(gè)采樣點(diǎn)的劃分區(qū)域.如公式(1)所示:

(1)

其中,vi屬于P,uj屬于Q;n表示待匹配點(diǎn)集P得到的采樣點(diǎn)集元素的個(gè)數(shù),即P的區(qū)域數(shù)目;m表示參考點(diǎn)集Q得到的采樣點(diǎn)集元素的個(gè)數(shù),即Q的區(qū)域數(shù)目.區(qū)域的數(shù)目由顱骨的五官和點(diǎn)云數(shù)目決定,在實(shí)際應(yīng)用中通常設(shè)置為1～9.當(dāng)顱骨點(diǎn)云數(shù)目比較多時(shí),區(qū)域數(shù)目需要設(shè)置大一些.當(dāng)區(qū)域數(shù)目取1時(shí),算法退化為全局配準(zhǔn)算法.

4 顱骨自動配準(zhǔn)算法

經(jīng)典ICP算法對兩個(gè)輸入顱骨數(shù)據(jù)的初始相對位置要求較高,若兩個(gè)數(shù)據(jù)相對位置較大則易陷入局部最優(yōu),得到錯(cuò)誤的配準(zhǔn)結(jié)果.算法時(shí)間復(fù)雜度較大,為O(NpNq),其中Np和Nq分別是待匹配點(diǎn)集和參考點(diǎn)集的點(diǎn)云個(gè)數(shù).該算法中使用Euclidean距離作為唯一約束,約束條件單一,配準(zhǔn)準(zhǔn)確度不高.

本文提出的算法在初始配準(zhǔn)中,運(yùn)用區(qū)域劃分思想,有效縮小配準(zhǔn)規(guī)模,極大緩解了時(shí)間復(fù)雜度大的問題,為精確配準(zhǔn)階段提供了有效的初始值,本文借助的RANSAC算法,對噪聲具有較高的魯棒性,提高了配準(zhǔn)的準(zhǔn)確度[16].在精確配準(zhǔn)中,采用kd樹(kd-tree)算法加快搜索速度,引入Destimation來剔除誤匹配點(diǎn)對,有效提高了配準(zhǔn)的準(zhǔn)確度,降低了配準(zhǔn)的時(shí)耗.

4.1 初始配準(zhǔn)

4.1.1 區(qū)域配準(zhǔn)

兩個(gè)區(qū)域的配準(zhǔn)是兩個(gè)更小規(guī)模的點(diǎn)云數(shù)據(jù)的配準(zhǔn).首先計(jì)算區(qū)域Vk和區(qū)域Ut的質(zhì)心,如公式(2)所示:

(2)

然后,根據(jù)Euclidean距離相似度和K近鄰(k-Nearest Neighbor,KNN)算法,對每一個(gè)區(qū)域Vk,在待匹配區(qū)域Ut中搜索與Vk的質(zhì)心之間Euclidean距離最近的K個(gè)相似區(qū)域Ua(a=1,2,3,…,K),如公式(3)所示:

‖o(pk)-o(qt)‖≤ξ1

(3)

其中,ξ1為大于0的實(shí)數(shù).

(4)

(5)

其中,ξ2為大于0的實(shí)數(shù).

(6)

圖1(a)、(b)中左圖為顱骨掃描數(shù)據(jù)集中待匹配點(diǎn)集的鼻子和下顎區(qū)域,圖1(a)、(b)中右圖表示鼻子和下顎區(qū)域與參考點(diǎn)集相似區(qū)域的配準(zhǔn)結(jié)果,區(qū)域配準(zhǔn)的匹配程度對應(yīng)φ的取值.可以看出,對于2個(gè)部分重疊區(qū)域進(jìn)行配準(zhǔn)能夠有效的得到精確的配準(zhǔn)結(jié)果.

4.1.2 求解全局剛體變換

(7)

(8)

圖1 區(qū)域配準(zhǔn)結(jié)果示意圖Fig.1 Result of Region registration

E(ω)=ωTAω
‖ω‖2=1,ω≥0

(9)

為了得到全局剛體變換矩陣T,首先假設(shè)將T分解為一個(gè)旋轉(zhuǎn)矩陣R和一個(gè)平移矩陣t,然后將線性變換矩陣T′分解為一個(gè)旋轉(zhuǎn)矩陣R′和平移矩陣t′.令t=t′,根據(jù)四元素法求得旋轉(zhuǎn)矩陣R.最后,用T對P和Q完成初始配準(zhǔn).

4.2 引入動態(tài)估計(jì)的ICP算法

在點(diǎn)云數(shù)據(jù)配準(zhǔn)的研究中,經(jīng)典ICP算法是最基礎(chǔ)、最廣泛使用的配準(zhǔn)算法[19,20].該算法針對點(diǎn)集中的每個(gè)點(diǎn),在目標(biāo)點(diǎn)集中全局搜索與其Euclidean距離最近的匹配點(diǎn),然后估計(jì)變換矩陣參數(shù),反復(fù)迭代,直到迭代的目標(biāo)函數(shù)值小于等于預(yù)設(shè)的閾值為止.

在經(jīng)典ICP算法中,利用Euclidean距離計(jì)算每個(gè)點(diǎn)的最近點(diǎn),容易出現(xiàn)誤匹配點(diǎn)對,并且匹配點(diǎn)對數(shù)量大,這是配準(zhǔn)準(zhǔn)確度低和配準(zhǔn)時(shí)耗高的瓶頸所在.如果提前刪除誤配的匹配點(diǎn)對,在迭代過程中將顯著的減少算法的計(jì)算量,有效提高配準(zhǔn)的準(zhǔn)確度,降低了配準(zhǔn)的時(shí)間復(fù)雜度.

在統(tǒng)計(jì)學(xué)中,均差和標(biāo)準(zhǔn)差是判定樣本的有效性和無偏性的標(biāo)準(zhǔn),因此本文提出Destimation思想,匹配點(diǎn)對之間距離在一定范圍內(nèi)與均差越近似,正確率就越大,由此標(biāo)準(zhǔn)來有效剔除誤配或錯(cuò)配的匹配點(diǎn)對.計(jì)算均差μd和標(biāo)準(zhǔn)差δd如公式(10)所示:

(10)

其中,di表示匹配點(diǎn)對之間的距離‖pi-qj‖,Np表示根據(jù)kd-tree算法求得的匹配點(diǎn)對數(shù).

根據(jù)(10)式計(jì)算Destimation,如公式(11)所示:

(11)

其中,λ為動態(tài)參數(shù),為大于0的實(shí)數(shù).

根據(jù)(11)對每個(gè)匹配點(diǎn)對設(shè)置一個(gè)動態(tài)估計(jì)Destimation.根據(jù)給定閾值ξ3,如果匹配點(diǎn)對的Destimation<ξ3,則刪除該匹配點(diǎn)對,ξ3為大于0的實(shí)數(shù).

本文采用均方根誤差RMSE作為配準(zhǔn)誤差[21].此度量方法直觀、嚴(yán)謹(jǐn),有效的彌補(bǔ)了經(jīng)典ICP誤差函數(shù)的缺陷.如公式(12)所示:

(12)

其中Np表示待匹配點(diǎn)集P的元素個(gè)數(shù),即匹配點(diǎn)對數(shù).pi為待匹配點(diǎn)集P中的點(diǎn),qj為參考點(diǎn)集Q中與pi相匹配的匹配點(diǎn).

4.3 三維顱骨配準(zhǔn)算法

本文將顱骨配準(zhǔn)算法分兩階段進(jìn)行:

1)初始配準(zhǔn).給定待匹配點(diǎn)集P和參考點(diǎn)集Q,初始配準(zhǔn)步驟如下:

Step5.根據(jù)E(ω)求得最優(yōu)變換矩陣的線性組合T′,易得平移矩陣t,由四素法求得旋轉(zhuǎn)矩陣R,繼而得到剛性變換矩陣T,完成初始配準(zhǔn).

2)精確配準(zhǔn).輸入初始配準(zhǔn)的結(jié)果,初始化迭代次數(shù)η=0.步驟如下:

Step2.由公式(11)得到Destimation,刪除Dη(i)中Destimation<ξ3的匹配點(diǎn)對.

Step3.根據(jù)四元素法計(jì)算得到旋轉(zhuǎn)矩陣R和平移向量T.

Step4.根據(jù)R,T更新Dη(i),即Pη+1=RηPη+Tη,Dη+1(i)=qj,由公式(12)得到誤差函數(shù)RMSEη.

Step5.計(jì)算2次迭代之間的估計(jì)誤差,如果|RMSEη+1-RMSEη|≤χ,就終止迭代;否則返回步驟(1)進(jìn)行迭代.

5 實(shí)驗(yàn)結(jié)果與分析

本文實(shí)驗(yàn)采用的是西北大學(xué)可視化研究所采集的顱面樣本數(shù)據(jù).為了驗(yàn)證本文算法對屬性相似顱骨的配準(zhǔn)準(zhǔn)確度以及效率,本文用文獻(xiàn)[22]提到的方法,在顱骨數(shù)據(jù)庫中檢索到38套屬性相近的顱骨組成測試數(shù)據(jù)集.算法中基本參數(shù)設(shè)置如下:ξ1=0.05,ξ2=0.02,ξ3=0.01.

實(shí)驗(yàn)采用的PC配置為3.50GHz、內(nèi)存為8GB,操作系統(tǒng)為Windows 7.在Matlab R2014a平臺上實(shí)現(xiàn)了結(jié)合顱骨分區(qū)的初始配準(zhǔn)算法和改進(jìn)的ICP算法.

5.1 初始配準(zhǔn)結(jié)果分析

本文對基于區(qū)域的自動配準(zhǔn)、基于曲率圖的配準(zhǔn)和本文算法三種算法進(jìn)行初始配準(zhǔn)對比實(shí)驗(yàn). 其中2組不同點(diǎn)云數(shù)據(jù)實(shí)驗(yàn)結(jié)果如圖2和圖3所示. 圖2為第一組顱骨點(diǎn)云配準(zhǔn)結(jié)果,圖2(a)左側(cè)為參考顱骨(13201個(gè)點(diǎn)),右側(cè)為待匹配顱骨(11282個(gè)點(diǎn)).圖3為第二組顱骨點(diǎn)云配準(zhǔn)結(jié)果,圖3(a)左側(cè)為參考顱骨(102905個(gè)點(diǎn)),右側(cè)為待匹配顱骨(92837個(gè)點(diǎn)).實(shí)驗(yàn)分析可知,本文算法與其他兩種算法相比,本文初始配準(zhǔn)獲得的初始結(jié)果比較好.

圖2 第一組數(shù)據(jù)初始配準(zhǔn)結(jié)果Fig.2 Initial registration results of the first set of data

圖3 第二組數(shù)據(jù)初始配準(zhǔn)結(jié)果Fig.3 Initial registration results of the second set of data

將本文中所提的RMSE作為配準(zhǔn)準(zhǔn)確度的評判標(biāo)準(zhǔn).表1為圖2初始配準(zhǔn)結(jié)果的配準(zhǔn)時(shí)耗及準(zhǔn)確度分析,由表1和表2可知,基于區(qū)域的自配準(zhǔn)算法和基于曲率圖配準(zhǔn)算法得到的初始配準(zhǔn)誤差比較大,第一組數(shù)據(jù)RMSE分別為0.3564mm和0.4372mm,第二組數(shù)據(jù)誤差分別為0.3950mm和0.4028mm.而本文算法的到的初始配準(zhǔn)結(jié)果比較好,第一組和第二組RMSE僅為0.1942mm和0.2184mm.基于區(qū)域的自動配準(zhǔn)算法初始配準(zhǔn)階段,需要將每一個(gè)區(qū)域與剩余所有區(qū)域進(jìn)行初始配準(zhǔn);基于曲率圖的配準(zhǔn)算法初始配準(zhǔn)階段,在應(yīng)用曲率相似度求解匹配點(diǎn)對的基礎(chǔ)上,再計(jì)算每個(gè)點(diǎn)的曲率圖;而本文算法初始配準(zhǔn)階段應(yīng)用Euclidean距離相似度和KNN算法求解相似區(qū)域與匹配點(diǎn)對,只需要將一個(gè)區(qū)域與相似的幾個(gè)區(qū)域進(jìn)行配準(zhǔn)即可.因此,本文算法與其他兩種算法相比,有效的降低了初始配準(zhǔn)階段的時(shí)間復(fù)雜度.

表1 第一組數(shù)據(jù)初始配準(zhǔn)算法對比
Table 1 Comparison of the first set of data registration algorithm

算法 配準(zhǔn)時(shí)耗/sRMSE/mm本文初始配準(zhǔn)算法13.3050.1942基于區(qū)域初始配準(zhǔn)36.8270.3564基于曲率圖初始配準(zhǔn)53.2930.4372

表2 第二組數(shù)據(jù)初始配準(zhǔn)算法對比
Table 2 Comparison of the second set of data registration algorithm

算法 配準(zhǔn)時(shí)耗/sRMSE/mm本文初始配準(zhǔn)算法31.9380.2184基于區(qū)域初始配準(zhǔn)69.5710.3950基于曲率圖初始配準(zhǔn)74.9050.4028

5.2 精確配準(zhǔn)結(jié)果分析

在初始配準(zhǔn)算法的基礎(chǔ)上,對基于區(qū)域中的稀疏ICP算法、基于曲率圖中的經(jīng)典ICP算法和本文算法的配準(zhǔn)準(zhǔn)確度和配準(zhǔn)時(shí)耗,做了對比實(shí)驗(yàn),實(shí)驗(yàn)結(jié)果如圖4-圖5所示.

圖4 第一組數(shù)據(jù)精確配準(zhǔn)結(jié)果對比Fig.4 Comparison of accurate registration results from the first set of data

由圖4(a)和圖5(a)分析可知,基于曲率圖中的經(jīng)典ICP算法,可以大幅度提高配準(zhǔn)的準(zhǔn)確度,但是該算法的配準(zhǔn)結(jié)果出現(xiàn)了局部最優(yōu)的現(xiàn)象.由圖4(b)和圖5(b)顯示,基于區(qū)域中的稀疏ICP算法的精確配準(zhǔn),大部分區(qū)域有效的提高了配準(zhǔn)的準(zhǔn)確度,但是部分區(qū)域(下顎、臉頰)出現(xiàn)了誤配或錯(cuò)配的現(xiàn)象.

圖4(c)顯示,在Destimation不同取值的情況下,當(dāng)λ=3時(shí),配準(zhǔn)的準(zhǔn)確度平均達(dá)到了89.83%;當(dāng)λ=2時(shí),配準(zhǔn)的準(zhǔn)確度平均達(dá)到了82.71%,額頭和臉頰部分出現(xiàn)錯(cuò)配;當(dāng)λ=1時(shí),配準(zhǔn)的準(zhǔn)確度平均達(dá)到了73.63%,但是基于區(qū)域中的稀疏ICP與本文算法相比,配準(zhǔn)效果更好,由圖可知,眼睛和下顎部分出現(xiàn)錯(cuò)配以及誤配.由圖4(c)可知,在λ=3時(shí),本文算法與基于區(qū)域中的稀疏ICP和基于曲率圖中的經(jīng)典ICP算法相比,配準(zhǔn)的準(zhǔn)確度平均提高了15.26%和19.13%.圖5(c)顯示,在Destimation不同取值的情況下,當(dāng)λ=3時(shí),配準(zhǔn)的準(zhǔn)確度平均達(dá)到了91.70%;當(dāng)λ=2時(shí),配準(zhǔn)的準(zhǔn)確度平均達(dá)到了80.95%,下顎區(qū)域出現(xiàn)錯(cuò)配;當(dāng)λ=1時(shí),配準(zhǔn)的準(zhǔn)確度平均達(dá)到了72.43%,右眼及下顎部分出現(xiàn)錯(cuò)配以及誤配.由圖可知,在λ=3時(shí),本文算法與基于區(qū)域中的稀疏ICP和基于曲率圖中的經(jīng)典ICP算法相比,配準(zhǔn)的準(zhǔn)確度平均提高了16.82%和20.01%.由此說明本文提出的算法是有效可行的,在Destimation不同取值情況下,靈活有效的提高了配準(zhǔn)的準(zhǔn)確度,降低了配準(zhǔn)時(shí)間復(fù)雜度.

圖5 第二組數(shù)據(jù)精確配準(zhǔn)結(jié)果對比Fig.5 Comparison of accurate registration results from the second set of data

表3和表4分別為圖4和圖5的配準(zhǔn)時(shí)耗和配準(zhǔn)的RMSE,由表3和表4可知,基于曲率圖中的經(jīng)典ICP的配準(zhǔn)算法耗時(shí)較大,配準(zhǔn)的RMSE分別為0.0829mm和0.0764mm.基于區(qū)域中的稀疏ICP,配準(zhǔn)時(shí)耗與本文算法相比,配準(zhǔn)時(shí)耗比較大,配準(zhǔn)準(zhǔn)確度比本文算法λ=1時(shí)好,但當(dāng)λ=2和λ=3時(shí),本文算法比基于區(qū)域中的稀疏ICP好.因此可知,與前兩種算法相比,本文配準(zhǔn)效率整體得到了進(jìn)一步提高,本文提出的算法有效的提高配準(zhǔn)的準(zhǔn)確度.

為了驗(yàn)證本文算法的有效性和可行性,在初始配準(zhǔn)的基礎(chǔ)上,對基于曲率圖中的經(jīng)典ICP算法、基于區(qū)域中的稀疏ICP和本文改進(jìn)ICP算法(λ=3)三種算法,在不同迭代次數(shù)下進(jìn)行了誤差統(tǒng)計(jì)分析.實(shí)驗(yàn)結(jié)果如圖6所示.

表3 第一組精確配準(zhǔn)算法對比
Table 3 Comparison of the first set of data Accurate registration algorithms

算法RMSE/mm配準(zhǔn)時(shí)耗/s基于曲率圖經(jīng)典ICP0.0829135.574基于區(qū)域稀疏ICP0.0453105.570λ=30.0237本文算法λ=20.042671.968λ=10.0489

表4 第二組精確配準(zhǔn)算法對比
Table 4 Comparison of the second set of data Accurate registration algorithms

算法RMSE/mm配準(zhǔn)時(shí)耗/s基于曲率圖經(jīng)典ICP0.0764173.293基于區(qū)域稀疏ICP0.0372146.827λ=30.0261本文算法λ=20.0354λ=10.0493109.736

由圖6看出,當(dāng)iteration>35時(shí),本文算法的RMSE趨于穩(wěn)定,而基于區(qū)域中的稀疏ICP算法和基于曲率圖中的經(jīng)典ICP算法,當(dāng)iteration≥45時(shí)才趨于穩(wěn)定.因此,本文改進(jìn)的ICP算法與稀疏ICP算法和基于曲率圖中的經(jīng)典ICP算法相比,收斂速度更快,配準(zhǔn)的時(shí)間復(fù)雜度大大減小.當(dāng)iteration=45時(shí),本文改進(jìn)的ICP算法的配準(zhǔn)誤差比基于區(qū)域中的稀疏ICP算法提高了很多,與基于曲率圖中的經(jīng)典ICP算法相比,配準(zhǔn)的準(zhǔn)確度大大提高.說明本文用提出的Destimation來剔除誤匹配點(diǎn)對的方法有效可行.本文算法不僅有效降低了配準(zhǔn)的時(shí)間復(fù)雜度,還提高了配準(zhǔn)的準(zhǔn)確度.

圖6 不同迭代次數(shù)下的RMSE分析Fig.6 RMSE analysis of different iteration times

為了更加直觀的對比這三種方法,表5顯示,在38套屬性相近的顱骨測試數(shù)據(jù)集上,分別采用基于區(qū)域配準(zhǔn)算法、基于曲率圖的配準(zhǔn)算法和本文算法進(jìn)行配準(zhǔn),并根據(jù)實(shí)驗(yàn)結(jié)果,統(tǒng)計(jì)三種算法的平均配準(zhǔn)誤差以及平均配準(zhǔn)時(shí)耗.

分析圖4-圖5和表5可知,基于曲率圖的的配準(zhǔn)算法尋找匹配點(diǎn)對的過程比較簡單,此方法主要是計(jì)算曲率,根據(jù)每個(gè)點(diǎn)的曲率值構(gòu)造曲率圖,并根據(jù)曲率圖尋找匹配點(diǎn),即需要計(jì)算顱骨模型中的所有點(diǎn)的曲率構(gòu)造曲率圖以及待匹配顱骨中點(diǎn)到參考顱骨中所有點(diǎn)的曲率圖的加權(quán)距離,然后取距離最短的m個(gè)點(diǎn)作為潛在匹配點(diǎn),這個(gè)計(jì)算過程十分費(fèi)時(shí),大大增加了算法的時(shí)間復(fù)雜度,又由于幾何約束條件太少,導(dǎo)致此方法配準(zhǔn)的準(zhǔn)確度降低,該方法精確配準(zhǔn)部分應(yīng)用經(jīng)典ICP算法,即采用基于Euclidean距離的最近點(diǎn)對應(yīng)方法尋找匹配點(diǎn),而且對匹配點(diǎn)對沒有進(jìn)行篩選,容易出現(xiàn)錯(cuò)配現(xiàn)象,導(dǎo)致配準(zhǔn)準(zhǔn)確度大大降低.采用基于區(qū)域的自動配準(zhǔn)算法,該方法首先需要將顱骨樣本進(jìn)行區(qū)域劃分,然后對待匹配顱骨

表5 不同配準(zhǔn)算法的平均誤差
Table 5 Average error of different registration algorithms

算法平均RMSE/mm平均時(shí)耗/s基于曲率圖0.0619116.071基于區(qū)域0.0398142.619λ=30.0196本文算法λ=20.0366109.736λ=10.0417

中的區(qū)域與參考顱骨中所有區(qū)域進(jìn)行配準(zhǔn),這個(gè)過程計(jì)算復(fù)雜,使得配準(zhǔn)的時(shí)間復(fù)雜度大大增加,精確配準(zhǔn)階段采用稀疏ICP算法,該算法采用稀疏誘導(dǎo)范式代替經(jīng)典ICP算法中的二范式,雖然使得匹配點(diǎn)對的準(zhǔn)確度增加,但是依舊沒有進(jìn)行誤匹配點(diǎn)對的篩選,而且還需要計(jì)算待匹配顱骨中的點(diǎn)到參考顱骨中所有點(diǎn)的稀疏誘導(dǎo)范式,這個(gè)計(jì)算過程,大大增加了配準(zhǔn)的時(shí)間復(fù)雜度,降低了配準(zhǔn)的準(zhǔn)確度.本文算法采用Euclidean距離相似度等幾何特征作為約束條件,根據(jù)KNN算法以區(qū)域質(zhì)心作為配準(zhǔn)單位減少待匹配顱骨中區(qū)域與參考顱骨區(qū)域配準(zhǔn)的區(qū)域數(shù)目,減少了計(jì)算量,大大提高配準(zhǔn)效率,精確配準(zhǔn)階段,采用kd-tree算法,有效的縮小了待匹配顱骨在參考顱骨中尋找匹配點(diǎn)的范圍,并且提出Destimation對匹配點(diǎn)對進(jìn)行篩選,有效剔除了誤匹配點(diǎn)對,減少了計(jì)算量,降低了配準(zhǔn)的時(shí)間復(fù)雜度,提高了配準(zhǔn)的準(zhǔn)確度.實(shí)驗(yàn)表明,分別采用這三種方法進(jìn)顱骨配準(zhǔn),本文算法的配準(zhǔn)時(shí)耗是基于曲率圖配準(zhǔn)算法的20%,是基于區(qū)域配準(zhǔn)算法的16%,配準(zhǔn)配準(zhǔn)的準(zhǔn)確度是其他兩種方法的2倍之余.

6 總 結(jié)

本文提出了結(jié)合分區(qū)和改進(jìn)ICP的顱骨自動配準(zhǔn)算法.首先在初始配準(zhǔn)階段,進(jìn)行顱骨分區(qū),有效避免了特征點(diǎn)的標(biāo)定以及曲率等幾何特征的計(jì)算,并運(yùn)用基于質(zhì)心的區(qū)域匹配得到相似的區(qū)域?qū)?運(yùn)用RANSAC算法得到最優(yōu)變換矩陣.然后進(jìn)行精確配準(zhǔn),初始配準(zhǔn)為此階段提供了良好的初始值,避免了算法陷入局部最優(yōu),運(yùn)用kd-tree算法尋找匹配點(diǎn),并提出Destimation來剔除誤配的匹配點(diǎn)對.實(shí)驗(yàn)表明,本文算法能有效緩解顱骨配準(zhǔn)過程中出現(xiàn)的效率低和準(zhǔn)確度不高的問題,滿足顱面復(fù)原中的顱骨形態(tài)統(tǒng)計(jì)要求.

[1] Shui Wu-yang,Zhou Ming-quan,Wu Zhong-ke,et al.An approach of craniofacial reconstruction based on registration [J].Journal of Computer-Aided Design & Computer Graphics(CAD&CG),2011,23(4):607-614.

[2] Besl P J,Mc Kay N D.A method for registration of 3-D shapes[J].IEEE Transactions on Pattern Analysis and Machine Intelligence(PAMI),1992,14(2):239-256.

[3] Medioni G,Chen Y.Object modeling by registration of multiple range images[J].Image & Vision Computing,1992,10(3):145-155.

[4] Wang Xin,Zhang Ming-ming,Yu Xiao,et al.Point cloud registration based on improved iterative closest point method[J].Opticsand and Precision Engineering,2012,20(9):2068-2077.

[5] Liu Chun-guo,Liu Chang,An Bai-ling.Three dimensional surface registration on the basis of genetic algorithm[J].China Metal Forming Equipment & Manufacturing Technology,(CMET),2009,44(4):110-113.

[6] Ko K H,Maekawa T,Patrikalakis N M.An algorithm for optimal free-form object matching [J].Computer-Aided Design,2003,35(10):913 -923.

[7] Zhu Yan-juan,Zhou Lai-shui,Zhang Li-yan.Registration of scattered cloud data[J].Journal of Computer-Aided Design & Computer Graphics(CAD&CG),2006,18(4):475-481.

[8] Ge Bao-zhen,Peng Bo,Tian Qing-guo.Registration of three-dimensional point-cloud data based on curvature map [J].Journal of Tianjin University,2013,46(2):174-180.

[9] Zhao Fu-qun,Zhou Ming-quan,Geng Guo-hua.Point cloud registration algorithm based on bounded rotation angle[J].Microelectronics & Computer,2017,34(3):46-49.

[10] Zhang C,Du S,Liu J,et al.Robust 3D point set egistration using iterative closest point algorithm with bounded rotation angle[J].Signal Processing,2016,120(C):777-788.

[11] Du H,Geng G,Li K,et al.3D skull registration based on registration points automatic correspondence[C].International Conference on Virtual Reality and Visualization(ICVRV),IEEE,2013,8768 (12):293-296.

[12] He Yi-yue,Ma Zi-ping,Gao Ni,et al.Local craniofacial morphological relationship modeling based on physiological registration of vertices[J].Computer Engineering and Applications(CEA),2016,52(13):19-24.

[13] Berar M,Desvignes M,Bailly G,et al.3D meshes registration:application to statistical skull model[J].Lecture Notes in Computer Science(LNCS),2004,3212(1):100-107.

[14] Sahilliog?Lu Y,Kavan L.Skuller:a volumetric shape registration algorithm for modeling skull deformities[J].Medical Image Analysis,2015,23(1):15-27.

[15] Song Wei-yan.RANSAC algorithm and its application in remote sensing image processing[D].Beijing:North China Electric Power University,2011.

[16] Papazov C,Burschka D.An efficient RANSAC for 3D object recognition in noisy and occluded scenes[C].Asian Conference on Computer Vision(ACCV).Springer-Verlag,2010,6492:135-148.

[17] Han Bao-chang,Cao Jun-jie,Su Zhi-xun.Automatic point clouds registration based on regions [J].Journal of Computer-Aided Design & Computer Graphics(CAD&CG),2015,(2):313-319.

[18] Brandes U,Cornelsen S,Pampel B,et al.Blocks of hypergraphs:applied to hypergraphs and outerplanarity[C].International Conference on Combinatorial Algorithms,Springer-Verlag,2010:201-211.

[19] Bouaziz S,Tagliasacchi A,Pauly M.Sparse iterative closest point[J].Computer Graphics Forum,2013,32(5):113-123.

[20] Li W,Song P.A modified ICP algorithm based on dynamic adjustment factor for registration of point cloud and CAD model[J].Pattern Recognition Letters,2015,65(1):88-94.

[21] Senin N,Colosimo B M,Pacella M.Point set augmentation through fitting for enhanced ICP registration of point clouds in multisensor coordinate metrology[J].Robotics and Computer-Integrated Manufacturing,2013,29(1):39-52.

[22] Wang Meng-yang.Research and implementation of partition-interpolation craniofacial reconstruction method based on knowledge base[D].Xi′an:Northwest University,2010.

附中文參考文獻(xiàn)：

[1] 稅午陽,周明全,武仲科,等.數(shù)據(jù)配準(zhǔn)的顱骨面貌復(fù)原方法[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2011,23(4):607-614.

[4] 王 欣,張明明,于 曉,等.應(yīng)用改進(jìn)迭代最近點(diǎn)方法的點(diǎn)云數(shù)據(jù)配準(zhǔn)[J].光學(xué)精密工程,2012,20(9):2068-2077.

[5] 劉純國,劉 暢,安百玲.基于遺傳算法的三維曲面配準(zhǔn)[J].鍛壓裝備與制造技術(shù),2009,44(4):110-113.

[7] 朱延娟,周來水,張麗艷.散亂點(diǎn)云數(shù)據(jù)配準(zhǔn)算法[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2006,18(4):475-481.

[8] 葛寶臻,彭 博,田慶國.基于曲率圖的三維點(diǎn)云數(shù)據(jù)配準(zhǔn)[J].天津大學(xué)學(xué)報(bào),2013,46(2):174-180.

[9] 趙夫群,周明全,耿國華.基于有界旋轉(zhuǎn)角的點(diǎn)云配準(zhǔn)算法[J].微電子學(xué)與計(jì)算機(jī),2017,34(3):46-49.

[12] 賀毅岳,馬自萍,高 妮,等.生理點(diǎn)配準(zhǔn)的顱面局部形態(tài)關(guān)系建模[J].計(jì)算機(jī)工程與應(yīng)用,2016,52(13):19-24.

[15] 宋衛(wèi)艷.RANSAC算法及其在遙感圖像處理中的應(yīng)用[D].北京：華北電力大學(xué),2011.

[17] 韓寶昌,曹俊杰,蘇志勛.一種區(qū)域?qū)哟紊系淖詣狱c(diǎn)云配準(zhǔn)算法[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),2015，(2):313-319.

[22] 王孟陽.基于知識庫的分區(qū)插值顱面復(fù)原方法研究與實(shí)現(xiàn)[D].西安：西北大學(xué),2010.