点集配准技术（ICP、RPM、KC、CPD）

点集配准技术（ICP、RPM、KC、CPD）
在计算机视觉和模式识别中，点集配准技术是查找将两个点集对齐的空间变换过程。

寻找这种变换的⽬的主要包括：1、将多个数据集合并为⼀个全局统⼀的模型；2、将未知的数据集映射到已知的数据集上以识别其特征或估计其姿态。

点集的获取可以是来⾃于3D扫描仪或测距仪的原始数据，在图像处理和图像配准中，点集也可以是通过从图像中提取获得的⼀组特征（例如⾓点检测）。

点集配准研究的问题可以概括如下：假设{M，S}是空间R d中的两个点集，我们要寻找⼀种变换T，或者说是⼀种从R d空间到R d空间的映射，将其作⽤于点集M后，可以使得变换后的点集M和点集S之间的差异最⼩。

将变换后的点集M记为T(M)，那么转换后的点集T(M)与点集S的差异可以由某种距离函数来定义，⼀种最简单的⽅法是对配对点集取欧式距离的平⽅：
点集配准⽅法⼀般分为刚性配准和⾮刚性配准。

刚性配准：给定两个点集，刚性配准产⽣⼀个刚性变换，该变换将⼀个点集映射到另⼀个点集。

刚性变换定义为不改变任何两点之间距离的变换，⼀般这种转换只包括平移和旋转。

⾮刚性配准：给定两个点集，⾮刚性配准产⽣⼀个⾮刚性变换，该变换将⼀个点集映射到另⼀个点集。

⾮刚性变换包括仿射变换，例如缩放和剪切等，也可以涉及其他⾮线性变换。

下⾯我们来具体介绍⼏种点集配准技术。

⼀. Iterative closest point（ICP）
ICP算法是⼀种迭代⽅式的刚性配准算法，它为点集M中每个点m i寻找在点集S中的最近点s j，然后利⽤最⼩⼆乘⽅式得到变换T，算法伪代码如下：
ICP算法对于待配对点集的初始位置⽐较敏感，当点集M的初始位置与点集S⽐较接近时，配准效果会⽐较好。

另外ICP算法在每次计算迭代过程中都会改变最近点对，所以其实很难证明ICP算法能准确收敛到局部最优值，但是由于ICP算法直观易懂且易于实现，因此它⽬前仍然是最常⽤的点集配准算法。

在计算变换T时可以使⽤SVD⽅式，过程如下：
优化⽬标：
1.将点集去中⼼化：
2.计算协⽅差矩阵并SVD分解：
3.计算旋转矩阵和平移向量：
⼆. Robust point matching（RPM）
RPM算法是⼀种使⽤退⽕和软对应⽅式的配准算法。

ICP算法在迭代计算时利⽤距离最近原则来产⽣待配准点对，⽽RPM算法利⽤软对应⽅式为任意点对赋予0到1之间的值，并最终收敛到0或1，如果是1则代表这两个点是配准点对。

RPM算法最终计算得到的配准点对是⼀⼀映射的，⽽ICP算法通常不是。

RPM算法定义配准的损失函数如下：
其中：
t是平移向量。

A是变换矩阵，如果是2D情况，矩阵A可以分解成4个变量{a，θ，b，c}，分别代表缩放系数，旋转⾓度，垂直和⽔平剪切量。

g(A)代表变换矩阵的正则项，其作⽤是防⽌变换矩阵中出现较⼤值。

µ是配准矩阵，定义如下：
α项作⽤是是配准矩阵中有更多的配准点对数量
算法伪代码如下：
三. Kernel correlation（KC）
KC算法是⼀种相似性测量⽅法，它将点集配准问题转化为寻找待配准点集之间最⼤相关性的过程。

对于两个点x i和x j，它们之间的核相关性（kernel correlation）定义为：
其中K(x, x i)是中⼼位置为x i的核函数，核函数通常为⾮负对称函数，常⽤的核函数有⾼斯核函数等。

⾼斯核函数形式如下：
由⾼斯核函数可以得到的核相关性如下：
可以看到核相关性是与距离||x i- x j||相关的函数，值得注意的是，由其他核函数得到的核相关性也是距离||x i- x j||的函数。

KC算法定义点集配准的损失函数如下：
定义核密度估计（kernel density estimates）：
损失函数可以表⽰为两个核密度估计之间的相关性：
四. Coherent point drift（CPD）
CPD算法将点集配准问题转换为概率密度估计问题，其将点集M的分布表⽰成混合⾼斯模型，当点集M与点集S完成配准后，对应的似然函数达到最⼤。

