推导四对对应点单应矩阵的计算公式?

本文主要是介绍推导四对对应点单应矩阵的计算公式?，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

转载：http://www.zhihu.com/question/23310855

王小龙数学话题优秀回答者数学，计算机视觉，图形图像处理

26 人赞同

题主说的是二维图片之间的单应变换吧。咱一步一步来：

一、单应和齐次坐标
一个单应矩阵是大小为3*3的矩阵 $H$ ，满足给定一个点 $p_1=[x_1,y_1,w_1]^T$ ， $H$ 把点 $p_1$ 变成一个新的点 $p_2=[x_2,y_2.w_2]^T=Hp_1$ 。注意由于它们都是齐次坐标，对应的图像上的两个点分别是 $[\frac{x_1}{w_1},\frac{y_1}{w_1}]^T$ 和 $[\frac{x_2}{w_2},\frac{y_2}{w_2}]^T$ 。

二、单应的自由度
如果给定一个单应 $H=\{h_{ij}\}$ ，给它的元素乘上同一个数 $a$ ，得到的的单应 $aH$ 和 $H$ 作用相同，因为新单应无非把齐次点 $p_1$ 变成了齐次点 $ap_2$ ， $ap_2$ 和 $p_2$ 对应的图像上的点相同。所以一个单应中只有8个自由元素，一般令右下角的那个元素 $h_{33}=1$ 来归一化。

三、求解单应
8个未知数，需要8个方程来求解，之所以四对点能够求解，是因为一对点提供两个方程。我们假设有两个图像上的点 $[x_1,y_1]^T$ 和 $[x_2,y_2]^T$ ,它们的齐次坐标为： $[x_1, y_1, 1]^T$ 和 $[x_2, y_2, 1]^T$ 带到上面的推倒里可以得到：
$x_2=\frac{x_1h_{11}+y_1h_{12}+h_{13}}{x_1h_{31}+y_1h_{32}+1}$
$y_2=\frac{x_1h_{21}+y_1h_{22}+h_{23}}{x_1h_{31}+y_1h_{32}+1}$
把这两个式子重新组织一下，得到等价的矩阵形式：
$Au=v$
其中：
$A=[\begin{array}{cccccccc}x_1 & y_1 & 1 & 0 & 0 & 0 & -x_1x_2 & -x_2y_1\\0 & 0 & 0 & x_1 & y_1 & 1 & -x_1y_2 & -y_1y_2\end{array} ]$
$u=[\begin{array}{cccccccc}h_{11} & h_{12}&h_{13}&h_{21}&h_{22}&h_{23}&h_{31}&h_{32}\end{array} ]^T$
$v=[\begin{array}{cc}x_2&y_2\end{array} ]^T$
如果有四对不共线匹配点对，这个方程组就能够垒到8行，存在唯一解，如果多于四对点，比如有n对点，方程就垒到2n行，用最小二乘法或SVD分解就可以求解 $H$ 。

四、程序实现
由于点对中可能存在不少错误匹配，因此往往需要使用RANSAC算法剔除错误匹配点对，整个流程在OpenCV中已经集成为函数 findhomography，输入点对坐标，输出单应和一个标记哪些点是正确匹配哪些点是错误匹配的数组。

编辑于 2014-04-26 6 条评论感谢收藏 • 没有帮助 • 举报 • 作者保留权利

胡知知谢邀

3 人赞同

谢邀请。
目的：寻找两组对应3D点之间的Rotation（R）和Translation（T）。
我用的是Absolute Orientation Quaternion
Horn 的《Closed-form solution of absolute orientation using unit quaternions》
方法比较老，如果有其他推荐请忽略。。。。

Paper大致逻辑是：
1. 两组中各取3个点，建立两组新坐标系（3个点：一个原点。另两个和原点可以定义一个面。两个中的一个和原点构成一个轴。垂直于面并过原点的方向为另一个轴。垂直于两个轴为第三个轴。三个相互垂直的轴就是一个新坐标系啦～）

2. 原目的是寻找两组3D点之间的Rotation（R）和Translation（T）. 现在则是寻找新定义的两个坐标系间的R，T。因为我们有第四个点，第四个点可以通过三个坐标系表示：原有的Global Coordinate, 以及1中新定义的两个坐标系。

于是Paper开始推公式找关系。。。。详情请见paper：（ Index of /bkph/papers）

