基于点云能量计算的半刚性配准算法

林宝尉　王法胜　田　勇　肖　乐　赵方达

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