CN108257163B - 一种已知扫描点位置下的两点点云配准方法 - Google Patents

一种已知扫描点位置下的两点点云配准方法 Download PDF

Info

Publication number
CN108257163B
CN108257163B CN201711291269.XA CN201711291269A CN108257163B CN 108257163 B CN108257163 B CN 108257163B CN 201711291269 A CN201711291269 A CN 201711291269A CN 108257163 B CN108257163 B CN 108257163B
Authority
CN
China
Prior art keywords
point
registered
station
point cloud
points
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201711291269.XA
Other languages
English (en)
Other versions
CN108257163A (zh
Inventor
姜光
赵晓娜
才长帅
贾静
彭亲利
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jiangsu Haohan Information Technology Co ltd
Original Assignee
Xidian University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Priority to CN201711291269.XA priority Critical patent/CN108257163B/zh
Publication of CN108257163A publication Critical patent/CN108257163A/zh
Application granted granted Critical
Publication of CN108257163B publication Critical patent/CN108257163B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10028—Range image; Depth image; 3D point clouds

Abstract

本发明公开了一种已知扫描点位置下的两点点云配准方法，其实现步骤为：(1)获取待扫描物体的点云数据；(2)预处理点云数据；(3)计算迭代次数；(4)判断当前迭代次数是否达到阈值，若是，则执行步骤(12)，否则，执行步骤(5)；(5)计算第一次旋转矩阵；(6)计算第二次旋转矩阵；(7)计算待配准点云的旋转矩阵；(8)计算当前正确对应点个数；(9)计算当前代价评估值；(10)更新参数；(11)判断是否达到迭代终止条件，若是，执行步骤(12)，否则，执行步骤(4)；(12)配准两站待配准点云。本发明需要的样本个数仅为2，采样过程中取得误匹配概率降低，配准结果可靠性高，具有配准结果精度高的优点。

Description

一种已知扫描点位置下的两点点云配准方法
技术领域
本发明属于图像处理技术领域，更进一步涉及计算机视觉技术领域中的一种已知扫描点位置下的两点点云配准方法。本发明可以应用于3D重建等具体应用场景，针对扫描点空间位置确知的相邻点云，实现对两站三维数据同时旋转配准，从而构建完整的三维点云模型。
背景技术
点云配准是将多视角扫描到的多站点云转换到统一的坐标系中，得到完整的三维模型。随机抽样一致性算法RANSAC在点云配准中是一种从样本中准确拟合出数学模型的方法。Chen首次将RANSAC思想应用于点云配准中，之后Aiger又提出了基于R ANSAC的四点法。目前提出的基于RANSAC思想的点云配准方法至少需要三组对应点，需要的样本点越少，取得误匹配点的概率就越低。
北京建筑大学在其拥有的专利技术“一种基于球标靶探测的激光点云数据自动配准方法”(专利申请号：201310746219.1，授权公告号：103646156B)中公开了一种基于球标靶探测的激光点云数据自动配准方法，包括以下步骤：分别从相邻两个站点的激光点云数据中提取球标靶特征；通过寻找相邻点云中由任意三个非共线球标靶球心构成的同名三角形，计算待配准站和参考站的变换关系。该方法通过在场景中设置标靶计算相邻点云的变换坐标，破坏了场景信息。并且该方法通过球形拟合的方式获得球标靶的球心，因此该方法存在的不足之处是，具有一定误差，无法保证配准结果的精确度。
王新、刘永山、朱代春、刘新、张小潮在其发表的论文“三维点云配准算法的研究”(燕山大学学报，2016,40(6):524-531)中提出了一种基于几何特征和RANSAC思想的粗配准算法以及基于点的领域几何特征的迭代配准方法。该方法在初始配准过程中采用的三点RANSAC，每次需要采样3组对应点作为样本点，通过判断由样本点构成的两个三角形是否满足空间几何中全等三角形的性质，将满足该性质的样本点视为合理的初始对应关系，并对合理的样本点进行刚体变换矩阵的求解。由于该方法采样需要的样本点越多，样本点之间的约束关系越复杂，存在的误差就越多，并且使得在采样过程中取得误匹配的概率变高，因此该方法存在的不足之处是，所用的三点RAN SAC在采样过程中取得误匹配点的概率高，存在的误差多，配准效果难以达到理想目标，初始配准精度差。
发明内容
本发明的目的在于克服上述已有技术的不足，提出了一种已知扫描点位置下的两点点云配准方法。
本发明实现的思路是，利用地面三维激光扫描仪获得待扫描物体的点云数据，并使用地面三维激光扫描仪配备的高精度全球定位系统GPS设备确定对应站点的扫描点位置，取相邻两站点云作为待配准点云，通过局部特征描述匹配的方法获取待两站待配准点云的对应点集，通过不断在对应点集中随机采样两对对应点作为样本，将样本绕扫描点位置旋转配准，直到满足迭代条件为止，将获得的最优旋转矩阵作用于两站待配准点云，完成配准。
本发明实现的具体步骤包括如下：
(1)获取待扫描物体的点云数据：
用地面三维激光扫描仪获取待扫描物体的点云数据；
(2)预处理点云数据：
(2a)从点云数据中任意选取一个站点的点云数据，将选取的点云数据的原点平移至所选取站点的扫描点位置，将平移后的点云数据作为第一站待配准点云；
(2b)从第一站待配准点云平移前的点云数据附近任意选取一个站点的点云数据，将选取的点云数据的原点平移至所选取站点的扫描点位置，将平移后的点云数据作为第二站待配准点云；
(2c)利用局部特征匹配方法，提取第一站待配准点云中与第二站待配准点云中相匹配的点对，组成两站待配准点云的对应点集，对应点集中包含正确对应点和错误对应点；
(3)利用迭代次数计算函数，计算迭代次数，确定迭代次数的最大值和最小值；
(4)判断当前迭代次数是否小于最小迭代次数或小于上次最大迭代次数，若是，则执行步骤(5)；否则，执行步骤(12)；
(5)计算两站待配准点云的第一次旋转矩阵：
(5a)从两站待配准点云的对应点集中，任意选取一组对应点作为样本A；
(5b)以第一站待配准点云的扫描点位置为球心O1，第一站待配准点云的扫描点位置到样本A中第一个点的欧式距离为半径，做球面S1；
(5c)以第二站待配准点云的扫描点位置为球心O2，第二站待配准点云的扫描点位置到样本A中第二个点的欧式距离为半径，做球面S2；
(5d)判断两个球面的半径之和是否大于两个球心之间的欧式距离，若是，则执行步骤(5e)；否则，执行步骤(5a)；
(5e)在两个球面的相交圆环的边缘上任取一点作为配准点；
(5f)按照下式，计算样本A中的第一个点、球心O1、配准点三点所在平面的法向量：
Figure GDA0002222439640000031
其中，n1表示样本A中的第一个点、球心O1、配准点三点所在平面的法向量，→表示向量符号，×表示向量叉乘操作，|·|表示向量的单位化操作，P1表示样本A中的第一个点，O1表示球面S1的球心，C表示配准点；
(5g)按照下式，计算球心O1和样本A中的第一个点确定的直线与球心O1和配准点确定的直线间的夹角：
Figure GDA0002222439640000032
其中，α1表示球心O1和样本A中的第一个点确定的直线与球心O1和配准点确定的直线间的夹角，arccos表示反余弦操作，r1表示球面S1的半径，d1表示样本A中第一个点到配准点的欧式距离；
(5h)将法向量n1和夹角α1代入绕任意轴旋转矩阵，得到第一站待配准点云的第一次旋转矩阵；
(5i)按照下式，计算样本A中的第二个点、球心O2、配准点三点所在平面的法向量：
Figure GDA0002222439640000041
其中，n2表示样本A中的第二个点、球心O2、配准点三点所在平面的法向量，→表示向量符号，×表示向量叉乘操作，|·|表示向量的单位化操作，Q1表示样本A中的第二个点，O2表示球面S2的球心；
(5j)通过下式，计算球心O2和样本A中第二个点确定的直线与球心O2和配准点确定的直线间的夹角：
Figure GDA0002222439640000042
其中，α2表示球心O2和样本A中第二个点确定的直线与球心O2和配准点确定的直线间的夹角，r2表示球面S2的半径，d2表示样本A中第一个点到配准点的欧式距离；
(5k)将法向量n2和夹角α2代入绕任意轴旋转矩阵，得到第二站待配准点云的第一次旋转矩阵；
(6)计算两站待配准点云的第二次旋转矩阵：
(6a)从两站待配准点云的对应点集中，任意选取除样本A之外的一组对应点，作为样本B；
(6b)以球心O1到匹配点的向量为旋转轴，将样本B中的第一个点旋转一个
Figure GDA0002222439640000043
角度，其中，arctan表示反正切操作，F表示点云数据的精度，d表示球心O1到球心O2的欧式距离；
(6c)以球心O2到匹配点的向量为旋转轴，将样本B中的第二个点旋转
Figure GDA0002222439640000044
角度，记录本次样本B中两个点的欧式距离及其对应的两个点的旋转角度；
(6d)判断样本B中的第二个点的旋转角度是否大于360°，若是，则将样本B中第二个点的旋转角度置为0后执行步骤(6e)；否则，执行步骤(6c)；
(6e)判断样本B中的第一个点的旋转角度是否大于360°，若是，则执行步骤(6f)；否则，执行步骤(6b)；
(6f)在所有旋转记录中，寻找样本B中两个点的欧式距离最小时对应的样本B中第一个点的旋转角度β1和第二个点的旋转角度β2；
(6g)将球心O1到匹配点的向量和旋转角度β1，代入绕任意轴旋转矩阵中，得到第一站待配准点云的第二次旋转矩阵；
(6h)将球心O2到匹配点的向量和旋转角度β2，代入绕任意轴旋转矩阵中，得到第二站待配准点云的第二次旋转矩阵；
(7)计算当前的两站待配准点云的旋转矩阵：
(7a)将第一站待配准点云的第二次旋转矩阵与第一次旋转矩阵相乘，得到第一站待配准点云的旋转矩阵；
(7b)将第二站待配准点云的第二次旋转矩阵与第一次旋转矩阵相乘，得到第二站待配准点云的旋转矩阵；
(8)计算当前正确对应点个数：
(8a)对两站待配准点云的对应点集进行配准，得到配准后的对应点集：
(8b)将配准后的对应点集中所有满足欧式距离阈值的对应点数，作为当前正确对应点的个数；
(9)利用代价函数，计算当前代价评估值；
(10)更新参数：
(10a)判断当前正确对应点个数是否等于上次正确对应点个数，若是，执行步骤(10b)；否则，执行步骤(10c)；
(10b)判断当前代价评估值是否小于上次代价评估值，若是，用当前代价评估值更新上次代价评估值后执行步骤(10e)；否则，执行步骤(10c)；
(10c)判断当前正确对应点个数是否大于上次正确对应点个数，若是，执行步骤(10d)；否则，执行步骤(11)；
(10d)用当前正确对应点个数更新上次正确对应点个数，利用迭代次数计算函数，计算当前最大迭代次数，用当前最大迭代次数更新上次最大迭代次数；
(10e)用当前两站待配准点云的旋转矩阵更新上次两站待配准点云的旋转矩阵；
(11)判断当前正确对应点个数是否等于两站待配准点云的对应点集包含的点对数，若是，执行步骤(12)；否则，将当前迭代次数加1后执行步骤(4)；
(12)配准两站待配准点云：
按照下式，对两站待配准点云中的每个点进行配准，完成两站待配准点云的配准：
Xi＝R1*(Pi-O1)+O1
Yi＝R2*(Qi-O2)+O2
其中，Xi表示对第一站待配准点云中第i个点进行配准后得到的点，R1表示第一站待配准点云的旋转矩阵，Pi表示第一站待配准点云中第i个点，Yi表示对第二站待配准点云中第i个点进行配准后得到的点，R2表示第二站待配准点云的旋转矩阵，Qi表示第二站待配准点云中第i个点。
本发明与现有技术相比具有如下优点：
第一，由于本发明在获取待扫描物体的点云数据时，仅利用地面三维激光扫描仪进行采集，克服了现有技术点云配准时，通过在场景中设置的标靶计算相邻点云的变换坐标，破坏场景信息的缺点，使得本发明点云配准方法可以更好的保护场景信息，提高了配准结果真实性。
第二，由于本发明在计算两站待配准点云的旋转矩阵过程中，仅获取了两组样本，克服了现有技术点云配准时，采集的样本点较多，取得误匹配点的概率较大的缺点，使得本发明在点云配准过程中取得误匹配点的概率更小，提高了配准结果可靠性。
第三，由于本发明在计算待配准点云旋转矩阵过程中，利用高精度全球定位系统GPS设备确定地面三维激光扫描仪扫描点中心位置作为约束条件，克服了现有技术点云配准时，约束关系复杂的缺点，使得本发明在点云配准过程中约束关系简单精确，提高了配准结果的精确性。
附图说明
图1是本发明的流程图；
图2是本发明两站待配准点云的第一次旋转矩阵的示意图；
图3是本发明两站待配准点云的第二次旋转矩阵的示意图；
图4是本发明计算两站待配准点云的第二次旋转矩阵步骤的流程图。
具体实施方式
下面结合附图对本发明作进一步的描述。
参照图1，本发明的具体实现步骤如下。
步骤1，获取待扫描物体的点云数据。
用地面三维激光扫描仪获取待扫描物体的点云数据。
步骤2，预处理点云数据。
从点云数据中任意选取一个站点的点云数据，将选取的点云数据的原点平移至所选取站点的扫描点位置，将平移后的点云数据作为第一站待配准点云，扫描点位置是指，利用地面三维激光扫描仪配备的高精度全球定位系统GPS设备所确定的站点对应的地面三维激光扫描仪的扫描点位置。
从第一站待配准点云平移前的点云数据附近任意选取一个站点的点云数据，将选取的点云数据的原点平移至所选取站点的扫描点位置，将平移后的点云数据作为第二站待配准点云。
利用局部特征匹配方法，提取第一站待配准点云中与第二站待配准点云中相匹配的点对，组成两站待配准点云的对应点集，对应点集中包含正确对应点和错误对应点，局部特征匹配方法是指，使用局部特征描述子对两站待配准点云中的点做局部特征描述，用局部特征对两站待配准点云之间的点进行相似性匹配，将所有的匹配点的集合作为两站待配准点云的对应点集。
步骤3，利用下式，计算迭代次数，确定迭代次数的最大值和最小值：
Figure GDA0002222439640000071
其中，K表示迭代次数，
Figure GDA0002222439640000072
表示向上取整操作，log表示以10为底的对数操作，η0表示两站待配准点云的对应点集中的所有点均为正确对应点的概率，其取值范围为[0.95,0.99]，
Figure GDA0002222439640000073
表示从m个不同元素中取出n个元素的组合数操作，N表示两站待配准点云的对应点集包含的对应点个数，Nin在计算最小迭代次数时表示预估计的正确对应点个数，在计算最大迭代次数时表示当前正确对应点个数。
步骤4，判断当前迭代次数是否小于最小迭代次数或小于上次最大迭代次数，若是，则执行步骤5；否则，执行步骤12。
上次最大迭代次数是指，在第一次迭代时的最大迭代次数为106，其余迭代过程中的最大迭代次数为当前迭代之前更新的最大迭代次数。
步骤5，计算两站待配准点云的第一次旋转矩阵。
为了不破坏已知扫描点位置信息，在配准第一组对应点时，将对应点绕各自的扫描位置旋转进行配准。
参照图2，对本发明计算两站待配准点云的第一次旋转矩阵过程作进一步的描述。
(5.1)从两站待配准点云的对应点集中，任意选取一组对应点作为样本A。
(5.2)如图2所示，以第一站待配准点云的扫描点位置为球心O1，第一站待配准点云的扫描点位置到样本A中第一个点的欧式距离为半径，做球面S1。
(5.3)如图2所示，以第二站待配准点云的扫描点位置为球心O2，第二站待配准点云的扫描点位置到样本A中第二个点的欧式距离为半径，做球面S2。
(5.4)判断两个球面的半径之和是否大于两个球心之间的欧式距离，若是，则执行本步骤的(5.5)；否则，执行本步骤的(5.1)。
(5.5)当两个球面的半径之和大于两个球心之间的欧式距离时，如图2中所示，两个球面会相交于一个圆环，在两个球面的相交圆环的边缘上任取一点作为配准点。
(5.6)按照下式，计算样本A中的第一个点、球心O1、配准点三点所在平面的法向量：
Figure GDA0002222439640000081
其中，n1表示样本A中的第一个点、球心O1、配准点三点所在平面的法向量，→表示向量符号，×表示向量叉乘操作，|·|表示向量的单位化操作，P1表示样本A中的第一个点，O1表示球面S1的球心，C表示配准点。
(5.7)按照下式，计算球心O1和样本A中的第一个点确定的直线与球心O1和配准点确定的直线间的夹角：
Figure GDA0002222439640000091
其中，α1表示球心O1和样本A中的第一个点确定的直线与球心O1和配准点确定的直线间的夹角，arccos表示反余弦操作，r1表示球面S1的半径，d1表示样本A中第一个点到配准点的欧式距离。
(5.8)将法向量n1和夹角α1代入绕任意轴旋转矩阵，得到第一站待配准点云的第一次旋转矩阵，任意旋转矩阵如下：
Figure GDA0002222439640000092
nx 2+ny 2+nz 2＝1
其中，R表示绕任意轴旋转任意角度的旋转矩阵，nx表示旋转轴在x轴方向的投影，ny表示旋转轴在y轴方向的投影，nz表示旋转轴在z轴方向的投影，
Figure GDA0002222439640000095
表示旋转角度。
(5.9)按照下式，计算样本A中的第二个点、球心O2、配准点三点所在平面的法向量：
Figure GDA0002222439640000093
其中，n2表示样本A中的第二个点、球心O2、配准点三点所在平面的法向量，→表示向量符号，×表示向量叉乘操作，|·|表示向量的单位化操作，Q1表示样本A中的第二个点，O2表示球面S2的球心。
(5.10)通过下式，计算球心O2和样本A中第二个点确定的直线与球心O2和配准点确定的直线间的夹角：
Figure GDA0002222439640000094
其中，α2表示球心O2和样本A中第二个点确定的直线与球心O2和配准点确定的直线间的夹角，r2表示球面S2的半径，d2表示样本A中第一个点到配准点的欧式距离。
(5.11)将法向量n2和夹角α2代入绕任意轴旋转矩阵，得到第二站待配准点云的第一次旋转矩阵。
步骤6，计算两站待配准点云的第二次旋转矩阵。
为了不破坏已知扫描点位置信息和步骤5中已经配准好的第一组对应点，在配准第二组对应点时，以扫描点位置到匹配点的向量为旋转轴，将第二组对应点配准到一起。参照图3的示意图和图4的流程图，对本发明计算两站待配准点云的第二次旋转矩阵过程作进一步的描述。
(6.1)从两站待配准点云的对应点集中，任意选取除样本A之外的一组对应点，作为样本B。
(6.2)如图3所示，以球心O1到匹配点的向量为旋转轴，将样本B中的第一个点旋转一个
Figure GDA0002222439640000101
角度，其中，arctan表示反正切操作，F表示点云数据的精度，d表示球心O1到球心O2的欧式距离，图中虚线所指的方向为对应点的旋转方向。
(6.3)如图3所示，以球心O2到匹配点的向量为旋转轴，将样本B中的第二个点旋转
Figure GDA0002222439640000102
角度，记录本次样本B中两个点的欧式距离及其对应的两个点的旋转角度，图中虚线所指的方向为对应点的旋转方向。
(6.4)判断样本B中的第二个点的旋转角度是否大于360°，若是，则将样本B中第二个点的旋转角度置为0后执行本步骤的(6.5)；否则，执行本步骤的(6.3)。
(6.5)判断样本B中的第一个点的旋转角度是否大于360°，若是，则执行本步骤的(6.6)；否则，执行本步骤的(6.2)。
(6.6)在所有旋转记录中，寻找样本B中两个点的欧式距离最小时对应的样本B中第一个点的旋转角度β1和第二个点的旋转角度β2。
(6.7)将球心O1到匹配点的向量和旋转角度β1，代入绕任意轴旋转矩阵中，得到第一站待配准点云的第二次旋转矩阵。
(6.8)将球心O2到匹配点的向量和旋转角度β2，代入绕任意轴旋转矩阵中，得到第二站待配准点云的第二次旋转矩阵。
步骤7，计算当前的两站待配准点云的旋转矩阵。
将第一站待配准点云的第二次旋转矩阵与第一次旋转矩阵相乘，得到第一站待配准点云的旋转矩阵。将第二站待配准点云的第二次旋转矩阵与第一次旋转矩阵相乘，得到第二站待配准点云的旋转矩阵。
步骤8，计算当前正确对应点个数。
按照下式，对两站待配准点云的对应点集进行配准，得到配准后的对应点集：
M＝{(R1*(Pi-O1)+O1,R2*(Qi-O2)+O2)}
其中，M表示两站待配准点云的对应点集进行配准后的对应点集，{}表示集合符号，(·,·)表示对应点对符号，R1表示当前第一站待配准点云的旋转矩阵，Pi表示在两站待配准点云的对应点集中属于第一站待配准点云的第i个点，R2表示当前第二站待配准点云的旋转矩阵，Qi表示在两站待配准点云的对应点集中属于第二站待配准点云的第i个点。
将配准后的对应点集中所有满足欧式距离阈值的对应点数，作为当前正确对应点的个数。
步骤9，利用代价函数，计算当前代价评估值，代价函数如下：
Figure GDA0002222439640000111
其中，J表示两站待配准点云的完成一次配准的代价评估值，N表示两站待配准点云的对应点集包含的点对数，∑表示求和操作，Ri在配准后的对应点集中第i个对应点的欧式距离大于对应点的欧式距离阈值时，Ri表示对应点的欧式距离阈值，否则，Ri表示配准后的对应点集中第i个对应点的欧式距离；Ij在配准后的对应点集中第i个对应点的欧式距离大于对应点的欧式距离阈值时，Ij表示0，否则，Ij表示1。
步骤10，更新参数。
(10.1)判断当前正确对应点个数是否等于上次正确对应点个数，若是，执行本步骤的(10.2)；否则，执行本步骤的(10.3)。上次正确对应点个数是指，在第一次迭代时的正确对应点个数为2，其余迭代过程中的正确对应点个数为当前迭代之前更新的正确对应点个数。
(10.2)判断当前代价评估值是否小于上次代价评估值，若是，用当前代价评估值更新上次代价评估值后执行本步骤的(10.5)；否则，执行本步骤的(10.3)。上次代价评估值是指，第一次迭代时的代价评估值为106，其余迭代过程中的代价评估值为当前迭代之前更新的代价评估值。
(10.3)判断当前正确对应点个数是否大于上次正确对应点个数，若是，执行本步骤的(10.4)；否则，执行步骤11。
上次正确对应点个数是指，在第一次迭代时的正确对应点个数为2，其余迭代过程中的正确对应点个数为当前迭代之前更新的正确对应点个数。
(10.4)用当前正确对应点个数更新上次正确对应点个数，利用迭代次数计算函数，计算迭代次数，用计算得到的迭代次数更新上次最大迭代次数。
上次正确对应点个数是指，在第一次迭代时的正确对应点个数为2，其余迭代过程中的正确对应点个数为当前迭代之前更新的正确对应点个数。上次最大迭代次数是指，在第一次迭代时的最大迭代次数为106，其余迭代过程中的最大迭代次数为当前迭代之前更新的最大迭代次数。
迭代次数计算函数如下：
Figure GDA0002222439640000121
其中，K表示迭代次数，
Figure GDA0002222439640000122
表示向上取整操作，log表示以10为底的对数操作，η0表示两站待配准点云的对应点集中的所有点均为正确对应点的概率，其取值范围为[0.95,0.99]，
Figure GDA0002222439640000123
表示从m个不同元素中取出n个元素的组合数操作，N表示两站待配准点云的对应点集包含的对应点个数，Nin在计算最小迭代次数时表示预估计的正确对应点个数，在计算最大迭代次数时表示当前正确对应点个数。
(10.5)用当前两站待配准点云的旋转矩阵更新上次两站待配准点云的旋转矩阵。上次两站待配准点云的旋转矩阵是指，在第一次迭代时的两站待配准点云的旋转矩阵为步骤5计算出的两站待配准点云的第一次旋转矩阵，其余迭代过程中的两站待配准点云的旋转矩阵为当前迭代之前更新的两站待配准点云的旋转矩阵。
步骤11，判断当前正确对应点个数是否等于两站待配准点云的对应点集包含的点对数，若是，执行步骤12；否则，将当前迭代次数加1后执行步骤4。
步骤12，配准两站待配准点云。
按照下式，对两站待配准点云中的每个点进行配准，完成两站待配准点云的配准：
Xi＝R1*(Pi-O1)+O1
Yi＝R2*(Qi-O2)+O2
其中，Xi表示对第一站待配准点云中第i个点进行配准后得到的点，R1表示第一站待配准点云的旋转矩阵，Pi表示第一站待配准点云中第i个点，Yi表示对第二站待配准点云中第i个点进行配准后得到的点，R2表示第二站待配准点云的旋转矩阵，Qi表示第二站待配准点云中第i个点。

Claims (10)

1.一种已知扫描点位置下的两点点云配准方法，其特征在于，包括如下步骤：
(1)获取待扫描物体的点云数据：
用地面三维激光扫描仪获取待扫描物体的点云数据；
(2)预处理点云数据：
(2a)从点云数据中任意选取一个站点的点云数据，将选取的点云数据的原点平移至所选取站点的扫描点位置，将平移后的点云数据作为第一站待配准点云；
(2b)从第一站待配准点云平移前的点云数据附近任意选取一个站点的点云数据，将选取的点云数据的原点平移至所选取站点的扫描点位置，将平移后的点云数据作为第二站待配准点云；
(2c)利用局部特征匹配方法，提取第一站待配准点云中与第二站待配准点云中相匹配的点对，组成两站待配准点云的对应点集，对应点集中包含正确对应点和错误对应点；
(3)利用迭代次数计算函数，计算迭代次数，确定迭代次数的最大值和最小值；
(4)判断当前迭代次数是否小于最小迭代次数或小于上次最大迭代次数，若是，则执行步骤(5)；否则，执行步骤(12)；
(5)计算两站待配准点云的第一次旋转矩阵：
(5a)从两站待配准点云的对应点集中，任意选取一组对应点作为样本A；
(5b)以第一站待配准点云的扫描点位置为球心O1，第一站待配准点云的扫描点位置到样本A中第一个点的欧式距离为半径，做球面S1；
(5c)以第二站待配准点云的扫描点位置为球心O2，第二站待配准点云的扫描点位置到样本A中第二个点的欧式距离为半径，做球面S2；
(5d)判断两个球面的半径之和是否大于两个球心之间的欧式距离，若是，则执行步骤(5e)；否则，执行步骤(5a)；
(5e)在两个球面的相交圆环的边缘上任取一点作为配准点；
(5f)按照下式，计算样本A中的第一个点、球心O1、配准点三点所在平面的法向量：
Figure FDA0002314160510000021
其中，n1表示样本A中的第一个点、球心O1、配准点三点所在平面的法向量，
Figure FDA0002314160510000022
表示向量符号，×表示向量叉乘操作，|·|表示向量的单位化操作，P1表示样本A中的第一个点，O1表示球面S1的球心，C表示配准点；
(5g)按照下式，计算球心O1和样本A中的第一个点确定的直线与球心O1和配准点确定的直线间的夹角：
Figure FDA0002314160510000023
其中，α1表示球心O1和样本A中的第一个点确定的直线与球心O1和配准点确定的直线间的夹角，arccos表示反余弦操作，r1表示球面S1的半径，d1表示样本A中第一个点到配准点的欧式距离；
(5h)将法向量n1和夹角α1代入绕任意轴旋转矩阵，得到第一站待配准点云的第一次旋转矩阵；
(5i)按照下式，计算样本A中的第二个点、球心O2、配准点三点所在平面的法向量：
Figure FDA0002314160510000024
其中，n2表示样本A中的第二个点、球心O2、配准点三点所在平面的法向量，
Figure FDA0002314160510000025
表示向量符号，×表示向量叉乘操作，|·|表示向量的单位化操作，Q1表示样本A中的第二个点，O2表示球面S2的球心；
(5j)通过下式，计算球心O2和样本A中第二个点确定的直线与球心O2和配准点确定的直线间的夹角：
Figure FDA0002314160510000026
其中，α2表示球心O2和样本A中第二个点确定的直线与球心O2和配准点确定的直线间的夹角，r2表示球面S2的半径，d2表示样本A中第一个点到配准点的欧式距离；
(5k)将法向量n2和夹角α2代入绕任意轴旋转矩阵，得到第二站待配准点云的第一次旋转矩阵；
(6)计算两站待配准点云的第二次旋转矩阵：
(6a)从两站待配准点云的对应点集中，任意选取除样本A之外的一组对应点，作为样本B；
(6b)以球心O1到匹配点的向量为旋转轴，将样本B中的第一个点旋转一个
Figure FDA0002314160510000031
角度，其中，arctan表示反正切操作，F表示点云数据的精度，d表示球心O1到球心O2的欧式距离；
(6c)以球心O2到匹配点的向量为旋转轴，将样本B中的第二个点旋转
Figure FDA0002314160510000032
角度，记录本次样本B中两个点的欧式距离及其对应的两个点的旋转角度；
(6d)判断样本B中的第二个点的旋转角度是否大于360°，若是，则将样本B中第二个点的旋转角度置为0后执行步骤(6e)；否则，执行步骤(6c)；
(6e)判断样本B中的第一个点的旋转角度是否大于360°，若是，则执行步骤(6f)；否则，执行步骤(6b)；
(6f)在所有旋转记录中，寻找样本B中两个点的欧式距离最小时对应的样本B中第一个点的旋转角度β1和第二个点的旋转角度β2；
(6g)将球心O1到匹配点的向量和旋转角度β1，代入绕任意轴旋转矩阵中，得到第一站待配准点云的第二次旋转矩阵；
(6h)将球心O2到匹配点的向量和旋转角度β2，代入绕任意轴旋转矩阵中，得到第二站待配准点云的第二次旋转矩阵；
(7)计算当前的两站待配准点云的旋转矩阵：
(7a)将第一站待配准点云的第二次旋转矩阵与第一次旋转矩阵相乘，得到第一站待配准点云的旋转矩阵；
(7b)将第二站待配准点云的第二次旋转矩阵与第一次旋转矩阵相乘，得到第二站待配准点云的旋转矩阵；
(8)计算当前正确对应点个数：
(8a)对两站待配准点云的对应点集进行配准，得到配准后的对应点集；
(8b)将配准后的对应点集中所有满足欧式距离阈值的对应点数，作为当前正确对应点的个数；
(9)利用代价函数，计算当前代价评估值；
(10)更新参数：
(10a)判断当前正确对应点个数是否等于上次正确对应点个数，若是，执行步骤(10b)；否则，执行步骤(10c)；
(10b)判断当前代价评估值是否小于上次代价评估值，若是，用当前代价评估值更新上次代价评估值后执行步骤(10e)；否则，执行步骤(10c)；
(10c)判断当前正确对应点个数是否大于上次正确对应点个数，若是，执行步骤(10d)；否则，执行步骤(11)；
(10d)用当前正确对应点个数更新上次正确对应点个数，利用迭代次数计算函数，计算当前最大迭代次数，用当前最大迭代次数更新上次最大迭代次数；
(10e)用当前两站待配准点云的旋转矩阵更新上次两站待配准点云的旋转矩阵；
(11)判断当前正确对应点个数是否等于两站待配准点云的对应点集包含的点对数，若是，执行步骤(12)；否则，将当前迭代次数加1后执行步骤(4)；
(12)配准两站待配准点云：
按照下式，对两站待配准点云中的每个点进行配准，完成两站待配准点云的配准：
Xi＝R1*(Pi-O1)+O1
Yi＝R2*(Qi-O2)+O2
其中，Xi表示对第一站待配准点云中第i个点进行配准后得到的点，R1表示第一站待配准点云的旋转矩阵，Pi表示第一站待配准点云中第i个点，Yi表示对第二站待配准点云中第i个点进行配准后得到的点，R2表示第二站待配准点云的旋转矩阵，Qi表示第二站待配准点云中第i个点。
2.根据权利要求1所述的一种已知扫描点位置下的两点点云配准方法，其特征在于，步骤(2a)、步骤(2b)、步骤(5b)、步骤(5c)中所述站点的扫描点位置是指，利用地面三维激光扫描仪配备的高精度全球定位系统GPS设备所确定的站点对应的地面三维激光扫描仪的扫描点位置。
3.根据权利要求1所述的一种已知扫描点位置下的两点点云配准方法，其特征在于，步骤(2c)中所述的局部特征匹配方法是指，使用局部特征描述子对两站待配准点云中的点做局部特征描述，用局部特征对两站待配准点云之间的点进行相似性匹配，将所有的匹配点的集合作为两站待配准点云的对应点集。
4.根据权利要求1所述的一种已知扫描点位置下的两点点云配准方法，其特征在于，步骤(3)、步骤(10d)中所述的迭代次数计算函数如下：
Figure FDA0002314160510000051
其中，K表示迭代次数，
Figure FDA0002314160510000052
表示向上取整操作，log表示以10为底的对数操作，η0表示两站待配准点云的对应点集中的所有点均为正确对应点的概率，其取值范围为[0.95,0.99]，
Figure FDA0002314160510000053
表示从m个不同元素中取出n个元素的组合数操作，N表示两站待配准点云的对应点集包含的对应点个数，Nin在计算最小迭代次数时表示预估计的正确对应点个数，在计算最大迭代次数时表示当前正确对应点个数。
5.根据权利要求1所述的一种已知扫描点位置下的两点点云配准方法，其特征在于，步骤(4)、步骤(10d)中所述的上次最大迭代次数是指，在第一次迭代时的最大迭代次数为106，其余迭代过程中的最大迭代次数为当前迭代之前更新的最大迭代次数。
6.根据权利要求1所述的一种已知扫描点位置下的两点点云配准方法，其特征在于，步骤(8a)中所述的对两站待配准点云的对应点集进行配准是按照下式实现的：
M＝{(R1*(Pi-O1)+O1,R2*(Qi-O2)+O2)}
其中，M表示两站待配准点云的对应点集进行配准后的对应点集，{}表示集合符号，(·,·)表示对应点对符号，R1表示当前第一站待配准点云的旋转矩阵，Pi表示在两站待配准点云的对应点集中属于第一站待配准点云的第i个点，R2表示当前第二站待配准点云的旋转矩阵，Qi表示在两站待配准点云的对应点集中属于第二站待配准点云的第i个点。
7.根据权利要求1所述的一种已知扫描点位置下的两点点云配准方法，其特征在于，步骤(9)中所述的代价函数如下：
Figure FDA0002314160510000061
其中，J表示两站待配准点云的完成一次配准的代价评估值，N表示两站待配准点云的对应点集包含的点对数，∑表示求和操作，Ri在配准后的对应点集中第i个对应点的欧式距离大于对应点的欧式距离阈值时，Ri表示对应点的欧式距离阈值，否则，Ri表示配准后的对应点集中第i个对应点的欧式距离；Ij在配准后的对应点集中第i个对应点的欧式距离大于对应点的欧式距离阈值时，Ij表示0，否则，Ij表示1。
8.根据权利要求1所述的一种已知扫描点位置下的两点点云配准方法，其特征在于，步骤(10a)、步骤(10c)、步骤(10d)中所述的上次正确对应点个数是指，在第一次迭代时的正确对应点个数为2，其余迭代过程中的正确对应点个数为当前迭代之前更新的正确对应点个数。
9.根据权利要求1所述的一种已知扫描点位置下的两点点云配准方法，其特征在于，步骤(10b)中所述的上次代价评估值是指，第一次迭代时的代价评估值为106，其余迭代过程中的代价评估值为当前迭代之前更新的代价评估值。
10.根据权利要求1所述的一种已知扫描点位置下的两点点云配准方法，其特征在于，步骤(10e)中所述的上次两站待配准点云的旋转矩阵是指，在第一次迭代时的两站待配准点云的旋转矩阵为步骤5计算出的两站待配准点云的第一次旋转矩阵，其余迭代过程中的两站待配准点云的旋转矩阵为当前迭代之前更新的两站待配准点云的旋转矩阵。
CN201711291269.XA 2017-12-08 2017-12-08 一种已知扫描点位置下的两点点云配准方法 Active CN108257163B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711291269.XA CN108257163B (zh) 2017-12-08 2017-12-08 一种已知扫描点位置下的两点点云配准方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711291269.XA CN108257163B (zh) 2017-12-08 2017-12-08 一种已知扫描点位置下的两点点云配准方法

Publications (2)

Publication Number Publication Date
CN108257163A CN108257163A (zh) 2018-07-06
CN108257163B true CN108257163B (zh) 2020-04-07

Family

ID=62721039

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711291269.XA Active CN108257163B (zh) 2017-12-08 2017-12-08 一种已知扫描点位置下的两点点云配准方法

Country Status (1)

Country Link
CN (1) CN108257163B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109712174B (zh) * 2018-12-25 2020-12-15 湖南大学 一种复杂异形曲面机器人三维测量的点云误配准滤除方法及系统
CN110415361B (zh) * 2019-07-26 2020-05-15 北京罗森博特科技有限公司 断裂物体拼接方法及装置
CN110342153B (zh) * 2019-07-29 2022-03-15 齐霄强 一种基于三维点云的垃圾桶识别抓取方法
CN111398936B (zh) * 2020-03-11 2021-04-06 山东大学 一种多路侧激光雷达点云配准装置及其使用方法
CN111402308B (zh) * 2020-03-17 2023-08-04 阿波罗智能技术(北京)有限公司 障碍物速度的确定方法、装置、设备和介质
CN112525161B (zh) * 2021-02-09 2021-05-14 南京景曜智能科技有限公司 一种旋转轴标定方法
CN115100258B (zh) * 2022-08-29 2023-02-07 杭州三坛医疗科技有限公司 一种髋关节图像配准方法、装置、设备以及存储介质

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102393183A (zh) * 2011-04-19 2012-03-28 程效军 基于控制网的海量点云快速配准方法
CN103955939A (zh) * 2014-05-16 2014-07-30 重庆理工大学 三维扫描系统中点云拼接用边界特征点配准方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10121223B2 (en) * 2015-03-02 2018-11-06 Aerial Sphere, Llc Post capture imagery processing and deployment systems

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102393183A (zh) * 2011-04-19 2012-03-28 程效军 基于控制网的海量点云快速配准方法
CN103955939A (zh) * 2014-05-16 2014-07-30 重庆理工大学 三维扫描系统中点云拼接用边界特征点配准方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
基于三维点云的机场通道面指标检测研究;才长帅等;《中国水运》;20170930;第17卷(第9期);全文 *
大规模城市三维重建中点云配准及平面提取研究;王瑞岩等;《信号处理》;20150630;第31卷(第6期);全文 *
应用摄像机位估计的点云初始配准;郭清达等;《光学精密工程》;20170630;第25卷(第6期);全文 *

Also Published As

Publication number Publication date
CN108257163A (zh) 2018-07-06

Similar Documents

Publication Publication Date Title
CN108257163B (zh) 一种已知扫描点位置下的两点点云配准方法
CN111145232A (zh) 一种基于特征信息变化度的三维点云自动配准方法
EP1319216A2 (en) Sar and flir image registration method
Yu et al. Robust robot pose estimation for challenging scenes with an RGB-D camera
CN108022262A (zh) 一种基于点的邻域重心向量特征的点云配准方法
CN109559340A (zh) 一种并行的三维点云数据自动化配准方法
CN114241018B (zh) 一种牙齿点云配准方法、系统及可读存储介质
CN109615645A (zh) 基于视觉的特征点提取方法
CN109965979A (zh) 一种无需标志点的稳健的神经导航自动注册方法
CN115100258B (zh) 一种髋关节图像配准方法、装置、设备以及存储介质
CN111028345B (zh) 一种港口场景下圆形管道的自动识别与对接方法
CN111243008A (zh) 一种用于高精度工件的圆弧数据拟合方法
CN117173227A (zh) 基于平面拟合的点云配准方法、装置和电子设备
CN113298838B (zh) 一种物体轮廓线提取方法及系统
Fretes et al. A review of existing evaluation methods for point clouds quality
CN115131433A (zh) 一种非合作目标位姿的处理方法、装置及电子设备
CN111639691B (zh) 一种基于特征匹配和贪婪搜索的图像数据采样方法
CN114545412A (zh) Isar图像序列等效雷达视线拟合的空间目标姿态估计方法
CN110197509B (zh) 一种基于彩色人工标识的相机位姿求解法
Huang et al. Pose determination of a cylinder using reprojection transformation
Tsai et al. Detection of vanishing points using Hough transform for single view 3D reconstruction
Wu et al. Point cloud registration algorithm based on the volume constraint
CN112183596A (zh) 结合局部网格约束和几何约束的直线段匹配方法与系统
CN111223139A (zh) 目标定位方法及终端设备
CN117076704B (zh) 一种图像质量分析评价方法及系统

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right

Effective date of registration: 20201201

Address after: 226300 Jianghai Yuanmeng Valley, No. 998, Century Avenue, high tech Zone, Nantong City, Jiangsu Province

Patentee after: JIANGSU HAOHAN INFORMATION TECHNOLOGY Co.,Ltd.

Address before: 710071, No. 2 Taibai South Road, Yanta District, Shaanxi, Xi'an

Patentee before: XIDIAN University

TR01 Transfer of patent right
PE01 Entry into force of the registration of the contract for pledge of patent right
PE01 Entry into force of the registration of the contract for pledge of patent right

Denomination of invention: A two-point point cloud registration method with known scanning point positions

Effective date of registration: 20211207

Granted publication date: 20200407

Pledgee: Nantong Jiangsu rural commercial bank Limited by Share Ltd.

Pledgor: JIANGSU HAOHAN INFORMATION TECHNOLOGY Co.,Ltd.

Registration number: Y2021980014300

PC01 Cancellation of the registration of the contract for pledge of patent right

Date of cancellation: 20221102

Granted publication date: 20200407

Pledgee: Nantong Jiangsu rural commercial bank Limited by Share Ltd.

Pledgor: JIANGSU HAOHAN INFORMATION TECHNOLOGY Co.,Ltd.

Registration number: Y2021980014300

PC01 Cancellation of the registration of the contract for pledge of patent right
PE01 Entry into force of the registration of the contract for pledge of patent right

Denomination of invention: A registration method of two-point point cloud with known scanning point position

Effective date of registration: 20221207

Granted publication date: 20200407

Pledgee: Jiangsu Nantong Rural Commercial Bank Co.,Ltd. Si'an Sub branch

Pledgor: JIANGSU HAOHAN INFORMATION TECHNOLOGY Co.,Ltd.

Registration number: Y2022980025338

PE01 Entry into force of the registration of the contract for pledge of patent right
PC01 Cancellation of the registration of the contract for pledge of patent right

Date of cancellation: 20231017

Granted publication date: 20200407

Pledgee: Jiangsu Nantong Rural Commercial Bank Co.,Ltd. Si'an Sub branch

Pledgor: JIANGSU HAOHAN INFORMATION TECHNOLOGY Co.,Ltd.

Registration number: Y2022980025338

PC01 Cancellation of the registration of the contract for pledge of patent right
PE01 Entry into force of the registration of the contract for pledge of patent right

Denomination of invention: A Two Point Point Cloud Registration Method with Known Scan Point Positions

Effective date of registration: 20231020

Granted publication date: 20200407

Pledgee: Jiangsu Nantong Rural Commercial Bank Co.,Ltd. Si'an Sub branch

Pledgor: JIANGSU HAOHAN INFORMATION TECHNOLOGY Co.,Ltd.

Registration number: Y2023980061751

PE01 Entry into force of the registration of the contract for pledge of patent right