基于點云能量計算的半剛性配準算法

林寶尉　王法勝　田　勇　肖　樂　趙方達

1(大連東軟信息學院　遼寧 大連 116023)2(大連理工大學　遼寧 大連 116023)3(廣島大學　日本 7390041)

?

基于點云能量計算的半剛性配準算法

林寶尉1王法勝1田勇1肖樂2趙方達3

1(大連東軟信息學院遼寧 大連 116023)2(大連理工大學遼寧 大連 116023)3(廣島大學日本 7390041)

針對3D復雜點云模型，提出一種基于能量值的半剛性配準算法。該算法首先對點云進行基于法向量夾角距離的分割處理從而生成多個剛性子點云。通過整數規劃法計算子點云間的配準能量從而推測配準的初始3D變換矩陣。基于初始矩陣，點云的最終配準結果由子點云間剛性配準方法ICP實現。實驗結果證明該方法可有效應用于模擬場景以及真實復雜場景的配準應用中。

半剛性配準點云分割點云能量整數規劃法

0　引　言

近年來環境監測技術的發展對常規場景觀察起到了有效的推動作用。例如海洋、山脈或高速公路檢測技術的應用，對地質災害的防治起到了積極的作用。然而傳統的監控系統存在著容易受損、監控范圍窄、設備需求量高以及成本大的缺點，不適合在野外廣闊場景中使用。那么如何智能進行野外場景變化檢測就成為了要解決的重點問題。雖然已經有一些方法被用來進行相關檢測，但其往往是基于2D-2D[1]圖片的靜態檢測方法，或基于2D-3D[9]的SfM(Structure-from-Motion)方法。此類方法并不能針對同一場景不同時期的3D點云實現高精度配準。在本文中，被檢測的對象為海岸防浪堤的石塊堆(如圖7所示)，每個石塊的尺度約為3×3×3m3。普通的點云生成產品，如Kinect，并不能有效地對其進行測量；在僅有基于SfM技術重建的點云情況下，并沒有一種可行2D方法能夠對其進行精確檢測。不同于上述方法，3D-3D配準法往往會得到比較精確的結果。那么需要解決的問題即：如何針對基于SfM生成的復雜點云進行基于3D方法的配準，并實現變化檢測。

本文提出了一種基于點云能量值的算法來實現點云剛性配準的方法。3D點云配準技術已廣泛應用于3D場景重建、可視化、醫學圖像處理、機器人視覺以及增強現實等研究領域。所謂3D點云配準即對兩組3D坐標點進行匹配以達到3D場景重疊的操作。根據需求的目的不同，點云的配準操作可以有效地應用于廣域場景的點云拼接、老舊場景的變化檢測以及更新等。在傳統點云配準任務中主要分為剛性配準和非剛性配準兩種方法。

剛性配準主要應用于靜態場景3D點云的匹配。通過推測點云間的變換矩陣從而找出兩個3D場景中的重疊部分。該方法經常被應用于房屋、建筑、家具以及山體石塊等外形不易改變的場景對象。剛性匹配針對于剛性場景的對應關系來計算場景間的變換矩陣。其中Besl[2]等人提出的ICP(IterativeClosestPoint)方法以及ICP的衍生方法如Softassign[16]、EMICP[17]是應用最廣泛的剛性匹配算法。該類型方法在匹配過程中主要有兩個處理步驟：在每一個迭代中尋找近鄰點；根據近鄰點計算剛性變換矩陣使兩個場景的L2誤差達到最優值。該類型方法也稱做正交Procrestes問題。其缺點是需要提供配準的初始位置以及尺度猜測，并且對點云數據的質量有一定要求。不同于ICP，基于特征描述的匹配不需要提供初始變換矩陣。例如Spinimages[3],NARF[4],Shapecontext[5]等算法，通過計算3D點云關鍵點的局部特征描述向量來尋找關鍵點之間的對應從而對點云進行配準。這種方法往往需要一個預先設置的搜索臨界范圍來計算局部描述符。然而對于尺度未知的場景點云，臨界范圍的取值工作是非常困難的。為了解決尺度變化問題，一些研究提出了尺度不變的特征描述方法，如3DSIFT[6]、nDSIFT[7]、2D與3D結合法[8-9]等。雖然這些方法能夠在某種程度上解決尺度不變的問題，但其更側重整體3D場景的描述，或者需要完整的2D輔助信息，對于不是均勻分布、重建質量不高的點云并不能有效的進行配準。

