CN110874849B - 一种基于局部变换一致的非刚性点集配准方法 - Google Patents

Abstract

本发明提供一种基于局部变换一致的非刚性点集配准方法，用以解决变形较大的非刚性物体，在配准过程中关节点周围出现的配准结果不理想的问题。为解决该问题，基于非刚性变换相邻节点之间的变换是比较一致的思想，在配准过程中通过计算每一个点的空间变换分别与其K个相邻点的空间变换的差异，令差异最小化实现对非刚性变换的局部约束，保证点集在配准过程中保持着局部结构，提高了关节点周围配准的效果。本发明在一些数据集上取得了较好的配准效果。

Description

一种基于局部变换一致的非刚性点集配准方法

技术领域

本发明属于三维视觉处理技术领域，尤其涉及点云变形较大情况下导致关节点周围配准结果不够理想时的点云配准，具体为一种基于局部变换一致的非刚性点集配准方法。

背景技术

在计算机信息技术和三维传感技术飞速发展的时代，三维重建技术在很多领域得到了广泛的应用，比如医学图像、3D打印、3D电影和动画。该技术是将从扫描设备中获得的物体点云信息建立成物体三维模型。由于扫描设备的范围有限，扫描一次是不能把物体的所有三维点云数据得到，所以需要多次扫描才能获取一个完整的模型。三维点云配准技术就是通过计算不同扫描视觉下的点云的变换关系，将点云转换到一个统一的坐标系下。根据物体类型的不同，点云配准可以分为刚性配准和非刚性配准，刚性配准针对的是刚体，刚体的变换形式是全局旋转和平移，称为刚性变换。非刚性配准针对的是非刚体，非刚体是随着时间和任意运动而变形的物体，对应的变换形式较复杂，称为非刚性变换。作为三维重建中的核心部分，三维点云配准技术一直是人们关注的重点内容。可以说点云配准技术的进步能够带动三维重建技术的发展。但点云配准技术不仅是科研学者的研究热点，同时也是研究难点。点云配准的目标是找出点云之间精确的对应关系，计算出点云之间最佳的变换矩阵，然而其实现过程中受到了很多限制。比如由扫描设备和扫描物体所处环境的噪声对数据本身产生的影响，不同扫描视角下物体的变形过大等，这些影响无疑给点云配准技术带来了巨大的挑战。

近些年来，国内外学者已经对点云配准技术进行了大量的研究，配准效果得到不断提高。点集配准方法可以大致分为：基于距离的方法和基于概率的方法。基于距离的方法分为两步，第一步计算两个点集之间的距离并找出对应关系，第二步是使具有确定的对应关系的两个点集之间的距离最小化。其中迭代最近点(Iterative Closest Point，简称ICP)是两个点集刚性配准领域中众所周知的方法。然而，ICP方法对初始条件比较敏感，容易陷入局部最小值。一致性点漂移(Coherent Point Drift，简称CPD)是基于概率的点集配准的流行方法，该方法将两个点集的配准公式化为概率密度估计问题。一个点集表示高斯混合模型(Gaussian Mixture Model，简称GMM)质心，通过将该点集一致地移动拟合到另一个点集。该方法在噪声和异常值情况下，取得了较好的配准效果。但其仍存在一些缺点，比如，该方法中GMM的分量的隶属概率相等，事实上，从GMM的一个分量分得到的样本点越多，该组件的先验概率应该越大。该方法只考虑了点集之间的全局结构信息，然而，当物体发生大的变形时，点集的局部结构一般是很强也很稳定的。后来基于这些问题，提出了许多扩展CPD算法。比如，重新分配隶属概率：提出了形状上下文特征，通过计算两个点集形状上下文的相似度来分配GMM的分量的隶属概率，从而提高了点集配准的性能。提取点集特征：通过合并点集的颜色信息，提出了一种新的概率方法进行点集配准。利用与每个点相关的颜色特征信息，合理地确定了点集之间点与点的对应概率，提高了点集配准的准确性和鲁棒性。以上算法由于没有考虑到点集的局部结构信息，对变形较大的物体进行配准时，在关节点周围容易产生错误的对应点对，从而导致错误的配准结果。构建局部结构信息：CPD算法中，仅考虑了两个点集的全局结构信息，为保持点集的局部结构信息，提出了对点集本身的局部结构进行约束，即一个点可以由其相邻点加权线性表示，这种表示关系经过非刚性变换后依然可以保持。该算法虽然考虑了点集的局部结构信息，然而当点集稀疏或者不均匀时，该局部结构不够稳定，会在点集部分区域产生一些不匹配的对应点对，从而导致错误的配准结果。

发明内容

本发明提供了一种基于局部变换一致的非刚性点集配准方法。针对点云变形较大情况下，点云关节点周围出现配准效果不理想的问题，采用对相邻点空间变换进行局部约束的方法，在配准过程中使点云保持着局部结构，从而提高了点云关节点周围配准的效果。

发明目的：本发明所要解决的是点云变形较大情况下导致关节点周围配准结果不够理想的问题，提出一种基于局部变换一致的非刚性点集配准方法。

技术方案：

一种基于局部变换一致的非刚性点集配准方法，其方法为：

