趙夫群， 周明全， 耿國華

(1.西北大學(xué) 信息科學(xué)與技術(shù)學(xué)院,陜西,西安 710127；2.咸陽師范學(xué)院 教育科學(xué)學(xué)院,陜西,咸陽 712000；3.北京師范大學(xué) 信息科學(xué)與技術(shù)學(xué)院，北京 100875)

厚度不可忽略的剛體碎塊的匹配可以通過對碎塊斷裂面的匹配來實(shí)現(xiàn)，即為兩個不同坐標(biāo)系下的斷裂面尋找一個合適的變換，使得兩個曲面能夠完全配準(zhǔn)或重合. 目前，曲面匹配方法已經(jīng)在三維建模[1]、文物復(fù)原[2-3]、航空影像處理[4]以及醫(yī)學(xué)研究[5]等領(lǐng)域得到了廣泛的應(yīng)用.

目前有很多剛體碎塊斷裂面的匹配方法，可以利用點(diǎn)特征進(jìn)行匹配，如Pokrass 等[6]通過對碎塊斷裂面進(jìn)行離散采樣得到的若干采樣點(diǎn)進(jìn)行匹配；也可以利用斷裂面的輪廓曲線特征進(jìn)行匹配，如李群輝等[7]提出了一種基于輪廓曲線特征的剛體碎塊斷裂面匹配算法，王洋等[8]提出用傅里葉變換邊緣曲線法來實(shí)現(xiàn)W型輪廓曲線的快速匹配方法；此外還可以利用區(qū)域特征進(jìn)行匹配[9]，如Papaioannou等[10]通過獲取斷裂面上的局部特征區(qū)域?qū)崿F(xiàn)了碎塊的匹配. 以上這些算法都在一定程度上實(shí)現(xiàn)了剛體碎塊的匹配問題，但是它們大多只考慮匹配中的剛體變換問題，在實(shí)際的匹配應(yīng)用中往往還存在著尺度因素.

本文充分考慮尺度因素，針對剛體碎塊的三角網(wǎng)格數(shù)據(jù)模型，采用一種改進(jìn)的ICP算法進(jìn)行碎塊的匹配，即在ICP算法中添加尺度矩陣、旋轉(zhuǎn)角約束和動態(tài)迭代系數(shù)，來解決碎塊的尺度匹配問題、因剛體變換中旋轉(zhuǎn)角變換過大引起的匹配失敗的問題以及匹配速度的問題.

1 斷裂面的提取

為了提高計(jì)算速度、減少存儲空間、突出建模特征，首先采用基于平均曲率的精簡算法[11]對初始數(shù)據(jù)模型進(jìn)行精簡. 然后進(jìn)行碎塊外面表面的分割，即將碎塊的外表面以棱邊為界限分割成多張曲面. 這里采用區(qū)域生長算法[12]實(shí)現(xiàn)曲面分割，具體思想為：首先判斷相鄰三角面片法矢的夾角是否大于給定閾值，若大于給定閾值則將這兩個三角面片分割在兩個不同的曲面上，否則分割在同一曲面上，重復(fù)該過程直到所有三角面片分割完成為止；然后逐一檢查分割曲面，若曲面特別小，則將其合并到相鄰的曲面中，否則判斷該曲面與其相鄰曲面的平均法矢夾角，若小于給定閾值則進(jìn)行合并，否則不合并.

針對剛體碎塊的斷裂面較原始面更加粗糙的情況，可以根據(jù)斷裂面上所有頂點(diǎn)的法向量的變化情況來區(qū)分?jǐn)嗔衙婧驮济?，表現(xiàn)在幾何特征上就是斷裂面上頂點(diǎn)的法向量變化較大.

定義曲面φ上頂點(diǎn)p的彎能量ek(p)為

(1)

式中：qj表示頂點(diǎn)p在曲面φ上k近鄰的頂點(diǎn)；np和nqj分別表示p和qj的法向量.

設(shè)置一個以p為球心r為半徑的包圍球，則點(diǎn)p在鄰域Br(p)內(nèi)的局部彎曲能力定義為

