CN101887525A - 基于分级的正反互逆的三维稠密点集快速配准方法 - Google Patents

Abstract

本发明提供一种基于分级的正反互逆的三维稠密点集快速配准方法，解决的技术问题是：给定两个点集，其中一个为参考点集，另一个为目标点集，现需求解两个变换函数；正变换函数h：从参考点集到目标点集的变换函数，反变换函数g：从目标点集到参考点集的变换函数，并要求h≈g^-1，g≈h^-1，其中h^-1和g^-1分别表示h和g的逆函数。最终通过这两个变换函数，得到两个点集的精确的对应关系。与现有技术相比，本方法的优点在于基于全局配准和局部配准的分级方式进行配准，并利用正方向和反方向上对应关系的差异来定义逆一致性误差和用于局部配准的特征点，避免了现有正反互逆配准算法中复杂耗时的函数求逆运算。实现简单、效率高，同时配准精度也高。

Description

基于分级的正反互逆的三维稠密点集快速配准方法

技术领域

本发明属于计算机应用领域，涉及计算机视觉和模式识别的一种三维稠密点集对点集的配准方法。

背景技术

配准是一个基础问题，它源自多个领域的很多实际问题，如不同传感器获得的信息融合；不同时间、条件获得图像或模型的差异检测；成像系统和物体场景变化情况下获得的图像的三维信息获取；图像中的模式或目标识别等。虽然参与配准的物体可以有多种形式，如三维物体的几何模型、二维或三维图像，但其实质是两个二维或三维点集的配准问题，即找到两个点集之间的对应关系。

迭代最近点算法(Iterated Closest Point，ICP)是最常用的点集配准算法，它通过迭代优化矩阵，在每次迭代过程中，对参考点集上的每个点，在目标点集中寻找最近点，并利用这样的对应点，计算相应的旋转矩阵和平移向量，将其用于参考点集上，得到新的参考点集并进入下次迭代过程，最终得到最优的转换矩阵，实现两点集的配准。但最初提出的ICP算法在寻找最近点时很耗时，因此很多学者提出了改进的ICP算法以提高配准效率。

ICP算法虽然使用广泛，但它属于刚性配准，在很多应用中不能满足需要，因为很多形变的性质是非线性的。因此必须采用非刚性配准算法，包括基于B样条或薄板样条(Thin PlateSpline，TPS)的弹性配准算法、基于互信息的弹性配准方法等。

通常一种配准算法得到的从参考点集到目标点集的变形函数(称之为正变换)和从目标点集到参考点集的变形函数(称之为反变换)是不互逆的。以基于薄板样条(TPS)的配准方法为例，虽然它能通过全局能量最低的限制得到从一个点集到另一个点集的平滑的变形函数，但它只在特征点上满足互逆的条件，在其他顶点不满足，因此会在一些局部区域出现大的配准误差。而正反变换互逆是建立两个点集元素之间一一对应关系的必要条件。为此一些学者展开了对正反互逆的配准算法的研究。提出了一些有效的方法，可实现两个点集的精确配准。但现有方法几乎都涉及逆函数求解，由于处理对象的离散性，求解过程耗时长，在很多具有实时性要求的领域和行业中得不到应用。虽然已有一些关于如何提高正反互逆配准效率的研究，但目前结果还不能满足实际应用要求。

发明内容

为克服上述现有技术存在的不足，本发明提供一种快速而又准确度高的三维点集对点集的配准方法，因为二维是三维的特例，所以也支持二维点集间的配准，主要用于三维人脸识别、CAD产品质量检测、医学图像融合和分析、三维模型统计等。

本发明欲解决的技术问题是：给定两个点集，其中一个点集为参考点集，表示成P，包含N_p个顶点，即{p_i，i＝1，...，N_p}；另一个点集为目标点集，表示成Q，包含N_q个顶点，即{q_j，j＝1，...，N_q}。现需快速求解两个变换函数；正变换函数h：从参考点集到目标点集的变换函数和反变换函数g：从目标点集到参考点集的变换函数，并要求h≈g^-1，g≈h^-1，其中h^-1和g^-1分别表示h和g的逆函数。最终通过这两个变换函数，得到两个点集的精确的对应关系。