非剛性配準作為配準問題中另外一個非常重要的應用，非剛性配準方法同樣運用在多種計算機視覺配準的研究中。如Sclaroff[10]等人提出的基于模式分析的2D圖片和3D點云對應問題研究、Huang[11-13]等研究提出的基于變形模型的3D外形推測法等。通過能量最優化、統計學等概率優化計算，非剛性匹配已經有效應用于人體器官、皮膚、衣服等配準檢測。然而在上述各方法中，完整的3D模型或3D變形對于非剛性的配準計算是至關重要的，在實際應用時，并不是所有的3D點云都擁有完整的點云信息。

半剛性匹配除了以上兩種情況，我們將現實存在的另外一種剛性場景變化定義為半剛性變化——在目標場景中存在的眾多剛性對象中的某個或某些對象發生的變化。對該半剛性變化進行的配準，稱之為半剛性匹配。例如本文中被檢測的防浪石塊對象，雖然該場景中的所有對象都是剛性的，然而傳統的剛性變化方法并不能準確對半剛性的兩個場景進行完全配準。同樣的，非剛性的方法也不能有效的作用于全部由剛性對象組成的3D場景。

本文提出了一種基于能量的3D點云場景配準方法。通過基于能量的整數規劃法，該方法可以有效地推測出場景中獨立對象間的一一對應關系，從而計算點云間的對應關系并計算配準矩陣。其創新處在于不僅對傳統的剛性配準有效，同時可以應用于上文中提出的半剛性配準，該方法適用于各種不同類型3D模型，對于不同尺度值的模型也是魯棒的。

1　點云能量計算

本文提出基于點云能量計算進行剛性及半剛性的匹配算法。該算法首先利用獨立的剛性塊對能量進行計算，獲取獨立剛性塊的3D分割技術介紹如下。

1.13D分割

算法第一步需要對3D場景進行分割。我們使用3D場景臨近點間的距離以及法向量夾角來進行分簇。每個簇代表場景中一個剛性物體對象。根據文獻[14]中的方法，場景中地面將作為潛在噪音被去除。如圖1所示，(a)中2D石塊照片被用來進行3D重建得到(b)圖結果。經過3D分割得到圖1(c)的結果。其中不同的圖像深度代表被分割出來的獨立的剛性塊石塊。該步驟作為處理的中間結果，其精度上的不足會在之后能量計算步驟中得到彌補。我們只需要得到初始的分割結果即可，即使圖中兩個簇被合并為一類，作為一個類同樣可以找到其相對應的概率最大配準對象。

圖1　3D場景及分割結果圖

1.2分割對應

完成分割后，需要計算剛性塊之間的對應關系。我們提出的是基于3D剛性塊間能量的整數規劃法來計算對應關系。這里定義的能量存在于一對或者一組剛性塊中，如果兩個點云場景中的某一對剛性塊的對應關系是正確的，那么這兩只剛性塊間擁有最大的相似度。計算場景能量的過程就是計算各個剛性塊間的相似度的過程。在接下來的小節中，我們將詳細介紹能量公式的定義以及使用。該公式的構建使用了與研究[24]中相似的構建方法，但是更有針對性。

(1) 問題公式化

設有兩個3D點云P和Q。P被分割成一組剛性塊P={Pα}，其中α=1,2,…,M，同樣的Q={Qβ}，β=1,2,…,N，并將每一個點云中的剛性塊的對應域表達為R=[1,2,…,M]×[1,2,…,N]。

我們使用一組二值向量x∈{0,1}MN來表達剛性塊的對應關系。每一組對應的下標使用a∈R來表示，xa∈x代表一對剛性塊的對應關系。如果一組a對應是有效的，那么xa=1，否則xa=0。

為了確保一個剛性塊只有一個對應關系，設定限制條件如下：

(1)