(2)

式中：Nr(p)=Br(p)∩φ.

那么，曲面φ的平均粗糙度為

(3)

若ρ(φ)大于給定閾值ρ′，則認(rèn)為該曲面為斷裂面，否則認(rèn)為該曲面為原始面.

2 斷裂面的匹配

2.1 ICP算法

設(shè)兩個待匹配的斷裂面為φ1和φ2，其三角網(wǎng)格數(shù)據(jù)模型所對應(yīng)的三維點(diǎn)集分別為D={di}和M={mj}，i=1,2,…,Nd，j=1,2,…,Nm，Nd和Nm分別表示點(diǎn)集D和M的大小. 那么斷裂面φ1和φ2的匹配問題實(shí)際上就是點(diǎn)集D和M的配準(zhǔn)問題，即尋找變換T使得D和M能夠最佳匹配，可以采用最小二乘(least square,LS)法計(jì)算，計(jì)算式為

(4)

由Besl等[13]提出的ICP算法解決的就是兩個點(diǎn)集的配準(zhǔn)問題，即尋找從點(diǎn)集D到M的旋轉(zhuǎn)和平移變換，使得它們能夠達(dá)到最佳配準(zhǔn)，即令式 (4) 中的T為剛體變換(包括旋轉(zhuǎn)和平移變換)，計(jì)算式為

s.t.RR=Im,det(R)=1.

(5)

式中:R表示旋轉(zhuǎn)矩陣;t表示平移變換.

ICP算法的具體實(shí)現(xiàn)步驟如下.

① 建立點(diǎn)集D和M的相關(guān)性，計(jì)算式為

(6)

式中i=1,2,…,Nd.

② 計(jì)算點(diǎn)集D和M的新旋轉(zhuǎn)和平移變換，計(jì)算式為

(R*,t*)=

(7)

然后，更新Rk+1和tk+1，計(jì)算式為

(8)

重復(fù)步驟①和②，直到滿足終止條件為止.

雖然ICP算法是一種精度較高的點(diǎn)集配準(zhǔn)算法，但是它沒有考慮尺度因素. 在實(shí)際應(yīng)用中，經(jīng)常要考慮尺度配準(zhǔn)的問題，也就是既有尺度變換又有剛體變換的情況.

2.2 SICP算法

既有尺度變換又有剛體變換的點(diǎn)集配準(zhǔn)可描述為[14]

(9)

式中:S=diag[s1s2…sm]為一個尺度矩陣；R為一個正交矩陣；t為一個平移矢量.

將式 (9) 代入到式 (4) 中，則有

(10)

那么式 (10) 就描述了尺度配準(zhǔn)問題，式中正交矩陣R表示旋轉(zhuǎn)變換. 為了避免一個點(diǎn)集收斂到另外一個點(diǎn)集中的一個很小的子集情況，也就是尺度矩陣為0矩陣的情況，就必須為尺度矩陣S設(shè)置邊界. 于是，該尺度配準(zhǔn)問題又進(jìn)一步轉(zhuǎn)換為

(11)

式中:S為一個帶有邊界的尺度矩陣；R為一個旋轉(zhuǎn)矩陣；t為一個平移矢量.

那么SICP算法的實(shí)現(xiàn)步驟如下.

① 通過當(dāng)前的變換 (Sk,Rk,tk)建立點(diǎn)集D和M的相關(guān)性，計(jì)算式為

(12)

② 令S=diag[s1s2…sm]，計(jì)算新的變換(Sk+1,Rk+1,tk+1)，計(jì)算式為

(Sk+1,Rk+1,tk+1)=

(13)

重復(fù)步驟①到②，如果S的變化量 ΔS=|Sk+1-Sk|小于閾值ε或達(dá)到最大迭代次數(shù)εmax，則停止該迭代過程.

2.3 改進(jìn)的ICP算法