利用三维点云获取设备，获取两个点云，分别记为源点云和目标点云，分别对这两个点云进行降采样，得到它们各自的关键点集，分别记为源关键点集S和目标关键点集O，迭代计算这两个关键点集之间的对应关系，源关键点集S的邻域索引矩阵和空间变换矩阵。在迭代过程中，为保持点集的局部结构，基于相邻点的空间变换一致思想，对非刚性变换进行局部约束，最终得到最优的空间变换矩阵，经空间变换矩阵变换后的源关键点集S'，目标关键点集O作为配准结果并输出。补充地说，配准的目的就是将两个关键点集统一到一个坐标系下，配准过程中点集O是没有变化的，只是将点集S不断进行变换，使点集S逐渐与点集O拟合。将点集O和变换后的点集S'一起输出才可以定性地展示配准效果。

进一步说，本发明所述的一种基于局部变换一致的非刚性点集配准方法，是通过计算机，并按如下步骤进行：

步骤1，向计算机输入由三维点云获取设备获取的源点云和目标点云，获取这两个点云的关键点集：源关键点集S和目标关键点集O。

步骤2，由步骤1中获取的两个关键点集，转换获取关键点集之间的对应关系。

步骤3，建立源关键点集S的邻域索引矩阵。

步骤4，获得空间变换矩阵，并以此获得变换后的源关键点集S'。

步骤5，设定迭代的参数值，迭代执行步骤2到步骤4，当达到参数值时停止迭代，将此时的源关键点集S'和目标关键点集O作为配准结果并输出。

进一步说，步骤1，获取两个不同视角的点云，对这些点云进行关键点提取，获取各自的关键点集，即源关键点集S和目标关键点集O。此处的点云是指一些只包含三维坐标信息的点的集合，关键点提取是指对点云进行简化，获得保留点云基本形状特征和空间结构信息的一些点，这些点也就是关键点。

步骤2，对步骤1中获取的两个关键点集，计算源关键点集S中的每个点在目标关键点集O中的对应点，获取这两个关键点集之间的对应关系。

步骤3，对源关键点集S中的每一个点进行K近邻搜索，即对关键点集S中每一个点，在整个源关键点集S中找到与其欧几里得距离(Euclidean distance，简称欧式距离)最近的K个点，这K个点便称为它的相邻点，并根据这些相邻点在源关键点集S中的位置信息，建立源关键点集S的邻域索引矩阵。此处的欧式距离是指两个三维点之间的实际距离，K取值不小于2，优选地取值为3。

步骤4，利用相邻点的空间变换一致，对非刚性变换进行局部约束，即对源关键点集S中的每一个点的变换都进行局部约束，计算每一个点的空间变换分别与其K个相邻点的空间变换的差异，在保证差异较小的情况下获得较优的空间变换矩阵，以此获得变换后的源关键点集S'。

步骤5，设定一个最大迭代次数和一个参数阈值，迭代执行步骤2到步骤4，参数值超过阈值或者迭代次数超过最大迭代次数，停止迭代，此时经空间变换矩阵变换后的源关键点集S'及目标关键点集O作为配准结果并输出。

优选的技术方案是：本发明所述的一种基于局部变换一致的非刚性点集配准方法的详细步骤如下：

步骤1，点云获取与提取点云关键点：获取仅包含三维坐标信息的两个点云，即源点云和目标点云，分别对这两个点云采用降采样的方式提取关键点，得到源关键点集S和目标关键点集O。所述的降采样是指对点云进行稀疏化处理，源关键点集S和目标关键点集O分别对应的是源点云和目标点云进行降采样后的数量不大于1万个点的两个稀疏点云，这两个点集分别保留了源点云和目标点云的形状特征及空间结构信息。

步骤2，计算源关键点集S与目标关键点集O之间的对应关系：在步骤1获取关键点集的基础上，对源关键点集S中的每一个点，计算它与目标关键点集O中的每一个点的欧式距离，依据对应的欧式距离确定两个关键点集中每个点与点之间的对应概率。

步骤3，建立源关键点集S的邻域索引矩阵：对步骤1获取的源关键点集S中的每一个点，在整个源关键点集S范围内计算与其欧式距离最近的K个相邻点，并记录这些相邻点在源关键点集S中的位置索引，根据已记录的索引建立源关键点集S的邻域索引矩阵。点集中的所有点都有标号的，比如第一个点，第二个点，......索引指的就是这样一个标号，用来标识相邻点属于源关键点集中的第几个点。

步骤4，计算源关键点集S的空间变换矩阵：根据源关键点集S及其邻域索引矩阵，定义源关键点集S中每个点的空间变换与其K个相邻点的空间变换的差异程度的函数式，即约束项。该约束项由源关键点集S中每个点的空间变换分别与其K个相邻点的空间变换的差值之和组成，用来约束相邻点之间的变换保持一致。数据项表示的是源关键点集S和目标关键点集O中每一对对应点之间的距离。约束项与数据项组成目标函数式的作用是最小化目标函数式，获取此时对应的空间变换矩阵，并将此空间变换矩阵作用到源关键点集S，获取新的源关键点集S'。此处的空间变换指一个点集位置发生变化后所经历的变换形式，该变换形式用一个矩阵表示，相应地，源关键点集S中的点和其相邻点的空间变换指的是一个点位置发生变化时所经历的变换形式，该变换形式用一个向量表示。

步骤5，迭代执行步骤2到步骤4：设置一个最大迭代次数和参数σ²的阈值，每循环执行一次步骤2到步骤4，便重新计算一下参数σ²的值，迭代次数加一。

当参数值大于阈值或者迭代次数达到最大迭代次数，停止迭代，此时将空间变换矩阵作用在源关键点集S上，得到新的源关键点集S'，将其和目标关键点集O，作为配准结果并输出。

