基于点云配准的曲面三维差异检测*

罗锦鸿，陈新度，吴 磊

(广东工业大学机电工程学院，广州 510006)

0 引言

在视觉引导机器人的智能化自动化加工任务中，工件的形状误差分析和三维差异数据提取是一个关键研究问题。本文针对陶瓷素胚工件的机器人自动化打磨任务，需要对陶瓷素胚工件的部分表面进行差异比对以识别定位加工区域。为了获取待加工工件的点云数据，如图1a所示，利用了三维扫描仪来扫描待加工工件的表面，如图1b所示，有标准模型曲面点云(白色)和通过扫描得到的待检测曲面点云(蓝色)。

近年来，研究人员通过对三维视觉方法获取的实测曲面的三维数据进行分析，来实现曲面差异检测，而使用点云配准的方法来让标准模型数据与实测数据对齐用来实现差异检测逐渐成为趋势，比如李宇萌等[1]采用了ICP精确配准的方法进行模型与实测数据对齐，在102 s内将2万左右点云的精确对齐，实现了转子表面的缺陷检测。马维康等[2]通过手持式智能激光扫描仪对试验件不同工序段进行扫描测量,将扫描的三维点云数据与设计的理论数模进行对比分析,得出型面精度的偏差情况,以此来检验试验件是否合格。

但是前人的方法大部分是针对全局配准或者位姿容易确定的情形，可以使用ICP算法直接进行数据对齐。在面向待配准点云数据是模型的局部实测数据的情形时，这些方法并不适用的。本文针对这种局部曲面配准情况，测试了两种基于局部特征的配准方法，实现了快速粗配准。在精确配准阶段，由于ICP配准算法在面对大型点云数据时耗时长，满足不了本文的实时性的要求。本文实现了一种由CUDA加速的ICP算法，实现了大型曲面点云的快速对齐。而且在曲面对齐后，通过本文设计的引入Octree预检测的差异提取方法可以快速地提取曲面的三维差异。

通过本文设计的实验证明，在面对要求实时性比对检测目标曲面与标准模型曲面的任务时，本文的方法具有精度高、可靠、快速的特点。

(a) 得到工件单帧扫描点云 (b) 扫描点云(蓝)和模型点云(白)

1 基于Hough投票的粗配准

点云配准是指利用两幅点云的匹配点，利用匹配关系来估算两幅点云的旋转平移关系，从而使两幅点云正确地重合在一起。首先要解决的是精确配准的初值问题，即解决粗配准问题。

点云的粗配准方法有很多，研究人员针对工件点云的粗配准问题上进行过很多研究，比如陈巍等[3]、周欣康[4]和项辉宇等[5]分别在校正叶片工件的装夹误差研究中、航空叶片型面检测中以及板金属冲压零件表面模型与实测点云配准研究中，都使用了PCA(主成分分析)的全局配准算法，对实测点云与模型点云进行粗对齐。但是由于PCA算法要求实测点云形状与模型相近，这就需要在获取数据时，完整地扫描工件。

本文要对工件的待加工处的部分扫描点云与整体模型点云进行配准，显然PCA全局配准方法不合适，文章考虑了基于局部特征的粗配准算法，基于局部特征的配准方法主要有SAC-IA[6](Sample Consensus Initial Alignment)、RANSAC[7]和三维霍夫投票[8](Hough voting)。本文分别对三维霍夫投票法和SAC-IA两种算法进行了测试。

本文总算法流程如图2所示。有模型点云M={m1,…,mi}和通过实测扫描点云S={s1,…,sj}，分别对M和S进行滤波得到稀疏的下采样点云M′={m1′,…,mk′ }和S′={s1′,…,sl′ }。计算M′和S′中的每个点的局部特征描述子F(M′)={fm1′,…,fmk′ }，F(S′)={fs1′,…,fsl′}和局部参考坐标系LRF(M′)={rm1′,…,rmk′ }，FRF(S′)={rs1′,…,rsl′}。使用Hough投票法或采样一致性初始对齐(SAC-IA)算法得到粗对齐矩阵T1，使用T1变换S中的每个点，得到粗对齐点云S_initial，然后用CUDA加速的ICP算法完成S_initial与M的精确对齐，得到最终配准完成的扫描点云S_dest，在提取差异点云步骤中，遍历S_dest点云中的所有点，当查询点与其在M中最近的点之间的距离大于设定阈值时，将该查询点看作是差异点，加入到差异点集合S_d，遍历完成后，最终的集合S_d是差异点云。

图2 三维差异检测算法流程图

1.1 体素栅格采样

体素栅格采样是点云数据精简的一种非常实用的方法，它首先创建一个三维的体素栅格，输入点云将包含在栅格中的三维立方体内，然后用三维立方体中所有的点计算平均值得到重心点，用来近似表示体素中其他点，统计重心点集合将得到下采样点云。

