%  dip07_1.m        gr23mar15  nov16
%     geometric transformations
clear all;

f = double(imread('.\img\cameraman.tif'));
f = f(1:2:end,1:2:end); %%%% beware
[nr,nc] = size(f);  % in the slides, 'x'=columns, 'y'=rows

%        | a11  a12  b1 |
%    A = | a21  a22  b2 |
%        |  0    0    1 |

 A = [1  0  0;  0  1  0;  0  0  1];  %   nil
%  A = [1.5  0  0;  0  .8  0;  0  0  1];  % anisotropic scaling
%  th = pi/3;   A = [cos(th)  -sin(th)  0;  sin(th)  cos(th)  0;  0  0  1];  % rotation
%  A = [1  .7  0;  0  1  0;  0  0  1];  % shearing,  parallel to rows
%  A = [1  0  90;  0  1  -110;  0  0  1];  % translation

Am1 = inv(A);                  % perform transformation
g = zeros(3*size(f));
for r=-nr+1:2*nr,  for c=-nc+1:2*nc,
	pos = round(Am1 * [c; r; 1]);     %  position of nearest grid-point
	if (sum(pos>[0;0;0])==3  && sum(pos<=[nc;nr;1])==3);  % avoid negative indices
				g(r+nr,c+nc) = f(pos(2),pos(1));
	end;
end; end;
figure(1); imshow(g,[]);
return

%%  ----  determine transformation based on 3 pairs of tiepoints
%legend  [c1i,r1i,1;   c2i,r2i,1;   c3i,r3i,1]';          [c1f,r1f,1;   c2f,r2f,1;   c3f,r3f,1]'; 
 p_ini = [25, 25, 1;  100, 35, 1;  70, 100, 1]';   p_fin = [ 5,  5, 1;  111, 20, 1;  70, 120, 1]';
% p_ini = [8, 20, 1;  100, 10, 1;  85, 120, 1]';   p_fin = [25, 18, 1;  111, 20, 1;  70, 120, 1]';
% p_ini = [8,  8, 1;  100, 20, 1;  85, 120, 1]';   p_fin = [25,  8, 1;  111, 10, 1;  60, 120, 1]';

% or give coordinates:
% figure(2); imshow(f,[]);  
% [c1i, r1i] = ginput(1);   [c1f, r1f] = ginput(1);
% [c2i, r2i] = ginput(1);   [c2f, r2f] = ginput(1);
% [c3i, r3i] = ginput(1);   [c3f, r3f] = ginput(1);
% p_ini = [c1i,r1i,1; c2i,r2i,1; c3i,r3i,1]';   p_fin = [c1f,r1f,1; c2f,r2f,1; c3f,r3f,1]';

figure(2); imshow(f,[]);   text(p_ini(1,:),p_ini(2,:),'\bf\color{blue}x'); 
    text(p_fin(1,:),p_fin(2,:),'\bf\color{blue}o');    pause();

A = p_fin * inv(p_ini);
Am1 = inv(A);
g = zeros(3*size(f));
for r=-nr+1:2*nr,  for c=-nc+1:2*nc,
		pos = round(Am1 * [c; r; 1]);
		if (sum(pos>[0;0;0])==3  && sum(pos<=[nc;nr;1])==3);
					g(r+nr,c+nc) = f(pos(2),pos(1));
		end;
end; end;
figure(1); imshow(g,[]); figure(2);

%% projective (bilinear) transform
% compact 'maketform' and 'imtransform' commands exists:

q = imread('.\img\P101.jpg');
figure(1); imshow(q); pause;
[nr,nc,dum] = size(q);
U = [0 0; nc 0; 0 nr; nc nr];
X = [16 0; nc-11 3; 0 nr; nc nr+8];
T = maketform('projective',U,X);
    % maps the coord. in each row of U to the corresponding ones in X.  
    % U and X define the corners of input and output quadrilaterals
for k=1:3,
    qq = q(:,:,k);
    qqt = imtransform(qq,T);
%    imshow(qqt); drawnow;
    qt(:,:,k) = qqt;
    end;
imshow(qt);