当参数值不大于阈值、迭代次数未达到最大迭代次数时，说明对应点之间的距离并没有接近于0，没有实现将两个点集统一到一个坐标系下，即没有完成两个关键点集的配准，返回步骤2。

更进一步说，步骤1的实现过程具体如下：

步骤1.1，点云获取：获取两个不同视角的点云，利用三维几何处理软件，如MeshLab，将这两个点云保存为只包含三维坐标信息的点云，得到源点云和目标点云。

步骤1.2，关键点提取：根据步骤1.1得到的源点云和目标点云，在源点云数据范围内建立一个网格，所述网格由三维体素组成(三维体素的大小是固定的)，通过将源点云包围在三维体素网格中，每个体素中包含源点云的多个三维点，为简化点的数量，则在每个三维体素中选一个点，这样在该三维体素内所有点就可以用这一个点表示，完成了对点云数据的简化。对网格中所有体素处理后得到的点云即为源点云对应的源关键点集S。同理，目标点云对应一个目标关键点集O。

更进一步说，步骤2具体如下：

步骤2.1，在步骤1.2得到的源关键点集S和目标关键点集O的基础上，计算源关键点集S中的每个点与目标关键点集O中每个点之间的欧几里得距离：假设y_m为源关键点集S中的一个点，x_n为目标关键点集O中的一个点，则两点之间的距离公式如下：

d(x_n,y_m)＝||x_n-y_m||²     (2-1)

步骤2.2，计算高斯混合模型(GMM)概率密度函数:将源关键点集S_M×D＝(y₁,…,y_M)^T表示为高斯混合模型的质心，将目标关键点集O_N×D＝(x₁,…,x_N)^T作为由高斯混合模型生成的数据点，在步骤2.1得到的点与点之间的欧几里得距离基础上，建立高斯混合模型概率密度函数公式如下：

其中，M和N分别表示源关键点集S和目标关键点集O的点的个数，m表示高斯混合模型的第m个高斯分量，n表示目标关键点集O中点的下标，D表示两个关键点集的维数，取值为3，exp表示以自然常数e为底的指数函数。σ²表示每个高斯分量的协方差，初始值为：

假设GMM的所有分量都是独立同分布的，则联合高斯混合模型概率密度函数公式如下：

其中，M表示的是源关键点集S的点的个数，N表示的是目标关键点集O的点的个数，P(m)＝1/M表示每个高斯分量的隶属概率。

步骤2.3，计算源关键点集S中的每个点与目标关键点集O中每个点之间的对应关系：基于高斯混合模型的点集配准方法，当源关键点集S和目标关键点集O对齐时，对于一个已知点x_n，其与点y_m之间的对应关系利用最大化高斯混合模型的后验概率获得。在步骤2.2的基础上，计算高斯混合模型的后验概率，获得两个点集之间的对应概率。高斯混合模型的后验概率公式如下：

P(m|x_n)＝P(m)p(x_n|m)/p(x_n)    (2-5)

更进一步说，步骤3具体如下：

源关键点集S邻域索引矩阵的建立：对步骤2中的源关键点集S,采用被广泛应用的kd-tree最近邻搜索算法，将源关键点集S中的每一个点作为查询点，检索在kd-tree树中与查询点距离最近的K个相邻点，计算这K个相邻点在源关键点集S中的位置索引，根据索引建立源关键点集S的邻域索引矩阵Idy＝[idy₁,…,idy_i,…,idy_M]^T。

更进一步说，步骤4具体如下：

步骤4.1，基于GMM点集配准的目标函数的获取：基于GMM点集配准将点集配准问题转化为概率密度估计问题，通过重新参数化GMM质心位置，使质心逐渐拟合数据点。由步骤2.2获得的联合GMM概率密度函数，则其负对数似然函数公式如下：

一般分别对相应的参数求导，令导数为零，可以得到GMM新的参数。然而，负对数似然函数，对数里面还有求和，实际上没有办法通过求导的方法来求负对数似然函数的最小值。于是，采用EM算法进行GMM参数估计问题。EM算法分为两步，第一步寻找目标函数，通过计算E的改变量，可以得到目标函数的公式如下：

其中，θ表示一组变换参数，p'(m|x_n)表示点x_n和点y_m的初始对应概率，T(y_m,θ)表示应用于源关键点集S的变换函数，

刚性变换与非刚性变换通过确定上式的T函数来区分，刚性变换定义为T(R,t)＝RS+t，其中R表示旋转矩阵，t表示位移向量，在被大家公用的CPD算法中将非刚性变换定义为一个基于高斯径向基函数的位移函数，公式如下：

T(S,W)＝S+GW      (4-3)

其中，G是一个M×M的高斯核矩阵，其元素为，

W是一个M×D的高斯核权重矩阵，通过规范权重矩阵W以强制运动一致性，使点集在配准期间保持整体空间的连通性。全局约束项表示如下：

E_g(W)＝Tr(W^TGW)   (4-4)

将全局约束项添加到目标函数中，则目标函数重新表示如下：

其中，G(m,·)对应矩阵G的第m行。

步骤4.2，构建局部约束项：在非刚性配准的背景下，相邻点的空间变换是一致的。基于此条件，建立了局部约束，旨在配准期间能够保持点集的局部结构。基于步骤3获得的源关键点集S的邻域索引矩阵，计算源关键点集S中每个点的空间变换分别与其K个相邻点空间变换的差值，并计算所有差值的和，使和的值尽可能的小来实现局部约束，源关键点集S空间变换的局部约束项公式如下：

