The computation of homography, essential and fundamental matrix-阿里云开发者社区

本次打算梳理下最基本的几个矩阵之间的关系以及计算，总结大体内容：

1. 单应性矩阵的基本概念

什么是单应性矩阵？

单应性变换包含什么样的射影组合(projective transformation)？

单应性关系的前提条件？

单应性与极几何的联系？

2. 单应性矩阵的计算

本质矩阵和基础矩阵的性质，上一篇博文有详细介绍，所以此处只讲计算方法了。

3. 基础矩阵的计算及计算前提条件

4.本质矩阵的计算

5.给定本质矩阵，如何提取相机朝向R和相对位移关系(基线方向T)

1. 单应性矩阵的基本概念

什么是单应性矩阵？

单应性是几何中的一个概念。单应性是一个从实射影平面到射影平面的可逆变换，直线在该变换下仍映射为直线。它是一对透视投影的组合射影变换群。它描述了当观察者视角改变时，被观察物体的感知位置会发生何种变化。射影变换并不保持大小和角度，但会保持重合关系和交比。

在计算机视觉领域中，空间中同一平面的任意两幅图像通过单应性关联在一起（假定是针孔相机）。比如，一个物体可以通过旋转相机镜头获取两张不同的照片(这两张照片的内容不一定要完全对应,部分对应即可),我们可以把单应性设为一个二维矩阵M,那么照片1乘以M就是照片2. 这有着很多实际应用，比如图像校正、图像对齐或两幅图像之间的相机运动计算（旋转和平移）等。一旦旋转和平移从所估计的单应性矩阵中提取出来，那么该信息将可被用来导航或是把3D物体模型插入到图像或视频中，使其可根据正确的透视来渲染，并且成为原始场景的一部分。

假设图A中有一个点x的坐标为(u,v,1)，图B中对应匹配点x’坐标为(u’,v’,1)，那么单应性矩阵可以将两个坐标联系起来：

注意，这种单应性联系不依赖于观察场景结构，即无论相机拍摄什么画面，两个画面的单应性关系保持不变。

单应性变换包含什么样的射影组合(projective transformation)？

首先看等距变换Isometries transformation，它相当于是平移变换和旋转变换的复合。

上式也可以写作：

上式中ε等于正负1，当ε=1时，等距变换与原结构保持相同方向，即是欧氏变换，表达刚性物体运动；当ε=-1时，等距变换是逆向的。另外，R是正交矩阵。

等距变换有三个自由度，一个旋转，两个平移。等距变换后目标结构的长度、角度、面积与原结构相同。

再看相似变换(Similarity Transformation)，它相当于是等距变换和均匀缩放的一个复合。

相似变换可以从两组匹配点计算出来，它有四个自由度，一个尺度，一个旋转，一个平移。相似变换后目标结构内部的角度、长度比、面积比保持不变。R是正交矩阵

再看仿射变换(Affine Transformation)，它相当于一个平移变换和一个非均匀变换的复合。

仿射矩阵A是2阶非奇异矩阵，A通过奇异值分解为旋转rotation和非均匀缩放deformation：

上式中U和V都是正交矩阵，对A简单描述下：λ1和λ2代表旋转后的目标在x和y方向的缩放比例，所以是非均匀的。Φ代表缩放的方向，θ代表旋转的方向。

仿射变换中A有六个自由度，2个尺度，2个旋转，2个平移。经过仿射变换后，原来平行的线还是平行的，并且平行线段的长度之比保持不变，面积之比也是不变的。另外仿射变换是否保向，取决于A的行列式值的正负。

最后看射影变换projectivetransformation，也就是单应性矩阵，他表达齐次坐标的非奇异线性变换。

v是射影变换后平行直线可能不再平行的原因。

