CN104835153A - 基于稀疏表示的非刚性表面对齐方法 - Google Patents

Abstract

本发明属于计算机应用领域、三维测量、配准，为在三维测量中更好的识别正常的数据和野值。为此，本发明采取的技术方案是，基于稀疏表示的非刚性表面对齐方法，利用表面与纹理的特征描述SHOT(Unique signatures of histograms for surface and texture description)为模板点集中的每个点找到在目标点集中的对应点，用tv-l₁建模对非刚性配准的数据模型进行约束，来将具有形变的三维表面数据进行非刚性的配准；在此过程中，寻找对应点，利用全变差-1范数tv-l₁建模与对方程组求最优解。本发明主要应用于三维测量。

Description

基于稀疏表示的非刚性表面对齐方法

技术领域

本发明属于计算机应用领域、三维测量、配准，具体讲，涉及基于稀疏表示的非刚性表面对齐方法。

背景技术

在计算机应用领域，三维曲面配准是一个非常重要的中间步骤，它在表面重建、三维物体识别、机器导航等问题中有着极其重要的应用。许多扫描系统只能提供一部分的表面数据，必须经过配准融合之后才可以获取物体完整的数字模型。在获取数字模型时，由于自遮挡和传感器范围的限制，需要从多个角度对物体进行扫描。或者，在有了一个完整的模型之后，只需要扫描动态物体的一部分，将部分数据与完整的模型进行配准，便可以得到一个新的完整的模型。近日，一些三维扫描设备使得大家对鲁棒的刚性配准算法更加感兴趣，一些采集设备使用退化的数据来降低成本，这就需要可以处理大量噪声和离群值的配准算法。

目前为止，多数的曲面配准算法都主要针对刚性配准，也就是说从源点到目标点的运动时刚性的。针对刚性配准，研究者提出了很多解决方案，如点标记法、自旋图像、主曲率方法、遗传算法、随机采样一致性算法等等，这些算法各有特色，在许多特定的情况下能够解决配准的问题。但是应用最广泛的，影响最大的还是由Besl等在1992年提出的ICP(IterativeClosest Point，迭代最近点)算法，它是基于纯粹几何模型的三维物体对准算法，由于它的强大功能以及高的精确度，简单易解的优点，很快就成为了曲面配准中的主流算法。随着迭代最近点算法的广泛应用，许多研究者对迭代最近点算法做了详细的研究，分析了该算法的缺陷和特点，提出了许多有价值的改进，推动了这一重要算法的发展。

由于近期刚性配准取得的进展，加上快速发展的扫描设备逐渐可以捕捉随时间变化的三维曲面，非刚性配准逐渐受到了关注。与刚性配准相比，非刚性配准仍然处于起步阶段。与刚性配准技术相似，传统的非刚性配准技术同样是使用迭代最近点与最小二乘法的传统ICP方法。这种方法简单易懂，复杂度相对较低。所以有很多有用的技术已经被提出，改变了数字几何处理的现状，慢慢着眼于动态场景和动态动作。梯度下降法，高斯牛顿算法是曲面配准中最常用的优化技术。近年来，非刚性匹配的应用已经越来越广泛。但是这种利用传统ICP算法来解决非刚性配准问题时，同样对野值比较敏感，数据丢失的问题也是一个大的挑战。

发明内容

为克服现有技术的不足，在三维测量中更好的的识别正常的数据和野值。为此，本发明采取的技术方案是，基于稀疏表示的非刚性表面对齐方法，利用表面与纹理的特征描述SHOT(Unique signatures of histograms for surface and texture description)为模板点集中的每个点找到在目标点集中的对应点，用tv-l₁建模对非刚性配准的数据模型进行约束，来将具有形变的三维表面数据进行非刚性的配准；在此过程中，寻找对应点，利用全变差-1范数tv-l₁建模与对方程组求最优解。

为模板点集中的每个点找到在目标点集中的对应点具体步骤是，将高维特征之间的相似性进行排序，只取前面一部分更为可靠的对应点，找到特殊点的对应关系之后，再根据特殊点的对应关系，用k-d树最近邻搜索的方法来找目标点集所有剩余点的对应点。

将非刚性配准问题利用tv-l₁建模具体步骤是：

