王帥帥, 柏艷紅, 王 銀, 孫志毅

(太原科技大學(xué)電子信息工程學(xué)院, 太原 030024)

三維激光掃描技術(shù)廣泛應(yīng)用于醫(yī)學(xué)研究[1]、文物保護(hù)[2]和三維重建[3]等領(lǐng)域．在三維數(shù)據(jù)采集過程中, 由于受激光掃描儀測量角度或外界障礙物遮擋等因素的影響, 被掃描物體的部分點(diǎn)云信息易于遺漏．為了獲取目標(biāo)物體表面完整的點(diǎn)云數(shù)據(jù), 須采用點(diǎn)云配準(zhǔn)[4-5]技術(shù)將多視角下的點(diǎn)云數(shù)據(jù)轉(zhuǎn)換到同一坐標(biāo)系下, 為后續(xù)的三維重建奠定基礎(chǔ)．由于點(diǎn)云配準(zhǔn)是三維重建的關(guān)鍵技術(shù), 故快速且精確的點(diǎn)云配準(zhǔn)技術(shù)具有重要的應(yīng)用價(jià)值．目前,應(yīng)用最廣泛且最經(jīng)典的點(diǎn)云配準(zhǔn)方法是Besl等[6]提出的迭代最近點(diǎn)(iterative closest point, ICP)算法, 該方法以四元數(shù)運(yùn)算模型為基礎(chǔ), 不斷迭代獲取使得兩片點(diǎn)云間歐式距離最小的變換矩陣, 但對兩片點(diǎn)云的初始位姿依賴性較強(qiáng),在迭代時容易陷入局部極值, 導(dǎo)致配準(zhǔn)結(jié)果較差．針對ICP算法的缺陷, 研究者大多通過有效的粗配準(zhǔn)獲得良好的初始配準(zhǔn)位姿, 再采用ICP算法完成精配準(zhǔn)[7]．Lei等[8]針對重疊小、噪聲大且位姿差異大的點(diǎn)云, 提出了一種基于均勻采樣點(diǎn)的全局點(diǎn)云配準(zhǔn)方法, 利用由距離誤差表示的質(zhì)量函數(shù)識別最優(yōu)變換并實(shí)現(xiàn)點(diǎn)云的粗配準(zhǔn), 進(jìn)而采用改進(jìn)的ICP算法完成精配準(zhǔn), 顯著提升了配準(zhǔn)效率和精度; 荊路等[9]基于尺度不變特征變換(scale invariant feature transform, SIFT)關(guān)鍵點(diǎn)結(jié)合快速點(diǎn)特征直方圖(fast point feature histogram, FPFH)描述子, 利用采樣一致性初始配準(zhǔn) (sample consensus initial alignment, SAC-IA)算法完成初始配準(zhǔn), 該算法雖提高了配準(zhǔn)精度, 但配準(zhǔn)時間較長; Zaganidis等[10]采用點(diǎn)云的語義信息作為搜索對應(yīng)關(guān)系的先驗(yàn)信息, 輔助正態(tài)分布變換(normal distribution transform, NDT)配準(zhǔn)點(diǎn)云, 改善了配準(zhǔn)精度和效率,但適用場景有限; 李宇翔等[11]基于內(nèi)部形態(tài)描述子(intrinsic shape signatures, ISS)關(guān)鍵點(diǎn)[12]利用余弦相似度匹配對應(yīng)點(diǎn)對獲得良好的初始位姿, 然后采用點(diǎn)到平面的ICP算法進(jìn)行精配準(zhǔn), 該算法對低重疊率的點(diǎn)云具有較好的魯棒性, 但方向直方圖(signature of histograms of orientations, SHOT)描述子匹配效率較低; Li等[13]提出了一種基于對應(yīng)關(guān)系的靈活配準(zhǔn)方法, 該方法配準(zhǔn)精度高、速度快且魯棒性較強(qiáng), 但其配準(zhǔn)精度在一定程度上依賴于關(guān)鍵點(diǎn)檢測器的定位精度; Li等[14]采用高斯混合模型擬合點(diǎn)云模型, 配準(zhǔn)精度雖得以提高, 但配準(zhǔn)效率較低．本文擬提出一種結(jié)合內(nèi)部形態(tài)描述子關(guān)鍵點(diǎn)及二進(jìn)制方向直方圖描述子(binary signature of histograms of orientations, BSHOT)的點(diǎn)云配準(zhǔn)方法, 通過粗配準(zhǔn)求解初始變換矩陣, 利用該矩陣將源點(diǎn)云轉(zhuǎn)換至目標(biāo)點(diǎn)云坐標(biāo)系, 再采用ICP算法完成精配準(zhǔn)．