为实现上述目的，本发明采用了基于全局和局部配准的分级的方法。全局配准为预配准，能使参考点集和目标点集大体匹配上，但在一些局部区域会因存在较大的逆一致性误差而出现错配现象，局部配准即对这些局部区域进行调整，逐步减小逆一致性误差，实现精确配准。根据这一思想。函数h和g被分解成一系列函数的复合，即h＝h_lnо...h_l1оh_g，g＝g_lnо...g_l1оg_g，其中h_g和g_g为全局的正反变换函数；h_li(i＝1，...，n)和g_li(i＝1，...，n)为局部的正反变换函数。此外，为了提高配准的效率，采用了简化的公式求取逆一致性误差，避免了耗时的函数求逆运算以及大矩阵运算。

本发明算法的主要步骤包括：

1)输入两个点集：参考点集P和目标点集Q；

2)全局配准；

全局配准可采用两种方法：薄板样条函数(TPS)算法或迭代最近点(ICP)算法。

如采用TPS方法，包含两个子步骤：

a)从两个点集中采用人机交互的方式选取两组对应特征点；

b)利用TPS算法和特征点，计算从参考点集到目标点集的变换函数(正变换h_g)，以及从目标点集到参考点集的变换函数(反变换g_g)，并把正反变换分别作用于参考点集和目标点集，得到变形后的两个点集P’和Q’。

如采用ICP方法，也包含两个子步骤：

a)判断待配准的两个点集中的顶点数量是否超过用户设定的阈值，如果是从两个点集中随机选取两组顶点，组成两个新的点集，新点集中各顶点间的距离不小于用户设定的距离阈值，且顶点数量满足数量阈值要求；

b)利用ICP算法计算两个新点集的正变换h_g和反变换g_g，并把正反变换分别作用于原始参考点集和目标点集，得到变形后的两个点集P’和Q’。

3)局部配准；包括局部特征点的自动定义和局部变形两个步骤。

3.1)局部特征点的自动定义，包含5个子步骤：

a)求变形后的参考点集上的每一点到目标点集上的最近点，如果两点之间的距离小于用户设定的阈值，则以该最近点作为它的对应点，最终得到一个对应点对列表；

b)求变形后的目标点集上的每一点到参考点集上的最近点，如果两点之间的距离小于用户设定的阈值，则以该最近点作为它的对应点，最终得到另一个对应点对列表；

c)根据前两个步骤得到的结果，计算参考点集和目标点集中存在对应点的每一个顶点的逆一致性误差，并对其对应点进行调整；

d)把两个对应点对列表中的经过调整操作的对应点对挑选出来，存成一个新列表，并根据逆一致性误差从大到小的顺序对新列表进行排序；

e)从排序后的列表中提取用于局部配准的两组特征点。

3.2)局部变形：采用具有紧支撑的径向基函数(CSRBF)，并根据3.1)步骤中得到的特征点计算正变换h_li和反变换g_li，并把正反变换分别作用于已变形的参考点集P’和目标点集Q’，得到两个新的变形点集，依然采用P’和Q’表示。

4)判断局部配准的次数是否超过用户设定的阈值，如果是，则配准结束，保存两个点集中顶点的对应关系；否则返回步骤3)继续执行局部配准。

与现有技术相比，本方法的特点在于采用基于全局配准和局部配准的分级方式进行配准，并利用正方向和反方向上对应关系的差异来定义逆一致性误差和用于局部配准的特征点，避免了现有正反互逆配准算法中复杂耗时的函数求逆运算。实现简单、效率高，同时配准精度也高。

附图说明

图1为本发明所述一种基于分级的正反互逆的三维稠密点集快速配准方法的算法流程图；

图2为本发明实施例1的两套颅骨点云模型的初始位置；

