% est_pchirp2w                                         Steve Mann, 1992-1993
%
% estimate 9 pchirp2w (8coord. transf, and 1 ag) parameters from pair of images
%
% The model is:   u = W((ax+by+c)/(gx+hy+1))
%                 v = W((dx+ey+f)/(gx+hy+1))
%
% Use: p = est_pchirp2w(A,B);
%  or: [p1,p2,p3,p4,p5,p6,p7,p8] = est_pchirp2(A,B);
%  or : est_pchirp2(A,B);         % nicely formatted screen output
%
% Example application: p=est_pchirp2w(A,B1);
%                      B2=pchirp2w(p,B1);      % now B2 will look alot like A
%                      p=est_pchirp2w(A,B2);   % residual difference
%                      B3=pchirp2w(p,B);       % B3 almost exactly like A by now
%
% sign convention: m (or X) down, n (or Y) to right
% properly handles situation when one or both images contain some NaN points
%
%
%
% Calls: est_pseudow, pseudo2p, which calls corners2r, etc

function [p1,p2,p3,p4,p5,p6,p7,p8,p9]=f(E1,E2);  % two images in, 0,1,or 9 outs

[M, N]  = size(E1);
[M2,N2] = size(E2);
if (M2~=M) | (N2~=N)
  error(...
    sprintf('pchirp2: images must be same size; %gby%g and %gby%g differ',...
            M,N,M2,N2)...
       )
end%if

%q = est_bilinear(E1,E2);
%parameters = bilinear2p(q);

%q = est_qchirp2(E1,E2);
%parameters = qchirp22p(q);

q = est_pseudow(E1,E2);
parameters = pseudo2p(q(1:8)); parameters(9)=q(9);

if (nargout == 9) | (nargout == 0)
  % p1 = parameters(1);   % if called with no semicolon, do not want to see
                          % any ans=[], etc, so make no reference to p1
  p2 = parameters(2);
  p3 = parameters(3);
  p4 = parameters(4);
  p5 = parameters(5);
  p6 = parameters(6);
  p7 = parameters(7);
  p8 = parameters(8);
  p9 = parameters(9);
end%if
if nargout == 9
  p1 = parameters(1);
end%if

if nargout == 1
  p1 = parameters(:);
end%if

if nargout == 0
  disp('est_pchirp2: 0 output args given; displaying result in p and correspondence:')
  disp(sprintf('p_vec: %+10.6f %+10.6f %+10.6f    %+10.6f %+10.6f %+10.6f',...
               parameters(1),p2,p3,p4,p5,p6))
  disp(sprintf('       %+10.6f                          %+10.6f   ag = %+9.3f',p7,p8,p9))
  x=[0; 0; 1; 1];
  y=[0; 1; 0; 1];
  one=ones(4,1);
  xp=(parameters(1)*x+p2*y+p3)./(p7*x+p8+one);
  yp=           (p4*x+p5*y+p6)./(p7*x+p8+one);
  disp1=sprintf('4 cor: %+9.1f/%g %+9.1f/%g %+9.1f/%g %+9.1f/%g',...
                M*xp(1),M,M*xp(2),M,M*xp(3),M,M*xp(4),M);
  disp([disp1])
  disp1=sprintf('       %+9.1f/%g %+9.1f/%g %+9.1f/%g %+9.1f/%g',...
                N*yp(1),N,N*yp(2),N,N*yp(3),N,N*yp(4),N);
  disp([disp1])
end%if

if (nargout ~=0) & (nargout ~= 1) & (nargout ~= 9)
  error('est_pchirp2w: number of output arguments must be 0, 1, or 9')
end%if