接下來將介紹代表獨立剛性塊間能量值的能量函數E(x)。整個能量函數由如下4組能量項組成。

(2) 能量項描述

E1(x)=∑a∈R(Cαβ(∑β=1Cαβ)-1)-1xa

(2)

能量項E2代表3D關鍵點間的L2誤差值。不同于E1中計算Pα與Qβ中對應點的個數，在這個能量項中將計算對應點間描述向量的歐式距離。具體的能量公式如下：

E2(x)=∑a∈RD1(Pα,Qβ)xa

(3)

其中D1代表Pα與Qβ中對應的3D關鍵點歐式距離。

能量項E3使用每一對或者一組3D剛性塊的剛性配準誤差來表達。如果對應關系正確，那么一個點云中的剛性塊必然有另一個點云中的剛性塊與其對應。換句話說，P中的Pα1和Pα2會有Q中的Qβ1和Qβ2與其對應。不同于單個剛性塊的匹配，將剛性塊間的關系變為連通圖，按圖的形式尋找配準將有助于提高配準精度，盡管預設需配準的兩個點云中某部分會發生變化，但對整體場景配準影響是可以被接受的。E3的定義如下:

E3(x)=∑a=(α1,β1)∈Rb=(α2,β2)∈RD2(Pα1,Pα2;Qβ1,Qβ2)xaxb

(4)

其中D2代表由ICP計算的剛性配準誤差。

為了避免當不存在對應關系時，算法卻仍然做配準計算，最后一個能量項E4引入了懲罰機制。換句話說，x中的所有元素為0但是能量值也變成了0，這種情況不是我們要得到的結果。為了避免該情況的發生，定義能量項E4為：

E4(x)=∑a∈Rε(1-xa)

(5)

其中ε代表一個懲罰系數。如果xa為0的比例較大，那么點云中的大部分剛性塊找不到其對應，則E4的值會增大，從而說明該組取值不是最優解。

通過合并式(2)-式(5)，最終得到的能量公式為：

E(x)=E1(x)+E2(x)+E3(x)+E4(x)

(6)

則成本函數為：

minE(x)=∑a∈REaxa+∑a,b∈REa,bxaxb

(7)

其中式(1)為其限制條件。由于式(7)屬于離散最優化問題，我們使用0-1整數規劃法來計算最優解。

圖2　人工合成的3D場景(虛線內石塊被移動)

圖3　初始匹配的3D結果(虛線內為被移動石塊)

圖4　模擬數據實驗結果

(3) 整數規劃法

為了求解上述離散問題，我們將dualdecomposition[15]的方法應用在能量最優化中。dualdecomposition方法將原問題分解為眾多subproblem，通過迭代求解subproblem的最優解來趨近原問題最優解。在此過程中，所有subproblem的上界與下界距離最小化即原問題的最優解。定義subproblemσ∈I,其中I為所有的subproblem域。則式(7)等價為:

(8)

對于某個subproblemσ(σ∈I)有：

(9)

式中，ρσ為能量系數。則各個σ下界為:

Φσ(Eσ)≤minE(x|Eσ)

(10)

(11)

通過式(8)-式(11)可以得出原問題的解的下界為：

(12)

式(9)-式(11)中的ρσ應與式(12)相同。

其中x*為原問題的最優解，通過重復的求σsubproblem的上界與下界最小距離Φσ(Eσ)即可求出x*。

此時的x*對于能量成本公式是最優解，然而對于整體的3D場景是初始解，即已經找到了各個剛性塊的對應關系。下一步要對每一個剛性塊進行精確剛性配準。

1.3剛性配準

通過初始對應關系x，可以得到相對準確的配準結果如圖3所示。接下來，將要進行精確的剛性配準。這一步中，每一個P中的Pα都可以找到一個Q中的對應Qβ。通過3D配準方法如RANSAC、Spinimages、3DSIFT或者ICP等算法即可對其進行精確配準。

2　實驗結果及分析

我們將通過實驗證明所提出的算法可以有效應用于人工合成以及真實場景數據。

2.1人工合成場景數據

