% IMTRANS
%
% Applies a geometric transform to an image
%
%  newim = imTrans(im, T, region, sze);
%
%  Where: 
%        im is the image to be transformed.
%        T is a 3x3 homogeneous transformation matrix.
%        region is an optinal 4 element vector specifying 
%               [minrow maxrow mincol maxcol] to transform.
%        sze is an optional desired size of the transformed image
%               (this is the maximum No of rows or columns).
%
%  The region argument is used when one is inverting a perspective
%  transformation and the vanishing line of the plane lies within
%  the image.  Attempts to transform any part of the vanishing line
%  will position you at infinity.  Accordingly one should specify 
%  a region that excludes any part of the vanishing line.
%
%  The sze parameter is optionally used to control the size of the
%  output image.  When inverting a perpective or affine transformation
%  the scale parameter is unknown/arbitrary, and without specifying
%  it explicitly the transformed image can end up being very small 
%  or very large.

%  Peter Kovesi  April 2000
%  Department of Computer Science
%  The University of Western Australia

function newim = imTrans(im, T, region, sze);

im = double(im)/255;  % Make sure image is double 

threeD = (ndims(im)==3);  % A colour image

if threeD    % Transform red, green, blue components separately
  r = transformImage(im(:,:,1), T, region, sze);
  g = transformImage(im(:,:,2), T, region, sze);
  b = transformImage(im(:,:,3), T, region, sze);

  % Fudge to correct the image - for some reason
  % it comes out mirrored left-to right
  %  ... I have to find what the real problem is
  r = fliplr(r);
  g = fliplr(g);
  b = fliplr(b);

  newim = repmat(uint8(0),[size(r),3]);
  newim(:,:,1) = uint8(round(r*255));
  newim(:,:,2) = uint8(round(g*255));
  newim(:,:,3) = uint8(round(b*255));

else                % Assume the image is greyscale
  newim = transformImage(im, T, region, sze);
  % Fixing fudge
  newim = fliplr(newim);      
end

%------------------------------------------------------------

% The internal function that does all the work

function newim = transformImage(im, T, region, sze);

[rows, cols] = size(im);

% Determine default parameters if needed
if nargin == 2
  region = [0 rows 0 cols];
  sze = max(rows,cols);
elseif nargin == 3
  sze = max(rows,cols);
elseif nargin ~= 4
  error('Incorrect arguments to imtrans');
end

% Find where corners go - this sets the bounds on the final image
B = bounds(T,region);
nrows = B(2) - B(1);
ncols = B(4) - B(3);

% Determine any rescaling needed
s = sze/max(nrows,ncols);

S = [s 0 0        % Scaling matrix
     0 s 0
     0 0 1];

T = S*T;
Tinv = inv(T);

% Recalculate the bounds of the new (scaled) image to be generated
B = bounds(T,region);
nrows = B(2) - B(1);
ncols = B(4) - B(3);

newim = zeros(nrows,ncols);

[x,y] = meshgrid(1:ncols,1:nrows);    % All possible xy coords in the image.

% Transform these xy coords.
sxy = homoTrans(Tinv, [x(:)'+B(3) ; y(:)'+B(1) ; ones(1,ncols*nrows)]);

xi = reshape(sxy(1,:),nrows,ncols);
yi = reshape(sxy(2,:),nrows,ncols);
[x,y] = meshgrid(1:cols,1:rows);

newim = interp2(x,y,double(im),xi,yi); % Interpolate values from source image


%---------------------------------------------------------------------
%
% Internal function to find where the corners of a region, R
% defined by [minrow maxrow mincol maxcol] are transformed to 
% by transform T and returns the bounds, B in the form 
% [minrow maxrow mincol maxcol]

function B = bounds(T, R)

P = [R(3) R(4) R(4) R(3)      % homogeneous coords of region corners
     R(1) R(1) R(2) R(2)
      1    1    1    1   ];
     
PT = round(homoTrans(T,P)); 

B = [min(PT(2,:)) max(PT(2,:)) min(PT(1,:)) max(PT(1,:))];
%      minrow          maxrow      mincol       maxcol