1 配準(zhǔn)方案

本文配準(zhǔn)過程由粗配準(zhǔn)和精配準(zhǔn)兩部分組成: 1) 粗配準(zhǔn)．由于BSHOT特征描述子結(jié)合了點(diǎn)云的密度和噪聲等幾何分布信息, 所以首先計(jì)算點(diǎn)云的分辨率, 利用ISS算法提取兩片點(diǎn)云的關(guān)鍵點(diǎn)并在關(guān)鍵點(diǎn)的鄰域構(gòu)建局部坐標(biāo)系, 統(tǒng)計(jì)鄰域點(diǎn)和查詢點(diǎn)法向量的夾角余弦值并計(jì)入直方圖中; 然后將高維直方圖特征轉(zhuǎn)化為BSHOT描述子, 利用漢明距離匹配對應(yīng)點(diǎn)對,再采用隨機(jī)采樣一致性(random sampling consistency, RSC)算法去除誤匹配點(diǎn)對, 得到點(diǎn)云配準(zhǔn)的初始變換矩陣; 2) 精配準(zhǔn)．通過點(diǎn)對點(diǎn)的ICP算法估算出最優(yōu)變換矩陣．具體配準(zhǔn)流程如圖1所示．

圖1 本文配準(zhǔn)算法流程Fig.1 Registration algorithm process of this paper

2 特征點(diǎn)提取算法及對應(yīng)點(diǎn)匹配

2.1 ISS關(guān)鍵點(diǎn)提取

當(dāng)傳感器視角發(fā)生稍許變化時,點(diǎn)云中的關(guān)鍵點(diǎn)依然保持不變,其穩(wěn)定性強(qiáng)、辨識度高且信息量豐富,關(guān)鍵點(diǎn)的數(shù)量遠(yuǎn)小于原始點(diǎn)云數(shù)量,故本文采用ISS算法提取關(guān)鍵點(diǎn)．一個內(nèi)部形態(tài)描述子ISS由一個局部參考坐標(biāo)系和一個編碼三維形態(tài)特征的高分辨率特征向量組成．

ISS算法通過計(jì)算各個查詢點(diǎn)與鄰域范圍內(nèi)所有點(diǎn)的協(xié)方差矩陣的特征值, 并設(shè)定特征值間的比例閾值選取關(guān)鍵點(diǎn)．設(shè)點(diǎn)云P中共有m個點(diǎn), 任意一點(diǎn)pi的原始坐標(biāo)為(xi,yi,zi),i=1,2,…m, ISS關(guān)鍵點(diǎn)提取過程如下:

1) 對點(diǎn)云數(shù)據(jù)P中的任一查詢點(diǎn)pi構(gòu)建局部參考坐標(biāo)系, 并將rf視為搜索半徑;

2) 統(tǒng)計(jì)點(diǎn)云P中每個查詢點(diǎn)pi在半徑為rf的球形鄰域內(nèi)的所有點(diǎn)pj, 并計(jì)算其權(quán)值

(1)

3) 利用rf范圍內(nèi)的所有鄰域點(diǎn)pj計(jì)算查詢點(diǎn)pi的加權(quán)協(xié)方差矩陣

(2)

5) 設(shè)定兩個小于1的閾值ε1,ε2,若滿足條件:

(3)

即視查詢點(diǎn)pi為ISS關(guān)鍵點(diǎn)．

2.2 BSHOT特征描述子

為使特征描述子保持旋轉(zhuǎn)、平移和尺度不變性,同時減少描述子匹配的時間, 受文獻(xiàn)[15]啟發(fā), 采用BSHOT算法對所提取關(guān)鍵點(diǎn)的鄰域特征進(jìn)行描述．由于BSHOT描述子是對SHOT描述子的二進(jìn)制編碼,故須計(jì)算SHOT描述子, 步驟如下:

1) 求解鄰域協(xié)方差矩陣

(4)

式中r為鄰域半徑,di為鄰域點(diǎn)qi到關(guān)鍵點(diǎn)q的歐氏距離, 并根據(jù)所求特征向量構(gòu)建局部參考坐標(biāo)系;

2) 求解協(xié)方差矩陣按降序排列的特征值(λ1>λ2>λ3)及相互正交的單位特征向量(v1,v2,v3),分別對應(yīng)于局部參考坐標(biāo)系的x+,y+,z+方向．以局部參考坐標(biāo)系為基準(zhǔn)建立關(guān)鍵點(diǎn)的球形鄰域,沿球的徑向、經(jīng)度、緯度方向分別劃分為2, 8, 2個部分,則整個球形鄰域被劃分為如圖2所示的32個子空間;