图3为本发明实施例1基于TPS的正方向和反方向的配准结果；

图4为本发明实施例1基于TPS和一次CSRBF的正方向和反方向的配准结果；

图5为本发明实施例1基于TPS和五次CSRBF的正方向和反方向的配准结果；

图6为本发明实施例1经TPS和若干次局部配准后目标点集的逆一致性误差；

图7为本发明实施例1经TPS和若干次局部配准后参考点集的逆一致性误差；

图8为本发明实施例1基于ICP和若干次CSRBF的正方向和反方向的配准结果；

图9为本发明实施例1经ICP和若干次局部配准后目标点集的逆一致性误差；

图10为本发明实施例1经ICP和若干次局部配准后参考点集的逆一致性误差；

图11为本发明实施例1的两套人脸点云模型的初始位置；

图12为本发明实施例1人脸在正反两个方向上的全局配准结果和最终结果。

具体实施方式

下面结合具体实施方式和优选的实施例作为具体的说明，但不作为对发明的限制。

实施例1：一种基于分级的正反互逆的三维稠密点集快速配准方法，其方法步骤如图1所示：

步骤1S1：给定两个点集；给定两个点集，其中一个点集为参考点集，表示成P，另一个点集为目标点集，表示成Q；

步骤2S2：进行全局配准；

全局配准可采用薄板样条算法(Thin-Plate-Spline，TPS)S200或者迭代最近点算法(IteratedClosest Point，ICP)S20。其中前者在配准过程中需要人机交互，后者是自动的。

1)基于TPS的全局配S200

薄板样条配准算法(TPS)属于非刚性配准算法，其函数f由两部分组成，第一部分由径向基函数表示的弹性变换构成；第二部分为全局仿射变换200。具体公式如下：

$f(x,y,z)=R_s(x,y,z)+{\mathrm{\φ}}_s(x,y,z)$

$={\mathrm{\Σ}}_{i=1}^n{\mathrm{\α}}_{i}U\left(\right||p_i^t-(x,y,z)|\left|\right)+{\mathrm{\β}}_1+{\mathrm{\β}}_{2}x+{\mathrm{\β}}_{3}y+{\mathrm{\β}}_{4}z---\left(1\right)$

其中n为特征点的个数，

为特征点

和顶点(x，y，z)之间的欧式距离，α_i(i＝1，...n)，β_j(j＝1，2，3，4)为待求的权重。

对于弹性变换部分，还有四个附加的边界条件，分别表示如下：

TPS在配准过程中能使变形模型的全局弯曲能量最低(公式(3))，所以TPS被认为是光滑性最好的一种配准算法之一。

为求解TPS函数f中的未知变量α_i(i＝1，...n)和β_j(j＝1，2，3，4)，需要以人机交互的方式在参考点集和目标点集中确定两组对应的特征点

和

其中n由用户设定，一般不超过100；然后把特征点

一一映射到对应特征点

如下公式所示：

$f\left(p_i^t\right)=q_i^t,i=\mathrm{1,2}...n---\left(4\right)$

这一求解过程通常表示成矩阵形式：

$\left(\begin{array}{cc}K& P\\ P^T& 0\end{array}\right)\left(\begin{array}{c}\mathrm{\α}\\ \mathrm{\β}\end{array}\right)=\left(\begin{array}{c}Q\\ 0\end{array}\right)---\left(5\right)$

其中K为n×n的矩阵，矩阵中的元素

P为n×4的矩阵，矩阵中的元素

α＝(α₁，α₂...α_n)^T，β＝(β₁，β₂，β₃，β₄)^T，

一旦确定了权重α_i(i＝1，...n)和β_j(j＝1，2，3，4)的值，就得到了全局的正变换函数h_g201，并可根据公式(1)对参考点集进行变形，使之匹配到目标点集上202。反变换函数g_g的求解过程类似。

2)基于ICP的全局配准S20

ICP算法的目标是求解两组点集P和Q的坐标系之间的3×3旋转矩阵R和平移向量S，使其满足以下最优化条件20：

