基于特征向量的點(diǎn)云配準(zhǔn)方法研究

俞 浩, 高 飛

(合肥工業(yè)大學(xué) 土木與水利工程學(xué)院,安徽 合肥 230009)

0 引 言

三維激光掃描技術(shù)是20世紀(jì)90年代中葉出現(xiàn)的新興技術(shù),它可以不接觸被測(cè)物體得到被測(cè)物體表面大量的密度高的三維點(diǎn)信息以及反射率等信息,由于三維激光掃描的快速性、準(zhǔn)確性,因此該技術(shù)又被稱為“實(shí)景復(fù)制技術(shù)”。該技術(shù)被廣泛運(yùn)用在逆向工程、虛擬現(xiàn)實(shí)、古建筑保護(hù)與修復(fù)等領(lǐng)域。

通常在外業(yè)采集數(shù)據(jù)時(shí),往往對(duì)某一物體不能通過一次掃描得到完整的三維點(diǎn)云數(shù)據(jù),因此在實(shí)際作業(yè)時(shí),需要從不同角度對(duì)物體進(jìn)行數(shù)據(jù)采集,因?yàn)樵诓煌嵌认逻M(jìn)行測(cè)量時(shí)坐標(biāo)系不同,所以必須要通過一定的方法對(duì)將不同角度下的點(diǎn)云數(shù)據(jù)進(jìn)行點(diǎn)云配準(zhǔn),即將不同坐標(biāo)系下的點(diǎn)云數(shù)據(jù)轉(zhuǎn)換到同一個(gè)坐標(biāo)系之下,從而得到物體完整的三維點(diǎn)云數(shù)據(jù)。常見的粗配準(zhǔn)方法有七參數(shù)轉(zhuǎn)換法和四元數(shù)法。

七參數(shù)轉(zhuǎn)換法就是要求得兩點(diǎn)云所在坐標(biāo)系之間的3個(gè)平移量、3個(gè)旋轉(zhuǎn)量以及1個(gè)尺度變換參數(shù)。轉(zhuǎn)換公式如下:

X=ΔX+(1+k)RP

(1)

R=R(α)R(β)R(γ)

(2)

(3)

(4)

(5)

由(1)式、(5)式可得:

(6)

其中,X、Y、Z為目標(biāo)點(diǎn)云坐標(biāo);Δx、Δy、Δz為3個(gè)平移參數(shù);k為尺度變換參數(shù);α、β、γ為3個(gè)旋轉(zhuǎn)參數(shù);px、py、pz為源點(diǎn)云坐標(biāo)。

由(6)式可知，得到7個(gè)轉(zhuǎn)換參數(shù)就可以實(shí)現(xiàn)源點(diǎn)云坐標(biāo)到目標(biāo)點(diǎn)云的轉(zhuǎn)化。

四元數(shù)是Hamilton在1843年找到的一種擴(kuò)展的復(fù)數(shù),于1986年由Faugeras和Hebert用四元數(shù)進(jìn)行配準(zhǔn)研究[1]。一個(gè)四元數(shù)h擁有1個(gè)實(shí)部和3個(gè)虛部,即

h=q0+q1i+q2j+q3k

(7)

其中,i、j、k為四元數(shù)的3個(gè)虛部單位。四元數(shù)是一個(gè)包含4個(gè)矢量的陣列:

(8)

(9)

(10)

定義源點(diǎn)云與目標(biāo)點(diǎn)云的質(zhì)心為:

(11)

兩點(diǎn)云之間的協(xié)方差矩陣為:

(12)

由協(xié)方差矩陣構(gòu)成的4×4對(duì)稱矩陣為:

(13)

計(jì)算Q(R)的特征向量與特征值,其最大特征值對(duì)應(yīng)的特征向量為qR=[q0q1q2q3]T。

qT=μx-R(qR)μp

(14)

本文通過引入特征向量的歐氏距離來解決兩片點(diǎn)云間的特征匹配問題,通過七參數(shù)轉(zhuǎn)換法對(duì)點(diǎn)云坐標(biāo)進(jìn)行轉(zhuǎn)換,確保在進(jìn)行最近點(diǎn)迭代法前兩點(diǎn)云有著比較好的初始位置。本文通過應(yīng)用實(shí)例驗(yàn)證該方法的可行性。

