Opencv2.4.9源码分析——Stitching（二）

2、计算单应矩阵

2.1 原理

在得到了图像特征点以后，我们就可以根据这些特征点，实现图像匹配，即得到重叠区域。而要把多幅图像拼接成一幅图像，就需要以某幅图像为基准，把其他图像映射到该图像所在的平面。映射的关键所在就是根据重叠区域的特征点计算图像间的单应矩阵。

我们可以通过最近邻算法（K-NN，这里的K表示最近邻特征点的数量）得到与图像A中的某一特征点a最相似的其他任意一幅图像B的K个特征点，对这K个特征点进行距离比较来最终判断是否找到了a的最相似特征点，即配对成功与否。例如我们这里选择K为2，即2-NN算法，图像B中的与a最相似的两个特征点为b和c，设L_ab和L_ac分别表示a与b和c的距离（距离是通过特征点描述符的比较实现的），并且L_ab＜L_ac，如果满足下式，则我们认为a的匹配特征点为b：

（1）

式中，ρ表示匹配阈值。

那么如何得到这K个特征点呢？最简单、也是最粗暴的方法就是计算a与图像B中的所有特征点的距离，这种方法适用于特征点不多的情况。当特征点较多，并且需要匹配的图像也较多时，应用K-D树（K-维树）的算法会更有效。

K-D树是一棵二叉树，树中存储的是一些K维数据，这正好对应于特征点的K维描述符。K-D树本质上是对该K维数据集合的一种K维空间的划分，即树中的每个叶节点都对应于一个K维超矩形区域。当K-D树构建好后，可以很迅速的得到某一数据所属的超矩形区域。所以K-D树非常适合范围搜索和最近邻搜索，但它不适合维数太多的数据搜索。另外，需要注意的是K-NN中的K和K-D树中的K表示的含义是不同的。

首先我们需要利用已知数据构建K-D树（这里的数据指的是图像B中的所有特征点，而特征点是用描述符表示的）。要想构建K-D树，需要解决两个问题：第一是选择哪个维度进行分割；第二是如何把数据分割为左分支和右分支。对于第一个问题，通常的做法是交替循环的依次选择所有的维度，例如在3-D树中，“父节点”用的x维，则“子节点”用y维，“孙节点”用z维，“重孙节点”又回到x维，“重重孙节点”用y维，以此类推。对于第二问题，一般的做法是选择该维度的中值进行分割，小于该值的数据归为左分支，大于该值的数据归为右分支。以上方法可以保证K-D树是平衡的，即每个叶节点与根节点的距离基本相同。

当K-D树构建好了以后，我们就可以搜索待查询数据（这里指的是图像A中的特征点a）的最近邻数据，它的过程是：从根节点出发，依据分叉规则（即选择分割维度和比较中值），沿着搜索路径达到最佳叶节点；然后再从该叶节点出发，回溯搜索路径，找到是否有比最佳叶节点更好的节点，用以替代它。从概念上来说，回溯的过程就是以待查询点为中心，以待查询点到最佳叶节点的距离为半径，画一个超球体，如果该超球体与分割维度所在的超平面相交，则需要搜索该超平面的另一侧（即另一个分支），如果不相交，则无需再搜索另一侧。

特征点的匹配点对得到后，我们就可以利用这些匹配点对评估计算图像间的单应矩阵H。

如果平面到平面之间的映射关系是投影映射的话，那么一定存在一个非奇异的3×3的矩阵H，当p表示一个平面的任意一点时，那么该点映射到另一平面的点就为Hp。当相机围绕一点进行旋转而得到不同的两幅图像时（如图3所示），这两幅图像之间就是投影映射，它们之间的关系就可以用矩阵H来进行描述，即图像A中的点通过矩阵H就可以映射到图像B中，从而实现两幅图像的拼接。

图3 通过旋转相机，得到不同的图像

需要指出的是，当矩阵H乘以任意一个不为零的常数时，并不改变两个平面之间的投影映射关系，因此H被认为是只有8个自由度的单应矩阵。

设图像A中的某一点的齐次坐标p为(u, v,1)，该点在图像B中的匹配点的齐次坐标p’为(u’, v’,1)，则：

（2）

即

（3）

图像坐标(x, y)所对应的齐次坐标为(x, y, 1)，而齐次坐标(x, y, z)所对应的图像坐标为(x/z, y/z)。

我们把式3展开，则：

（4）

把式4写成矩阵的形式：

（5）

我们观察式5等号左侧的第一个矩阵会发现，第三行是前两行的线性组合，即

（6）

因此式5应改为

（7）

一对匹配点对可以列写两个方程，那么要想求解只有8个自由度的矩阵H只需要四个匹配点对即可，但要求这四对点每三个都不能共线：

（8）

式中，h是矩阵H的向量形式。Ah＝0被称为齐次线性最小二乘问题，因此求解单应矩阵H就转换为齐次线性最小二乘问题。该方法也称为直接线性变换（DLT）。

式8表示的是一个齐次线性方程组，式中A为矩阵，h为相量，我们需要求非零的h。而h属于A的零空间，有时也称为A的零（右）相量，因此h可以由A的奇异值分解（SVD）得到，也就是说并不需要解方程组就可以得到h。A的奇异值分解形式为：

（9）

则h为V的最后一列，再把h重新整理，就可得到H。

通过式9得到h的方法非常依赖于图像的坐标原点及尺度，这会使算法不稳定，产生数值误差。但我们可以通过归一化的方法使数值收敛于正确的结果。归一化的方法是首先平移图像坐标，使图像的坐标原点为匹配点对的重心；然后是改变尺度，使匹配点到平移后的坐标原点的距离为1。匹配点对所在的两幅图像，以及它们的横、纵坐标都必须分别独立的进行上述归一化处理。归一化的具体计算过程为：

