基于全局与局部相似性测度的非刚性点集配准

彭磊 杨秀云 张裕飞 李光耀

摘 要：非刚性点集配准算法中，能否找到正确的对应关系对配准结果起着至关重要的作用，而通常两个点集中的对应点除了距离比较接近之外还具有相似的邻域结构，因此提出基于全局与局部相似性测度的非刚性点集配准算法。首先，使用一致性点漂移（CPD）算法作为配准框架，采用高斯混合模型对点集进行建模。然后，对全局局部混合距离进行改进，形成全局与局部相似性测度准则。最后，采用期望最大化（EM）算法迭代地求解对应关系和变换公式：在迭代初期局部相似性所占比重较大，从而能够尽快地找到正确的对应关系;随着迭代的进展全局相似性比重逐渐增大，从而确保得到较小的配准误差。实验结果表明，与薄板样条鲁棒点匹配（TPS-RPM）算法、高斯混合模型点集配准（GMMREG）算法、基于L2E估计的鲁棒点匹配算法（RPM-L2E）、基于全局局部混合距离与薄板样条的点集配准算法（GLMDTPS）和CPD算法相比，所提算法的均方根误差（RMSE）分别下降了39.93%、42.45%、32.51%、22.36%和11.76%，说明该算法具有较好的配准效果。

关键词：点集配准;图像配准;非刚性配准;高斯混合模型;相似性测度

中图分类号：TP391.41

文献标志码：A

Abstract： In the non-rigid point set registration algorithm， whether the correct correspondence can be found plays an important role. Generally the corresponding points in two point sets have similar neighborhood structures besides the close distance. Therefore， a non-rigid point set registration algorithm based on global and local similarity measurement was proposed. Firstly， the Coherent Point Drift （CPD） algorithm was used as the registration framework， and the Gaussian mixture model was used to model the point sets.Secondly， the global and local mixture distancewas improved to form the global and local similarity measurement criterion.Finally， the correspondence and the transformation formulawere solved by the Expectation Maximization （EM） algorithm. In the initial stage of the iteration， the proportion of local similarity was larger so that the correct correspondence was able to be found rapidly; with the progress of the iteration， the proportion of global similarity was increased to ensure the smaller registration error. Experimental  results show that compared with the Thin Plate Spline Robust Point Matching （TPS-RPM） algorithm，the Gaussian Mixture Models point set REGistration （GMMREG） algorithm，the Robust Point Matching algorithm based on L2E estimation （RPM-L2E）， the Global and Local Mixture Distance and Thin Plate Spline based point set registration algorithm （GLMDTPS） and the CPD algorithm， the proposed algorithm has the Root Mean Squared Error （RMSE） decreased by 39.93%， 42.45%， 32.51%， 22.36% and 11.76% respectively， indicating the proposed algorithm has better registration performance.Key words：  point set registration; image registration; non-rigid registration; Gaussian mixture model; similarity measurement

0 引言

點集配准是一个基础而关键的问题，广泛应用在计算机视觉[1]、医学图像分析[2-4]、遥感图像处理[5]和模式识别[6]等领域。其目标是寻找两组给定点集之间的对应关系，通过空间变换将一个点集对齐到另一个点集。根据变换方式的不同可分为刚性配准和非刚性配准。刚性点集配准的变换方式包括放缩、平移、旋转，仅涉及少量的参数。非刚性点集配准的变换方式往往是未知的、复杂的、难以建模的，通常包含较多的变换参数。

迭代最近点算法[7]是著名的点集配准方法之一，主要用于自由曲线和曲面的配准，它假设两个点集间距离最近的点为对应点，通过迭代地最小化均方距离来达到目标函数的收敛。Chui等[8]提出了一个非刚性点匹配的通用框架，使用薄板样条作为空间变换模型，并通过软分配和确定性退火技术联合求解优化问题。一致性点漂移（Coherent Point Drift， CPD）算法[9]将配准两个点集考虑为概率密度估计问题，使用运动一致性理论进行约束，通过迭代方法将高斯混合模型（Gaussian Mixture Model， GMM）[9-10]的质心拟合到数据点。Jian等[11]使用两个GMM来表示点集，并通过两个GMM之间的差异最小化来解决点集配准问题。文献[12-13]采用高斯域的方法实现了异常点移除和点集配准。这些方法主要使用点之间的全局距离来确定对应关系，基本没有考虑点的邻域结构信息。

