基于邻域旋转体积的点云配准方法与流程

技术特征：

1.一种基于邻域旋转体积的点云配准方法，其特征在于，所述点云配准方法包括：

步骤1，载入视角相邻的源点云source和目的点云target，并分别计算源点云的表面分辨率mr_src和目的点云的表面分辨率mr_tgt；

步骤2，对源点云source和目的点云target，分别以n*mr_src和n*mr_tgt为邻域，计算源点云的法向量normal_src和目的点云的法向量normal_tgt；

步骤3，利用Harris3D算法，根据源点云的法向量normal_src计算源点云的源关键点集合keypoints_src，根据目的点云的法向量normal_tgt计算目的点云的目的关键点集合keypoints_tgt；

步骤4，对源关键点集合keypoints_src和目的关键点集合keypoints_tgt中的每个关键点，分别以n倍表面分辨率为邻域计算出局部坐标系，并根据局部坐标系计算出邻域内点的旋转体积作为描述向量，得到源关键点描述向量集features_src和目的关键点描述向量集features_tgt；

步骤5，根据源关键点描述向量集features_src和目的关键点描述向量集features_tgt，在源关键点集合keypoints_src和目的关键点集合keypoints_tgt中确定多个匹配关键点对；根据匹配关键点对计算出从源关键点描述向量集features_src到目的关键点描述向量集features_tgt的旋转平移矩阵Rt；利用旋转平移矩阵Rt对源关键点描述向量集features_src进行旋转和平移，得到中间关键点描述向量集features_src'；计算中间关键点描述向量集features_src'中的每个描述向量与目的关键点描述向量集features_tgt中的对应描述向量之间的平均欧式距离；多次执行该步骤，取平均欧式距离最小时的旋转平移矩阵Rt作为最终的旋转平移矩阵；

步骤6，根据最终的旋转平移矩阵将源点云source变换到目的点云target所在的坐标系下，得到中间点云source′，以完成根据关键点对源点云source和目的点云target的粗配准；

步骤7，利用ICP算法对中间点云source′和目的点云target进行精确配准。

2.根据权利要求1所述的点云配准方法，其特征在于，所述步骤1中，分别计算源点云的表面分辨率mr_src和目的点云的表面分辨率mr_tgt，包括：

步骤1.1，对于源点云source和目的点云target中的每一个点，通过KD-tree算法找到其最近点；

步骤1.2，计算每一个点与其最近点之间的距离d_i；

步骤1.3，通过如下公式(1)计算点云的表面分辨率mr：

$mr=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{d}_i---\left(1\right)$

其中，N为点云中点的个数。

3.根据权利要求1所述的点云配准方法，其特征在于，所述步骤4中，对源关键点集合keypoints_src和目的关键点集合keypoints_tgt中的每个关键点，分别以n倍表面分辨率为邻域计算出局部坐标系，并根据局部坐标系计算出邻域内点的旋转体积作为描述向量，得到源关键点描述向量集features_src和目的关键点描述向量集features_tgt，包括：

对于源关键点集合keypoints_src和目的关键点集合keypoints_tgt中的任一关键点p，计算p的描述向量的过程通过如下步骤4.1至步骤4.4实现：

步骤4.1，建立p点及其邻域点Nbhd(p)的局部坐标系，具体为：

首先，以p的法向量为z轴；然后，在Nbhd(p)内找到法向量与p的法向量之间的夹角最小的关键点q，并连接p与q在p的切平面上的投影点q'，将该方向作为x轴；最后，根据z×x计算出y轴，形成p处的局部坐标系其中p为该局部坐标系的原点，代表x轴，代表y轴，代表z轴；

步骤4.2，对邻域球S进行空间几何区间划分，形成m个区间，其中，邻域球S是指Nbhd(p)形成的空间，邻域球S的圆心为p，半径为R，R＝n×mr，mr为源点云的表面分辨率mr_src或目的点云的表面分辨率mr_tgt，具体为：

