一种基于差分进化算法和TrimmedICP算法的点云配准方法与流程

本发明属于计算机图形学、测绘学、逆向工程等领域，涉及一种基于差分进化算法和trimmedicp算法的点云配准方法。

背景技术：

在物体三维重建领域获取目标表面深度信息常见的有两种方式：一是使用三维激光扫描仪进行全方位扫描，二是通过双目立体视觉方法利用图像匹配获取深度信息。对于前者来说由于物体的复杂性不可能一次扫描获取目标的整体点云数据，对于后者则更是需要从不同角度拍摄双目图片，以获取不同角度的点云信息，这就面临着要处理多测站点云配准问题。多测站点云配准问题就是要将不同位置获得的点云数据对齐合并在一个统一的坐标系下，这类问题可以分解为两两点云的配准问题。目前解决点云配准的方法可以分为两大类：公共点法和迭代最近点(icp)算法。公共点法要求两个点云间至少有3对的人工标志点或者特征点，可以通过最小二乘法求得转换矩阵参数，但是这种方法需要去单独识别特定的某些标志点并编号记录带来很大的数据复杂度；icp方法以最近距离标准为基础迭代地分配对应关系，并且获得关于两个点云的刚性变换最小二乘，然后再重新决定对应关系并继续迭代直到到达最小值，这种方式简单计算复杂度较低，但是icp算法前提要求点云重叠率大于50％，并且要给它设置一个较好的初值才不会使之陷入局部最优解。因此，如何获得一个可以解决点云重叠率较小情况下的高效、精确的点云配准方法是目前三维重建过程中亟待解决的问题。

技术实现要素：

本发明的目的是提供可以解决点云重叠率较小情况下的一种高效、精确的基于差分进化算法和trimmedicp(裁剪迭代最近点)算法的点云配准方法。

为了实现上述发明目的，本发明采用如下技术方案：

一种基于差分进化算法和trimmedicp算法的点云配准方法，包括如下步骤：

(1)对于参考点云d和配准点云m，设置差分进化算法的参数，包括种群大小pop、交叉概率cr、最大进化代数ng、缩放因子f、终止阈值t，以及trimmedicp算法参数包括的修剪阈值下界trmin和修剪阈值上界trmax；

(2)初始化种群，随机生成式中xj,i(0)表示第0代的第i个个体的第j个基因，表示第i个个体的第j个基因的上界，表示第i个个体的第j个基因的下界，rand(0,1)表示在(0,1)区间随机数；

(3)利用trimmedicp算法由初始种群得到精确变换参数r0′和t0′，其中r0′表示由初始种群得到的配准点云m进行变换的旋转矩阵，t0′表示由初始种群得到的配准点云m进行变换的平移矩阵；

(4)将配准点云m中的点进行一次空间变换d′＝r0′m+t0′，然后求出此时的最优配准代价cost(0)，其中d′表示点云m经过变换得到的点云，cost(0)表示由初始种群得到的最优配准代价；

(5)对此时的种群进行变异操作，随机选取种群中不同的两个个体，将其向量差缩放后与待变异个体进行向量合成：

vi(g+1)＝xr1(g)+f·(xr2(g)-xr3(g))，其中i≠r1≠r2≠r3；式中，vi(g+1)表示经过变异产生的第g+1代中第i个变异中间体，f表示缩放因子，xr1(g)表示第g代种群中第r1个个体，xr2(g)表示第g代种群中第r2个个体，xr3(g)表示第g代种群中第r3个个体；

(6)对此时的种群和变异中间体进行交叉操作：

式中，uj,i(g+1)表示经过交叉操作产生的第g+1代中第i个交叉中间体的第j个基因，vj,i(g+1)表示第g+1代中第i个变异中间体的第j个基因，xj,i(g)表示第j代种群中第i个个体的第j个基因，cr表示交叉概率，jrand表示[1,2，…，pop]的随机整数。

(7)对此时的种群和交叉中间体进行选择操作：

首先求取输入为ui(g+1)时的最优配准代价和输入为xi(g)时的最优配准代价，根据两个最优配准代价的大小关系选择进入下一代的个体，式中，xi(g+1)表示经过选择之后得到的g+1代种群第i个个体，ui(g+1)表示第g+1代第i个交叉中间体，xi(g)表示第g代种群第i个个体，f表示求取最优配准代价的函数，f(ui(g+1))表示输入为ui(g+1)时的最优配准代价，f(xi(g))表示输入为xi(g)时的最优配准代价；