Belongie等[6]引入了形状上下文（Shape Context，SC）的概念，用来测量两个形状之间的相似度。Ma等[14-15]首先使用SC描述符建立粗糙的对应关系，然后用L2E估计器来估算变换。Yang等[16]采用鲁棒的全局和局部混合距离（Global and Local Mixture Distance，GLMD）测度准则进行非刚性点集配准。虽然这些算法使用点集的结构信息作为相似性测量，但是它们在查找对应关系和估计变换模型时是分别处理的。

为了充分利用点之间的距离和点的邻域结构信息，快速找到正确的对应关系，本文以CPD算法[9]为框架，采用GMM表示点集，对GLMD[16]进行改进，并融入到GMM中，形成了全局与局部相似性测度准则。GMM的高斯分量反映了全局相似性测度准则，而每个高斯分量使用不同的比例系数，则反映了点的局部相似性测度准则。通过权重系数调节二者的比重。采用期望最大化（Expectation Maximization， EM）算法迭代地求解对应关系和变换公式。在迭代初期，局部相似性所占比重较大，从而尽快找到正确的对应关系;迭代过程中逐渐提高全局相似性比重，最终得到更小的配准误差。该算法在保证快速找到点的正确对应关系的同时，能够得到较好的配准效果。

1 一致性点漂移算法

设D维（通常D=2或3）空间的两个点集X和Y，浮动点集为X={xm|m=1，2，…，M}，参考点集为Y={yn|n=1，2，…，N}。非刚性点集配准算法的目标是根据相似性测度准则找到X和Y之间点的对应关系，并估算一个非刚性变换公式T，使X通过变换后对齐到Y。

CPD算法将两个点集的配准看成是概率密度估计问题。使用GMM对问题进行建模[9]，把X中的点看成GMM的质心，Y中的点看成是由GMM产生的数据点。方法的核心是保持点集的拓扑结构，迫使GMM的质心作为一个整体一致地向参考点集移动。在优化过程中，浮动点集通过变换T与参考点集逐渐对齐，并可以根据后验概率获得点的对应关系。GMM的概率密度函数为：

式（1）中的均匀分布p（yn|M+1）=1/N，用来处理点集中存在的噪声和异常点，所占比重为ω（0≤ω≤1）。其他所有高斯组成部分所占比重为1-ω，其中每个高斯分量使用相同的比例系数ηmn=1/M。

为了使變换比较平滑，采用Tikhonov正则化框架[9]。定义再生核希尔伯特空间的正则项λφ（T）/2，其中φ（T）为一个平滑函数，λ为正则项的权重系数。使用θ表示参数集合，通过极大似然估计来求得参数。将负对数似然函数作为CPD算法的能量函数，如式（3）所示：

2 全局和局部混合距离

2.1 全局距离

浮动点集上一点xm与参考点集上一点yn之间的全局距离[16]定义为其欧氏距离的平方，如式（4）所示：

Gmn可以描述两个点集之间的全局性结构差异。通常两个点集基本对齐后，距离较近的点很可能就是对应点。

2.2 局部距离

某个点q的邻域定义为其所在点集中距离q最近的K个点。用φ（q）k来表示第k个最近点。局部距离[16]用于度量两个点xm与yn的邻域结构相似性，定义为：

其中ψ是位移函数，表示点xm及其邻域点整体向点yn移动直至xm与yn点重合，其定义为：

3 本文算法

3.1 局部相似性的改进

GLMD算法中局部距离是通过计算两个点的邻域中具有相同编号k的点对之间距离的和得到的。本文将式（5）中φ（yn）的下标由k改为f（k），即寻找两个邻域中两两点对之间距离之和最小的值作为局部相似性。式（5）改写为：

f为两个邻域中点的某种一一映射关系，使得两个邻域（分别包含K个点）之间K个点对的距离之和达到最小，即能够使Lmn得到最小值。相当于二分图的最大匹配问题，可以通过匈牙利算法来求解f。