$\begin{array}{c}{E}_x={\mathrm{\α}}_{\mathrm{rigid}}{E}_{\mathrm{rigid}}+{\mathrm{\α}}_{\mathrm{smooth}}{E}_{\mathrm{smooth}}\\ ={\mathrm{\Σ}}_{i=1}^{N_1}{\left|\right|q_i-(R_i\×p_i+t_i)\left|\right|}_F^2+{\mathrm{\Σ}}_{i,j\∈\mathrm{\ϵ}}{\left|\right|R_i-R_j\left|\right|}_1\end{array}---\left(1\right)$

设是两个点集，点的个数分别为N₁、N₂，代表同一个三维物体扫描结果的两个有重叠的表面，设P为模板(source)点集，Q为目标(target)点集，ε表示点集中三角面片的边，p_i是模板点集中的第i个点，q_i是目标点集中与之相对的对应点，R_i、R_j表示第i、j个点的旋转矩阵，表示矩阵的F范数。||·||₁表示矩阵的1范数。t_i表示第i个点的平移量，在获得了较好的初值情况下，迭代最近点算法可以得到很好的收敛性，其中α_rigid与α_smooth分别表示两项所具有的权重，E_rigid表示数据项的目标函数能量，用来保证经旋转平移变换后对应点之间的距离总和最小，E_smooth表示平滑项的目标函数能量，保证临近的点拥有尽量相似的变换，E_x是目标方程总的能量。将上述方程重新描述，

$\min \mathrm{\α}{\left|\right|V\left|\right|}_1+{\left|\right|W(\mathrm{DX}-U)\left|\right|}_F^{2}s.t.V=M\⊗G---\left(2\right)$

s.t.表示使满足，表示克罗内克积，设每个模板点集中的点的变换参数是一个3×4的变换矩阵X_i，模板中找到对应点的总数为n，则最终要求的未知数为一个4n×3的矩阵X：[X₁ … X_n]^T，每个模板中顶点的坐标为v_i，设v_i＝[x，y，z，1]^T，与之相对应的目标点集中的点为u_i，用w_i表示匹配的可靠性，如果模板中的点v_i在目标点集中没有找到与之相对应的点，则将w_i设为0，否则设为1，W＝[w₁，…，w_n]，U＝[u₁ … u_n]^T，上标T表示转置，α表示数据项的权重。

定义稀疏矩阵D、M，得到要解的最优化方程具体步骤是：

1)

则其中U＝[u₁ … u_n]^T，I₃是3*3的单位阵；

2)定义稀疏矩阵M，M的行数为模板点集中边的个数，M的列数为模板点集中定点的个数如果第r条边连接第i个和第j个定点，则M_r，i＝1，M_r，j＝-1，G＝diag(1，1，1，1)。此时， $E_{\mathrm{smooth}}={\left|\right|(M\⊗G)X\left|\right|}_1.$

对方程组求最优解是利用增广拉格朗日方法进行最终求解，具体步骤是：

1)用表示系数矩阵的库SuiteSparse库来表示稀疏矩阵、进行克列斯基Cholesky分解、并求解稀疏矩阵的线性方程。

2)总式： $E_{\mathrm{\&mu;}}(V,X,Y)=\mathrm{\&alpha;}{\left|\right|V\left|\right|}_1+{\left|\right|\mathrm{DX}-U\left|\right|}_F^2+<Y,V-\mathrm{BX}>+\frac{\mathrm{\&mu;}}{2}{\left|\right|V-\mathrm{BX}\left|\right|}_F^2,$ 其中，V＝BX，E是目标方程的总能量，μ是一个正数，选择一个恰当的μ使算法收敛，Y是拉格朗日乘子：

$V^{k+1}=\arg {\min F}_{{\mathrm{\&mu;}}^k}(V,X^k,Y^k),X^{k+1}=\arg \min F_{{\mathrm{\&mu;}}^k}(V^{k+1},X,Y^k);$

对变量V与X进行迭代最小化，更新拉格朗日乘子Y，最终得到需要的变换矩阵X。

与已有技术相比，本发明的技术特点与效果：

本发明方法针对具有形变的非刚性数据进行配准，在迭代最近点的基础上，提出了基于稀疏表示的非刚性配准，1范数比2范数更好的描述有噪声的扫描模型，可以更好的的识别正常的数据和野值，具有以下特点：

1、程序简单，易于实现。