圖2 BSHOT空間結(jié)構(gòu)Fig.2 Space structure of BSHOT

3) 統(tǒng)計(jì)子空間內(nèi)各鄰域點(diǎn)qi與關(guān)鍵點(diǎn)q法向量間的夾角余弦值cosθ并計(jì)入11維直方圖, 再組合32個子空間的直方圖, 構(gòu)成352維的SHOT高維直方圖特征．

一個352維的SHOT描述子需要1 408 B的存儲空間, 將其編碼為352位(即44 B)的二進(jìn)制描述子,不僅可降低描述子所需存儲空間,而且能提高特征匹配的效率．對于一個SHOT描述子Si,i={0,1,2,…,351},Si中的每個值都是0和1之間的十進(jìn)制數(shù)．假設(shè)以Bi表示BSHOT描述子,i={0,1,2,…,351}, 且Bi的每個值都為二進(jìn)制數(shù)0或1．從Si的初始端連續(xù)取4個值{S0,S1,S2,S3},對應(yīng)位被編碼為BSHOT描述子{B0,B1,B2,B3}．

設(shè)Ssum=S0+S1+S2+S3, 若Si的4個值都為0, 則{B0,B1,B2,B3}為{0,0,0,0}; 若有一個值Si∈{S0,S1,S2,S3}超過Ssum的90%, 則對應(yīng)的Bi編碼為1,如S2值大于Ssum的90%, 則{B0,B1,B2,B3}編碼為{0,0,1,0}．若上述情況不成立,再檢查任意兩個值的和是否大于Ssum的90%, 如若S0+S1超過Ssum的90%, 則將{B0,B1,B2,B3}編碼為{1,1,0,0}．以此類推, BSHOT描述子通過上述迭代方式將SHOT描述子由十進(jìn)制編碼為二進(jìn)制．

2.3 對應(yīng)關(guān)鍵點(diǎn)匹配

經(jīng)過計(jì)算點(diǎn)云分辨率、提取ISS關(guān)鍵點(diǎn)及計(jì)算BSHOT描述子后,可得到源點(diǎn)云和目標(biāo)點(diǎn)云中各片點(diǎn)云的關(guān)鍵點(diǎn)及其特征向量．由于點(diǎn)云密度不一、拓?fù)浣Y(jié)構(gòu)復(fù)雜且存在噪聲點(diǎn),所以不同點(diǎn)云中的同一個關(guān)鍵點(diǎn)的描述子會出現(xiàn)偏差,當(dāng)利用雙向漢明距離進(jìn)行對應(yīng)關(guān)鍵點(diǎn)匹配時易產(chǎn)生錯誤的匹配點(diǎn)對,這將會影響后續(xù)初始變換矩陣的求解．

為了解決可能出現(xiàn)的誤匹配問題,現(xiàn)采用RSC算法剔除誤匹配點(diǎn)對,并通過隨機(jī)取樣剔除局外點(diǎn),建立一個僅由局內(nèi)點(diǎn)數(shù)據(jù)組成的點(diǎn)集．下面對輸入點(diǎn)云數(shù)據(jù)隨機(jī)選取點(diǎn)集計(jì)算預(yù)設(shè)的模型參數(shù)．假設(shè)模型的初始狀態(tài)至少含3個局內(nèi)點(diǎn),局內(nèi)點(diǎn)在數(shù)據(jù)中的占比

ε=ni/N,

(5)

式中ni為局內(nèi)點(diǎn)數(shù),N為局內(nèi)點(diǎn)與局外點(diǎn)之和．若以(1-ε)k表示經(jīng)k次迭代計(jì)算后模型至少采樣到一個局外點(diǎn)的概率, 則局內(nèi)點(diǎn)的采樣概率P=1-(1-ε)k,P越接近1表示模型中局內(nèi)點(diǎn)數(shù)量越多．在保證一定置信度下迭代結(jié)束后估計(jì)出最優(yōu)的模型參數(shù)．

剔除錯誤的對應(yīng)點(diǎn)對后,利用剩余的正確對應(yīng)點(diǎn)關(guān)系, 以目標(biāo)點(diǎn)云作為參考, 求解兩片點(diǎn)云間的初始變換矩陣M, 將源點(diǎn)云轉(zhuǎn)換到以目標(biāo)點(diǎn)云為標(biāo)準(zhǔn)的坐標(biāo)系中．兩片點(diǎn)云間的變換關(guān)系為

(6)

式中(x′i,y′i,z′i)為目標(biāo)點(diǎn)云上的點(diǎn),(xi,yi,zi)為源點(diǎn)云上與目標(biāo)點(diǎn)云對應(yīng)的點(diǎn)．由于變換矩陣M中含有12個與旋轉(zhuǎn)或平移相關(guān)的參數(shù),所以至少需要6組對應(yīng)點(diǎn)通過求解線性方程組得出變換矩陣M．

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