如图1所示，两个点p和q，各自的邻域内分别有编号为1、2、3、4的点，虚线表示两个邻域之间点的对应关系。通过计算对应点间距离之和得到局部距离。GLMD按照编号确定对应点，有可能产生错误的结果，从而得到错误的局部距离，如图1（a）所示。本文算法能够找到一个最佳匹配来确定对应点，因此能够求得正确的局部距离，同时所求局部距离也是最小的，如图1（b）所示。

3.2 全局局部相似性测度

CPD算法利用了点之间的全局距离，即认为全局距离较近的点很可能为对应点。因此，式（2）可以写为：

CPD算法中每个高斯分量使用相同的比例系数ηmn=1/M，没有考虑点的邻域结构。然而两个点的邻域结构越相似，即局部相似性越高，它们成为对应点的可能性也越大。为了充分利用点的局部相似性，对GMM中的每个高斯分量使用不同的比例系数：

其中：S=∑Mi=1 exp（-βLin）， β为权重系数。两个点的局部相似性越大时，对应的ηmn越大，即该高斯分量在GMM中所占的比重越大，其成为对应点的概率就越大。

式（8）和（9）共同体现了全局与局部相似性结合的测度准则。权重系数β控制全局相似性和局部相似性之间的比例关系。当β非常大时，相当于最小化局部相似性Lmn;当β很小时，趋向于最小化全局相似性Gmn。全局与局部相似性测度提供了一种灵活的方式，通过最小化两个点集之间的全局相似性和局部相似性上的差异来估算对应关系。

3.3 全局与局部测度引导的配准过程

在非刚性点集配准中，能否找到正确的对应关系对配准效果有很大的影响。当两个点集形态近似时，距离较近的点，往往是对应点;而当两个点集形态差别较大时，具有相似邻域结构的点，很有可能是对应点。因此采用全局与局部相似性结合的测度准则引导配准过程。在优化算法迭代的初期局部相似性起较大的作用，从而尽快找到正确的对应关系;在迭代后期，两个点集基本对齐，此时全局相似性测度起的作用变大，从而确保算法收敛到较小的配准误差。

配准过程首先以浮动点集X和参考点集Y为输入;其次使用GMM进行建模，把X看成GMM的质心，Y看成由GMM产生的数据点;再次构造目标函数并进行参数设置;然后使用EM优化策略对目标函数中的参数进行估计，交替地执行期望步（Expectation step，E-step）和最大化步（Maximization step， M-step）两步，直至算法收敛;最后输出浮动点集的变换公式和两个点集之间的对应关系。本文算法流程如图2所示。

3.3.2 参数设置

优化算法最大迭代次数设为100。邻域中点的个数设为K=4。初始时权重系数β取一个较大的值，一般设β=K2，每一次迭代均乘以一个系数r=0.95，类似于退火率，即β←rβ。这样能够保证优化迭代的前期局部相似性起主导作用，后期全局相似性起主导作用，从而在尽快找到正确对应关系的同时得到较小的配准误差。

3.3.3 E-step阶段

在E-step阶段，首先计算Y与T（X）之间的局部相似性，更新每个高斯分量的比例系数ηmn;然后根据贝叶斯定理，使用当前参数计算后验概率分布Pmn。

3.3.4 M-step阶段

在M-step阶段，计算并更新参数。对式（10）求导，并使其等于0，可以求得相应参数。

非刚性变换公式T定义为点的初始位置加上一个位移函数v，即T（X）=X+v（X）[9]。使用高斯矩阵核定义再生核希尔伯特空间，可以得到v的函数形式：

其中：wm是系数矩阵W=[w1 w2 … wM]T中的元素;Γ（x，xm）是核矩阵。正则项φ（T）在这里的具体形式为WΓWT[14]。则变换公式可表示为T=X+ΓW。同时求得参数ω和σ2的值。

3.3.5 算法输出

EM算法交替地执行E-step和M-step，直至最大迭代次数或配准误差小于给定阈值。最终浮动点集的变换公式由T给出，点的对应关系由后验概率矩阵P给出。