2、利用SHOT方法采用高维特征进行特征点的寻找，这样找到的特征点更加可靠。

3、由于0范数的非凸特性，使得求解变得非常困难，所以我们采用L0范数的最优凸近似L1范数进行约束，1范数最小化是凸优化问题，可以进行线性方程的求解。

4、考虑到三维数据点集数据庞大，用SuiteSparse库来表示稀疏矩阵、进行Cholesky分解、并求解稀疏矩阵的线性方程。

附图说明

本发明上述的和/或附加的方面和优点从下面结合附图对实施例的描述中将变得明显和容易理解:

图1为本发明方法的总框图；

图2为利用标准数据库的数据，寻找对应点的结果,红色的点云代表目标点集，绿色的点云代表源点集；

图3为利用标准数据库的数据，得到的最终配准结果，红色的点云代表目标点集，绿色的点云代表源点集。

具体实施方式

本发明用SHOT(Unique signatures of histograms for surface and texture description，表面与纹理的特征描述)方法为模板点集中的每个点找到在目标点集中的对应点，用tv-l₁建模对非刚性配准的数据模型进行约束，来将具有形变的三维表面数据进行非刚性的配准；在此过程中，寻找对应点，利用tv-l₁建模与对方程组求最优解，具体方法包括以下步骤：

1)为模板点集中的每个点找到在目标点集中的对应点。用SHOT方法为模板点集中的每个点找到在目标点集中的对应点，为了确保准确性，我们将高维特征之间的相似性进行排序，只取前面一部分更为可靠的对应点。找到特殊点的对应关系之后，再根据特殊点的对应关系，用k-d树最近邻搜索的方法来找目标点集所有剩余点的对应点。

2)将非刚性配准问题利用tv-l₁建模：

$\begin{array}{c}{E}_x={\mathrm{\&alpha;}}_{\mathrm{rigid}}{E}_{\mathrm{rigid}}+{\mathrm{\&alpha;}}_{\mathrm{smooth}}{E}_{\mathrm{smooth}}\\ ={\mathrm{\&Sigma;}}_{i=1}^{N_1}{\left|\right|q_i-(R_i\&times;p_i+t_i)\left|\right|}_F^2+{\mathrm{\&Sigma;}}_{i,j\&Element;\mathrm{\&epsiv;}}{\left|\right|R_i-R_j\left|\right|}_1\end{array}---\left(1\right)$

设是两个点集，点的个数分别为N₁、N₂，代表同一个三维物体扫描结果的两个有重叠的表面。设P为模板(source)点集，Q为目标(target)点集。ε表示点集中三角面片的边，p_i是模板点集中的第i个点，q_i是目标点集中与之相对的对应点。R_i表示第i个点的旋转矩阵，t_i表示第i个点的平移量。表示矩阵的F范数。||·｜|₁表示矩阵的1范数。在获得了较好的初值情况下，迭代最近点算法可以得到很好的收敛性，其中α_rigid与α_smooth分别表示两项所具有的权重。E_rigid表示数据项的目标函数能量，用来保证经旋转平移变换后对应点之间的距离总和最小，E_smooth表示平滑项的目标函数能量，保证临近的点拥有尽量相似的变换，E_x是目标方程总的能量。将上述方程重新描述，

$\min \mathrm{\&alpha;}{\left|\right|V\left|\right|}_1+{\left|\right|W(\mathrm{DX}-U)\left|\right|}_F^{2}s.t.V=M\&CircleTimes;G---\left(2\right)$

其中矩阵的构建，方法见4)，α表示数据项的权重。

3)为模板点集中的点在目标点集中寻找对应点。

31)用SHOT方法，找到前面一部分的可靠对应点，此时所找到点w_i设为1，其余设为0，在第一次迭代以后，除了前面的可靠对应点，还寻找每个点的最近点，如果两个对应点之间的距离小于某个阈值，则认为是可靠的对应点，所找到点w_i设为1，其余设为0，W＝[w₁，…，w_n]。

4)定义稀疏矩阵D、M，得到要解的最优化方程：

$\min \mathrm{\&alpha;}{\left|\right|V\left|\right|}_1+{\left|\right|W(\mathrm{DX}-U)\left|\right|}_F^{2}s.t.V=M\&CircleTimes;G$