(8)把(7)中生成的第g+1代种群作为初值，利用trimmedicp算法得到精确变换参数rg+1′和tg+1′，其中rg+1′表示得到的配准点云m进行变换的精确旋转矩阵，tg+1′表示得到的配准点云m进行变换的精确平移矩阵；

(9)将配准点云m中的点进行一次空间变换d′＝rg+1′m+tg+1′，然后求出g+1代的最优配准代价cost(g+1)，其中d′表示点云m经过变换得到的点云。

(10)判断是否满足终止条件：终止条件为第g代和第g+1代种群最优配准代价差值小于设置的阈值t或者代数达到设置的最大代数ng；若满足终止条件则进入(11)，若不满足则返回步骤(5)；

(11)输出g+1代最优配准代价以及精确变换参数rg+1′和tg+1′。

所述步骤(4)、步骤(7)和步骤(9)中求取最优配准代价的具体步骤是：

1)输入旋转矩阵初值ri和平移矩阵初值ti，经过trimmedicp算法处理获得了配准点云m进行精确配准所需要的精确旋转矩阵ri′和精确平移矩阵ti′；

2)建立变换方程：d′＝mri′+ti′，式中，d′为m经过变换得到的点云；

3)对于d′中的每一个点寻找距离参考点云d中最近的点，以欧式距离作为这个点的配准代价，则整体点云的最优配准代价：式中，cost表示点云最优配准代价，z表示点云d′中的任意一点，zclose表示点云d中距离点z空间最近的点，||z-zclose||2表示点z和点zclose之间欧氏距离。

本发明具有以下有益效果：

本发明应用差分进化算法和trimmedicp算法相结合的方式去解决点云配准的问题。trimmedicp算法可以解决点云配准中重叠率较低的问题，但是需要首先设置很好的旋转矩阵初值和平移矩阵初值，否则会陷入局部最优；差分进化算法可以通过随机方式生成初始种群，并且应用变异、交叉、选择操作使种群分布在所有的范围。两者相结合既解决了选取旋转矩阵和平移矩阵初值的困难问题，又可以获得全局的最优解，以获得精确的配准结果。

附图说明

图1为基于差分进化算法和trimmedicp算法的点云配准方法流程图；

图2为trimmedicp算法流程图；

图3为差分进化算法流程图；

图4为差分进化算法中交叉操作过程示意图；

图5(a)(b)为重叠率较低的两个点集在配准后的不同视角结果图；

图6(a)为本专利方法和基于粒子滤波的点云配准方法以及基于遗传算法的点云配准方法在平均进化代数上的比较；

图6(b)为本专利方法和基于粒子滤波的点云配准方法以及基于遗传算法的点云配准方法在平均耗费时间上的比较；

图6(c)为本专利方法和基于粒子滤波的点云配准方法以及基于遗传算法的点云配准方法在平均最优配准代价上的比较；

具体实施方式

下面结合附图对本发明做进一步详细描述：

如图1所示，对于给定的有一定重叠度的参考点云d和配准点云m，通过初始化操作生成旋转矩阵r0和平移矩阵t0；将初始值r0和t0输入到trimmedicp算法可以求出在此初值下的精确的旋转矩阵r0′和精确的平移矩阵t0′，但r0′和t0′有可能是局部最优解，求出此时的最优配准代价cost(0)；经过一代差分进化操作后获得更新后的初值：旋转矩阵r1和平移矩阵t1；经过trimmedicp算法得到精确配准所需的旋转矩阵r1′和平移矩阵t1′，并求出此时的最优配准代价cost(1)；比较cost(0)与cost(1)差值和阈值大小关系以确定是否满足终止条件，如果满足则输出cost(1)和对应的r1′、t1′，如果不满足继续进行差分进化操作获得更新后的初值：旋转矩阵r2和平移矩阵t2，然后通过trimmedicp算法求精确配准所需的旋转矩阵r2′和平移矩阵t2′，求出此时的最优配准代价cost(2)，判断是否满足终止条件，以此循环操作，直至满足终止条件。

本发明的一种基于差分进化算法和trimmedicp算法的点云配准方法，具体包括如下步骤：

s100：对于参考点云d和配准点云m，设置差分进化算法的参数，包括种群大小pop、交叉概率cr、最大进化代数ng、缩放因子f、终止阈值t，以及trimmedicp算法参数包括修剪阈值下界trmin和修剪阈值上界trmax。