該數據取自模擬的海岸防浪堤，人工水泥澆灌5×5×5cm3的石塊。如圖2所示，通過SfM技術重建出左邊3D場景后，虛線中石塊被挪動了位置，再次重建出右圖顯示的3D場景。兩個場景點云點數分別為15 804和25 178，點云的初始尺度比明顯不同。該模擬數據的信息重復性較高，對于傳統方法來說很難進行有效配準。通過該論文提出的配準算法，初始配準結果如圖3所示。可以看出經過移動的石塊能夠找到各自對應。最終的精確配準結果在圖4(a)中顯示，圖中兩點云完全重疊在一起。

作為對比，ICP以及2D-3D[9]算法也被應用于該模擬模型數據。失敗的配準結果如圖4(b)(c)所示，雖然ICP方法對于尺度匹配有一定效果，但在不提供初始位置猜測的情況下，結果明顯是錯誤；同樣，由于該點云的信息重復性較高，2D-3D的方法也無法提高準確的配準。傳統的配準算法不能提供初始的對應以及兩個點云的大小尺度比，往往會導致配準失敗。我們將200次的迭代所計算的L2能量顯示在配準能量圖5中，從圖中可以看出，ICP方法達到的局部最優解(能量值約為100)就趨于平緩，然而本文方法在通過前期能量計算的幫助下，可以在最終剛性配準階段迅速達到全局最優解并停止變化。這說明該配準方法對于模擬數據是有效的。

圖5　200次迭代的配準能量變化

圖7　真實場景復原三維圖

圖8　真實場景配準結果(圖7中兩點云重疊結果圖)

2.2真實場景數據

真實場景數據同樣選擇為海岸防浪堤的真實消波塊，在海岸線按規則堆放的消波塊對于海浪的侵蝕有很大的抵消作用。在相同的地點，圖6(a)拍攝于2009年10月，圖片的解析度為1600×1200，(b)拍攝于2013年12月，解析度為2256×1504。進行3D復原重建之后的場景如圖7所示。不同點云的點數分別為42 162和59 717，不同點云的尺度值同樣也為未知。由于真實數據的復雜性，2D-3D算法和ICP算法都很難匹配得到準確結果。如圖8所示，本文提出的配準方法可以有效地對真實數據進行配準。

3　結　語

本文中，我們提出了一種基于點云能量值整數規劃的半剛性配準算法來實現點云的重疊區域查找。為了提高配準精度，除了有針對性地提出能量成本公式以外，還借鑒了3D分割以及對偶分解的算法來實現配準。實驗證明，我們的方法不僅對模擬數據有效，同時對真實的復雜數據也是有效的。

對于真實數據，雖然整數規劃過程中可以降低3D分割結果對最終配準的影響，然而提高3D分割的準確性將對提高最終結果的處理速度與精度起到有效作用。另外，本文中使用的對偶分解法是否是最高效的求解方法并沒有得到驗證，如何選擇求解方法來提高運行速度將是我們今后的研究重點。

[1]BaoweiLin,YujiUeno,KSakai,etal.ImagebasedDetectionof3DSceneChange[J].IEEJTransactionsonElectronicsInformationandSystems,2013,133(1):103-110.

[2]BeslPJ,McKayND.AMethodforRegistrationof3-DShapes[J].IEEETransactionsonPatternAnalysisandMachineIntelligence(PAMI),1991,14(2):239-256.

[3]JohnsonAE,HebertM.Usingspinimagesforefficientobjectrecognitionincluttered3Dscene[J].IEEETransactionsonPatternAnalysisandMachineIntelligence(PAMI),1999,21(5):433-449.

[4]StederB,RusuRB,KonoligeK,etal.:NARF: 3Drangeimagefeaturesforobjectrecognition[C]//ProceedingsofWorkshoponDefiningandSolvingRealisticPerceptionProblemsinPersonalRoboticsatIROS,2010.

[5]BelongieS,MalikJ，PuzichaJ.Shapematchingandobjectrecognitionusingshapecontexts[J].IEEETransactionsonPatternAnalysisandMachineIntelligence(PAMI),2002,24(4):509-522.

[6]ScovannerP,AliS,ShahM.A3-dimensionalsiftdescriptoranditsapplicationtoactionrecognition[C]//Proceedingsofthe15thinternationalconferenceonMultimedia,ACM,2007:357-360.