41)设每个模板点集中的点的变换参数是一个3×4的变换矩阵X_i，模板中找到对应点的总数为n，则我们最终要求的未知数为一个4n×3的矩阵X：[X₁ … X_n]^T。每个模板中顶点的坐标为v_i，设v_i＝[x，y，z，1]^T，与之相对应的目标点集中的点为u_i，用w_i表示匹配的可靠性，如果模板中的点v_i在目标点集中没有找到与之相对应的点，则将w_i设为0，否则设为1。

则其中U＝[u₁ … u_n]^T。

42)定义稀疏矩阵M，M的行数为模板点集中边的个数，M的列数为模板点集中定点的个数如果第r条边连接第i个和第j个定点，则M_r，i＝1，M_r，j＝-1，G＝diag(1，1，1，1)。此时， $E_{\mathrm{smooth}}={\left|\right|(M\&CircleTimes;G)X\left|\right|}_1.$

5)利用增广拉格朗日方法进行最终求解。

51)用SuiteSparse库来表示稀疏矩阵、进行Cholesky分解、并求解稀疏矩阵的线性方程。

52)总式： $E_{\mathrm{\&mu;}}(V,X,Y)=\mathrm{\&alpha;}{\left|\right|V\left|\right|}_1+{\left|\right|\mathrm{DX}-U\left|\right|}_F^2+<Y,V-\mathrm{BX}>+\frac{\mathrm{\&mu;}}{2}{\left|\right|V-\mathrm{BX}\left|\right|}_F^2,$ 其中，Y是拉格朗日乘子，E是目标方程的总能量，μ是一个正数，选择一个恰当的μ可以使我们的算法收敛，这里μ取1.2。V-problem:X-problem: $X^{k+1}=\arg \min F_{{\mathrm{\&mu;}}^k}(V^{k+1},X,Y^k).$ V＝BX。

对变量V与X进行迭代最小化，更新拉格朗日乘子Y，最终可以得到我们需要的变换矩阵X。

下面结合附图和具体实例，进一步详细说明本发明。

本发明针对具有形变的非刚性数据进行配准，在迭代最近点的基础上，提出了基于稀疏表示的非刚性配准，用SHOT方法为模板点集中的每个点找到在目标点集中的对应点，用tv-l₁建模对非刚性配准的数据模型进行约束，来将具有形变的三维表面数据进行非刚性的配准。

本发明的主要框架迭代最近点有两个步骤，首先为模板点集中的每个点找到在目标点集中的对应点。用SHOT方法与最近邻搜索的方法来找目标点集的对应点，然后解优化方程，得出最优解X。将X代入第一步骤中重新找对应点，这样迭代多次之后便可以得到一个最优解。

1)为模板点集中的每个点找到在目标点集中的对应点。用SHOT方法为模板点集中的每个点找到在目标点集中的对应点，为了确保准确性，我们将高维特征之间的相似性进行排序，只取前面一部分更为可靠的对应点。找到特殊点的对应关系之后，再根据特殊点的对应关系，用k-d树最近邻搜索的方法来找目标点集所有剩余点的对应点。

2)将非刚性配准问题利用TV-l₁建模：

$\begin{array}{c}{E}_x={\mathrm{\&alpha;}}_{\mathrm{rigid}}{E}_{\mathrm{rigid}}+{\mathrm{\&alpha;}}_{\mathrm{smooth}}{E}_{\mathrm{smooth}}\\ ={\mathrm{\&Sigma;}}_{i=1}^{N_1}{\left|\right|q_i-(R_i\&times;p_i+t_i)\left|\right|}_F^2+{\mathrm{\&Sigma;}}_{i,j\&Element;\mathrm{\&epsiv;}}{\left|\right|R_i-R_j\left|\right|}_1\end{array}---\left(1\right)$

是两个点集，代表同一个三维物体扫描结果的两个有重叠的表面。称P为模板点集，Q为目标点集。ε表示点集中三角面片的边。其中α_rigid与α_smooth分别表示两项所具有的权重。E_rigid表示数据项的目标函数能量，用来保证经旋转平移变换后对应点之间的距离总和最小，E_smooth表示平滑项的目标函数能量，保证临近的点拥有尽量相似的变换，E_x是目标方程总的能量，R_i、R_j表示第i、j个点的旋转矩阵。将上述方程重新描述，

