% function i2 = warpimage(i1, A, rangex, rangey)
%
% warp image i1 by transformation A into i2, pixels with i2-coordinates
% rangex, rangey will be rendered; if rangex, rangey is not specified,
% bounding box of warped i1 will be the range

function i2 = warpimage(i1, A, rangex, rangey)

d = size(i1);

[ x1, y1 ] = meshgrid(1:d(2), 1:d(1));

if nargin < 4
	C1 = [ x1(:)'; y1(:)'; ones(1,d(1)*d(2)) ];
	% transform them by A
	C2 = A * C1;
	C2 = C2 ./ (ones(3,1) * C2(3,:));

	% get coordinates of corners and size of the new image
	corners = [ floor(min(C2')); ceil(max(C2')) ]'
	nd = corners(:,2) - corners(:,1) + 1;
	rangex = corners(1,1):corners(1,2);
	rangey = corners(2,1):corners(2,2);
else
	nd(1) = rangex(end) - rangex(1) + 1;
	nd(2) = rangey(end) - rangey(1) + 1;
end

% where the coordinates in the new image are looking in the old one
[ x2, y2 ] = meshgrid(rangex, rangey);
N1 = [ x2(:)'; y2(:)'; ones(1,nd(1)*nd(2)) ];
N2 = inv(A) * N1;
N2 = N2 ./ (ones(3,1) * N2(3,:));

% use interp2 (interpolating value at float coords between points on 2D grid)
for c=1:size(i1,3)
	i2(:,:,c) = uint8(round(interp2(x1, y1, double(i1(:,:,c)), reshape(N2(1,:),size(x2)), reshape(N2(2,:), size(y2)))));
end
