一种基于割线特征相关的曲面点云配准方法

1.本发明涉及点云配准技术领域，具体涉及一种基于割线特征相关的曲面点云配准方法。


背景技术：

2.点云配准技术广泛应用于三维重建、目标跟踪与识别、文物修复等领域，存在很高的应用价值和广阔的应用前景。如何提高跟踪技术的准确性、实时性，以适应实际应用领域的各种需求是目前的点云配准领域关注的关键问题。
3.中国专利公开号为“cn115861397a”，名称为“一种基于改进fpfh-icp的点云配准方法”，该方法首先对点云进行体素滤波进行下采样，通过iss3d算法提取特征点，引入点云重心对特征点加以约束，剔除错误点对；然后用改进的快速点特征直方图(fpfh)对特征点构建特征描述子，通过sca-ia算法对点云进行粗配准，最后再对粗配准后点云利用改进的icp算法进行精确配准。该方法在特征点提取需全局搜索计算导致计算量大，而对特征点构建特征描述子时，fpfh算法计算复杂，计算量会随着特征点数量的增加而提高，该方法在最后使用了icp的改进方法进行全局配准，这对粗配准后两点云的位姿有一定要求，而且存在计算量大效率低的问题。


技术实现要素：

4.本发明为了解决现有技术的方法存在计算量大配准效率低的问题，提出一种基于割线特征相关的曲面点云配准计算方法。
5.本发明解决技术问题的技术方案是：
6.一种基于割线特征相关的曲面点云配准计算方法，其包括如下步骤：
7.步骤一，给定三角网格化目标点云p和待配准点云q；
8.步骤二，获取目标点云p和待配准点云q的割线；
9.首先，对目标点云p进行主成分分析(pca)得到三个主成分分别为n1(x1,y1,z1)、n2(x2,y2,z2)、n3(x3,y3,z3)，通过式(1)求出均值点c(c
x
,cy,cz)；
[0010][0011]
其中：pi为目标点云p中第i个点，n为pi的个数；
[0012]
然后，过点c分别以n1和n2为法向量建立a，b平面方程如式(2)；
[0013][0014]
过点c以n3为方向向量建立直线z方程如下：
[0015][0016]
以直线z为轴线，以r为半径建立圆柱面l方程如下：
[0017]
(x-xm)2 (y-ym)2 (z-zm)2＝r2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0018]
其中：(xm,ym,zm)为圆柱面上点(x,y,z)在轴线z上的投影点的坐标，
[0019]
最后，设与平面有交点的三角网格化点云中的边为ei，e
i0
和e
i1
为ei边的两个端点，为点云与平面的交点，则通过式(5)与式(2)联立得到交点坐标，将所有交点求出后便可得到平面a，b与点云p，q的交线；
[0020][0021]
其中：0≤λ≤1，设与圆柱面又交点的三角网格化点云中的边为be，be边上两个端点为ea和eb，pe为点云与圆柱面的交点，则通过式(6)与式(4)联立得到交点坐标，将所有交点求出后便可得到圆柱面l与点云p，q的交线；
[0022][0023]
其中：0≤k≤1，此处所得到的交线即为割线；
[0024]
步骤三，对割线进行重采样；
[0025]
以a平面所获得的成对割线a
p
、aq为例，分别对割线a
p
、aq用矩阵g
p
和矩阵gq来表示，并对矩阵g
p
和gq整体进行主成分分析并以重心点为原点，以三个主成分为坐标轴建立局部坐标系，通过式(7)将割线由世界坐标系转到局部坐标系得到矩阵j
p
，jq；
[0026][0027]
其中：ra为旋转矩阵，ta为平移矩阵；割线都是在平面内的，局部坐标系下的割线点的z值均为0，只需针对重合部分在第一主成分方向即x轴方向通过式(8)对x坐标以步长λ进行重新采样得到y值，最终得到局部坐标系下的重采样的割线点；
[0028]
yi＝f(xi i
×
λ)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0029]
其中：i＝(1,2,...,s)，s为重采样点的个数，f(x)为重采样函数；
[0030]
步骤四，对重采样割线构建割线特征模型；
[0031]
设重采样的割线点集合为d，通过式(9)获得集合d中元素的特征；
[0032][0033]
其中：di为d中第i个元素，ei为第i个元素特征值，d
pre
和d
next
为以di为圆心rc为半径圆与割线相交的点，得到割线的特征描述子为e＝{e1,e2,...,e
l
}，同一面所得的成对割线特征描述子则为该平面的割线特征模型g＝{e
ap
,e
aq
}；
[0034]
步骤五，对割线特征模型进行相关性分析；
[0035]
a平面的割线特征模型中目标点云的线特征描述子为e
ap
，待配准点云描述子为e
aq
，首先，通过式(10)对模型进行去中心化处理得到
[0036]
[0037][0038]
然后，通过式(11)对和进行相关性分析得到相关系数ε，其中0≤ε≤1；接下来对e
ap
和e
aq
按照λ步长做相对相移，每移动一步λ均计算相关系数，得到相关系数最大时的运动步长s
λ
，以及相关系数ε，此时的相关系数最大，两个割线是一一对应的，求出重合部分并分别从每个割线选取三个对应点对，得到对应点对集合t，o；
[0039]
步骤六，对对应点对通过奇异值分解求解配准的旋转矩阵和平移矩阵；
[0040]
首先，根据式(12)求出两个对应点对集合t，o的质心，并去中心化处理；
[0041][0042]
其中：ti为集合t中第i个点，oi为集合o中第i个点，t和o分别为集合t，o的质心，t
′
和o
′
分别为集合t和o去中心化后结果；
[0043]
然后，构建h矩阵：并对h矩阵进行奇异值分解获得u矩阵和v矩阵；
[0044]
最后，得到点云配准的旋转矩阵为r＝uv-1
，平移矩阵为
[0045]
步骤七，判断是否达到迭代次数或相邻误差差值以及误差均小于阈值；
[0046]
若没有达到迭代次数δ或相邻误差差值以及误差τ均小于阈值条件时，则根据式(13)对点云q进行坐标更新并回到步骤二，不断迭代步骤二至步骤七直到达到迭代条件时完成点云配准。
[0047]
q＝r
×
p t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0048]
本发明的有益效果在于：本发明利用主成分方向建立标准平面和圆柱面并求出各面与曲面点云的割线，计算割线上各点曲率并以此为特征建立割线特征模型进行相关性分析来进行配准，仅需要割线上的部分点参与配准，而不是所有点均参与计算，大大节省了计算量提升了配准速度，而且还能够得到较高的配准精度。
附图说明
[0049]
图1为本发明一种基于割线特征相关的曲面点云配准方法的流程图；
[0050]
图2为本发明所述配准前点云初始状态示意图；
[0051]
图3为本发明所述各面割线示意图；
[0052]
图4为本发明所述特征示意图；
[0053]
图5为本发明所述相关分析前a平面特征描述子示意图；
[0054]
图6为本发明所述相关分析后a平面特征描述子示意图；
[0055]
图7为本发明所述a平面配准结果示意图；
[0056]
图8为本发明所述点云配准误差变化示意图；
[0057]
图9为本发明所述点云配准结果示意图。
具体实施方式
[0058]
以下将结合附图对本发明的具体实施方式做进一步的详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。
[0059]
如图1所示，一种基于割线特征相关的曲面点云配准方法，包括如下步骤：
[0060]
步骤一，给定三角网格化目标点云p和待配准点云q。
[0061]
通过二维正态分布概率密度函数方程建立一个200
×
300大小的目标点云p，其中曲面方程如下式:
[0062][0063]
其中：k＝30，ρ＝0，μ1＝0，μ2＝0，σ1＝50，σ2＝100。再将点云p在沿x方向平移10个单位，沿y方向沿5个单位，绕z轴旋转15度，得到待配准点云q，p和q初始状态，如图2所示。
[0064]
步骤二，获取目标点云p和待配准点云q的割线。
[0065]
首先，对目标点云p进行主成分分析(pca)得到三个主成分分别为n1(x1,y1,z1)、n2(x2,y2,z2)、n3(x3,y3,z3)，通过式(1)求出均值点c(c
x
,cy,cz)。
[0066][0067]
其中：pi为点云p中第i个点，n为点云p的大小，然后，过点c分别以n1和n2为法向量建立a，b平面方程如式(2)。
[0068][0069]
过点c以n3为方向向量建立直线z方程如下：
[0070][0071]
以直线z为轴线，以r为半径建立圆柱面l方程如下：
[0072]
(x-xm)2 (y-ym)2 (z-zm)2＝r2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0073]
其中：(xm,ym,zm)为圆柱面上点在轴线z上的投影点的坐标；最后，设与平面有交点的三角网格化点云中的边为ei，e
i0
和e
i1
为ei边的两个端点，为点云与平面的交点，则通过式(5)与式(2)联立得到交点坐标，将所有交点求出后便可得到平面a，b与点云p，q的交线。
[0074][0075]
其中：0≤λ≤1，设与圆柱面有交点的三角网格化点云中的边为be，be边上两个端点为ea和eb，pe为点云与圆柱面的交点，则通过式(6)与式(4)联立得到交点坐标，将所有交点求出后便可得到圆柱面l与点云p，q的交线。
[0076][0077]
其中：0≤k≤1，此处所得到的交线即为割线，如图3所示。
[0078]
步骤三，对割线进行重采样。
[0079]
以a平面所获得的成对割线a
p
、aq为例，分别对割线a
p
、aq用矩阵g
p
和矩阵gq来表示，并对矩阵g
p
和gq整体进行主成分分析并以重心点为原点，以三个主成分为坐标轴建立局部坐标系，通过式(7)将割线由世界坐标系转到局部坐标系得到矩阵j
p
，jq。
[0080][0081]
其中：ra为旋转矩阵，ta为平移矩阵。因为割线都是在平面内的，所以局部坐标系下的割线点的z值均为0，因此，只需针对重合部分在第一主成分方向即x轴方向通过式(8)对x坐标以步长λ进行重新采样得到y值，最终得到局部坐标系下的重采样的割线点。
[0082]
yi＝f(xi i
×
λ)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
[0083]
其中：i＝(1,2,...,s)，s为重采样点的个数，f(x)为重采样函数。
[0084]
步骤四，对重采样割线构建割线特征模型。
[0085]
设重采样的割线点集合为d，通过式(9)获得集合d中元素的特征。
[0086][0087]
其中：di为d中第i个元素，ei为第i个元素特征值，d
pre
和d
next
为以di为圆心rc为半径圆与割线相交的点，如图4所示，因此，可得割线的特征描述子为e＝{e1,e2,...,e
l
}，而同一面所得的成对割线特征描述子则为该平面的割线特征模型g＝{e
ap
,e
aq
}。
[0088]
步骤五，对割线特征模型进行相关性分析。
[0089]
a平面的割线特征模型中目标点云的割线特征描述子为e
ap
，待配准点云特征描述子为e
aq
，a平面的两个割线特征描述子如图5所示，横坐标代表割线点的位置，纵坐标代表每个位置上的特征值。首先，通过式(10)对模型进行去中心化处理得到
[0090][0091][0092]
然后，通过式(11)对和进行相关性分析得到相关系数ε，其中0≤ε≤1，相关分析后a平面两特征描述子如图6所示。接下来对e
ap
和e
aq
按照λ步长做相对相移，每移动一步λ均计算相关系数，得到相关系数最大时的运动步长s
λ
＝553，以及相关系数ε＝0.9997，此时的相关系数最大，两个割线是一一对应的，求出重合部分并分别从每个割线选取三个对应点对，得到对应点对集合t，o。
[0093]
步骤六，对对应点对通过奇异值分解求解配准的旋转矩阵和平移矩阵。
[0094]
首先，根据式(12)求出两个对应点对集合t，o的质心，并去中心化处理。
[0095]
[0096]
然后，构建h矩阵：并对h矩阵进行奇异值分解获得u矩阵和v矩阵；最后，得到点云配准的旋转矩阵为r＝uv-1
，平移矩阵为svd求解方法求出a平面配准的旋转矩阵为平移矩阵为
[0097]
步骤七，判断是否达到迭代次数或相邻误差差值以及误差均小于阈值；
[0098]
若没有达到迭代次数δ或相邻误差差值以及误差τ均小于阈值条件时，则根据式(13)对点云q进行坐标更新并回到步骤二，不断迭代步骤二至步骤七直到达到迭代条件时完成点云配准。
[0099]
q＝r
×
p t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)
[0100]
点云配准过程中的均方根误差变化结果图，如图8所示，最终的误差为0.0658mm，总用时为2.191s，且最终点云配准结果图，如图9所示。