$\min \mathrm{\&alpha;}{\left|\right|V\left|\right|}_1+{\left|\right|W(\mathrm{DX}-U)\left|\right|}_F^{2}s.t.V=M\&CircleTimes;G---\left(2\right)$

其中矩阵的构建，方法见2)、3)。

3)为模板点集中的点在目标点集中寻找对应点。

31)用SHOT方法，找到前20％的可靠对应点，此时所找到点w_i设为1，其余设为0，在第一次迭代以后，除了前20％的可靠对应点，还寻找每个点的最近点，如果两个对应点之间的距离小于数据总长度的1/20，则认为是可靠的对应点，所找到点w_i设为1，其余设为0。

4)定义稀疏矩阵D、矩阵M。

41)设每个模板点集中的点的变换参数是一个3×4的变换矩阵X_i，模板中找到对应点的总数为n，则我们最终要求的未知数为一个4n×3的矩阵X：[X₁ … X_n]^T。每个模板中顶点的坐标为v_i，设v_i＝[x，y，z，1]^T，与之相对应的目标点集中的点为u_i，用w_i表示匹配的可靠性，如果模板中的点v_i在目标点集中没有找到与之相对应的点，则将w_i设为0，否则设为1，W＝[w₁，…，w_n]。

则 $E_{\mathrm{rigid}}={\left|\right|W(\mathrm{DX}-U)\left|\right|}_F^2.$ 其中

U＝[u₁ … u_n]^T，I₃是3*3的单位阵；

42)定义稀疏矩阵M，M的行数为模板点集中边的个数，M的列数为模板点集中定点的个数如果第r条边连接第i个和第j个定点，则M_r，i＝1，M_r，j＝-1，G＝diag(1，1，1，1)。此时， $E_{\mathrm{smooth}}={\left|\right|(M\&CircleTimes;G)X\left|\right|}_1.$

综合以上两式，要解的最优化方程为：

$\min \mathrm{\&alpha;}{\left|\right|V\left|\right|}_1+{\left|\right|W(\mathrm{DX}-U)\left|\right|}_F^{2}s.t.V=M\&CircleTimes;G$

α表示数据项的权重。

5)利用增广拉格朗日方法进行最终求解。

51)用SuiteSparse库来表示稀疏矩阵、进行Cholesky分解、并求解稀疏矩阵的线性方程。52)总式： $E_{\mathrm{\&mu;}}(V,X,Y)=\mathrm{\&alpha;}{\left|\right|V\left|\right|}_1+{\left|\right|\mathrm{DX}-U\left|\right|}_F^2+<Y,V-\mathrm{BX}>+\frac{\mathrm{\&mu;}}{2}{\left|\right|V-\mathrm{BX}\left|\right|}_F^2,$ 其中，Y是拉格朗日乘子。第k+1次迭代时：

$X^{k+1}=\arg \min F_{{\mathrm{\&mu;}}^k}(V^{k+1},X,Y^k).$

对变量V与X进行迭代最小化，更新拉格朗日乘子Y，最终可以得到我们需要的变换矩阵X。α根据每次迭代中E_rigid项与E_smooth项的值来选取，使得每次E_rigid项与E_smooth项所占比重相同。

Claims (5)

1.一种基于稀疏表示的非刚性表面对齐方法，其特征是，利用表面与纹理的特征描述SHOT(Unique signatures of histograms for surface and texture description)为模板点集中的每个点找到在目标点集中的对应点，用tv-l₁建模对非刚性配准的数据模型进行约束，来将具有形变的三维表面数据进行非刚性的配准；在此过程中，寻找对应点，利用全变差-1范数tv-l₁建模与对方程组求最优解。

2.如权利要求1所述的基于稀疏表示的非刚性表面对齐方法，其特征是，为模板点集中的每个点找到在目标点集中的对应点具体步骤是，将高维特征之间的相似性进行排序，只取前面一部分更为可靠的对应点，找到特殊点的对应关系之后，再根据特殊点的对应关系，用k-d树最近邻搜索的方法来找目标点集所有剩余点的对应点。

3.如权利要求1所述的基于稀疏表示的非刚性表面对齐方法，其特征是，将非刚性配准问题利用tv-l₁建模具体步骤是：