其中，K取值为3，G(m,·)对应矩阵G的第m行，Idy(m,k)表示源关键点集S中第m个点的第k个相邻点的索引。

步骤4.3，建立最终的目标函数：基于步骤4.1获取的目标函数，步骤4.2得到的局部约束项，将局部约束项添加到步骤4.1获得的目标函数中，本发明的目标函数公式如下：

Q(W,σ²)＝Q_d(W)+λE_l(W)        (4-7)

其中，λ表示局部约束的权重系数，取值为50000。

然后执行EM算法的第二步，最小化目标函数Q，分别对W和σ²求导，

对W求导：

对σ²求导：

其中，P是一个M×N的矩阵，其元素值表示源关键点集S中的每个点与目标关键点集O中每个点之间的对应概率。P1是P与值全为1的列向量的乘积，d(P1)表示由向量P1组成的对角矩阵。

令式(4-8)与式(4-9)的和为0，得到W的值，令式(4-10)为0，得到σ²的值，通过S+GW变换源关键点集S得到变换后的源关键点集S'，将变换后的源关键集S'作为下次迭代的源关键点集S。

更进一步说，步骤5具体如下：

设定一个最大迭代次数和参数σ²的阈值，循环执行步骤2到步骤4，每循环执行一次，先对源关键点集S进行变换，即S'＝S+GW，然后计算此时的σ²值，迭代次数相应加一，判断迭代次数是否超过最大迭代次数，比较σ²的值是否大于设定的阈值，当迭代次数超过最大迭代次数或者σ²的值大于给定的阈值时，迭代终止。若迭代未终止，将变换后的源关键点集S'作为下次迭代的源关键点集S。最终得到的源关键点集S'和目标关键点集O作为配准结果并输出。

有益的技术效果

本发明所提供的非刚性点集配准方法，是用来解决点云变形较大情况下关节点周围配准效果不佳的问题。采用的是对每个点的空间变换进行局部约束的方法。步骤1，获取两个不同视角的点云，对这两个点云进行关键点提取，获取这两个点云的关键点集，即源关键点集S和目标关键点集O。步骤2，对步骤1中获取的两个关键点集，计算源关键点集S中的每个点在目标关键点集O中的对应点，获取关键点集之间的对应关系。步骤3，对源关键点集S中的每一个点进行K近邻搜索，确定每一个点的K个相邻点，并根据这些相邻点在源关键点集S中的位置信息，建立源关键点集S的邻域索引矩阵。步骤4，利用相邻点的空间变换一致，对非刚性变换进行局部约束，即对源关键点集S中的每一个点的变换都进行局部约束，计算每一个点的空间变换分别与其K个相邻点的空间变换的差异，在保证差异较小的情况下获得较优的空间变换矩阵，以此获得变换后的源关键点集S'。步骤5，设定一个最大迭代次数和一个参数阈值，迭代执行步骤2到步骤4，参数值超过阈值或者迭代次数超过最大迭代次数，停止迭代，此时得到的源关键点集S'和目标关键点集O作为配准结果并输出。

本发明是一种基于局部变换一致的非刚性点集配准方法，针对点云变形较大情况下关节点周围配准效果不佳的问题，设计有效的解决方法。与现有技术不同的是，本发明没有改变GMM的分量的隶属概率，也没有采用丰富的点集特征，仅仅只使用了点的位置信息。没有对点集本身的局部结构进行约束，而是对点集的空间变换进行局部约束。本发明主要建立在相邻点空间变换通常是一致这一思想基础上，利用点集中点的变换与其相邻点的变换之间的关系，计算点的空间变换与其相邻点空间变换的差值，再对这些差值求和，令和尽可能小达到对点集的空间变换的局部约束，使点集在配准过程中保持着局部结构，有效地提高了点云变形较大情况下关节点周围的配准效果。

附图说明

图1为本发明方法的基本流程图。

图2为实验配准误差图，对3个样例分别进行配准，CPD算法和本发明的配准误差对比图。

图3为实验配准效果图，对源关键点集S及目标关键点集O进行配准，CPD算法和本发明的配准效果对比图。

图4为实验配准效果图，对源关键点集S及目标关键点集O进行配准，CPD算法和本发明的配准效果对比图。

具体实施方法

现结合附图详细说明本发明的特点。

参见图1，一种基于局部变换一致的非刚性点集配准方法，利用三维点云获取设备，获取两个点云，分别记为源点云和目标点云，分别对这两个点云进行降采样，得到它们各自的关键点集，分别记为源关键点集S和目标关键点集O，迭代计算这两个关键点集之间的对应关系，源关键点集S的邻域索引矩阵和空间变换矩阵。在迭代过程中，为保持点集的局部结构，基于相邻点的空间变换一致思想，对非刚性变换进行局部约束，最终得到最优的空间变换矩阵，经空间变换矩阵变换后的源关键点集S'，目标关键点集O作为配准结果并输出。

本发明所述的一种基于局部变换一致的非刚性点集配准方法，是通过计算机，并按如下步骤进行：

步骤1，向计算机输入由三维点云获取设备获取的源点云和目标点云，获取这两个点云的关键点集：源关键点集S和目标关键点集O。

步骤2，由步骤1中获取的两个关键点集，转换获取关键点集之间的对应关系。

