clear
imbase = ['../text/images/'];
f = imread([imbase,'engineer.tif']); f = double(f);

[nr nc] = size(f);
x = linspace(-0.5,0.5,nc)'; 
y1 = linspace(-0.5,0.5,nr)';
y = flipud(y1);

figure(1);clf;colormap gray

imagesc(x,y,f);axis square;axis xy
print -djpeg im1.jpg


[X,Y] = meshgrid(x,y);

% Rotation
th = pi/4;
P = X*cos(th) + Y*sin(th);
Q = -X*sin(th) + Y*cos(th);
g = interp2(X,Y,f,P,Q,'*linear',0);
imagesc(x,y,g);colormap gray
axis square ; axis xy
print -djpeg im_rot.jpg

% Translation
P = X-0.2;
Q = Y-0.4;
g = interp2(X,Y,f,P,Q,'*linear',0);
imagesc(x,y,g);colormap gray
axis square ; axis xy
print -djpeg im_tr.jpg

% Scaling
P = 2*X;
Q = 2*Y;

g = interp2(X,Y,f,P,Q,'*linear',0);
imagesc(x,y,g);colormap gray
axis square ; axis xy
print -djpeg im_sc.jpg

% Nonuniform scaling
R = sqrt(X.^2 + Y.^2);
S = exp(-1.5*R);
P = S.*X;
Q = S.*Y;
g = interp2(X,Y,f,P,Q,'*linear',0);
imagesc(x,y,g);colormap gray
axis square ; axis xy
print -djpeg im_nusc.jpg

% Nonuniform scaling2
R = sqrt(X.^2 + Y.^2);
S = R
P = S.*X;
Q = S.*Y;
g = interp2(X,Y,f,P,Q,'*linear',0);
imagesc(x,y,g);colormap gray
axis square ; axis xy
print -djpeg im_nusc2.jpg

% Nonuniform Rotation
R = sqrt(X.^2 + Y.^2);
th = 5*R;
P = X.*cos(th) + Y.*sin(th);
Q = -X.*sin(th) + Y.*cos(th);
g = interp2(X,Y,f,P,Q,'*linear',0);
imagesc(x,y,g);colormap gray
axis square ; axis xy
print -djpeg im_nurot.jpg


R = sqrt(X.^2 + Y.^2);
th = 5*exp(-10*R);
P = X.*cos(th) + Y.*sin(th);
Q = -X.*sin(th) + Y.*cos(th);
g = interp2(X,Y,f,P,Q,'*linear',0);
imagesc(x,y,g);colormap gray
axis square ; axis xy
print -djpeg im_nurot.jpg


return

% Movie of rotate and scale 
th = linspace(0,2*pi,250);
st = linspace(1,3,250);
for idx = 1:length(th)
  thi = th(idx);
  P = X*cos(thi) + Y*sin(thi);
  Q = -X*sin(thi) + Y*cos(thi);
  P = st(idx)*P;
  Q = st(idx)*Q;
  g = interp2(X,Y,f,P,Q,'*linear',0);
  imagesc(x,y,g);
  colormap gray; axis square ; axis xy
  M(idx) = getframe;
end

% Movie of nonuniform rotate and scale 
th = linspace(0,2*pi,100);
st = linspace(1,1.1,100);
R = 1./(1+sqrt(X.^2 + Y.^2));
R = R/max(max(R));
for idx = 1:length(th)
  thi = th(idx);
  P = X.*cos(thi*R) + Y.*sin(thi*R);
  Q = -X.*sin(thi*R) + Y.*cos(2*thi*R);
  P = st(idx)*P;
  Q = st(idx)*Q;
  g = interp2(X,Y,f,P,Q,'*linear',0);
  imagesc(x,y,g);
  colormap gray; axis square ; axis xy
  M(idx) = getframe;
end