$\begin{array}{c}{E}_x={\mathrm{\&alpha;}}_{\mathrm{rigid}}{E}_{\mathrm{rigid}}+{\mathrm{\&alpha;}}_{\mathrm{smooth}}{E}_{\mathrm{smooth}}\\ ={\mathrm{\&Sigma;}}_{i=1}^{N_1}{\left|\right|q_i-(R_i\&times;p_i+t_i)\left|\right|}_F^2+{\mathrm{\&Sigma;}}_{i,j\&Element;\mathrm{\&epsiv;}}{\left|\right|R_i-R_j\left|\right|}_1\end{array}---\left(1\right)$

设是两个点集，点的个数分别为N₁、N₂，代表同一个三维物体扫描结果的两个有重叠的表面，设P为模板(source)点集，Q为目标(target)点集，

ε表示点集中三角面片的边，p_i是模板点集中的第i个点，q_i是目标点集中与之相对的对应点，R_i、R_j表示第i、j个点的旋转矩阵，表示矩阵的F范数，||·||₁表示矩阵的1范数，t_i表示第i个点的平移量，在获得了较好的初值情况下，迭代最近点算法可以得到很好的收敛性，其中α_rigid与α_smooth分别表示两项所具有的权重，E_rigid表示数据项的目标函数能量，用来保证经旋转平移变换后对应点之间的距离总和最小，E_smooth表示平滑项的目标函数能量，保证临近的点拥有尽量相似的变换，E_x是目标方程总的能量，将上述方程重新描述，

$\begin{array}{cc}\min & \mathrm{\&alpha;}{\left|\right|V\left|\right|}_1+{\left|\right|W(\mathrm{DX}-U)\left|\right|}_F^{2}s.t.V=M\&CircleTimes;G---\left(2\right)\end{array}$

s.t.表示使满足，表示克罗内克积，设每个模板点集中的点的变换参数是一个3×4的变换矩阵X_i，模板中找到对应点的总数为n，则最终要求的未知数为一个4n×3的矩阵X：[X₁ … X_n]^T，每个模板中顶点的坐标为v_i，设v_i＝[x，y，z，1]^T，与之相对应的目标点集中的点为u_i，用w_i表示匹配的可靠性，如果模板中的点v_i在目标点集中没有找到与之相对应的点，则将w_i设为0，否则设为1，W＝[w₁，…，w_n]，U＝[u₁ … u_n]^T，上标T表示转置，α表示数据项的权重。

4.如权利要求3所述的基于稀疏表示的非刚性表面对齐方法，其特征是，定义稀疏矩阵D、M，得到要解的最优化方程具体步骤是：

1)

则其中U＝[u₁ … u_n]^T，I₃是3*3的单位阵；

2)定义稀疏矩阵M，M的行数为模板点集中边的个数，M的列数为模板点集中定点的个数如果第r条边连接第i个和第j个定点，则M_r，i＝1，M_r，j＝-1，G＝diag(1，1，1，1)，此时， $E_{\mathrm{smooth}}={\left|\right|(M\&CircleTimes;G)X\left|\right|}_1.$

5.如权利要求1所述的基于稀疏表示的非刚性表面对齐方法，其特征是，对方程组求最优解是利用增广拉格朗日方法进行最终求解，具体步骤是：

1)用表示系数矩阵的库SuiteSparse库来表示稀疏矩阵、进行克列斯基Cholesky分解、并求解稀疏矩阵的线性方程；

2)总式： $E_{\mathrm{\&mu;}}(V,X,Y)=\mathrm{\&alpha;}{\left|\right|V\left|\right|}_1+{\left|\right|\mathrm{DX}-U\left|\right|}_F^2+<Y,V-\mathrm{BX}>+\frac{\mathrm{\&mu;}}{2}{\left|\right|V-\mathrm{BX}\left|\right|}_F^2,$ 其中，V＝BX，E是目标方程的总能量，μ是一个正数，选择一个恰当的μ使算法收敛，Y是拉格朗日乘子：

$V^{k+1}=\arg \min F_{{\mathrm{\&mu;}}^k}(V,X^k,Y^k),X^{k+1}=\arg \min F_{{\mathrm{\&mu;}}^k}(V^{k+1},X,Y^k);$

对变量V与X进行迭代最小化，更新拉格朗日乘子Y，最终得到需要的变换矩阵X。