s200：随机初始化种群。每个个体设置为1×7(1行7列)的矩阵，前3列分别表示x、y、z方向的旋转角度，取值范围是[-π，π]；4到6列分别表示x、y、z方向的平移量，取值范围根据两点云d和m数据的x、y、z分量方向上的最大差值确定；第7列为修剪参数，范围为[-0.5(trmax-trmin)，0.5(trmax-trmin)]。以上数据均为随机生成。然后，将初始种群转化为trimmedicp算法输入所需的格式，具体就是旋转角度转换为旋转矩阵：

式中，r0表示初始旋转矩阵，γβθ分别表示绕x、y、z轴旋转的角度，对应初始种群的前三列数据；设置初始平移矩阵t0，对应初始种群的4到6列数据。

s300：根据s200生成的初值r0和t0，通过trimmedicp算法求取精确的旋转矩阵r0′和精确的平移矩阵t0′，详细求解流程如图2所示：

1)设定变量值，包括当前重叠率phipre，最小重叠率minphi，重叠率差值阈值trphi，最大步数movestep。对点云m进行变换d′＝mr0+t0，式中d′为点云m变换之后的点云；

2)寻找d′中每个点对应点云d中最近的点，对应关系保存在corr中，每对点的距离保存在td中；

3)对td每行值的平方按由小到大排序，保存到sorttd2；

4)根据trmin和trmax计算修剪点云对的下界和上界：下界mintdindex＝floor(trmin×n)，上界maxtdindex＝ceiling(trmax×n)，其中，n表示点云m中点个数，mintdindex表示修剪点云索引号下界，maxtdindex表示修剪点云索引号上界，floor表示向下取整函数，ceiling表示向上取整函数；

5)计算所有点的均方误差：将sorttd2前j行的和保存在cumtd2的第j行，然后在cumtd2中取出步骤4)中确定的下界到上界中的元素得到cumtd2′，均方误差求解公式：

其中，j表示下界到上界中的一个索引号，cumtd2(j)表示cumtd2的第j行元素，mse(j-mintdindex)表示索引号为j-mintdindex的点的均方误差值。

6)更新最小重叠率minphi的值：其中tr(j-mintdindex)表示索引为j-mintdindex的重叠参数；

其中，mse(j)表示索引为j的均方误差，tr(j)表示索引为j的重叠参数，phi(j)表示索引为j的平均重叠率；赋值操作：phipre＝minphi，其中minphi表示最小重叠率phipre表示当前重叠率；搜索phi中最小值并赋给minphi，对应的索引为nindex。

7)更新标志重叠率，triksi＝tr(nindex)，其中triksi表示标志重叠率，tr(nindex)表示索引为nindex的重叠率。

8)获取当前重叠参数下裁剪后的点云对：将corr和td组成n×3的矩阵corrtd，按第3列td从小到大排序得到sortcorrtd。trlength＝floor(tr×n)，其中，trlength表示修剪长度，选取sortcorrtd前trlength行前两列元素命名为tcorr，。

9)应用奇异值分解求取精确的旋转矩阵r0′和精确的平移矩阵t0′：分别求两个裁剪后的m点云和d点云坐标均值，记为mm和md；进行中心化处理，点云中每个点坐标减去坐标均值转化为指向坐标均值的向量；mshifted和sshifted分别保存两个点云生成的向量组；k＝sshifted×mshifted′，对k′进行奇异值分解[uav]＝svd(k′)，其中u、a、v为k′进行svd奇异值分解后得到的矩阵；则r0′＝v×u′，t0′＝mm-r0′×md，其中u′表示矩阵u进行转置。

10)对点云m进行变换d′＝mr0′+t0′；运行步数加1：currstep＝currstep+1，其中currstep表示当前运行步数。

11)如果currstep<movestep并且abs(minphi-phipre)>trphi，则进入步骤2)，否则输出r0′和t0′，其中abs(minphi-phipre)表示最小重叠率和当前重叠率差值的绝对值大小。

s400：根据以上获得的r0′和t0′，建立变换方程：d′＝mr0′+t0′，式中，d′为m经过变换得到的点云；对于d′中的每一个点寻找距离d中最近的点以欧式距离作为这个点的配准代价，则点云整体最优配准代价：式中，cost(0)表示由初始种群得到的点云最优配准代价，z表示点云d′中的任意一点，zclose表示点云d中距离点z空间最近的点，||z-zclose||2表示点z和点zclose之间欧氏距离。。