步骤3，建立源关键点集S的邻域索引矩阵。

步骤4，获得空间变换矩阵，并以此获得变换后的源关键点集S'。

步骤5，设定迭代的参数值，迭代执行步骤2到步骤4，当达到参数值时停止迭代，将此时的源关键点集S'和目标关键点集O作为配准结果并输出。

进一步说，本发明的详细步骤为：

步骤1，获取两个不同视角的点云，对这些点云进行关键点提取，获取各自的关键点集，即源关键点集S和目标关键点集O。此处的点云是指一些只包含三维坐标信息的点的集合，关键点提取是指对点云进行简化，获得保留点云基本形状特征和空间结构信息的一些点，这些点也就是关键点。

步骤2，对步骤1中获取的两个关键点集，计算源关键点集S中的每个点在目标关键点集O中的对应点，获取这两个关键点集之间的对应关系。

步骤3，对源关键点集S中的每一个点进行K近邻搜索，即对关键点集S中每一个点，在整个源关键点集S中找到与其欧几里得距离(Euclidean distance，简称欧式距离)最近的K个点，这K个点便称为它的相邻点，并根据这些相邻点在源关键点集S中的位置信息，建立源关键点集S的邻域索引矩阵。此处的欧式距离是指两个三维点之间的实际距离，K取值不小于2，优选地取值为3。

步骤4，利用相邻点的空间变换一致，对非刚性变换进行局部约束，即对源关键点集S中的每一个点的变换都进行局部约束，计算每一个点的空间变换分别与其K个相邻点的空间变换的差异，在保证差异较小的情况下获得较优的空间变换矩阵，以此获得变换后的源关键点集S'。

步骤5，设定一个最大迭代次数和一个参数阈值，迭代执行步骤2到步骤4，参数值超过阈值或者迭代次数超过最大迭代次数，停止迭代，此时经空间变换矩阵变换后的源关键点集S'及目标关键点集O作为配准结果并输出。

本发明的优选方法如下：

步骤1，点云获取与提取点云关键点：获取仅包含三维坐标信息的两个点云，即源点云和目标点云，分别对这两个点云采用降采样的方式提取关键点，得到源关键点集S和目标关键点集O。所述的降采样是指对点云进行稀疏化处理，源关键点集S和目标关键点集O分别对应的是源点云和目标点云进行降采样后的数量不大于1万个点的两个点集，这两个点集分别保留了源点云和目标点云的形状特征及空间结构信息。

步骤2，计算源关键点集S与目标关键点集O之间的对应关系：在步骤1获取关键点集的基础上，对源关键点集S中的每一个点，计算它与目标关键点集O中的每一个点的欧式距离，依据对应的欧式距离确定两个关键点集中每个点与点之间的对应概率。

步骤3，建立源关键点集S的邻域索引矩阵：对步骤1获取的源关键点集S中的每一个点，在整个源关键点集S范围内计算与其欧式距离最近的K个相邻点，并记录这些相邻点在源关键点集S中的位置索引，根据已记录的索引建立源关键点集S的邻域索引矩阵。点集中的所有点都有标号的，比如第一个点，第二个点，......索引指的就是这样一个标号，用来标识相邻点属于源关键点集中的第几个点。

步骤4，计算源关键点集S的空间变换矩阵：根据源关键点集S及其邻域索引矩阵，定义源关键点集S中每个点的空间变换与其K个相邻点的空间变换的差异程度的函数式，即约束项。该约束项由源关键点集S中每个点的空间变换分别与其K个相邻点的空间变换的差值之和组成，用来约束相邻点之间的变换保持一致。数据项表示的是源关键点集S和目标关键点集O中每一对对应点之间的距离。约束项与数据项组成目标函数式的作用是最小化目标函数式，获取此时对应的空间变换矩阵，并将此空间变换矩阵作用到源关键点集S，获取新的源关键点集S’。此处的空间变换指一个点集位置发生变化后所经历的变换形式，该变换形式用一个矩阵表示，相应地，源关键点集S中的点和其相邻点的空间变换指的是一个点位置发生变化时所经历的变换形式，该变换形式用一个向量表示。

步骤5，迭代执行步骤2到步骤4：设置一个最大迭代次数和参数σ²的阈值，每循环执行一次步骤2到步骤4，便重新计算一下参数σ²的值，迭代次数加一。

当参数值大于阈值或者迭代次数达到最大迭代次数，停止迭代，此时将空间变换矩阵作用在源关键点集S上，得到的新的源关键点集S'，将其和目标关键点集O作为配准结果并输出。

当参数值不大于阈值、迭代次数未达到最大迭代次数时，说明对应点之间的距离并没有接近于0，没有实现将两个点集统一到一个坐标系下，即没有完成两个关键点集的配准，返回步骤2。

更进一步说，步骤1的实现过程具体如下：

步骤1.1，点云获取：获取两个不同视角的点云，利用三维几何处理软件，如MeshLab将这两个点云保存为只包含三维坐标信息的点云，得到源点云和目标点云。

步骤1.2，关键点提取：根据步骤1.1得到的源点云和目标点云，在源点云数据范围内建立一个网格，所述网格由三维体素组成，通过将源点云包围在三维体素网格中，每个体素中包含源点云的多个三维点，为简化点的数量，则在每个三维体素中选一个点，这样在该三维体素内所有点就可以用这一个点表示，完成了对点云数据的简化。对网格中所有体素处理后得到的点云即为源点云对应的源关键点集S。同理，目标点云对应一个目标关键点集O。