4 实验与分析

为了评估算法的性能，设计了二维与三维合成点集配准实验和真实医学图像配准实验，并对实验结果进行了比较和定量分析。所有實验的软件环境为Matlab R2013a，硬件环境为16GB内存、i7-4770K Inter处理器的PC。

4.1 合成点集配准实验

采用二维的Chui-Rangarajan合成数据[8]和三维的人脸点集[9]数据进行实验。且实验结果与薄板样条鲁棒点匹配（Thin Plate Spline Robust Point Matching， TPS-RPM）算法[8]、高斯混合模型点集配准（Gaussian Mixture Models point set REGistration， GMMREG）算法[11]、基于L2E估计的鲁棒点匹配算法（Robust Point Matching algorithm based on L2E estimation， RPM-L2E）[14-15]、基于全局局部混合距离与薄板样条的点集配准算法（Global and Local Mixture Distance and Thin Plate Spline based point set registration algorithm， GLMDTPS）[16]和CPD算法[9]进行比较。性能评价采用与文献[8，10，12-15]相同的方法。由于合成数据中每组测试样例的真正对应点已经预先标识出，即真正的对应关系是已知的，因此使用配准后两个点集的真正对应点之间的平均欧氏距离，即均方根误差（Root Mean Squared Error， RMSE）进行定量分析。误差计算公式为：

4.1.1 二维点集配准实验

选取Chui-Rangarajan合成数据小鱼形状的变形、遮挡、异常点3组数据进行实验。每组包括由低到高5种不同的退化等级，每一等级包含100个测试样例。图3为6种算法的配准结果展示。

图4是配准性能的定量分析与比较，用误差棒表示每种退化等级样例的平均误差和标准差。

图4中横坐标表示5种不同的退化等级。

每种算法所对应的误差棒的中心表示平均误差，误差棒的长度表示标准差。从图4可以看出，当变形、遮挡、异常点等退化程度较小时，各种算法均有不错的性能，但随着退化程度的增大，本文算法的优势越来越明显。表1为6种算法3组实验的平均误差统计，可看出本文算法在各种情况下都有一定的优势。

4.1.2 三维点集配准实验

三维点集配准实验选取文献[9]中的人脸点集。变形、遮挡和异常点3组数据，每组选取50个测试样例。本文算法的配准效果如图5所示。

在二維和三维点集上的实验结果表明，本文算法与TPS-RPM、GMMREG、RPM-L2E、GLMDTPS和CPD算法相比，平均RMSE分别下降了39.93%、42.45%、32.51%、22.36%和11.76%。

4.2 真实图像配准实验

为了比较在CPD配准框架中引入全局与局部相似性测度前、后的配准效果，选用20组医学MRI图像进行实验。首先提取图像的轮廓，生成点集，对点集进行配准;然后根据变换公式生成配准后的图像。图6为一组人脑MRI图像的测试样例，展示了CPD算法与本文算法的配准结果。图6第一行为初始的浮动图像与参考图像，以及两种算法配准后浮动图像与参考图像的灰度差图像;第二行为对应的点集。

配准效果的定量分析采用轮廓点集配准误差和图像灰度误差两种方法。

首先采用类似文献[17]人工标注点的方法，在两个轮廓点集上查找具有明显特征的点对，人工标注50个标准对应点。计算两个轮廓对应标注点之间的RMSE。两种算法的平均误差如表3中人工标注点的平均误差所示。召回率定义为配准后误差小于给定阈值的对应点与全部标注的对应点在数量上的比值。图7展示了20组测试样例在各种阈值下CPD与本文算法召回率的变化曲线。实验结果表明，本文算法优于CPD算法。

其中：Y（u，v）与X（u，v）是两个图像在坐标点（u，v）上的灰度值;U×V是图像的分辨率。表3中像素灰度的均方根误差展示了20组测试样例使用两种算法配准后的图像灰度误差，同样可以看出本文算法优于CPD算法。

5 结语

