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

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

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

26 人赞同

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

一、单应和齐次坐标
一个单应矩阵是大小为3*3的矩阵 ，满足给定一个点 ， 把点 变成一个新的点 。注意由于它们都是齐次坐标，对应的图像上的两个点分别是 和 。

二、单应的自由度
如果给定一个单应 ，给它的元素乘上同一个数 ，得到的的单应 和 作用相同，因为新单应无非把齐次点 变成了齐次点 ， 和 对应的图像上的点相同。所以一个单应中只有8个自由元素，一般令右下角的那个元素 来归一化。

三、求解单应
8个未知数，需要8个方程来求解，之所以四对点能够求解，是因为一对点提供两个方程。我们假设有两个图像上的点 和 ,它们的齐次坐标为： 和 带到上面的推倒里可以得到：


把这两个式子重新组织一下，得到等价的矩阵形式：

其中：



如果有四对不共线匹配点对，这个方程组就能够垒到8行，存在唯一解，如果多于四对点，比如有n对点，方程就垒到2n行，用最小二乘法或SVD分解就可以求解 。

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

编辑于 2014-04-26  6 条评论  感谢 

分享

  收藏  •  没有帮助  •  举报  •  作者保留权利

收起

3赞同 反对，不会显示你的姓名

胡知知 谢邀

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(' ')

发布于 2014-04-12  5 条评论  感谢 

分享

  收藏  •  没有帮助  •  举报  •  作者保留权利

0赞同 反对，不会显示你的姓名

庞阿困 别人笑我太疯癫，我夸他人看得准。

泻药。@王小龙已经回答得很清楚了。可以看看《multiple view geometry》第一章，讲的挺清楚。