estimatePerspectiveTransform.m

function [t1 t2] = estimatePerspectiveTransform(x,y,X,Y)
% given two planar quadrilaterals, compute a perspective transformation
% that maps the first onto the second
%
% useful, for example, for inverting the perspective transformation that an
% image undergoes when photographed. this function is based on a simplified
% model of the perspective transformation that maps quadrilaterals on an
% arbitrary, unspecified plane in world-space, to the plane in image-space.
% simplified, relative to a full estimate of camera position, orientation,
% view angle, etc etc. it simply determines a transformation that maps one
% quadrilateral to another; it does NOT inherently provide an estimate of
% the relative position/orientation of the world-plane/camera, or the 
% camera view angle. i dont think these things are needed for basic 
% photography correction...
%
% usage:
% [t1 t2] = estimatePerspectiveTransform(x,y,X,Y);
% x,y is source quad
% X,Y is dest quad
%
% given (x1,y1) in source-space, compute (X1,Y1) in dest-space:
%
% [X1 Y1 ?] = t1*[x1;y1;1]/(t2*[x1;y1;1]);
%
% this function can be used to demonstrate that although lines transform to
% lines, linearly spaced points along lines transform to nonlinearly-spaced
% points along lines - one way of describing perspective transformation
% distortion
% 
% source: http://alumni.media.mit.edu/~cwren/interpolator/
%
% TODO: is there a good way to estimate the transform using sets of
% image-space lines that are known to be parallel in corrected-space?

% this is the transform:
% [XW]   [a b c] [x]
% [YW] = [d e f] [y]
% [ W]   [g h 1] [1]
%
% followed by the division: (XW,YW)/W (the nonlinear component)
% equivalently:
% X = (ax+by+c)/(gx+hy+1)
% Y = (dx+ey+f)/(gx+hy+1)
%
% with 8 unknown coefficients, we need at least 8 equations to estimate 
% their values. each point provides two equations, so four points, which
% define a quadrilateral, can be used to create this system, with simple
% algebra:
% [x1 y1 1  0  0 0 -X1x1 -X1y1] [a]   [X1]
% [ 0  0 0 x1 y1 1 -Y1x1 -Y1y1] [b]   [Y1]
% [x2 y2 1  0  0 0 -X2x2 -X2y2] [c]   [X2]
% [ 0  0 0 x2 y2 1 -Y2x2 -Y2y2] [d] = [Y2]
% [...                        ] [.]   [..]
%
% or A*l = B
% solve this:
% A.'*A*l = A.'*B
% l = (A.'*A)^-1*A.'*B
%
% l contains the transform coefficients
% these are separated into the linear and nonlinear coefficients, and
% returned as two matrices:
% t1 = [a b c; d e f; 0 0 1]
% t2 = [g h 1]


% create A matrix
A = [ x y ones(size(x)) zeros(4,3)        -x.*X -y.*X ...
      zeros(4,3)        x y ones(size(x)) -x.*Y -y.*Y ];
A = reshape (A', 8 , 8 )';

% create B vector
B = [ X , Y ];
B = reshape (B', 8 , 1 );

% estimate solution vector
l = inv(A' * A) * A' * B;

% separate into linear and nonlinear components
t1 = reshape([l(1:6)' 0 0 1 ],3,3)';
t2 = [l(7:8)' 1];


% original, with stupid variable names:
% B = [ x y ones(size(x)) zeros(4,3)        -x.*X -y.*X ...
%       zeros(4,3)        x y ones(size(x)) -x.*Y -y.*Y ];
% B = reshape (B', 8 , 8 )';
% D = [ Xp , Yp ];
% D = reshape (D', 8 , 1 );
% l = inv(B' * B) * B' * D;
% A = reshape([l(1:6)' 0 0 1 ],3,3)';
% C = [l(7:8)' 1];