实验表明，对于以本文中的工件为代表的典型平滑曲面点云，由于点云密集以及曲面光滑，三维局部特征十分相似，如果直接使用原始数据进行粗配准会导致计算量非常大而且失败。于是本文采用了栅格滤波算法对模型点云M和扫描点云S进行滤波。数据量在210 000左右的模型点云和数据量在34 000左右扫描点云下栅格滤波后的数据量为8 200和1 300左右。输入点云栅格下采样后的点云如图3所示。点云滤波后，在空间分布上具有规律性，在计算三维局部特征时，通过采用栅格滤波后的点作为关键点将有效地去冗余点云，以及减少了计算特征描述子要考虑的近邻点数，加快了粗配准的速度以及成功率。

图3 栅格采样后的模型点云和扫描点云的状态

1.2 3D-Hough投票

3D-Hough[8]中算法具体步骤如下：

(1)输入模型的下采样曲面点云M′={m1′,…,mk′ }和下采样扫描点云S′={s1′,…,sl′}。

(2)计算M′和S′中的每个点的局部三维特征描述子F(M′)={fm1′,…,fmk′ }，F(S′)={fs1′,…,fsl′}和局部参考坐标系LRF(M′)={rm1′,…,rmk′ }，FRF(S′)={rs1′,…,rsl′}。

(3)在本实验中，将重心点作为Hough投票的参考点。于是先计算M′的重心点c_m，再计算M′中每个点到重心点c_m的向量vmk=c_m-mk′，得到向量集合V_M={vm1,…,vmk}。

(4)通过局部三维特征描述子欧式距离作为匹配度来计算大小为n的匹配集。

(5)遍历匹配集合，匹配上的对应点mn和sn对应的两个参考坐标系：rmk′到rsl′之间有变换矩阵Tn，于是有Tn*vmn=vsn，即估计的扫描点云的重心点为csn=Tn*vmn+sl′。

(6)通过步骤(5)的计算就会得到一个重心点集合C_S={cs1,…,csn}。集合C_S即是三维Hough空间中的投票。

(7)为集合C_S设置区间步长划分投票空间，当投票数超过一定时，选取投票数最大的区间，该峰值区间中的所有投票重心求平均值即为扫描点云的参考点(扫描点云重心)。

(8)通过对峰值区间筛选的对应关系，可以计算出模型点云与扫描点云之间的旋转平移变换矩阵。

本文分别对6组在不同的位置拍摄的单帧扫描数据进行了测试，粗配准后状态如图4所示。

图4 粗配准后的模型点云和扫描点云的状态

如表1所示，在同等效果的情况下，使用基于局部特征的的3D-Hough投票算法更快，平均运行时间在0.3 s以内。这是由于SAC-IA方法是在匹配点集中不断抽取必要数目的匹配点去估算变换矩阵，所以要进行多次抽取来选取最小误差的估算变换参数，而3D-Hough投票算法，通过找出峰值的方式来剔除错误匹配点对，一次估计就可以求解出变换矩阵，所以速度更快。

表1 测试SAC-IA和Hough投票法

续表

2 CUDA加速的点云精确配准

2.1 ICP

在点云精确对齐领域，Besl P J等[9]提出了经典的ICP(Iterative Closest Point)配准算法，通过迭代最近点的方法实现了两个点云的精确配准。通过上面的粗配准步骤后，再通过ICP算法进行精确配准可以使得两个点云精确对齐。

ICP算法基本原理是：分别在待配准的目标点云P和源点云Q中，按照一定的约束条件比如邻近关系，找到对应点对最邻近点对，构建基于距离的误差函数式(1)。不断迭代计算出最优匹配参数R和t，使得误差函数最小，直至达到满意的误差要求，一般使用SVD或者非线性优化方法求解。

(1)

2.2 CUDA-ICP

直到最近，几乎所有的视觉实现都使用了GPU或FPGA等方法进行加速[10-11]。由于本文模型以及点云点数较多，而且面对曲率变化小，特征不明显的点云下采样后在精确配准时会损失精度甚至可能配准失败。本文利用了GPU的超高性能计算能力，以及CUDA架构的简易开发性，编写了CUDA程序来加速ICP算法，以实现实时配准。

在实验中，首先利用了K-D tree(k-dimensional树)数据结构实现了快速的最近点搜索，对一个模型点云构建一次K-D tree可以用于多次迭代中的最近邻搜索。其次，通过在对应估计的步骤中利用多线程来同时计算搜索近邻。再者，在刚体变换扫描点云S步骤中，开启多线程同时对每个点进行矩阵运算。最后，在需要规约计算的地方使用CUDA提供的Thrust库的接口thrust::reduce实现快速规约运算，这些方法都极大地加速了配准过程。

