多视几何 007：直接线性变换

理论

矩阵H会将一幅图像上的一个点的坐标a=(x,y,1)映射成另一幅图像上的点的坐标b=(x1,y1,1)，但H是未知的。通常需要根据在同一平面上已知的一些点对(比如a,b)来求H。 假设已知点对(a,b)，则有下面的公式：

$$H=\left[\begin{array}{ccc}{h}_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\end{array}\right]$$
$$\begin{gathered} \left[\begin{array}{c}{x}_1\\ y_1\\ 1\end{array}\right]=\left[\begin{array}{ccc}{h}_{11}& h_{12}& h_{13}\\ h_{21}& h_{22}& h_{23}\\ h_{31}& h_{32}& h_{33}\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right] \\ 即\{\begin{array}{l}{x}_1=h_{11}x+h_{12}y+h_{13}\\ y_1=h_{21}x+h_{22}y+h_{23}\\ 1=h_{31}x+h_{32}y+h_{33}\end{array} \end{gathered}$$

$$\begin{gathered} \because 1=h_{31}x+h_{32}y+h_{33}\Rightarrow \\ \{\begin{array}{l}{x}_1=\frac{x_1}{1}=\frac{h_{11}x+h_{12}y+h_{13}}{h_{31}x+h_{32}y+h_{33}}\\ y_1=\frac{y_1}{1}=\frac{h_{21}x+h_{22}y+h_{23}}{h_{31}x+h_{32}y+h_{33}}\end{array} \\ 令A=\left[\begin{array}{ccccccccc}-x& -y& -1& 0& 0& 0& x{x}_1& y{x}_1& x_1\\ 0& 0& 0& -x& -y& -1& x{y}_1& y{y}_1& y_1\end{array}\right]， \\ h=[h_{11},h_{12},h_{13},h_{21},h_{22},h_{23},h_{31},h_{32},h_{33}{]}^T,则两点的对应可写为： \\ {A}_{2\times 9}\ast h_{9\times 1}=0 \end{gathered}$$

因为齐次坐标（即(x,y,1)）来表示平面上的点，aH与H的作用是相同的，所以变换的自由度为8，一对点能产生两个方程，共需要4对对应点，即可求取H矩阵；如超过4对，通过最小二乘法或RANSAC求取最优参数。

$$A=U\ast \Sigma \ast V^T$$
然后取V的最后一列出来作为求解h。因为矩阵A是行满秩，即只有一个自由度。

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

实现

添加链接描述

function v = homography_solve(pin, pout)
% HOMOGRAPHY_SOLVE finds a homography from point pairs
%   V = HOMOGRAPHY_SOLVE(PIN, POUT) takes a 2xN matrix of input vectors and
%   a 2xN matrix of output vectors, and returns the homogeneous
%   transformation matrix that maps the inputs to the outputs, to some
%   approximation if there is noise.
%
%   This uses the SVD method of
%   http://www.robots.ox.ac.uk/%7Evgg/presentations/bmvc97/criminispaper/node3.html

% David Young, University of Sussex, February 2008

if ~isequal(size(pin), size(pout))
    error('Points matrices different sizes');
end
if size(pin, 1) ~= 2
    error('Points matrices must have two rows');
end
n = size(pin, 2);
if n < 4
    error('Need at least 4 matching points');
end

% Solve equations using SVD
x = pout(1, :); y = pout(2,:); X = pin(1,:); Y = pin(2,:);
rows0 = zeros(3, n);
rowsXY = -[X; Y; ones(1,n)];
hx = [rowsXY; rows0; x.*X; x.*Y; x];
hy = [rows0; rowsXY; y.*X; y.*Y; y];
h = [hx hy];
if n == 4
    [U, ~, ~] = svd(h);
else
    [U, ~, ~] = svd(h, 'econ');
end
v = (reshape(U(:,9), 3, 3)).';
end