更进一步说，步骤2具体如下：

步骤2.1，在步骤1.2得到的源关键点集S和目标关键点集O的基础上，计算源关键点集S中的每个点与目标关键点集O中每个点之间的欧几里得距离：假设y_m为源关键点集S中的一个点，x_n为目标关键点集O中的一个点，则两点之间的距离公式如下：

d(x_n,y_m)＝||x_n-y_m||²     (2-1)

步骤2.2，计算高斯混合模型(GMM)概率密度函数:将源关键点集S_M×D＝(y₁,…,y_M)^T表示为高斯混合模型的质心，将目标关键点集O_N×D＝(x₁,…,x_N)^T作为由高斯混合模型生成的数据点，在步骤2.1得到的点与点之间的欧几里得距离基础上，建立高斯混合模型概率密度函数公式如下：

其中，M和N分别表示源关键点集S和目标关键点集O的点的个数，m表示高斯混合模型的第m个高斯分量，n表示目标关键点集O中点的下标，D表示两个关键点集的维数，取值为3，exp表示以自然常数e为底的指数函数。σ²表示每个高斯分量的协方差，初始值为：

假设GMM的所有分量都是独立同分布的，则联合高斯混合模型概率密度函数公式如下：

其中，M表示的是源关键点集S的点的个数，N表示的是目标关键点集O的点的个数，P(m)＝1/M表示每个高斯分量的隶属概率。

步骤2.3，计算源关键点集S中的每个点与目标关键点集O中每个点之间的对应关系：基于高斯混合模型的点集配准方法，当源关键点集S和目标关键点集O对齐时，对于一个已知点x_n，其与点y_m之间的对应关系利用最大化高斯混合模型的后验概率获得。在步骤2.2的基础上，计算高斯混合模型的后验概率，获得两个点集之间的对应概率。高斯混合模型的后验概率公式如下：

P(m|x_n)＝P(m)p(x_n|m)/p(x_n)   (2-5)。

更进一步说，步骤3具体如下：

源关键点集S邻域索引矩阵的建立：对步骤2中的源关键点集S,采用被广泛应用的kd-tree最近邻搜索算法，将源关键点集S中的每一个点作为查询点，检索在kd-tree树中与查询点距离最近的K个相邻点，计算这K个相邻点在源关键点集S中的位置索引，根据索引建立源关键点集S的邻域索引矩阵Idy＝[idy₁,…,idy_i,…,idy_M]^T。

更进一步说，步骤4具体如下：

步骤4.1，基于GMM点集配准的目标函数的获取：基于GMM点集配准将点集配准问题转化为概率密度估计问题，通过重新参数化GMM质心位置，使质心逐渐拟合数据点。由步骤2.2获得的联合GMM概率密度函数，则其负对数似然函数公式如下：

一般分别对相应的参数求导，令导数为零，可以得到GMM新的参数。然而，负对数似然函数，对数里面还有求和。实际上没有办法通过求导的方法来求负对数似然函数的最小值。于是，采用EM算法进行GMM参数估计问题。EM算法分为两步，第一步寻找目标函数，通过计算E的改变量，可以得到目标函数的公式如下：

其中，θ表示一组变换参数，p'(m|x_n)表示点x_n和点y_m的初始对应概率，T(y_m,θ)表示应用于源关键点集S的变换函数，

刚性变换与非刚性变换通过确定上式的T函数来区分，刚性变换定义为T(R,t)＝RS+t，其中R表示旋转矩阵，t表示位移向量，在被大家公用的CPD算法中将非刚性变换定义为一个基于高斯径向基函数的位移函数，公式如下：

T(S,W)＝S+GW      (4-3)

其中，G是一个M×M的高斯核矩阵，其元素为，

W是一个M×D的高斯核权重矩阵，通过规范权重矩阵W以强制运动一致性，使点集在配准期间保持整体空间的连通性。全局约束项表示如下：

E_g(W)＝Tr(W^TGW)     (4-4)

将全局约束项添加到目标函数中，则目标函数重新表示如下：

其中，G(m,·)对应矩阵G的第m行。

步骤4.2，构建局部约束项：在非刚性配准的背景下，相邻点的空间变换是一致的。基于此条件，建立了局部约束，旨在配准期间能够保持点集的局部结构。基于步骤3获得的源关键点集S的邻域索引矩阵，计算源关键点集S中每个点的空间变换分别与其K个相邻点空间变换的差值，并计算所有差值的和，使和的值尽可能的小来实现局部约束，源关键点集S空间变换的局部约束项公式如下：

其中，K取值为3，G(m,·)对应矩阵G的第m行，Idy(m,k)表示源关键点集S中第m个点的第k个相邻点的索引。

步骤4.3，建立最终的目标函数：基于步骤4.1获取的目标函数，步骤4.2得到的局部约束项，将局部约束项添加到步骤4.1获得的目标函数中，本发明的目标函数公式如下：

Q(W,σ²)＝Q_d(W)+λE_l(W)        (4-7)

其中，λ表示局部约束的权重系数，取值为50000。

然后执行EM算法的第二步，最小化目标函数Q，分别对W和σ²求导，

对W求导：

对σ²求导：

其中，P是一个M×N的矩阵，其元素值表示源关键点集S中的每个点与目标关键点集O中每个点之间的对应概率。P1是P与值全为1的列向量的乘积，d(P1)表示由向量P1组成的对角矩阵。