图5是本文中的CUDA-ICP算法计算和优化的流程图以及详细步骤。

图5 CUDA-ICP算法流程图

(1)初始化：输入M={m1,…,mm}和经过粗对齐后的扫描点云S={s1,…,sn}。输入收敛阈值τ和连续收敛次数阈值t以及最大近邻阈值d；

(2)对模型点云M构建K-D tree，用于最近邻查询；

(3)开启大于n数量的线程并行：通过已经构建好的K-D tree实现最近邻搜索，寻找到S中每个点sn在M中的最近点mm，当对应点之间的距离小于d时加入对应点集合。最终得到对应点集合M′和S′；

(4)使用thrust库进行规约：计算集合M′和S′之间的误差et；

(5)计算迭代误差前后变化Δe=|et-et-1|，当Δe大于设定收敛阈值τ，则继续执行。当Δe小于τ时则认为本次迭代收敛，若连续收敛的次数大于设定的次数阈值t，就结束计算，本实验中t=5；

(6)使用thrust库进行规约：计算M′和S′的重心点M_c和S_c；

(7)开启大于n数量的线程并行：计算去重心点集合M_2和S_2；

(8)开启大于n数量的线程并行：计算去重心点集合M_2和S_2之间的协方差矩阵集合W；

(9)使用thrust库进行规约：对集合W进行规约得到协方差矩阵w。使用SVD分解w得到本次迭代的变换矩阵T；

(10)开启大于n数量的线程并行：使用变换矩阵T来旋转平移扫描点集合S中的每个点，变换后的扫描点集合作为S跳到步骤(3)；

(11)结束计算。

表2 测试CUDA-ICP

续表

本文对6组实验数据进行了测试，如表2结果所示，本文改进的CUDA-ICP算法将原本运行时间在7 s左右的基于K-D tree的标准ICP配准算法加速到了0.5 s左右。本文配准精度控制在RMS误差小于0.5 mm，每次迭代时间控制在配准耗时控制在1.5 ms以内，满足本文中的三维差异定位任务的准确性和实时性的要求。

3 三维差异提取

本文认为扫描点云S_dest={sd1,…,sdj}的当前查询点与其在模型点云M={m1,…,mi}上的最近点的距离大于某个设定的距离阈值时，就认为该查询点是差异点。这样遍历算法的时间复杂度为O(i*j)。通过八叉树(Octree)数据结构实现了三维差异的快速提取。

如图6所示，八叉树结构通过对三维空间的集合实体进行体元剖分，每个体元具有相同的时间和空间复杂度，通过循环递归的划分方法对大小为2n×2n×2n的三维空间的几何对象进行剖分，从而构成一个具有根节点的方向图。通过递归地比较模型点云构建的Octree的树结构和扫描点云构建的Octree结构，可以对比得到不同的差异节点，通过记录新增节点和缺少节点并记录其对应点云，从而实现检测差异。

本文通过Octree节点预先得到存在差异点云的体素，再通过有差异点的体素及其周围体素中的点云对标准模型点云进行近邻搜索，得到最终的差异点云。Octree数据结构示意图如图6所示。

图6 Octree数据结构示意图

(1)设置距离阈值d_min，则Octree的最小子节点步长度为V=d_min，输入模型点云M构建Octree，输入扫描点云S_dest构建Octree；

(2)通过比对两个Octree得到新增的Octree节点，提取新增节点及其邻近节点处的点云。可以快速得到初步差异点云S_d_initia={s1″,…,sn″}，如图7a所示；

(3)遍历S_d_initia，查询点sq″与在模型点云M上的最近点mq之间的距离为dq=distance(mq,sq″)，当dq>d_min时认为当前查询点sq″是差异点，将该点加入到差异点集合S_d中。如图7b～图7d所示；

(4)输出差异点云S_d，如图7d所示。

通过上述计算，可以正确获得需要提取的差异点云，差异提取步骤在CPU端中运行时间为0.15 s。使用CUDA加速后耗时在1 ms以内。

(a) 初步差异点云S_d_initia (b) 差异点云S_d距离示意图

(c) 最终差异点云S_d(红) (d) 点云S_d细节图7 差异点云提取

4 结束语

本文针对初始位姿不确定的平滑点云数据的差异提取问题进行了研究，利用曲面配准，以及近邻搜索等方法实现了三维差异检测。而且针对算法运行速度问题，测试了SAC-IA与3D-Hough投票两种点云的粗配准方法，使用了Octree数据结构实现了差异的快速定位以及K-D tree实现快速近邻搜索，并利用了CUDA来加速程序。最终在本课题马桶陶瓷素胚的曲面形状差异检测上，将点对误差控制在0.5 mm以内，运行时间在控制在1 s以内，实现了快速精确的工件三维差异检测与定位。