s500：差分进化算法流程如图3所示，首先对此时的种群进行变异操作：

随机选取种群中不同的两个个体，将其向量差缩放后与待变异个体进行向量合成：vi(g+1)＝xr1(g)+f·(xr2(g)-xr3(g))，其中i≠r1≠r2≠r3。式中，vi(g+1)表示经过变异产生的第g+1代中第i个变异中间体，f表示缩放因子，xr1(g)表示第g代种群中第r1个个体，xr2(g)表示第g代种群中第r2个个体，xr3(g)表示第g代种群中第r3个个体。此过程要保证生成的中间体满足各分量的阈值范围，如果不满足则重新随机生成。

s600：对此时的种群和变异中间体进行交叉操作：

式中，uj,i(g+1)表示经过交叉操作产生的第g+1代中第i个交叉中间体的第j个基因，vj,i(g+1)表示第g+1代中第i个变异中间体的第j个基因，xj,i(g)表示第j代种群中第i个个体的第j个基因，cr表示交叉概率，jrand表示[1,2，…，pop]的随机整数。

s700：对此时的种群和交叉中间体进行选择操作，如图4所示，来生成第g+1代种群：

首先求取输入为ui(g+1)时的最优配准代价和输入为xi(g)时的最优配准代价，根据两个最优配准代价的大小关系选择进入下一代的个体，式中，xi(g+1)表示经过选择之后得到的g+1代种群第i个个体，ui(g+1)表示第g+1代第i个交叉中间体，xi(g)表示第g代种群第i个个体，f表示求取最优配准代价的函数，f(ui(g+1))表示输入为ui(g+1)时的最优配准代价，f(xi(g))表示输入为xi(g)时的最优配准代价。

s800：把s700中生成的第g+1代种群作为初值，利用trimmedicp算法可以得到精确变换参数rg+1′和tg+1′，其中rg+1′表示得到的配准点云m进行变换的精确旋转矩阵，tg+1′表示得到的配准点云m进行变换的精确平移矩阵。

s900：将配准点云m中的点进行一次空间变换d′＝rg+1′m+tg+1′，然后求出g+1代的最优配准代价cost(g+1)，其中d′表示点云m经过变换得到的点云。

s1000：判断是否满足终止条件。终止条件为第g代和第g+1代种群最优配准代价差值小于设置的阈值t或者代数达到设置的最大代数ng。若满足终止条件则进入s1100，若不满足则返回步骤s500。

s1100：输出g+1代最优配准代价以及精确变换参数rg+1′和tg+1′。

所述步骤s400、s700和s900中求取最优配准代价的具体步骤是：

1)输入旋转矩阵初值ri和平移矩阵初值ti，经过trimmedicp算法处理获得了配准点云m进行精确配准所需要的精确旋转矩阵ri′和精确平移矩阵ti′；

2)建立变换方程：d′＝mri′+ti′，式中，d′为m经过变换得到的点云；

3)对于d′中的每一个点寻找距离参考点云d中最近的点，以欧式距离作为这个点的配准代价，则整体点云的最优配准代价：式中，cost表示点云最优配准代价，z表示点云d′中的任意一点，zclose表示点云d中距离点z空间最近的点，||z-zclose||2表示点z和点zclose之间欧氏距离。

如图5(a)(b)所示，是利用所述的基于差分进化算法和trimmedicp算法的点云配准方法进行实验，所得的其中一组实验结果不同视角的示意图：两个点云在经过配准之后，重叠的区域经过空间变换移动到了一起。

图6(a)(b)(c)是对比实验结果，对比实验是将本专利方法和基于粒子滤波的点云配准方法以及基于遗传算法的点云配准方法在相同的10个数据集下进行实验，每个数据集用不同的方法分别运行10次，然后求三种方法在不同评价标准下评价量的平均值，评价标准分别为平均进化代数、平均耗费时间、平均最优配准代价。由图6(a)可以得到本专利方法相比另外两种方法需要较少的进化代数就能得到最优配准结果；由图6(b)可以得到本专利方法与基于遗传算法的配准方法平均耗费时间在不同数据集上互有优劣，但平均耗费时间都要比基于粒子滤波算法的配准方法少；由图6(c)可以得到本专利方法相比另外两种方法得到的最优配准代价最小，结果最准确。