$\min {\mathrm{\Σ}}_{k=1}^N{\left|\right|q_k-(R\·p_k+S)\left|\right|}^2$

当P和Q的顶点数量相同，并且一一对应时，可采用四元数奇异值分解(Singular ValueDecomposition，SVD)得到R和S，并可根据R和S构造4×4的齐次变换矩阵T。但在实际应用中无法保证P和Q满足上述条件。因此，ICP方法常用迭代方式寻找最优配准矩阵。具体步骤如下：

(1)粗配准：通常采用主元分析法(Principal Component Analysis，PCA)求得P和Q两组点集的3个特征向量，结合各自的重心，得到4个对应的特征点对，按照前述对应点集的配准方法求得初始配准矩阵T₀。

(2)精配准，包含4个步骤：

设k为迭代次数，N_p和N_q分别为点集P和Q的元素个数，ε为给定精度(例如ε＝10^-10)。

①令k：＝0，P⁰＝P，T⁰＝T⁰，对P⁰中的每一点

使用kd树在Q中寻找与

欧氏距离最近的点从而构成一个点集：

$\begin{array}{ccc}{D}^0\left(i\right)=q_j& \mathrm{st}& \underset{1\≤j{\≤N}_q}{\min }\left|\right|p_i^0-q_j\left|\right|\end{array}$

设置初始残差： $r^0={\mathrm{\Σ}}_{i=0}^{N_p}\left|\right|p_i^0-D^0\left(i\right)\left|\right|/N_p.$

②P^k+1＝T^k□P^k，D^k+1(i)＝q_j满足条件

更新平均残差：

③如果|r^k+1-r^k|≤ε或者k达到了迭代次数上限则算法结束，T_k·T_k-1…T₀即为所求的最终配准矩阵；否则继续④。

④对Q^k+1和D^k+1进行对应点集的SVD法求得齐次坐标下的配准矩阵T，令T^k+1＝T，令k：＝k+1，转步骤②。

ICP的配准效率直接与点集中的顶点数量相关，当顶点数量超过1万时，配准效率会很低。为了提高全局配准的速度，本发明设置了顶点数量上限：4000。当点集P和Q中顶点数量超过上限时，采用随机采样的方式从两个点集中选取4000个顶点，组成新的点集21。其中点集中的每两个顶点之间的距离不小于距离阈值，该阈值的初始值由用户设定，但随采样情况改变，利用新的点集执行ICP配准，得到旋转矩阵R和平移向量S后，即得到全局的正变换函数h_g22。以此函数作用于原始参考点集，完成全局配准23。反变换函数g_g的求解过程类似。

步骤3S3：经过全局配准后，两个点集在大部分区域能匹配上，但在一些局部区域存在较大配准误差。本发明采用基于紧支撑的径向基函数配准算法(Compact Support Radial BasisFunctions，CSRBF)对这些局部区域进行局部配准，最终达到精确配准的目的。

基于紧支撑的径向基函数配准算法CSRBF是一种基于特征点的配准算法，其特征点根据简化的逆一致性误差自动生成。

1)局部特征点的自动确定。局部配准中的特征点

和

根据逆一致性误差自动获得。通常，点集P中的顶点p_i的逆一致性误差定义为||h(p_i)-g^-1(p_i)||，同理，点集Q中的顶点q_i的逆一致性误差定义为||g(q_i)-h^-1(q_i)||。但这样的定义方式涉及到逆函数h^-1和g^-1的求解。由于点集的离散性，求解过程非常耗时。

为了加快配准速度，本发明采用了一种简单而又直观的方式求解顶点的逆一致性误差。理想状态下，点集P中的每一个顶点p_i通过正变换函数在映射到点集Q中的顶点q_i，同时顶点q_i也通过反变换函数映射到顶点p_i，反之亦然。这样就得到两个点集的一一对应关系，并且点集中任何一个顶点的逆一致性误差为0。但在实际应用中，两个点集的顶点数量不一定相同，顶点之间也不一定能一一对应。但逆一致性误差的本质不变。以点集P上的顶点为例。对于点集P中的一个顶点p_i，可通过查找最近点找到它在点集Q中的对应点q_i，同时点集Q中可能没有，也可能有几个顶点，记作点集Q_cor，在点集P上的对应点为p_i。这时q_i和Q_cor之间的差异就可用于衡量点p_i的逆一致性误差。它们之间的差异越大，则点p_i的逆一致性误差就越大。