1 特征點(diǎn)匹配

要實(shí)現(xiàn)不同坐標(biāo)系下的點(diǎn)云配準(zhǔn),就要求解出不同坐標(biāo)系之間的轉(zhuǎn)換參數(shù)(旋轉(zhuǎn)矩陣R和平移矩陣T),那么就必須要對(duì)不同坐標(biāo)系下的相同點(diǎn)進(jìn)行匹配。

目前特征點(diǎn)匹配主要分為2類:① 通過控制點(diǎn)來進(jìn)行坐標(biāo)系之間的轉(zhuǎn)換[2];② 根據(jù)點(diǎn)云重疊部分的點(diǎn)以及點(diǎn)與點(diǎn)之間的固有關(guān)系來進(jìn)行配準(zhǔn)。這種配準(zhǔn)方法主要有求解點(diǎn)與點(diǎn)之間的配準(zhǔn)算法[3-4]、求解點(diǎn)與平面之間距離的配準(zhǔn)算法、計(jì)算豪斯多夫距離來衡量點(diǎn)云相似度以及基于幾何特征的配準(zhǔn)算法[5-6]。著名的最近點(diǎn)迭代法就是通過求解點(diǎn)與點(diǎn)之間的歐式距離的最小值來進(jìn)行點(diǎn)對(duì)匹配,文獻(xiàn)[7]通過計(jì)算點(diǎn)鄰域的切平面與點(diǎn)之間的距離來進(jìn)行匹配,但是這種方法的魯棒性不是太好。文獻(xiàn)[8]進(jìn)行了一些豪斯多夫距離在點(diǎn)云配準(zhǔn)的理論研究,但未進(jìn)行實(shí)際配準(zhǔn)應(yīng)用。

總的來說,第2種配準(zhǔn)步驟主要分為尋找對(duì)應(yīng)點(diǎn)和求解對(duì)應(yīng)關(guān)系2步。

尺度不變特征變換(scale-invariant feature transform,SIFT)算法可以將不同坐標(biāo)系下的點(diǎn)云里的特征點(diǎn)提取出來。SIFT算法是David Lowe在1999年提出,2004年整理完善的[9],最初是一種根據(jù)電腦視覺的算法來偵測(cè)與描述影像的局部特征點(diǎn),這些局部特征點(diǎn)可以幫助識(shí)別物體。該算法同樣可以運(yùn)用到點(diǎn)云數(shù)據(jù)處理中來,其實(shí)質(zhì)是計(jì)算出特征點(diǎn)的方向。它所找出來的關(guān)鍵點(diǎn)其實(shí)是一些在點(diǎn)云數(shù)據(jù)中比較突出的點(diǎn),這些點(diǎn)不會(huì)因視角、光照等因素而改變。

在將特征點(diǎn)提取出來后,就可以用KDTree對(duì)這些特征點(diǎn)進(jìn)行其法向量的求解。通過對(duì)其鄰域點(diǎn)集求法向量來近似代替其法向量[10]。對(duì)鄰域點(diǎn)的三維坐標(biāo)進(jìn)行中心化,求得協(xié)方差矩陣并對(duì)角化,求得3個(gè)特征值、最小特征值對(duì)應(yīng)的特征向量即為法向量。

(15)

C·nj=λj·nj

(16)

通過對(duì)不同坐標(biāo)系下的點(diǎn)云數(shù)據(jù)P的所有特征向量ni(i=1,2,3,…,n)和點(diǎn)云數(shù)據(jù)Q的特征向量nj(j=1,2,3,…,n)求解歐式距離,若最小歐式距離小于一個(gè)閾值ε,則這2個(gè)特征向量所對(duì)應(yīng)的特征點(diǎn)即被認(rèn)為是在不同坐標(biāo)系下的相同點(diǎn)；若最小歐式距離大于這個(gè)閾值ε,則不進(jìn)行匹配,進(jìn)行下一個(gè)特征向量的最小歐式距離的求解。

運(yùn)用求解歐式距離對(duì)2個(gè)不同坐標(biāo)系下的特征點(diǎn)進(jìn)行匹配的具體流程如圖1所示。