令式(4-8)与式(4-9)的和为0，得到W的值，令式(4-10)为0，得到σ²的值，通过S+GW变换源关键点集S得到变换后的源关键点集S'，将变换后的源关键集S'作为下次迭代的源关键点集S。

更进一步说，其特征在于：

步骤5具体如下：设定一个最大迭代次数和参数σ²的阈值，循环执行步骤2到步骤4，每循环执行一次，先对源关键点集S进行变换，即S'＝S+GW，然后计算此时的σ²值，迭代次数相应加一，判断迭代次数是否超过最大迭代次数，比较σ²的值是否大于设定的阈值，当迭代次数超过最大迭代次数或者σ²的值大于给定的阈值时，迭代终止。若迭代未终止，将变换后的源关键点集S'作为下次迭代的源关键点集S。最终得到的源关键点集S'和目标关键点集O作为配准结果并输出。

实施例

本发明的实验硬件环境是：Intel(R)Xeon(R)CPU E3-1231 V3@3.40GHz3.40GHz,12G内存，Microsoft Windows10专业版，编程环境是MATLAB(R2017a)64位，测试图(详见图3、图4所示的实验配准效果图)来源于网上公开的点云数据。

实施例1，如图2所示，横坐标表示3个样例，纵坐标表示每一个样例的配准误差，配准误差指目标关键点集O和变换后的源关键点集S'对应点之间的平均欧式距离。虚线表示CPD算法配准误差，实线表示本发明配准误差。可以看出本发明的配准误差是小于CPD算法的。

实施例2，源关键点集S和目标关键点集O配准效果对比图。如图3所示，实验硬件环境是Intel(R)Xeon(R)CPU E3-1231 V3@3.40GHz 3.40GHz,12G内存，Microsoft Windows10专业版，编程环境是MATLAB(R2017a)64位。第一列第一行的两个点集分别是经过关键点提取后获得的源关键点集S和目标关键点集O，第二行和第三行的点集分别对应CPD算法和本发明的配准效果图。第二列为第一列的侧面效果展示。从对比图中，可以看出本发明在肘关节和膝关节周围的配准效果明显好于CPD算法。

实施例3，源关键点集S和目标关键点集O配准效果对比图。如图4所示，实验硬件环境是Intel(R)Xeon(R)CPU E3-1231 V3@3.40GHz 3.40GHz,12G内存，Microsoft Windows10专业版，编程环境是MATLAB(R2017a)64位。第一列第一行的两个点集分别是经过关键点提取后获得的源关键点集S和目标关键点集O，第二行和第三行的点集分别对应CPD算法和本发明的配准效果图。第二列为第一列的侧面效果展示。从对比图中，可以看出本发明在膝关节和脚周围的配准效果明显好于CPD算法。

通过与CPD算法定量和定性的分析，可以看出本发明的配准方法在点云变形较大情况下取得了较好的配准效果，由此证明了本发明的意义。

Claims (4)

1.一种基于局部变换一致的非刚性点集配准方法，其特征在于：通过计算机，并按如下步骤进行：

步骤1，向计算机输入由三维点云获取设备获取的源点云和目标点云，获取这两个点云的关键点集：源关键点集S和目标关键点集O；具体如下：

步骤1.1，点云获取：获取两个不同视角的点云，利用三维几何处理软件将这两个点云保存为只包含三维坐标信息的点云，得到源点云和目标点云；

步骤1.2，关键点提取：根据步骤1.1得到的源点云和目标点云，在源点云数据范围内建立一个网格，所述网格由三维体素组成，通过将源点云包围在三维体素网格中，每个体素中包含源点云的三维点，三维点在一个以上，为简化点的数量，则在每个三维体素中选一个点，这样在该三维体素内所有点就可以用这一个点表示，完成了对点云数据的简化；对网格中所有体素处理后得到的点云即为源点云对应的源关键点集S；同理，目标点云对应一个目标关键点集O

步骤2，由步骤1中获取的两个关键点集，转换获取关键点集之间的对应关系；

步骤3，建立源关键点集S的邻域索引矩阵；邻域索引矩阵的获取方法为：采用被广泛应用的kd-tree最近邻搜索算法，将源关键点集S中的每一个点作为查询点，检索在kd-tree树中与查询点距离最近的K个相邻点，计算这K个相邻点在源关键点集S中的位置索引，根据索引建立源关键点集S的邻域索引矩阵Idy＝[idy₁,···,idy_i,···,idy_M]^T；M表示的是源关键点集S的点的个数；

步骤4，根据源关键点集S及其邻域索引矩阵，定义源关键点集S中每个点的空间变换与其K个相邻点的空间变换的差异程度的函数式，即约束项；该约束项由源关键点集S中每个点的空间变换分别与其K个相邻点的空间变换的差值之和组成，用来约束相邻点之间的变换保持一致；数据项表示的是源关键点集S和目标关键点集O中每一对对应点之间的距离；约束项与数据项组成目标函数式的作用是最小化目标函数式，获取此时对应的空间变换矩阵，并将此空间变换矩阵作用到源关键点集S，获取新的源关键点集S'；

步骤5，设定迭代的参数值，迭代执行步骤2到步骤4，当达到参数值时停止迭代，将此时的源关键点集S'和目标关键点集O作为配准结果并输出。

