单应性矩阵的理解及求解

1. 齐次坐标（Homogeneous Coordinate）

一幅2D图像上的非齐次坐标为(x,y)，而齐次坐标为(x,y,1)，也可以写成(x/z,y/z,1)或(x,y,z)。齐次坐标有很多好处，比如可以很清楚的确定一个点在不在直线上：

T(x)*I=0, T表示转置矩阵；

还可以描述无穷远点：(x,y,0)；

还可以把平移和旋转写到一个矩阵里（也有不愿意这么干的，摄影测量里用的都是非齐次坐标，平移和旋转分开写）；

具体的可去看<计算机视觉中的多视图几何>,i.e.,<multiple view geometry in computer vision>。

2. 单应性矩阵（Homography Matrix）

两个不同视角的图像上的点对的homogeneous coordinate可以用一个射影变换（projective transformation）表述，即：x1 = H*x2

二维和三维的图示如下：

射影变换也叫“单应”--Homography，“Homo”前缀就是same的意思，表示“同”，homography就是用同一个源产生的graphy，中文译过来大概就是“单应”。

因此上面式子中的矩阵H就叫单应性矩阵。上式中的x1和x2都是3*1的齐次坐标，因此H是一个3*3的矩阵：{h00,h01,h02;h10,h11,h12;h20,h21,h22}。

然而，单应矩阵的自由度是多少呢？

如果给定一个单应H={h_ij}，给它的元素乘上同一个数a，得到的的单应a*H和H作用相同，因为新单应无非把齐次点x1变成了齐次点a*x1，都是一回事。因此我们可以把a换成1/h22，那么H就变成了只有8个自由元素的矩阵。

那么需要多少个点对求解这个H呢？如果需要唯一解的话，需要4个点对（对应8个方程，去解H中的8个未知数）。

3. 两个应用例子

例子1：如何进行两个图像的拼接？

解答：通过图像特征点匹配+RANSAC可以获得图像之间的单应性矩阵，然后把其中一个图像通过这个矩阵投影到另一个图像上就完成了基本的拼接。

例子2：在相机的内参数标定过程中会用到求解单应性矩阵吗？

对一个棋盘格拍照，棋盘格的世界坐标系是用户任意设定的，标定的时候，默认世界坐标系就是以标定板左上角点为原点，z轴垂直于标定板，xoy面与标定板重合的三维直角坐标系。棋盘格的格子长度已知，因此可以知道各个角点的世界坐标系坐标XYZ（Z=0）。由于Z=0，因此可以忽略掉Z这个维度，世界坐标系中某个坐标(X,Y,Z,1)到图像坐标(x,y,1)的变换就等价于(X,Y,1)到图像坐标(x,y,1)的变换。此时的变换矩阵就由3*4变为3*3，成为了单应性矩阵。

因此用4个角点就可以计算H的8个参数。以不同的角度对棋盘格拍摄3张就可以得到3个H，如果用张正友标定法，就可以得到6个约束方程，可以求解B矩阵（对称矩阵）的6个未知参数，进而通过Cholesky分解求解出内参矩阵A的参数。

4. 单应矩阵Homography求解

《从零开始学习「张氏相机标定法」（一）成像几何模型》中我们已经得到了像素坐标系和世界坐标系下的坐标映射关系：

其中，u、v表示像素坐标系中的坐标，s表示尺度因子，fx、fy、u0、v0、γ（由于制造误差产生的两个坐标轴偏斜参数，通常很小）表示5个相机内参，R，t表示相机外参，Xw、Yw、Zw（假设标定棋盘位于世界坐标系中Zw=0的平面）表示世界坐标系中的坐标。

单应性概念的引出

我们在这里引入一个新的概念：单应性（Homography）变换。可以简单的理解为它用来描述物体在世界坐标系和像素坐标系之间的位置映射关系。对应的变换矩阵称为单应性矩阵。在上述式子中，单应性矩阵定义为：

其中，M是内参矩阵。

从单应矩阵定义式子来看，它同时包含了相机内参和外参。在进一步介绍相机标定知识之前，我们重点来了解一下单应性，这有助于深入理解相机标定。因为在计算机视觉领域，单应性是一个非常重要的概念。

为了不让读者淹没在公式的汪洋大海中失去兴趣，我们颠倒一下顺序，先来看看单应性到底有什么用，然后再介绍单应矩阵的估计方法。

单应性在计算机视觉中的应用