在非刚性点集配准算法中，为了充分利用点的邻域结构信息来查找对应关系，将改进的全局和局部混合距离融入到CPD算法的高斯混合模型配准框架中，提出了基于全局与局部相似性的测度准则。使用权重系数调节全局相似性测度和局部相似性测度二者的比例。采用EM算法迭代地求解对应关系和变换公式，迭代过程中局部相似性所占比重由大变小，全局距离所占比重由小变大，从而在保证快速找到正确对应关系的前提下，提高配准的精确度。通过合成数据与真实数据的实验结果比对和分析，可以看出本文算法具有较好的配准性能。然而将点集邻域结构的计算引入到相似性测度中，加大了算法的计算量、增加了算法的运行时间，可能无法满足实时性要求较高的应用场合。因此未来的工作在进一步研究如何更好地利用点集结构信息的同时，还需要考虑如何提高算法的执行速度。

参考文献（References）

[1] MA J， ZHAO J， JIANG J， et al. Locality preserving matching[J]. International Journal of Computer Vision， 2019， 127（5）： 512-531.

[2] FERRANTE E， PARAGIOS N. Slice-to-volume medical image registration： a survey[J]. Medical Image Analysis， 2017， 39： 101-123.

[3] GUAN S， WANG T， MENG C， et al. A review of point feature based medical image registration[J]. Chinese Journal of Mechanical Engineering， 2018， 31： No.76.

[4] TANG H， PAN A， YANG Y， et al. Retinal image registration based on robust non-rigid point matching method[J]. Journal of Medical Imaging and Health Informatics， 2018， 8（2）： 240-249.

[5] WU Y， MA W， ZHANG J， et al. Point-matching algorithm based on local neighborhood information for remote sensing image registration[J]. Journal of Applied Remote Sensing， 2018， 12（1）： No.016002.

[6] BELONGIE S， MALIK J， PUZICHA J. Shape matching and object recognition using shape contexts[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2002， 24（4）： 509-522.

[7] BESL P J， McKAY N D. Method for registration of 3-D shapes[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 1992， 14（2）： 239-256.

[8] CHUI H， RANGARAJAN A. A new point matching algorithm for non-rigid registration[J]. Computer Vision and Image Understanding， 2003， 89（2/3）： 114-141.

[9] MYRONENKO A， SONG X. Point set registration： coherent point drift[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2010， 32（12）： 2262-2275.

[10] WANG G， CHEN Y. Fuzzy correspondences guided Gaussian mixture model for point set registration[J]. Knowledge-Based Systems， 2017， 136： 200-209.

[11] JIAN B， VEMURI B C. Robust point set registration using Gaussian mixture models[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2011， 33（8）： 1633-1645.

[12] WANG G， ZHOU Q， CHEN Y. Robust non-rigid point set registration using spatially constrained Gaussian fields[J]. IEEE Transactions on Image Processing， 2017， 26（4）： 1759-1769.

[13] WANG G， CHEN Y， ZHENG X. Gaussian field consensus： a robust nonparametric matching method for outlier rejection[J]. Pattern Recognition， 2018， 74： 305-316.

[14] MA J， QIU W， ZHAO J， et al. Robust L2E estimation of transformation for non-rigid registration[J]. IEEE Transactions on Signal Processing， 2015， 63（5）： 1115-1129.

[15] MA J， ZHAO J， TIAN J， et al. Robust estimation of nonrigid transformation for point set registration[C]// Proceedings of the IEEE 2013 Conference on Computer Vision and Pattern Recognition. Piscataway： IEEE， 2013： 2147-2154.

[16] YANG Y， ONG S H， FOONG K W C. A robust global and local mixture distance based non-rigid point set registration[J]. Pattern Recognition， 2015， 48（1）： 156-173.

[17] 王麗芳， 王雁丽， 蔺素珍， 等. 基于改进的Zernike矩的局部描述符与图割离散优化的非刚性多模态脑部图像配准[J]. 计算机应用， 2019， 39（2）： 582-588. （WANG L F， WANG Y L， LIN S Z， et al. Non-rigid multi-modal brain image registration based on improved Zernike moments based local descriptor and graph cuts discrete optimization[J]. Journal of Computer Applications， 2019， 39（2）： 582-588.）