由于齐次坐标对于尺度具有不变性，可以从H中提取该尺度(比如v)，因此单应性矩阵有8个自由度。两个平面之间的射影关系，可以通过4组匹配点计算得到(注意任意3点不可共线)。经过射影变换后，目标结构的交比cross ratio保持不变(两个平面上分别有4个点，平面1上线段AB与线段CD之比，等于平面2上线段A’B’与线段C’D’之比)。另外变换前后共点，共线，相切，拐点，切线的不连续性和岐点保持不变。

射影变换可以看作是一些列变换的组合：

上三角矩阵K的行列式值为1，v不等于0，如果s的正负确定，则H是唯一的。

HP (2 dof) moves theline at infinity; HA (2 dof) affects the affine properties, but does not movethe line at infinity; and finally, HS is a general similarity transformation (4dof) which does not affect the affine or projective properties.

总结：

射影变换后，面积的缩放与位置有关，直线变换方向也与面积有关。
射影变换是不可逆的。
仿射变换位于射影变换和相似变换之间，仿射变换推广相似变换使得夹角不再保持不变，造成物体形状在变换后产生歪斜，他在变换中，直线只与方向有关，与位置无关，不会像射影变换那样，越远越小

下图中参数数量与上面不符，因为他针对的是3D点，上面是2D点

单应性关系的前提条件？

如果两幅图像之间的相机运动只有旋转而没有平移，那么这两幅图像通过单应性关联在一起（假定是针孔相机）。也就是说(if and only if)：

a. Both images are viewing thesame plane from a different angle

b. Both images are taken from thesame camera but from a different angle.( Camera is rotated about its center ofprojection without any translation)

也就是说:

a.当两个相机所代表的视平面(image plane)相互平行，stereo normal case，单应性矩阵才是有效的。(baseline不为0，但是极点在无穷远)

b.当且仅当相机是绕相机中心旋转，没有偏移，baseline什么的都是零，单应性矩阵才是有效的。(baseline为0)

单应性与极几何的联系？

下图为极几何的经典视图，相机A的中心为O1，相机B的中心为O2，它们共同看到空间一点P。根据极几何的性质，空间中所有落在直线O1P上的点，最终都会落在在Image Plane B的极线上。

根据单应性的前提条件，单应性与极几何的区别在于两个相机主点之间的距离为0时，

总结起来：

当相机之间存在偏移，那么该用极几何关系描述点的对应关系，此时两个图片之间点的匹配依赖于物点的深度，即物点到相机的距离。

It maps a point in one image to a line in the other image. Actualmatching point depends on the depth of that point,that is its distance from thecamera

当相机之间不存在偏移时，极几何关系不再适用，该用单应性描述点的对应关系，此时两个图片之间点的匹配不依赖于物点的深度。(此时三角方程不成立)

It maps a point in one image to a point in the other image. Themapping does not depend on the depth (distance) of point to the camera. Thismakes sense because baseline z zero means ordinary stereo triangulationequations fail.

从下图看出相机拍摄时并未有偏移，故通过单应性做的图像拼接。

2. 单应性矩阵的算法

从单应性矩阵的定义出发：

除此之外，我们可以通过归一化坐标，使得h的结果是numerical stable。但是即使这样，这种做法对于噪声非常敏感，有几个outlier就完蛋了。

而第二种方法就很好，

为了使上式为0，咱只要对左侧ATA进行奇异值分解，挑选奇异值前九个最小的v的列向量作为h即可，SVD：

Let h be the column ofV (unit eigenvector) associated with the smallest eigenvalue in D.(if only 4points, that eigenvalue will be 0)

Matlab代码：

function H = find_H(P1, P2)
% H = find_H(P1, P2)
% 
%    P1 are coordinates from source frame ([n*2], n>= 4)
%    P2 are coordinates from destination frame ([n*2], n>= 4)
%    H is the 3*3 homography matrix from P1 to P2
%      [P2; 1] ~ H [P1; 1]
%  

x = P1(:, 1);
y = P1(:, 2);
X = P2(:, 1);
Y = P2(:, 2);

A = zeros(length(x(:))*2,9);