為了進(jìn)一步提高SICP算法的配準(zhǔn)精度和收斂速度，首先在SICP算法中引入旋轉(zhuǎn)角約束(即設(shè)置旋轉(zhuǎn)角的上界和下界)，解決由旋轉(zhuǎn)角變化過大引起的配準(zhǔn)失敗的問題，然后加入動態(tài)迭代系數(shù)h，以提高算法的收斂速度. 這里將這種改進(jìn)的SICP算法簡稱為改進(jìn)的ICP算法.

首先，針對兩個3D點(diǎn)集的配準(zhǔn)問題，定義旋轉(zhuǎn)矩陣R為R=RxRyRz.Rx,Ry,Rz分別為

式中:θx、θy和θz為旋轉(zhuǎn)角[15];θxb、θyb和θzb表示旋轉(zhuǎn)角的均值;Δθx、Δθy和Δθz為旋轉(zhuǎn)角的偏差;θxb-Δθx、θyb-Δθy和θzb-Δθz為旋轉(zhuǎn)角的下界;θxb+Δθx、θyb+Δθy和θzb+Δθz為旋轉(zhuǎn)角的上界.

旋轉(zhuǎn)角的主成分矢量采用主成分分析(princilalcomponentanalysis,PCA)方法來計(jì)算. 對于3D點(diǎn)集的剛體配準(zhǔn)，如果這兩個點(diǎn)集能夠正確配準(zhǔn)，那它們必然具有相似的點(diǎn)分布. 因此當(dāng)它們正確配準(zhǔn)時，旋轉(zhuǎn)角的3對主成分矢量都接近0°或180°. 所以旋轉(zhuǎn)角可以通過旋轉(zhuǎn)角的3對主成分矢量來估計(jì).

假設(shè)點(diǎn)集D和M的主成分分別為：VD={VD1,VD2,VD3}和VM={VM1,VM2,VM3}，VDi和VMi是3×1的矢量，i=1,2,3. 那么點(diǎn)集M中有8組主成分滿足條件，即

(14)

旋轉(zhuǎn)角的均值θxb、θyb和θzb可以通過式(14)的計(jì)算得到. 旋轉(zhuǎn)角的偏差Δθx、Δθy和Δθz的值取決于點(diǎn)集的相似程度，具體由實(shí)驗(yàn)情況來給定，為了計(jì)算方便可以設(shè)置為Δθx=Δθy=Δθz=10°.

那么點(diǎn)集D和M的配準(zhǔn)問題可進(jìn)一步描述為

動態(tài)迭代系數(shù)h是一個可以自動調(diào)整變換參數(shù)的整數(shù)值，它可以在不影響算法的配準(zhǔn)精度和收斂方向的情況下，大大減少算法的迭代次數(shù)、提高算法的收斂速度. 在SICP算法中加入動態(tài)迭代系數(shù)h的步驟如下：

① 計(jì)算剛體變換矢量q=[Rk|tk]T以及qk的相鄰兩次迭代的變化量Δqk；

② 用剛體變換矢量Δqk更新SICP算法步驟2)中的(Sk+1,Rk+1,tk+1)共h次，通常h是一個大于等于0的整數(shù).

那么改進(jìn)的ICP算法的具體實(shí)現(xiàn)步驟如下.

① 給定剛體變換初值q=[R0t0]T，R0為初始旋轉(zhuǎn)矩陣，t0為初始平移矢量，初始尺度值為S0，尺度的邊界為[ak,bk]；

② 令D0=R0S0D+t0，迭代次數(shù)k=0，動態(tài)迭代系數(shù)h=0；

③k=k+1；

④ 估計(jì)旋轉(zhuǎn)角度θx,θy,θz的邊界，即θx∈[θxb-Δθx,θxb+Δθx],θy∈[θyb-Δθy,θyb+Δθy]，θz∈[θzb-Δθz,θzb+Δθz]，Δθx=Δθy=Δθz=10°；

⑤ 計(jì)算最優(yōu)旋轉(zhuǎn)角θxo,θyo,θzo，最優(yōu)旋轉(zhuǎn)矩陣Ropt=RxoRyoRzo;