[7]CheungW,HamarnehG.n-sift:n-dimensionalscaleinvariantfeaturetransform[J].IEEETrans.ImageProcess,2009,18(9):2012-2021.

[8]BaatzG,KoserK,ChenD,etal.Leveraging3Dcitymodelsforrotationinvariantplace-of-interestrecognition[J].InternationalJournalofComputerVision(IJCV),2012,96(3):315-334.

[9]BaoweiLin,ToruTamaki,MalcosSlomp,etal. 3DKeypointsDetectionfroma3DPointCloudforReal-TimeCameraTracking[J].IEEJTransactionsonElectronics,InformationandSystems,2013,133(1):84-90.

[10]SclaroffS,PentlandA.ModalMatchingforCorrespondenceandRecognition[J].IEEETransactionsonPatternAnalysisandMachineIntelligence(PAMI),1995,17(6):545-561.

[11]HuangQX,AdamsB,WickeM,etal.Non-RigidRegistrationUnderIsometricDefomations[J].ComputerGraphicsForum,2008,27(5):1449-1457.

[12]LiH,SumnerRW，PaulyM.GlobalCorrespondenceOptimizationforNon-rigidRegistrationofDepthScans[J].ComputerGraphicsForum,2008,27(5):1421-1430.

[13]EcksteinI,PonsJP,TongY,etal.GeneralizedSurfaceFlowsforMeshProcessing[J].EurographicsSymposiumonGeometryProcessing,2007,11-20,28.

[14]HolzD,HolzerS,RusuRB,etal.Real-TimesegmentationusingRGB-DCameras[C]//ProceedingsofRoboCupInternationalSymposium,Springer,2012:306-317.

[15]TorresaniL,KolmogorovV,RotherC.ADualDecompositionApproachtoFeatureCorrespondence[J].IEEETransactionsonPatternAnalysisandMachineIntelligence(PAMI),2013,35(2):259-271.

[16]GoldS,RangarajanA,LuCP,etal.NewAlgorithmfor2Dand3DPointMatching:PoseEstimationandCorrespondence[J].PatterRecognition,1998,31(8):1019-1031.

[17]GrangerS,PennecX.Multi-scaleEM-ICP:AFastandRobustApproachforSurfaceRegistration[C]//Proceedingsof7thEuropeanConferenceonComputerVision(ECCV),2002(4):69-73.

SEMI-RIGIDREGISTRATIONMETHODBASEDONPOINTCLOUDENERGYCALCULATION

LinBaowei1WangFasheng1TianYong1XiaoLe2ZhaoFangda3

1(Dalian Neusoft University of Information,Dalian 116023,Liaoning,China)2(Dalian University of Technology,Dalian 116023,Liaoning,China)3(Hiroshima University,Higashi-Hiroshiam,7390041,Japan)

Weproposedanenergyvalue-basedsemi-rigidregistrationmethodfor3Dcomplexpointcloudmodel.Firstthemethodsegmentsthepointcloudbasedondistanceofnormalvectoranglesoastogenerateacoupleofrigidsub-pointclouds.Thenitcalculatesthealignmentenergybetweensub-pointcloudswithintegerprogrammingmethodsothattospeculateontheinitial3Dtransformationmatrixofregistration.Basedoninitialmatrix,thefinialregistrationresultisimplementedbyrigidalignmentsmethodICPamongsub-pointclouds.Experimentalresultsdemonstratedthatthemethodcanbeeffectivelyappliedtotheregistrationapplicationsinsimulativescenesandrealcomplexscenes.

Semi-rigidregistrationPointcloudsegmentationPointcloudenergyIntegerprogramming

2014-09-25。國家自然科學科學基金項目(613000 82)；遼寧省教育廳科學研究項目(L2014575)；大連東軟信息學院博士啟動基金項目(ZX2014KJ003)。林寶尉，講師，主研領域：計算機視覺，模式識別，機器視覺。王法勝，講師。田勇，副教授。肖樂,博士生。趙方達，博士生。

TP3

ADOI:10.3969/j.issn.1000-386x.2016.03.042