单应性在计算机视觉领域是一个非常重要的概念，它在图像校正、图像拼接、相机位姿估计、视觉SLAM等领域有非常重要的作用。

1. 图像校正

用单应矩阵进行图像矫正的例子如下图所示，最少需要四个对应点对（后面会给出原因）就可以实现。

2. 视角变换

单应矩阵用于视角变换的例子如下图所示，可以方便地将左边普通视图转换为右图的鸟瞰图。

3. 图像拼接

既然单应矩阵可以进行视角转换，那我们把不同角度拍摄的图像都转换到同样的视角下，就可以实现图像拼接了。如下图所示，通过单应矩阵H可以将image1和image2都变换到同一个平面。

单应矩阵用于图像拼接的例子如下所示。

4. 增强现实（AR）

平面二维标记图案（marker）经常用来做AR展示。根据marker不同视角下的图像可以方便的得到虚拟物体的位置姿态并进行显示，如下图所示。

如何估计单应矩阵？

了解了上述单应性的部分应用后，我们就有很大的动力来学习单应矩阵的推导和计算了。首先，我们假设两张图像中的对应点对齐次坐标为(x',y',1)和(x,y,1)，单应矩阵H定义为：

则有：

矩阵展开后有3个等式，将第3个等式代入前两个等式中可得：

也就是说，一个点对对应两个等式。在此插入一个讨论：单应矩阵H有几个自由度？

或许有人会说，9个啊，H矩阵不是9个参数吗？从h11到h33总共9个。真的是这样吗？实际上并不是，因为这里使用的是齐次坐标系，也就是说可以进行任意尺度的缩放。比如我们把hij乘以任意一个非零常数k并不改变等式结果：

所以实际上单应矩阵H只有8个自由度。8自由度下H计算过程有两种方法。

第一种方法：直接设置h33=1，那么上述等式变为：

第二种方法：将H添加约束条件，将H矩阵模变为1，如下：

以第2种方法（用第1种也类似）为例继续推导，我们将如下等式（包含||H||=1约束）:

乘以分母展开，得到：

整理，得到：

假如我们得到了两幅图片中对应的N个点对（特征点匹配对），那么可以得到如下线性方程组：

写成矩阵形式：

由于单应矩阵H包含了||H||=1约束，因此根据上图的线性方程组，8自由度的H我们至少需要4对对应的点才能计算出单应矩阵。这也回答了前面图像校正中提到的为何至少需要4个点对的根本原因。

但是，以上只是理论推导，在真实的应用场景中，我们计算的点对中都会包含噪声。比如点的位置偏差几个像素，甚至出现特征点对误匹配的现象，如果只使用4个点对来计算单应矩阵，那会出现很大的误差。因此，为了使得计算更精确，一般都会使用远大于4个点对来计算单应矩阵。另外上述方程组采用直接线性解法通常很难得到最优解，所以实际使用中一般会用其他优化方法，如奇异值分解、Levenberg-Marquarat（LM）算法（后续文章会介绍）等进行求解。

如何根据标定图得到单应矩阵？

经过前面一系列的介绍，我们应该大致明白如何根据打印的棋盘标定图和拍摄的照片来计算单应矩阵H。我们来总结一下大致过程。

1. 打印一张棋盘格标定图纸，将其贴在平面物体的表面。

2. 拍摄一组不同方向棋盘格的图片，可以通过移动相机来实现，也可以移动标定图片来实现。

3. 对于每张拍摄的棋盘图片，检测图片中所有棋盘格的特征点（角点，也就是下图中黑白棋盘交叉点，中间品红色的圆圈内就是一个角点）。我们定义打印的棋盘图纸位于世界坐标系Zw=0的平面上，世界坐标系的原点位于棋盘图纸的固定一角（比如下图中黄色点）。像素坐标系原点位于图片左上角。

4. 因为棋盘标定图纸中所有角点的空间坐标是已知的，这些角点对应在拍摄的标定图片中的角点的像素坐标也是已知的，如果我们得到这样的N>=4个匹配点对（越多计算结果越鲁棒），就可以根据LM等优化方法得到其单应矩阵H。当然计算单应矩阵一般不需要自己写函数实现，OpenCV中就有现成的函数可以调用，对应的C++函数是：

Mat findHomography(InputArray srcPoints, InputArray dstPoints, int method=0, double ransacReprojThreshold=3, OutputArray mask=noArray() )

从函数定义来看，只要输入匹配点对，指定具体计算方法即可输出结果。