⑥ 利用式 (12) 建立點(diǎn)集D和M的相關(guān)性ck(i)∈{1,2,…,ND}；

⑦ 利用奇異值分解法計(jì)算旋轉(zhuǎn)矩陣Rk、平移矢量tk和尺度Sk；

⑧ 計(jì)算Dk=RkSkD+tk及其相鄰兩次迭代的變化量ΔDk；

⑨ 判斷動態(tài)迭代系數(shù)h，若h>0，則利用式 (15) 計(jì)算新的變換(Sk+1,Rk+1,tk+1)，并執(zhí)行Δqk(Sk+1,Rk+1,tk+1)共h次來進(jìn)一步更新變換；

⑩ 判斷均方根誤差σRMS，若σRMS,k+1-σRMS,k>ε，則執(zhí)行h=h+1操作，否則執(zhí)行h=0操作，RMS的定義為

.

(16)

3 碎塊匹配實(shí)驗(yàn)

實(shí)驗(yàn)采用公共數(shù)據(jù)集石頭和灰泥碎塊[16]以及3D掃描儀獲取的兵馬俑碎塊完成匹配，所有碎塊均采用三角網(wǎng)格數(shù)據(jù)模型表示，如圖1所示，顯然匹配中存在尺度變換. 匹配過程為：首先提取碎塊的斷裂面，然后分別采用ICP算法、SICP算法和改進(jìn)的ICP算法對碎塊斷裂面進(jìn)行尺度配準(zhǔn). 石頭碎塊、灰泥碎塊和兵馬俑碎塊的匹配結(jié)果分別如圖2、圖3和圖4所示.

從圖2、圖3和圖4的匹配結(jié)果來看，ICP算法不能實(shí)現(xiàn)碎塊的尺度匹配，而SICP算法和改進(jìn)ICP算法均能實(shí)現(xiàn)碎塊的尺度匹配，并且改進(jìn)ICP算法的匹配誤差更小，效果更好. 實(shí)驗(yàn)中，ICP算法、

圖1 待匹配的碎塊 Fig.1 Blocks needed to be matched

圖2 石頭碎塊的匹配結(jié)果Fig.2 Matching results of stone blocks

圖3 灰泥碎塊的匹配結(jié)果Fig.3 Matching results of plaster blocks

圖4 兵馬俑碎塊的匹配結(jié)果Fig.4 Matching results of Terracotta Warriors blocks

SICP算法和改進(jìn)ICP算法的匹配參數(shù)(匹配誤差、迭代次數(shù)和耗時等)如表1所示.

表1 算法的匹配參數(shù)

從表1可見，ICP算法的匹配誤差較大，SICP算法和改進(jìn)ICP算法均有較高的匹配精度，并且改進(jìn)ICP算法的匹配速度和精度較SICP算法又有了大幅度的提高，是一種精度更高、收斂速度更快的碎塊斷裂面點(diǎn)集匹配算法.

4 結(jié)束語

針對厚度不可忽略的剛體碎塊的三角網(wǎng)格數(shù)據(jù)模型，本文提出了一種改進(jìn)的ICP算法來對兩個剛體碎塊的斷裂面進(jìn)行精確匹配. 該算法通過加入尺度矩陣、旋轉(zhuǎn)角度約束和動態(tài)迭代系數(shù)的方式來改進(jìn)ICP算法，使得算法在匹配精度、速度等方面的性能都得到了很大的提高，能夠?qū)蓚€不同尺度碎塊的斷裂面進(jìn)行精確、快速地對齊. 在今后的研究中，要進(jìn)一步考慮噪聲和外點(diǎn)對碎塊匹配精度和速度的影響，提出更加精確、快速和抗噪的碎塊匹配算法. 此外，還要進(jìn)一步研究特殊碎塊的匹配問題，如兵馬俑碎塊，分析其特點(diǎn)，尋找更加適合這種特殊復(fù)雜碎塊的匹配算法，以達(dá)到自動虛擬復(fù)原兵馬俑的目的.