2.根据权利要求1所述的一种基于局部变换一致的非刚性点集配准方法，其特征在于：在步骤1中，获取两个不同视角的点云，对这些点云进行关键点提取，获取各自的关键点集，即源关键点集S和目标关键点集O；

在步骤2中，对步骤1中获取的两个关键点集，计算源关键点集S中的每个点在目标关键点集O中的对应点，获取这两个关键点集之间的对应关系；

在步骤3中，对源关键点集S中的每一个点进行K近邻搜索，即对关键点集S中每一个点，在整个源关键点集S中找到与其欧几里得距离最近的K个点，这K个点便称为它的相邻点，并根据这些相邻点在源关键点集S中的位置信息，建立源关键点集S的邻域索引矩阵；

在步骤4中，利用相邻点的空间变换一致，对非刚性变换进行局部约束，即对源关键点集S中的每一个点的变换都进行局部约束，计算每一个点的空间变换分别与其K个相邻点的空间变换的差异，在保证差异较小的情况下获得较优的空间变换矩阵，以此获得变换后的源关键点集S'；

在步骤5中，设定一个最大迭代次数和一个参数阈值，迭代执行步骤2到步骤4，参数值超过阈值或者迭代次数超过最大迭代次数，停止迭代，此时经空间变换矩阵变换后的源关键点集S'及目标关键点集O作为配准结果并输出。

3.根据权利要求1或2所述的一种基于局部变换一致的非刚性点集配准方法，其特征在于：详细步骤如下：在步骤1中，点云获取与提取点云关键点：获取仅包含三维坐标信息的两个点云，即源点云和目标点云，分别对这两个点云采用降采样的方式提取关键点，得到源关键点集S和目标关键点集O；所述的降采样是指对点云进行稀疏化处理，源关键点集S和目标关键点集O分别对应的是源点云和目标点云进行降采样后的数量不大于1万个点的两个点集，这两个点集分别保留了源点云和目标点云的形状特征及空间结构信息；

在步骤2中，计算源关键点集S与目标关键点集O之间的对应关系：在步骤1获取关键点集的基础上，对源关键点集S中的每一个点，计算它与目标关键点集O中的每一个点的欧式距离，依据对应的欧式距离确定两个关键点集中每个点与点之间的对应概率；

在步骤3中，建立源关键点集S的邻域索引矩阵：对步骤1获取的源关键点集S中的每一个点，在整个源关键点集S范围内计算与其欧式距离最近的K个相邻点，并记录这些相邻点在源关键点集S中的位置索引，根据已记录的索引建立源关键点集S的邻域索引矩阵；

在步骤5中，迭代执行步骤2到步骤4：设置一个最大迭代次数和参数σ²的阈值，每循环执行一次步骤2到步骤4，便重新计算一下参数σ²的值，迭代次数加一；σ²表示每个高斯分量的协方差；σ表示标准差；

当参数值大于阈值或者迭代次数达到最大迭代次数，停止迭代，此时将空间变换矩阵作用在源关键点集S上，得到新的源关键点集S'，将其和目标关键点集O，作为配准结果并输出；

当参数值不大于阈值、迭代次数未达到最大迭代次数时，说明对应点之间的距离并没有接近于0，没有实现将两个点集统一到一个坐标系下，即没有完成两个关键点集的配准，返回步骤2。

4.根据权利要求1或2所述的一种基于局部变换一致的非刚性点集配准方法，其特征在于：步骤2具体如下：

步骤2.1，在步骤1.2得到的源关键点集S和目标关键点集O的基础上，计算源关键点集S中的每个点与目标关键点集O中每个点之间的欧几里得距离：假设y_m为源关键点集S中的一个点，x_n为目标关键点集O中的一个点，则两点之间的距离公式如下：

d(x_n,y_m)＝||x_n-y_m||²     (2-1)

步骤2.2，计算高斯混合模型概率密度函数:将源关键点集S_M×D＝(y₁,···,y_M)^T表示为高斯混合模型的质心，将目标关键点集O_N×D＝(x₁,···,x_N)^T作为由高斯混合模型生成的数据点，高斯混合模型简称GMM，在步骤2.1得到的点与点之间的欧几里得距离基础上，建立高斯混合模型概率密度函数公式如下：

其中，M和N分别表示源关键点集S和目标关键点集O的点的个数，m表示高斯混合模型的第m个高斯分量，n表示目标关键点集O中点的下标，D表示两个关键点集的维数，取值为3，exp表示以自然常数e为底的指数函数；σ²表示每个高斯分量的协方差，初始值为：

假设GMM的所有分量都是独立同分布的，则联合高斯混合模型概率密度函数公式如下：

其中，M表示的是源关键点集S的点的个数，N表示的是目标关键点集O的点的个数，P(m)＝1/M表示每个高斯分量的隶属概率；

步骤2.3，计算源关键点集S中的每个点与目标关键点集O中每个点之间的对应关系：基于高斯混合模型的点集配准方法，当源关键点集S和目标关键点集O对齐时，对于一个已知点x_n，其与点y_m之间的对应关系利用最大化高斯混合模型的后验概率获得；在步骤2.2的基础上，计算高斯混合模型的后验概率，获得两个点集之间的对应概率；高斯混合模型的后验概率公式如下：

P(m|x_n)＝P(m)p(x_n|m)/p(x_n)    (2-5)

P(m)表示每个高斯分量的隶属概率；m表示高斯混合模型的第m个高斯分量；x_n为已知点；p(x_n|m)高斯混合模型概率密度函数公式。

Images