圖1 計(jì)算歐式距離匹配點(diǎn)對(duì)流程

2 點(diǎn)云粗配準(zhǔn)

用匹配好的兩組不同坐標(biāo)系下的點(diǎn)云數(shù)據(jù)進(jìn)行七參數(shù)求解[11]，即

A=(1+k)RB+T

(17)

加上最小二乘的約束條件,得出一個(gè)最符合的七參數(shù)。其中

(17)式中,A為目標(biāo)點(diǎn)云;B為源點(diǎn)云;R為旋轉(zhuǎn)矩陣;T為平移矩陣;k為縮放比例。由于點(diǎn)云坐標(biāo)式基于掃描儀的極坐標(biāo)系求出,因此縮放比例k可以忽略。

(17)式可以變?yōu)?

A=RB+T

(18)

結(jié)合最小二乘的約束條件：

VTPV=min

(19)

可以得到一個(gè)最為準(zhǔn)確的旋轉(zhuǎn)矩陣R和平移矩陣T。

3 點(diǎn)云精配準(zhǔn)

運(yùn)用粗配準(zhǔn)求得的旋轉(zhuǎn)矩陣R和平移矩陣T可以將源點(diǎn)云數(shù)據(jù)轉(zhuǎn)換成目標(biāo)點(diǎn)云附近,實(shí)現(xiàn)2片點(diǎn)云的粗配準(zhǔn)。但是不同點(diǎn)云之間會(huì)有分叉的情況出現(xiàn),基于這種情況,引入最近點(diǎn)迭代法對(duì)點(diǎn)云位置進(jìn)行精配準(zhǔn)。

最近點(diǎn)迭代法(ICP算法)最先由Besl和Mckay提出[12-13],其本質(zhì)就是在最小二乘的約束條件下,求出2個(gè)點(diǎn)集之間的變換關(guān)系,假設(shè){Pi,i=1,2,3,…,n}為源點(diǎn)集,{Qi,i=1,2,3,…,m}為目標(biāo)點(diǎn)集,則構(gòu)建函數(shù)為：

(20)

通過迭代確定這2個(gè)點(diǎn)集的最優(yōu)匹配,直到達(dá)到設(shè)定的某個(gè)收斂準(zhǔn)則為止,也就是通過不停地求解這2個(gè)點(diǎn)集的變換關(guān)系，并且根據(jù)這個(gè)變換更新2個(gè)點(diǎn)集的位置,直到變換關(guān)系滿足這個(gè)收斂準(zhǔn)則,這時(shí)2個(gè)點(diǎn)集之間的變換關(guān)系就認(rèn)為是最優(yōu)解，即旋轉(zhuǎn)矩陣R與平移矩陣T最優(yōu),其具體步驟為:

(1) 首先從源點(diǎn)云及目標(biāo)點(diǎn)云中分別提取點(diǎn)集Pi∈P,Qi∈Q。

(2) 計(jì)算源點(diǎn)集中的點(diǎn)與目標(biāo)點(diǎn)集中的點(diǎn)的對(duì)應(yīng)關(guān)系,使得‖Qi-Pi‖2=min。

(3) 計(jì)算旋轉(zhuǎn)矩陣R與平移矩陣T,滿足關(guān)系如下：

(21)

(4) 運(yùn)用旋轉(zhuǎn)矩陣R與平移矩陣T變換P點(diǎn)云的位置,即Pi′={Pi′=RPi+T}。

(5) 計(jì)算均方差δ,具體公式如下：

(22)

(6) 比較δ與最初設(shè)定的閾值ε。若δ<ε,則停止計(jì)算；若δ>ε,則重復(fù)步驟(2)～步驟(6)。

這樣經(jīng)過每一次迭代,Pi就會(huì)更靠近點(diǎn)云Qi,基于這種思想,Besl證明了這種方法的收斂性。

4 應(yīng)用實(shí)例分析