首先，以的正方向和邻域球S的交点p_north为北极，的逆方向与邻域球S的交点p_south为南极，为中轴，沿着邻域球S球面的经线均匀地将邻域球S剖分为m个区间，该区间集合记为Region，每个独立区间记为Region_j(j∈[0m-1])；然后，将邻域球S正投影到upv平面后所得的圆沿圆心平行下移至以p_south为圆心处，记该圆为C_{prj_S}，该圆C_{prj_S}也即为圆心在p_south、半径为R的与邻域球S正切的圆；接下来，将Region_j分别投影到C_{prj_S}上，形成m个扇形，扇形集合记为Sector，每个扇形记为Sector_j(j∈[0m-1])，且Region_j与Sector_j一一对应；

步骤4.3，确定Nbhd(p)中每个点所属的区间，具体为：

首先，将Nbhd(p)中的任一点p_i，按照如下公式(2)投影到upv平面上，投影点记为q_i，并根据如下公式(3)计算出与坐标轴的逆时针夹角α_i；

然后，因为圆C_{prj_S}所在的平面与upv平面平行的关系，因此，可根据α_i的值计算q_i所属的扇形，从而确定出q_i所处的区间；接着，由p_i与q_i的一一对应性，确定出p_i所处的区间；其中，p_i所处的区间与q_i所处的区间相同；

步骤4.4，计算p的第k(k∈[0m-1])个维度的数值，具体为：

首先，对每个区间内的邻域点集中的每个点，计算其在upv平面上投影点与p的距离，将距离由小到大进行排序，得到每个区间的邻域点序列的各个元素，并将p点作为首元素加入到每个区间的邻域点序列中，形成每个区间的邻域点序列，该邻域点序列记为Nbhd(p)_k(k∈[0m-1])；

设任一区间Region_k的邻域点序列Nbhd(p)_k为{p,p₁,p₂,...,p_s}，按序依次取出p和p₁，并分别投影至C_{prj_S}平面上，投影点分为别为p_south和q₁，依次连接p_south、q₁、p₁和p形成一个梯形T₁；然后，将T₁绕着轴旋转360°，其扫过的体积记为V_k，1；同理，依次计算序列<p₁，p₂>、<p₂，p₃>，…，<p_s-1，p_s>扫过的体积V_k，2，V_k，3，…，V_k，s；接下来，利用如下公式(4)将所有体积累加从而计算出Region_k内所有邻域点序列扫过的体积和V_k，将其作为p的描述向量的第k个维度的数值；

其中，在计算Region_k中两相邻邻域点<p_i-1，p_i>及其投影点<q_i-1，q_i>所形成的梯形旋转后的模型的体积V_k,i时，可将该模型拆分为一个空心圆柱体Model_k,i,1和一个中间为空心圆柱的圆台Model_k,i,2；Model_k,i,1的体积V_k,i,1可通过如下公式(5)进行计算，其中，h_k,i,1代表Model_k,i,1的高，可通过如下公式(6)进行计算；r_k,i-1和r_k,i分别代表Model_k,i,1的内、外半径，均可通过如下公式(7)进行计算；Model_k,i,2的体积V_k,i,2可由圆台的体积V_k,i,3与圆柱的体积V_k,i,4相减而得，见如下公式(8)；如下公式(9)用于计算体积V_k,i,3，其中h_k,i,2代表圆台的高，r_k,i-1和r_k,i分别代表圆台的内、外半径，其也通过公式(7)进行计算；公式(10)用于计算体积V_k,i,4，其中h_k,i,2代表空心圆柱的高，其通过公式(11)进行计算，r_k,i-1代表圆柱的半径；另外，r_k,0＝0，h_k,0,1＝R；

$V_k=\underset{i=1}{\overset{S}{\&Sigma;}}{V}_{k,i}=\underset{i=1}{\overset{S}{\&Sigma;}}(V_{k,i,1}+V_{k,i,2})---\left(4\right)$

$V_{k,i,1}=\mathrm{\&pi;}*h_{k,i,1}*(r_{k,i}^2-r_{k,i-1}^2)---\left(5\right)$

$r_{k,i}=\sqrt{\left|\right|p_{south}{p}_i|{|}^2-h_{k,i,1}^2}---\left(7\right)$