for i = 1:length(x(:)),
  a = [x(i),y(i),1];
  b = [0 0 0];
  c = [X(i);Y(i)];
  d = -c*a;
  A((i-1)*2+1:(i-1)*2+2,1:9) = [[a b;b a] d]; %拼接矩阵A，自己对照截图，一次两行，总共2N行
end

[U S V] = svd(A);
h = V(:,9);
H = reshape(h,3,3)';

3. 基础矩阵F的计算及前提条件

分别取[xn’,yn’,1]，[xn’’,yn’’,1]与F的第一个元素，得到展开式

按照下图定义an和f，得到线性方程

我们同样使用奇异值分解去求Af=0这个方程。由于基础矩阵F是齐次的，矩阵A的秩不大于8，因此我们只要8对匹配点就可以解该方程。实际场景中有可能有不止8对匹配点，而且这些匹配点中存在噪声，矩阵A可能是满秩的(should not be)。此时，我们把A中最小的奇异值对应的奇异向量作为基础矩阵F的解。

通过下图可以清晰得看到SVD全貌，即V的最后一列(对应D的最小的九个奇异值)是F的解。但是这还不够，现在并不能保证F的秩等于2(否则Af全为0了，不符合噪声实际情况)。

现在的目标是使r(F)=2，并且F应该与我们的估计F’无限接近。怎么做？

我们对F本身进行第二次奇异值分解，从F’的奇异值中挑最小的，并将它设置为0.也就是说，根据上一步得到的F’，我们强制修改F’的奇异值，而不改变U’，V’，使得rank(F’)=2。

Matlab代码：

function F = F_from_point_pairs(xs, xss)
%xs, xss: Nx3 homologous point coordinates, N>7
%F: 3x3 fundamental matrix

%coefficient matrix
for n = 1 : size(xs, 1)
	A(n, :) = kron(xss(n, :), xs(n, :));
end

%singular value decomposition
[U, D, V] = svd(A);

%approximate F, possibly regular
Fa = reshape(V(:,9), 3, 3)';

% svd decomposition of F
[Ua, Da, Va] = svd(Fa);

%algebraically best F, singular
F = Ua * diag([Da(1,1), Da(2,2), 0]) * Va';

注：

在求单应性矩阵时扯到过坐标归一化，防止numerically unstable，这里详细讲下。

数值不稳定是指，比如对于一副3000*4000的图片，由于其图片像素坐标过大，使得计算结果不稳定。我们按照下图所示，将原点移动到图像中心，将像素点坐标约束到(-1,1)。

坐标归一化，会改变真实F的值，那么坐标归一化前后基础矩阵的对应关系为：

设存在变换T使得坐标x归一化

哪些情况下，基础矩阵的F会很不稳定?

a.如果所有匹配点对应的物点均在同一平面上，即rank(A)<8，此时Af=0无解。就好比很难对墙面上的画做三维重建

b.如果两个相机的光心重合，即相机没有移动，纯粹旋转，则skew-symmetric矩阵全为0，此时需要考虑使用单应性关系。

4. 本质矩阵E的计算

我们可以模仿之前求基础矩阵的方法求本质矩阵，不过此时E的约束关系除了秩必须为2，还要求E的奇异值均为1。

Matlab代码：

function E = E_from_point_pairs(xs, xss)
%xs, xss: Nx3 homologous point coordinates, N>7
%E: 3x3 fundamental matrix

%coefficient matrix
for n = 1 : size(xs, 1)
	A(n, :) = kron(xss(n, :), xs(n, :));
end

%singular value decomposition
[U, D, V] = svd(A);

%approximate F, possibly regular
Ea = reshape(V(:,9), 3, 3)';

% svd decomposition of E
[Ua, Da, Va] = svd(Ea);

%algebraically best F, singular
E = Ua * diag([1, 1, 0]) * Va';

注意:

上面代码函数的输入参数，不是图像坐标，而是乘以相机内参后的图像坐标，xs=K*xs

坐标归一化同样使用此处