根据逆一致性误差自动获取局部特征点的具体步骤如下：

(1)求变形后的参考点集上的每一点到目标点集上的最近点，如果两点之间的距离小于用户设定的阈值，则以该最近点作为它的对应点，得到一个对应点对列表30。

(2)求变形后的目标点集上的每一点到参考点集上的最近点，如果两点之间的距离小于用户设定的阈值，则以该最近点作为它的对应点，得到另一个对应点对列表31。

(3)求有对应点的顶点的逆一致性误差并对其对应点进行调整32。

以点集P中的顶点p_i为例(点集Q中的顶点q_i类似)，设顶点p_i的在目标点集中的对应点为q_i，同时在目标点集中寻找对应点为p_i的顶点，组成点集Q_cor，如果Q_cor为空，则p_i的逆一致性误差为无效值，同时也无需对p_i的对应点进行调整；否则，点p_i的逆一致性误差定义为q_i和Q_center之间的欧式距离，其中Q_center为点集Q_cor的重心。同时在这种情况下，把顶点p_i的对应点从点q_i调整为q_i和Q_cor组成的点集的重心。

(4)选出经过调整操作的对应点对，并按照逆一致性误差从大到小的顺序对其进行排序33。

(5)从排序后的点对列表中选出特征点

和特征点

在点集P中的分布满足下列条件：

$p_i^c-p_j^{c}0.5a,i\≠j$

其中a为紧支撑的径向基函数配准算法CSRBF的支撑集大小，定义为a＝3.66Δ，Δ为步骤(4)中所有对应点对在x、y、z轴上的坐标的最大差值。如在全局配准中采用的是TPS算法，则还有附加条件：

$\left|\right|p_i^c-{p^t}_j\left|\right|>0.5a/k,p_i^c\∈P^c\mathrm{and}{p^t}_j\∈P^t$

其中

是TPS配准时在点集P中手工定义的特征点，k是局部配准执行的次数。

根据步骤(2)定义的顶点p_i的逆一致性误差可近视等价于标准定义：||h(p_i)-g^-1(p_i)||，其中h(p_i)≈q_i，因g(q_j)≈p_i，q_j∈Q_cor，所以g^-1(p_i)≈Q_center。通过这种简化定义，避免了复杂的函数求逆运算，因此提高了配置效率。

从步骤(5)得知，局部特征点的数量m取决于Δ的值。Δ值越大，则CSRBF的支撑集a就越大，特征点的数量m就越小。在实验中，为避免大矩阵运算，m的最大值被限定为400。

此外，值得注意的是，虽然局部正变换和反变换中的特征点的序号和对应关系是相同的，但两种变换中的特征点的坐标是不同的。用于局部正变换中的特征点位于变形后的参考点集P’中，

位于原始的目标点集Q中；而用于局部反变换中的特征点

位于原始的参考点集P中，位于变形后的目标点集Q’中。

2)本发明具体采用了Wendland紧支撑径向基函数(Wendland Compact Support RadialBasis Functions，Wendland CSRBF)进行局部变形。其中径向基函数定义为Wendland函数。采用这一函数，每一个特征点在三维空间中的作用域为一个半径可调的球体，而不再是整个点集范围。因此在配准过程中可以只调整没有匹配上的区域，同时保持已经配好的区域。当给定空间维数d，平滑度C^2k(R)以及欧式距离r，Wendland函数ψ_d，k(r)可表示为：

其中：