（10）

式中，n表示匹配点对，这里n为4，该式得到了坐标平均值，即坐标平移量。

（11）

式11得到了坐标尺度，即匹配点到原点的距离为1，如果想使距离为 （许多文献上，此处的值都是 ，Opencv用的是1），则只需再乘以 即可。由式10和式11分别构成两个坐标变换矩阵：

（12）

则两幅图像的坐标分别变为：

（13）

即

（14）

这时，我们就需要用坐标变换后的4个匹配点对计算单应矩阵，得到 ，它与真实坐标下的H的关系为：

（15）

式中的T’^-1很容易由式12中的T’得到：

（16）

下面我们给出归一化直接线性变换的步骤：

●由式10和式11分别计算平移量和尺度；

●由式12计算坐标变换矩阵；

●由式13计算坐标变换后的特征点坐标值，并代入式8得到，再经过奇异值分解（式9）得到 ；

●由式15得到最终的单应矩阵H。

以上是由4个匹配点对计算出单应矩阵H，为了增加鲁棒性，可以用更多的匹配点对计算H。此时式8就是一个超定方程组，因为A是2n×9的矩阵，而n>4，所以2n>9。当匹配点对n的数量很大时，直接用SVD的方法就不是很方便，这时我们就可以直接应用最小二乘法求解。因为Ah＝0，它的误差平方函数f(h)可以写为：

（17）

很显然，我们应该使f(h)最小，此时的h即为求解。因此我们对f(h)求导，并使导数为0，则

（18）

A^TA是对称的9×9的方阵，我们对A^TA进行特征值分解，则h应该是A^TA的特征值为0所对应的特征向量，当然因为存在误差，有时会出现没有一个特征值为0的情况，这时h等于最接近0的特征值所对应的特征向量。可以很容易的证明，对A进行奇异值分解与对A^TA进行特征值分解得到的h是完全相等。因此特征值分解的方法也适用于n＝4的情况。

尽管使用了更多的匹配点对，但由于误差的存在，匹配点对的位置可能会出现偏差，即有好的匹配点对（称之为内点），也有不好的匹配点对（称之为外点）。如果用外点计算，则H肯定是不正确的，所以我们应该选择内点。那么如何区分内点和外点呢？通常的做法是应用RANSAC（随机抽样一致性，RANdom SAmple Consensus）算法。

RANSAC算法的核心思想是从样本集中抽取多少次才能保证至少有一次不是外点的概率为q。则据此，可以得到下列关系式：

（19）

式中，N表示抽取的次数，ε表示样本集中外点的比例，s表示每一次抽取样本的数量，在这里s为4，因为求单应矩阵H至少需要4个匹配点对，因此我们需要每一次抽取4个匹配点对。式19又可以表示为：

（20）

式中，q是需要事先确定的参数，q越大，得到的内点就越可信。ε由重映射误差的方法得到，即当得到了单应矩阵H后，可以把图像A中的特征点p重映射到图像B中，计算映射后的点与它的匹配点对的特征点p’之间的几何距离：

（21）

而重映射误差为d²(p’,Hp)。如果该距离小于事先确定好的阈值η，则认为这对匹配点对是内点，否则是外点。计算所有的匹配点对，这样就可以统计出内点和外点的数量，从而得到ε，进而由式20得到N。

图4 重映射误差

下面我们给出RANSAC计算单应矩阵H的步骤：

●给出参数q和η，以及迭代次数N；

●进入迭代，一共需迭代N次：

▲ 随机抽取4个匹配点对，要求任意3个不共线；

▲ 由归一化直接线性变换计算单应矩阵H；

▲ 由重映射误差方法（式21）计算内点的数量，并得到参数ε；

▲ 由式20更新最大迭代次数N；

●在迭代循环内已得到最大的内点数量，由这些内点计算单应矩阵，得到最终的H。

另外，我们还需要用下列的关系式来说明内点数n_i和匹配点对数量n_f之间的关系：

（22）

如果满足式22，则说明图像间有很好的匹配。或者我们用一个匹配置信度c来衡量：

（23）

c越大越好，当然如果c大于3（3为实验数据），则认为这两幅图像十分接近，可以被看成是一幅图像了；而如果c小于一定的值，则认为这两幅图像没有重叠的区域，所以就无法拼接在一起。

应用RANSAC算法的主要目的是得到内点，而用全体内点经过特征值分解计算得到的H误差较大，可以说比较粗糙，因此我们还需要用更加“精细”的方法得到精确的H。由式21可知，H应该使重映射误差（即匹配点对之间的几何距离）达到最小。因此这是一个参数优化估计求极值的问题。基于LM算法（Levenberg-Marquardt）的优化是目前应用较为广泛的一种无条件约束优化方法，该算法是在逼近某个极小点时平方收敛，因此它是专门用于误差平方和最小化的方法。该算法的特点是，它是高斯牛顿法和梯度下降法的结合，既有高斯牛顿法的快速收敛性，也有梯度下降法的全局搜索特性。

下面我们就给出应用LM算法计算H的方法。

首先给出误差指标函数，即所有内点的重映射误差的平方和：

（24）

式中，n表示内点的总数，h是单应矩阵H的向量形式，第i对匹配点对的重映射误差e_i为：

（25）

而e(h)则为：

（26）

LM算法用于迭代计算h的公式为：

（27）

式中，h的上标表示迭代的次数，步长∆h为：