趙夫群， 周明全， 耿國華

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

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

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

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

1 斷裂面的提取

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

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

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

(1)

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

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

(2)

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

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

(3)

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

2 斷裂面的匹配

2.1 ICP算法

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

(4)

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

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

(5)

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

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

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

(6)

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

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

(R*,t*)=

(7)

然后，更新Rk+1和tk+1，計算式為

(8)

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

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

2.2 SICP算法

既有尺度變換又有剛體變換的點集配準可描述為[14]

(9)

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

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

(10)

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

(11)

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

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

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

(12)

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

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

(13)

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

2.3 改進的ICP算法

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

首先，針對兩個3D點集的配準問題，定義旋轉(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)方法來計算. 對于3D點集的剛體配準，如果這兩個點集能夠正確配準，那它們必然具有相似的點分布. 因此當它們正確配準時，旋轉(zhuǎn)角的3對主成分矢量都接近0°或180°. 所以旋轉(zhuǎn)角可以通過旋轉(zhuǎn)角的3對主成分矢量來估計.

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

(14)

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

那么點集D和M的配準問題可進一步描述為

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

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

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

那么改進的ICP算法的具體實現(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；

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

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

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

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

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

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

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

.

(16)

3 碎塊匹配實驗

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

從圖2、圖3和圖4的匹配結(jié)果來看，ICP算法不能實現(xiàn)碎塊的尺度匹配，而SICP算法和改進ICP算法均能實現(xiàn)碎塊的尺度匹配，并且改進ICP算法的匹配誤差更小，效果更好. 實驗中，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算法和改進ICP算法的匹配參數(shù)(匹配誤差、迭代次數(shù)和耗時等)如表1所示.

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

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

4 結(jié)束語

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