${(1-r)}_{+}^v=\left\{\begin{array}{cc}{(1-r)}^v& 0\&le;r<1\\ 0& r\&GreaterEqual;1\end{array}\right.$

为截断多项式，而

$\mathrm{I\&psi;}\left(r\right)={\&Integral;}_r^{\&infin;}\mathrm{t\&psi;}\left(t\right)\mathrm{dt},r\&GreaterEqual;0$

为积分运算，在公式(6)中需执行k次。

从公式(6)可以看到，Wendland函数ψ_d，k(r)只在r≤1时有效。其有效范围可以缩放至a，缩放后函数的数学属性保持不变。

ψ_a(r)＝ψ(r/a)

对于3维空间和k＝0，1，2的情况，Wendland函数ψ_d，k(r)可分别表示如下：

${\mathrm{\&psi;}}_{\mathrm{3,0}}\left(r\right)={(1-r)}_{+}^2$

${\mathrm{\&psi;}}_{\mathrm{3,1}}\left(r\right)={(1-r)}_{+}^4(4r+1)$

${\mathrm{\&psi;}}_{\mathrm{3,2}}\left(r\right)={(1-r)}_{+}^6({35r}^2+18r+3)$

在实验中，我们采用ψ_3，1(r)作为CSRBF中的径向基函数，具体公式如下：

$u\left(\stackrel{\&OverBar;}{x}\right)=\stackrel{\&OverBar;}{x}+{\mathrm{\&Sigma;}}_{i=1}^m{\mathrm{\&alpha;}}_i{\mathrm{\&psi;}}_{a,\mathrm{3,1}}\left(\right||p_i^c-\stackrel{\&OverBar;}{x}|\left|\right)---\left(7\right)$

其中

是一个顶点，为顶点

到特征点

之间的欧式距离，m是特征点的个数，α_i(i＝1，...m)是未知权重。同TPS类似，这些权重可通过把参考点集上的特征点

一一映射到目标点集上的对应特征点

求解，如下所示：

$u\left(p_i^c\right)=q_i^c,i=\mathrm{1,2}...m---\left(8\right)$

但由于局部配准中的特征点

和

都是估计得到，可能并不完全对应。因此，本发明增加了一个误差阈值ξ，允许特征点在对应时存在一定的误差：

${\mathrm{\&Sigma;}}_{i=1}^m{\left|\right|u\left(p_i^c\right)-q_i^c\left|\right|}^2\&le;\mathrm{\&xi;}$

写成矩阵形式：

(K+λI)α＝ΔQ    (9)

其中K为m×m的矩阵，矩阵中的元素

I为单位矩阵，

在本发明中直接设置为λ＝0.001a，a是ψ_3，1(r)的支撑集大小，α＝(α₁，α₂...α_m)^T，

$\mathrm{\&Delta;Q}={(q_1^c-p_1^c,q_2^c-p_2^c...q_m^c-p_m^c)}^T.$

权重α_i(i＝1，...m)的值一旦确定，就得到了局部的正变换函数h_li，并可根据公式(7)对参考点集进行局部变形，使之更紧密地匹配到目标点集上35。局部的反变换函数g_li的求解过程类似。

步骤4S4：局部配准可执行多次，直到执行次数超过用户设定的阈值。通常对于顶点数量较小的点集，执行2次局部配准就能得到很好的结果，而对于复杂的点集，一般不超过5次。

步骤5S5：如果大于设定的阀值，配准结束，最近点即为对应点。

将此方法应用于两个颅骨点云模型和两个人脸点云模型的配准，参见图2和图11。这些模型都是从CT图像经分割得到。其中目标颅骨的顶点数为166,757，参考颅骨的顶点数为166,543；目标人脸的顶点数为62,589，参考人脸的定点数为59,901。所有实验都在一台HPxw4400工作站上进行，该工作站具有频率为1.86GHz的Intel(R)Xeon(R)CPU，3G内存。

图3是两个颅骨点云模型在正反两个方向经TPS全局配准后的结果。其中第一列为TPS中采用的两组特征点，第二列为正方向上的配准的结果，第三列是反方向上的配准结果。由于颅骨顶点数过大，如采用点绘制很多细节看不清楚，因此我们把点云模型进行了三角化操作，并采用面绘制的方法绘制其中一个颅骨，然后采用点绘制或面绘制的方法绘制另一个颅骨。从图3可以看出，经全局配准后，在正反两个方向上两个颅骨都能基本上匹配上，但在一些局部区域，如头顶和边缘区域，两个颅骨还存在一些差异。

图4显示了在TPS全局配准的基础上再执行一次CSRBF局部配准后的结果。其中第一列为125对自动生成的用于CSRBF配准的特征点，第二列为正方向上的配准的结果，第三列是反方向上的配准结果。经一次局部调整后，两个颅骨在顶部和背部区域的差异明显减小。

图5显示了最终的配准结果。可以看出，经过五次局部配准后，两个颅骨无论在正方向上，还是反方向上几乎都能完全匹配上。与图3相比，边缘区域部分的配准准确度有很大提高。

图6和7分别显示了经TPS和0次/1次/3次/5次局部配准后目标点集和参考点集的逆一致性误差。如按照3.3.2的逆一致性误差定义，有些顶点的误差会为无效值。因此，我们采用了另一种定义来求参考点集上顶点p_i的误差：||p_i-gоh(p_i)||；同理，目标点集上顶点q_i的误差为：||q_i-hоg(q_i)||。在实验中，误差的单位为毫米。从结果可以看出，随着局部配准的不断执行，误差不断减小，最后误差基本控制在1毫米范围内。

图8是全局配准改用ICP算法后，在正反两个方向，分别经0次、1次和5次CSRBF的配准结果。图9和10是对应的逆一致性误差结果。同样可以看出，误差随着局部配准的次数增加而减低。

图12是图11中显示的两套人脸点云模型经ICP配准后的结果，以及经ICP全局配准和3次CSRBF局部调整后的配准结果。后者的准确度明显比前者高。

此外，对配准耗时情况也有统计。其中对于颅骨模型的配准，采用TPS全局配准和5次CSRBF局部配准共耗时32.5秒；采用ICP全局配准和5次CSRBF局部配准共耗时34.9秒。人脸模型的配准采用ICP全局配准和3次CSRBF局部配准共耗时8.7秒。相对于现有算法由于涉及函数求逆操作，处理两幅分辨率为181×217需近半小时时间，我们算法的效率有很大提高。

本方法的特色在于采用基于全局配准和局部配准的分级方式进行配准，并利用正方向和反方向上对应关系的差异来定义逆一致性误差和用于局部配准的特征点，避免了现有正反互逆配准算法中复杂耗时的函数求逆运算。实现简单、效率高，同时配准精度也高。

上述实验结果表明本方法对拓扑结构复杂、顶点数量繁多的点云模型都能快速得到准确度很高的配准结果，可以用于计算机视觉和模式识别等各应用领域，具有快速准确、操作简单、应用前景广的特点。

最后所应说明的是，以上实施方式仅用以说明本发明的技术方案而非限制。尽管参照实施方式对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (7)

1.一种基于分级的正反互逆的三维稠密点集配准方法，其特征在于：包括下列步骤：

步骤一(S1)：输入两个点集，一个参考集，一个目标集；

步骤二(S2)：全局配准；

步骤三(S3)：局部配准；

步骤四(S4)：判断局部配准次数；

步骤五(S5)：配准结束，最近点即为对应点；

其中，步骤二(S2)所述的全局配准可采用两种方法，薄板样条函数算法(TPS)(S200)或迭代最近点算法(ICP)(S20)；

其中，步骤三(S3)所述的局部配准方法包含六个子步骤：

1)求变形后的参考点集上的每一点到目标点集上的最近点，两点之间的距离小于用户设定的阈值，则以该最近点作为它的对应点，最终得到一个对应点对列表(30)；

