% ===============================================================================
% Tutorial for FAIR: distances, hands, SSD versus rotation (extended version of E7_basic.m)
% (c) Jan Modersitzki 2009/03/25, see FAIR.2 and FAIRcopyright.m.
% note: '%$' is used so hide comments in the book                              
%
%   - data                 PETCT, Omega=(0,140)x(0,151), level=4:7, m=[128,128]
%   - viewer               viewImage2D
%   - interpolation        splineInter2D
%   - distance             'SSD','NCC','MI','NGF'
%   - transformation       rotation2D
% ===============================================================================

clear, close all, help(mfilename);

setupPETCTdata; level = 6; omega = MLdata{level}.omega; m = MLdata{level}.m;
viewImage('reset',viewOptn{:},'axis','off');
inter('reset','inter','splineInter2D','regularizer','moments','theta',1e0);
[T,R] = inter('coefficients',MLdata{level}.T,MLdata{level}.R,omega,'out',0);
distance('reset','distance','SSD');
center = (omega(2:2:end)-omega(1:2:end))'/2;
trafo('reset','trafo','rotation2D','c',center);
fprintf('%20s : %s\n','viewImage',viewImage);
fprintf('%20s : %s\n','inter',inter);
fprintf('%20s : %s\n','distance',distance);
fprintf('%20s : %s\n','trafo',trafo);

xc = getCenteredGrid(omega,m);
Rc = inter(R,omega,xc);

wc = linspace(-pi/2,pi/2,51);
Dc = zeros(size(wc));

distances = {'SSD','NCC','MI','NGF'};
diffImage = @(Tc) viewImage(128+(Tc-Rc)/2,omega,m);

for k=1:length(distances),

  distance('reset','distance',distances{k});
  fprintf('---- %s\n\n',distance);
  
  Dc = zeros(size(wc));
  % adapt visualization functions
  switch distance,
    case 'SSD',
      resStr = '|T(yc)-R|';
      res = @(rc) 128+rc/2;
      plotRes = @(rc) viewImage(res(rc),omega,m);
      setRes  = @(ph,rc) set(ph,'cdata',reshape(res(rc),m)');
    case 'NCC',
      resStr = '|T(yc)-R|';
      res = @(rc) 128+rc/2;
      plotRes = @(rc) viewImage(res(rc),omega,m);
      setRes  = @(ph,rc) set(ph,'cdata',reshape(res(rc),m)');
    case 'MI',
      %set-up default parameter
      distance('set','tol',1e-7,'minT',0,'maxT',256,'nT',60,'minR',0,'maxR',256,'nR',60);
      distance('disp');
      nT = distance('get','nT') + 4;
      nR = distance('get','nR') + 4;
      ax = [wc(1),wc(end),50,60];
      resStr = '\rho(|T(yc),R)';
      res = @(rc) reshape((1+min(rc,0.001)),nT,nR);
      plotRes = @(rc) imagesc(res(rc));
      setRes  = @(ph,rc) set(ph,'cdata',res(rc)');
    case 'NGF',
      %set-up default parameter
      distance('set','edge',100);
      distance('disp');
      ax = [wc(1),wc(end),370,400];
      resStr = 'NGF(T(yc),R)';
      res = @(rc) reshape(rc,m);
      plotRes = @(rc) viewImage2Dsc(res(rc),omega,m);
      setRes  = @(ph,rc) set(ph,'cdata',res(rc)');
  end;

  % run the loop over all rotations
  for j = 1:length(wc),
    yc = trafo(wc(j),xc);
    Tc = inter(T,omega,yc);
    [dc,rc] = distance(Tc,Rc,omega,m);
    Dc(j) = dc;

    if j == 1,
      FAIRfigure(k,'figname',mfilename); clf;
      subplot(2,2,1);  viewImage(Rc,omega,m);          th(1) = title('R');
      subplot(2,2,2);  vh = viewImage(Tc,omega,m);     th(2) = title('T(yc)');
      subplot(2,2,3);  rh = plotRes(rc); axis image;   th(3) = title(resStr);
      subplot(2,2,4);  ph = plot(wc(1),Dc(1),'r.','markersize',markersize);
      th(4) = title(sprintf('%s versus rotation',distance));
      axis([wc(1),wc(end),-inf,inf]); hold on;
      set(th,'fontsize',fontsize);
      FAIRpause;

    else
      figure(k);
      set(vh,'cdata',reshape(Tc,m)')
      setRes(rh,rc);
      subplot(2,2,4);
      set(ph,'visible','off');
      plot(wc(1:j),Dc(1:j),'k-','linewidth',2);
      ph = plot(wc(j),Dc(j),'r.','markersize',markersize);
      drawnow, FAIRpause(1/100)
    end;
    fprintf('.');
    if ~rem(j,50) || j == length(wc), fprintf('\n'); end;
  end;
  FAIRpause;
end;