将点集M中每个点m i作为混合⾼斯模型中每个成分的中⼼，并且假设每个成分概率相等，那么其分布模型可以表⽰为：
其中：
考虑噪声的影响，在分布模型中⼜加⼊了均匀分布函数，其权重为w，这样上式分布模型可以进⼀步表⽰为：
混合⾼斯模型的中⼼在配准过程中与变换参数θ相关，为了求得模型参数，需要极⼩化负对数似然函数：
求解θ和σ2过程可以使⽤期望最⼤算法，算法包含两步：
E-Step：通过旧的分布模型参数计算后验概率分布P old(i|s j)
M-Step：极⼩化损失函数来得到新的模型参数θ和σ2
算法伪代码如下：
在M-step中，对于刚性变换和⾮刚性变换，变换模型参数求解分别如下：
⽆异常点⽆噪声情况
有异常点⽆噪声情况
⽆异常点有噪声情况
有异常点有噪声情况
本⽂为原创，转载请注明出处：。

参考⽂献：
[1]
[2] https://igl.ethz.ch/projects/ARAP/svd_rot.pdf
[3] Gold S, Rangarajan A, Lu C, et al. New algorithms for 2D and 3D point matching : Pose estimation and correspondence. Pattern Recognition, 1998, 31(8): 1019-1031.
[4] Chui H, Rangarajan A. A new point matching algorithm for non-rigid registration. Computer Vision and Image Understanding, 2003, 89(2): 114-141.
[5] Tsin Y, Kanade T. A Correlation-Based Approach to Robust Point Set Registration. european conference on computer vision, 2004: 558-569
[6] Myronenko A, Song X B. Point Set Registration: Coherent Point Drift. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2010, 32(12): 2262-2275.
附录
薄板样条（Thin Plate Spline）
薄板样条可以表⽰成径向基函数的形式，给定⼀系列控制点{c i, i = 1, ..., n}，径向基函数定义了空间中从位置x到f(x)之间的映射关系：
薄板样条对应的径向基函数核为σ(r) = r2log r。

对于⾮刚性变换，如果使⽤薄板样条来表⽰，那么径向基函数的⾃变量和应变量可以分别认为是变换前的点坐标x和变换后的点坐标y。

下⾯具体介绍如何使⽤薄板样条表⽰⾮刚性变换。

假设⾃变量是2维空间中的⼀系列点集x（A图所⽰），应变量也是2维空间中的⼀系列点集y（B图所⽰），我们对应变量的每个维度进⾏薄板样条插值，插值函数形式如下：
其中：c是标量，a ∈ R2×1，w ∈ R n×1，s(x)表达式式如下：
插值函数表达式中c + a T x代表仿射变换，w T s(x)代表⾮仿射变换。

对于变换后每个点的第⼀个维度，我们可以建⽴⼀个约束条件：
然后额外添加如下约束条件：
最终根据以上约束条件就可以求得薄板样条插值函数的参数了。

期望最⼤算法（Expectation Maximization）
期望最⼤算法是解决含隐变量模型参数的⼀种⽅法。

给定数据样本X = {x(1), ... , x(n)}，每个观察样本x(i)还对应着⼀个隐变量z(i)。

假设样本为独⽴同分布，为了拟合模型p(x; θ)参数，我们需要极⼤化对数似然估计函数：
对于每个样本x(i)，我们引⼊隐变量z的某种分布Q i(z)，那么上式可以进⼀步变化：
上式最后⼀步利⽤了Jensen 不等式。

于是我们得到似然函数L(θ)的下界函数，那么我们可以不断最⼤化下界函数G(θ, Q)来使得L(θ)不断提⾼，算法过程如下：
1.（E-Step）固定θ，调整Q(z)，使得似然函数L(θ)与下界函数G(θ, Q)在θ处相等
2.（M-Step）固定Q(z)，调整θ，使得下界函数G(θ, Q)达到最⼤值。

3.迭代上述步骤，直⾄收敛到局部最优解。

对于第1步，只有当为常数时等号成⽴，⼜由于，所以得到了Q(z)的表达式：
对于第2步，下界函数G(θ, Q)可以继续推导，得到⼀个更简单的表达形式：
下⾯具体介绍以期望最⼤算法来求解混合⾼斯模型的参数。

如上图所⽰，给定⼀系列⼆维数据点x，假设这些数据是由3个正态分布函数随机取样得到，这些正态分布函数线性组合在⼀起就组成了混合⾼斯模型的概率密度函数：
那么得到对数似然估计函数如下：
我们为每个数据点x i引⼊⼀个隐变量z，这个隐变量代表这个数据点x i是由混合⾼斯模型中的哪个正态分布函数⽣成得到。

根据前⾯期望最⼤算法步骤：
1.计算隐变量的分布Q i(z k)
Q i(z k)它代表了数据点x i是由第k个正态分布函数⽣成的概率：
2.计算混合⾼斯模型参数πk，µk，Σk 下界函数进⼀步推导：
于是可以求出模型参数：。