V_k,i,2＝V_k,i,3-V_k,i,4 (8)

$V_{k,i,3}=\frac{1}{3}*\mathrm{\&pi;}*h_{k,i,2}*(r_{k,i}^2+r_{k,i-1}^2+r_{k,i}*r_{k,i-1})---\left(9\right)$

$V_{k,i,4}=\mathrm{\&pi;}*h_{k,i,2}*r_{k,i-1}^2---\left(10\right)$

h_k,i,2＝h_k,i-1,1-h_k,i,1 (11)；

步骤4.5，依据上述步骤4.4，计算p的其余m-1个维度的数值，最终组成p的m维的描述向量。

4.根据权利要求1所述的点云配准方法，其特征在于，所述步骤5中，根据源关键点描述向量集features_src和目的关键点描述向量集features_tgt，在源关键点集合keypoints_src和目的关键点集合keypoints_tgt中确定多个匹配关键点对；根据匹配关键点对计算出从源关键点描述向量集features_src到目的关键点描述向量集features_tgt的旋转平移矩阵Rt；利用旋转平移矩阵Rt对源关键点描述向量集features_src进行旋转和平移，得到中间关键点描述向量集features_src'；计算中间关键点描述向量集features_src'中的每个描述向量与目的关键点描述向量集features_tgt中的对应描述向量之间的平均欧式距离；多次执行该步骤，取平均欧式距离最小时的旋转平移矩阵Rt作为最终的旋转平移矩阵，包括：

步骤5.1，随机从源关键点描述向量集features_src中选取3组描述向量fs0、fs1和fs2，置入集合fs中，利用KD-Tree算法在目的关键点描述向量集features_tgt中为fs0搜索5个欧式距离最近的描述向量，并从5个欧式距离最近的描述向量中随机选取一个描述向量ft0作为fs0对应的描述向量；同理，从目的关键点描述向量集features_tgt中确定fs1对应的描述向量ft1和fs2对应的描述向量ft2，将ft0、ft1和ft2置入集合ft中；

步骤5.2，在源关键点集合keypoints_src中找到fs0、fs1和fs2分别对应的源关键点ks0、ks1和ks2，在目的关键点集合keypoints_tgt中找到ft0、ft1和ft2分别对应的目的关键点kt0、kt1和kt2，形成匹配关键点对<ks0，kt0>、<ks1，kt1>和<ks2，kt2>；依次连接ks0和ks1、ks1和ks2、ks2和ks0形成三条边并分别计算三条边的长度，记为d_ks01、d_ks12和d_ks20；同理，依次连接kt0和kt1、kt1和kt2、kt2和kt0形成三条边并分别计算三条边的长度，记为d_kt01、d_kt12和d_kt20；依次计算对应边长度相似比，即计算min(d_ks01/d_kt01，d_kt01/d_ks01)、min(d_ks12/d_kt12，d_kt12/d_ks12)和min(d_ks20/d_kt20，d_kt20/d_ks20)；如果相似比均大于预设阈值δ(δ∈[0 1])，则执行步骤5.3，否则返回步骤5.1；

步骤5.3，基于匹配关键点对<ks0，kt0>、<ks1，kt1>和<ks2，kt2>，通过SVD算法分解计算从源关键点描述向量集features_src到目的关键点描述向量集features_tgt的旋转平移矩阵Rt；

步骤5.4，利用旋转平移矩阵Rt对源关键点描述向量集features_src进行旋转和平移，得到中间关键点描述向量集features_src'；对中间关键点描述向量集features_src'中的每个描述向量，利用KD-Tree算法在目的关键点描述向量集features_tgt中搜索欧式距离最近的描述向量作为其对应的描述向量，并计算中间关键点描述向量集features_src'中的每个描述向量与目的关键点描述向量集features_tgt中的对应描述向量之间的平均欧式距离，将该平均欧式距离作为一个评价标准fitness；

步骤5.5，多次执行上述步骤5.1至5.4，选取最小fitness时的旋转平移矩阵作为最终的旋转平移矩阵。