2)求变形后的目标点集上的每一点到参考点集上的最近点，两点之间的距离小于用户设定的阈值，则以该最近点作为它的对应点，最终得到另一个对应点对列表(31)；

3)根据前两个步骤得到的结果，计算参考点集和目标点集中存在对应点的每一个顶点的逆一致性误差，并对其对应点进行调整(32)；

4)把两个对应点对列表中的经过调整操作的对应点对挑选出来，存成一个新列表，并根据逆一致性误差从大到小的顺序对新列表进行排序(33)；

5)从排序后的列表中提取用于局部配准的两组特征点(34)；

6)采用具有紧支撑的径向基函数CSRBF，并根据上一个步骤得到的特征点计算正变换h_li和反变换g_li，并把正反变换分别作用于已变形的参考点集P’和目标点集Q’，得到两个新的变形点集，依然采用P’和Q’表示(35)；

其中，步骤四(S4)所述的判断局部配准次数小于预先设定的阀值的，重复执行步骤三(S3)。

2.根据权利要求1所述一种基于分级的正反互逆的三维稠密点集配准方法，其特征在于：步骤二(S2)所述的薄板样条函数算法(TPS)(S200)全局配准法包含如下两个子步骤：

1)从两个点集中采用人机交互的方式选取两组对应特征点(200)；