注意： 虽然四个点就能解出来R,T. 但这是在对应点无误差的假设下。这个方法是至少提供4组对应点，实际使用中通常要提供远多于四个点来降低误差对结果的影响。

ETH Zurich 写了个matlab code

使用时请注意Copyright:
========================
% Computes the orientation and position (and optionally the uniform scale
% factor) for the transformation between two corresponding 3D point sets Ai
% and Bi such as they are related by:
%
% Bi = sR*Ai+T
%
% Implementation is based on the paper by Berthold K.P. Horn:
% "Closed-from solution of absolute orientation using unit quaternions"
% The paper can be downloaded here:
% http://people.csail.mit.edu/bkph/papers/Absolute_Orientation.pdf
%
% Authors: Dr. Christian Wengert, Dr. Gerald Bianchi
% Copyright: ETH Zurich, Computer Vision Laboratory, Switzerland
%
% Parameters: A 3xN matrix representing the N 3D points
% B 3xN matrix representing the N 3D points
% doScale Flag indicating whether to estimate the
% uniform scale factor as well [default=0]
%
% Return: s The scale factor
% R The 3x3 rotation matrix
% T The 3x1 translation vector
% err Residual error
%
% Notes: Minimum 3D point number is N > 4
========================

如下：

function [s R T err] = absoluteOrientationQuaternion( A, B, doScale)

%Argument check
if(nargin<3)
doScale=1;
end
%Return argument check
if(nargout<1)
usage()
error('Specify at least 1 return arguments.');
end
%Test size of point sets
[c1 r1] = size(A);
[c2 r2] = size(B);
if(r1~=r2)
usage()
error('Point sets need to have same size.');
end
if(c1~=3 | c2~=3)
usage()
error('Need points of dimension 3');
end
if(r1<4)
usage()
error('Need at least 4 point pairs');
end

%Number of points
Na = r1;

%Compute the centroid of each point set
Ca = mean(A,2);
Cb = mean(B,2);

%Remove the centroid
An = A - repmat(Ca,1,Na);
Bn = B - repmat(Cb,1,Na);

%Compute the quaternions
M = zeros(4,4);
for i=1:Na
%Shortcuts
a = [0;An(:,i)];
b = [0;Bn(:,i)];
%Crossproducts
Ma = [ a(1) -a(2) -a(3) -a(4) ;
a(2) a(1) a(4) -a(3) ;
a(3) -a(4) a(1) a(2) ;
a(4) a(3) -a(2) a(1) ];
Mb = [ b(1) -b(2) -b(3) -b(4) ;
b(2) b(1) -b(4) b(3) ;
b(3) b(4) b(1) -b(2) ;
b(4) -b(3) b(2) b(1) ];
%Add up
M = M + Ma'*Mb;
end

%Compute eigenvalues
[E D] = eig(M);

%Compute the rotation matrix
e = E(:,4);
M1 = [ e(1) -e(2) -e(3) -e(4) ;
e(2) e(1) e(4) -e(3) ;
e(3) -e(4) e(1) e(2) ;
e(4) e(3) -e(2) e(1) ];
M2 = [ e(1) -e(2) -e(3) -e(4) ;
e(2) e(1) -e(4) e(3) ;
e(3) e(4) e(1) -e(2) ;
e(4) -e(3) e(2) e(1) ];

R = M1'*M2;

%Retrieve the 3x3 rotation matrix
R = R(2:4,2:4);
%Compute the scale factor if necessary
if(doScale)
a =0; b=0;
for i=1:Na
a = a + Bn(:,i)'*R*An(:,i);
b = b + Bn(:,i)'*Bn(:,i);
end
s = b/a;
else
s = 1;
end

%Compute the final translation
T = Cb - s*R*Ca;

%Compute the residual error
if(nargout>3)
err =0;
for i=1:Na
d = (B(:,i) - (s*R*A(:,i) + T));
err = err + norm(d);
end
err = (err)/Na;
end

%Displayed if an error occurs
function usage()
disp('Usage:')
disp('[s R T error] = absoluteOrientationQuaternion( A, B, doScale)')
disp(' ')
disp('Return values:')
disp('s The scale factor')
disp('R The 3x3 rotation matrix')
disp('T The 3x1 translation vector')
disp('err Residual error (optional)')
disp(' ')
disp('Input arguments:')
disp('A 3xN matrix representing the N 3D points')
disp('B 3xN matrix representing the N 3D points')
disp('doScale Optional flag indicating whether to estimate the uniform scale factor as well [default=0]')
disp(' ')