實(shí)驗(yàn)數(shù)據(jù)為2個(gè)不同坐標(biāo)系下的某古文物鼎的點(diǎn)云數(shù)據(jù)，初始點(diǎn)云如圖2所示,左側(cè)點(diǎn)云為源點(diǎn)云,右側(cè)點(diǎn)云為目標(biāo)點(diǎn)云,目的是要將左側(cè)點(diǎn)云向右側(cè)點(diǎn)云轉(zhuǎn)化。先用SIFT算法提取出特征點(diǎn)后再對(duì)特征點(diǎn)進(jìn)行法向量求解，并進(jìn)行點(diǎn)對(duì)配對(duì),得到的粗配準(zhǔn)的旋轉(zhuǎn)矩陣R和平移矩陣T分別為:

圖2 兩點(diǎn)云初始位置

粗配準(zhǔn)的結(jié)果如圖3所示。用最近點(diǎn)迭代法進(jìn)行迭代50次實(shí)現(xiàn)精配準(zhǔn),再迭代的過程中每一次迭代都會(huì)生成一個(gè)旋轉(zhuǎn)矩陣及平移矩陣,相對(duì)應(yīng)的源點(diǎn)云也就一次一次地更加逼近目標(biāo)點(diǎn)云。其中最后一次的旋轉(zhuǎn)矩陣和平移矩陣分別為:

圖3 粗配準(zhǔn)后的點(diǎn)云位置

將初始點(diǎn)云直接進(jìn)行ICP算法陷入的局部最優(yōu)的情況如圖4所示。

圖4 直接運(yùn)用ICP算法陷入的局部最優(yōu)情況

直接進(jìn)行ICP算法每一次迭代的具體誤差值以及運(yùn)行時(shí)間如圖5所示,經(jīng)過粗配準(zhǔn)之后再運(yùn)用ICP算法進(jìn)行精配準(zhǔn)后的效果如圖6所示。

圖6 經(jīng)過ICP算法迭代后的兩片點(diǎn)云

從圖4可以看出,直接進(jìn)行ICP算法進(jìn)行點(diǎn)云配準(zhǔn)會(huì)導(dǎo)致點(diǎn)云并不能很好地拼接在一起,甚至?xí)霈F(xiàn)局部最優(yōu)的情況,如圖5所示,雖然迭代后的誤差值逐漸趨于收斂,但是實(shí)際位置并沒有很好地配準(zhǔn)。因此本文的粗配準(zhǔn)算法有一定的必要性,通過粗配準(zhǔn)將兩點(diǎn)云調(diào)整至大致重疊的位置,它們之間僅僅相差較小的旋轉(zhuǎn)和平移。

圖5 直接ICP算法運(yùn)行時(shí)間及誤差值

從圖6可以看出,經(jīng)過粗配準(zhǔn)之后再進(jìn)行點(diǎn)云數(shù)據(jù)精配準(zhǔn),兩片點(diǎn)云幾乎完全重疊,迭代50次后均方誤差達(dá)到0.000 111 194 m,在運(yùn)行時(shí)間方面,文中程序所用的計(jì)算機(jī)配置為Intel(R)Core(TM)i5-7500處理器,16 GB內(nèi)存,直接ICP算法花費(fèi)時(shí)間為63 289 ms,經(jīng)過粗配準(zhǔn)之后再進(jìn)行ICP算法的花費(fèi)時(shí)間為48 011 ms。具體差異見表1所列。

表1 2種方法差異

5 結(jié) 論

由于多站數(shù)據(jù)的坐標(biāo)系不同,初始點(diǎn)云的位置差異比較大,直接進(jìn)行ICP算法無法準(zhǔn)確配準(zhǔn)拼接,甚至?xí)萑刖植孔顑?yōu)的問題。而僅僅用四元數(shù)法或七參數(shù)轉(zhuǎn)換坐標(biāo)進(jìn)行粗配準(zhǔn)并不能達(dá)到所要的精度。本文針對(duì)這個(gè)問題提出新的解決方法,即先利用特征點(diǎn)對(duì)應(yīng)的特征向量尋找點(diǎn)對(duì)并求解七參數(shù)的方法進(jìn)行粗配準(zhǔn),將點(diǎn)云位置之間的剛性變換變小,再通過ICP算法對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行精確配準(zhǔn)。這樣既可以解決直接進(jìn)行ICP算法陷入的局部最優(yōu)問題,也可以達(dá)到所需要的精度。