
function warpedIm = affineTransf( inputim)

[ny,nx] = size(inputim);

clf;
imagesc(inputim);
autopick = input('Do you want Western North America autopick?  (1=yes, 0=no)  ');

if( autopick==0)
	disp('Click on four points in the input image.');	
	[xi,yi] = ginput(4);
	xi = round(xi);  yi = round(yi);	
	hold on;
	for n = 1:4
		plot( xi(n), yi(n), 'r*');
	end
	disp('Now click on the locations of \n those four points in the \n output image.');	
	[xo,yo] = ginput(4);
	xo = round(xo);  yo = round(yo);
	for n = 1:4
		plot( xo(n), yo(n), 'b*');
	end
	hold off;
else
	ileft = 1;
	while( inputim(1,ileft) == 0)  
		ileft = ileft+1;
	end
	iright = nx;
	while( inputim(ny,iright) == 0)  
		iright = iright-1;
	end
	yi = [ny;1;1;ny];  xi = [1;ileft;nx;iright];
	yo = [ny;1;1;ny];  xo = [1;1;nx;nx];
end

disp('size(inputim)');
size(inputim)
nx
ny
xi
yi

%
% compute INVERSE transform coefficients, i.e.,
% a vector p(8) so that 
%
%  x = p(1)x' + p(2)y' + p(3)
%      ----------------------
%      p(7)x' + p(8)y' + 1
%
%  y = p(4)x' + p(5)y' + p(6)
%      ----------------------
%      p(7)x' + p(8)y' + 1
%
% defines the transformation from known output
% to input space -- a backward map.
%

T = zeros(8,8);

T(1:4,1) = xo;  T(1:4,2) = yo;  T(1:4,3) = 1;
T(1:4,7) = -xo.*xi;  T(1:4,8) = -yo.*xi;		
T(5:8,4) = xo;  T(5:8,5) = yo;  T(5:8,6) = 1;
T(5:8,7) = -xo.*yi;  T(5:8,8) = -yo.*yi;		

d = [xi; yi];

p = T\d;

[x,y] = meshgrid( [1:nx], [1:ny]);

xmap = ( p(1)*x + p(2)*y + p(3) ) ./ ( p(7)*x + p(8)*y + 1 );
ymap = ( p(4)*x + p(5)*y + p(6) ) ./ ( p(7)*x + p(8)*y + 1 );

% Compute bilinear interpolation
%--------------------------------
warpedIm = interp2( x, y, inputim, xmap, ymap, 'linear');
[i2bad,j2bad] = find( isnan(warpedIm));
for n = 1:length(i2bad)
  warpedIm(i2bad(n),j2bad(n)) = 0.0;
end