硬件環(huán)境為Inter (R) Core i5-4200H@2.8 GHz處理器,8.0 GB運(yùn)行內(nèi)存, 操作系統(tǒng)為Windows10 64位．軟件環(huán)境為cmake-gui, Visual Studio 2017, 開源點(diǎn)云庫PCL1.8.1．實(shí)驗(yàn)數(shù)據(jù)采用斯坦福大學(xué)不同視角下的dragon和bunny點(diǎn)云模型, dragon000與dragon024的點(diǎn)云數(shù)量分別為41 841, 34 836, bunny000與bunny045的點(diǎn)云數(shù)量分別為40 256,40 097．首先加載源點(diǎn)云和目標(biāo)點(diǎn)云,由于BSHOT特征結(jié)合了點(diǎn)云的幾何分布信息,故須在計(jì)算點(diǎn)云分辨率后提取兩片點(diǎn)云中的ISS關(guān)鍵點(diǎn), 然后根據(jù)關(guān)鍵點(diǎn)的鄰域特征計(jì)算BSHOT描述子．相同關(guān)鍵點(diǎn)的特征描述子也近似相同,本文算法剔除誤匹配后的對應(yīng)點(diǎn)如圖3所示．

圖3 對應(yīng)點(diǎn)匹配結(jié)果Fig.3 Corresponding point matching results

在相同軟硬件條件下, 將本文ISS-BSHOT+點(diǎn)到點(diǎn)的ICP算法與傳統(tǒng)點(diǎn)到點(diǎn)的ICP算法[6]、SIFT-FPFH+點(diǎn)到點(diǎn)ICP算法[9]、ISS-SHOT+點(diǎn)到面的ICP算法[11]進(jìn)行對比實(shí)驗(yàn)．為了驗(yàn)證本文算法的穩(wěn)健性,將目標(biāo)點(diǎn)云沿x軸正方向平移0.25 m后進(jìn)行配準(zhǔn)實(shí)驗(yàn)．配準(zhǔn)前源點(diǎn)云和目標(biāo)點(diǎn)云的初始位姿如圖3所示．由圖3可見,配準(zhǔn)前源點(diǎn)云與目標(biāo)點(diǎn)云初始位姿相差較大且無重疊．本文采用均方根誤差(root mean square error, RMSE)作為配準(zhǔn)精度的評價(jià)指標(biāo),衡量兩片點(diǎn)云在空間上的接近程度,結(jié)果如表1所示．圖4為不同算法的粗配準(zhǔn)結(jié)果,其中紅色點(diǎn)云為目標(biāo)點(diǎn)云,藍(lán)色點(diǎn)云為經(jīng)過初始變換矩陣變換后的源點(diǎn)云．

由圖4和表1可知:本文粗配準(zhǔn)算法的精度略低于文獻(xiàn)[11],較文獻(xiàn)[9]高47.4%;雖然本文粗配準(zhǔn)算法在將SHOT特征描述子轉(zhuǎn)換為BSHOT特性描述子時會丟失部分信息,使得利用BSHOT特征描述子查找到的對應(yīng)點(diǎn)數(shù)量低于SHOT描述子,進(jìn)而導(dǎo)致粗配準(zhǔn)精度略低,但是仍能為后續(xù)精配準(zhǔn)提供良好的初始位姿．

為兼顧配準(zhǔn)效率和精度,本文采用點(diǎn)到點(diǎn)的ICP算法進(jìn)行精配準(zhǔn),設(shè)定ICP最大迭代次數(shù)為10,幾種不同算法的精配準(zhǔn)結(jié)果如圖5所示,其中紅色點(diǎn)云為目標(biāo)點(diǎn)云,藍(lán)色點(diǎn)云為經(jīng)過ICP所求變換矩陣變換后的源點(diǎn)云．由圖5和表1可知:傳統(tǒng)ICP配準(zhǔn)算法的配準(zhǔn)精度最低且用時較長;本文配準(zhǔn)算法用時最少,配準(zhǔn)效率較文獻(xiàn)[9]和文獻(xiàn)[11]分別高32.5%,51.6%,配準(zhǔn)精度較文獻(xiàn)[9]高15.4%,但略低于文獻(xiàn)[11]．其原因是采用二進(jìn)制的BSHOT特征描述子進(jìn)行關(guān)鍵點(diǎn)的特征匹配時速度較快,但粗配準(zhǔn)精度略低于文獻(xiàn)[11],進(jìn)而導(dǎo)致精配準(zhǔn)精度下降．