2)利用TPS算法和特征点，计算从参考点集到目标点集的变换，以及从目标点集到参考点集的变换，并把正反变换分别作用于参考点集和目标点集，得到变形后的两个点集(201-202)。

3.根据权利要求1所述一种基于分级的正反互逆的三维稠密点集配准方法，其特征在于：步骤二(S2)所述的迭代最近点算法(ICP)(S20)也包含如下两个子步骤：

1)判断参考点集和目标点集中的顶点数量，顶点数量大于设定的阀值，则从两个点集中随机选取两组顶点，组成两个新的点集，新的点集中各顶点间的距离大于用户设定的距离阈值，且顶点数量满足数量阈值要求(20-22)；

2)利用ICP算法计算两个新点集的正变换和反变换，并把正反变换分别作用于参考点集和目标点集，得到变形后的两个点集(23)。

4.根据权利要求1所述的一种基于分级的正反互逆的三维稠密点集配准方法，其特征在于：所述步骤三(S3)中的子步骤1)(30)和2)(31)中采用kd树结构加快最近点求解。

5.根据权利要求1所述的一种基于分级的正反互逆的三维稠密点集配准方法，其特征在于：所述步骤三(S3)中的子步骤3)(32)中顶点的逆一致性误差的定义分两种情况：

情况1：参考点集中顶点p_i在目标点集中的对应点为顶点q_i，同时目标点集中有一个或几个顶点，记作点集Q_cor，其在参考点集中的对应点为顶点p_i，此时顶点p_i的逆一致性误差定义为q_i和Q_center之间的欧式距离，其中Q_center为点集Q_cor的重心；

情况2：参考点集中顶点p_i在目标点集中的对应点为顶点q_i，但目标点集中没有顶点的对应点为顶点p_i，此时顶点p_i的逆一致性误差定义为无效值；

上述定义也适应目标点集上的顶点。

6.根据权利要求1所述的一种基于分级的正反互逆的三维稠密点集配准方法，其特征在于：所述步骤三(S3)中的子步骤3)(32)步骤中对应点的调整只在一种情况下执行：

即当参考点集中顶点p_i在目标点集中的对应点为顶点q_i，同时目标点集中有一个或几个顶点，记作点集Q_cor，其在参考点集中的对应点为顶点p_i，此时顶点p_i的对应点调整为q_i和Q_cor组成的点集的重心；

上述对应点的调整情况也适应目标点集。

7.根据权利要求1所述的一种基于分级的正反互逆的三维稠密点集配准方法，其特征在于：所述步骤三(S3)中的子步骤5)(34)中从参考点集中提取出来的特征点之间的距离大于0.5a，其中a为紧支撑径向基函数CSRBF的紧支撑集大小，定义为a＝3.66Δ，Δ为列表中对应点对在x、y、z轴上的坐标的最大差值。

Images