% est_psuedoperspective                                Steve Mann, 1992-1993
%
% estimate the 9 Wyckoff pseudoperspective parameters from a pair of images
%
% The model is:   u = a + bx + cy   +   gx^2 + hxy
%                 v = d + ex + fy   +   gxy  + hy^2
%
% Use: p = est_pseudo(A,B);
%  or: [p1,p2,p3,p4,p5,p6,p7,p8] = est_pseudo(A,B);
%  or : est_pseudo(A,B);      % nicely formatted screen output
%
% Example application: pb=est_pseudow(A,B);
%                      pp=pseudow2pw(pb);
%                      B2=pchirp2w(pp,B);      % now B2 will look alot like A
%                      pb=est_pseudow(A,B2); % residual difference
%                      pp=pseudow2pw(pb);
%                      B3=pchirp2w(pp,B);      % 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
%
% See also: pseudo2pw, pchirp2w

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

%%%%mask = isnan(E1).*isnan(E2);  % NaN mask: regions where either is undefined
% doesnt work because nanmask makes it still nan

[El,Em,En] = derivatives(E1,E2);
[M,N]= size(El);                    % size of any of the 3 derivative matrices
%%% this might be wrong, but basing on Em and En might also be wrong ...   MNnotnan = M*N - sum(sum(isnan(El));      % number of pixels that are not nan
MNnotnan = M*N - sum(sum(isnan(Em)));      % number of pixels that are not nan
[n m]=meshdom(1:N,M:-1:1);          % variables ``x'' and ``y''
%disp('est_pseudo: converting from pixels to units on [0,1)')      doesnt work
%m=m/M;    doesnt work
%n=n/N;    doesnt work

D1=[sum2nan((m.*m.*Em+m.*n.*En) .*(m.*m.*Em+m.*n.*En)) sum2nan(m.*Em .*(m.*m.*Em+m.*n.*En)) sum2nan(n.*Em .*(m.*m.*Em+m.*n.*En)) sum2nan(Em .*(m.*m.*Em+m.*n.*En))...
    sum2nan((m.*n.*Em+n.*n.*En) .*(m.*m.*Em+m.*n.*En)) sum2nan(m.*En .*(m.*m.*Em+m.*n.*En)) sum2nan(n.*En .*(m.*m.*Em+m.*n.*En)) sum2nan(En .*(m.*m.*Em+m.*n.*En)) MNnotnan];
D2=[sum2nan((m.*m.*Em+m.*n.*En) .*m.*Em)               sum2nan(m.*Em .*m.*Em)               sum2nan(n.*Em .*m.*Em)               sum2nan(Em .*m.*Em)...
    sum2nan((m.*n.*Em+n.*n.*En) .*m.*Em)               sum2nan(m.*En .*m.*Em)               sum2nan(n.*En .*m.*Em)               sum2nan(En .*m.*Em)               MNnotnan];
D3=[sum2nan((m.*m.*Em+m.*n.*En) .*n.*Em)               sum2nan(m.*Em .*n.*Em)               sum2nan(n.*Em .*n.*Em)               sum2nan(Em .*n.*Em)...
    sum2nan((m.*n.*Em+n.*n.*En) .*n.*Em)               sum2nan(m.*En .*n.*Em)               sum2nan(n.*En .*n.*Em)               sum2nan(En .*n.*Em)               MNnotnan];
D4=[sum2nan((m.*m.*Em+m.*n.*En) .*Em)                  sum2nan(m.*Em .*Em)                  sum2nan(n.*Em .*Em)                  sum2nan(Em .*Em)...
    sum2nan((m.*n.*Em+n.*n.*En) .*Em)                  sum2nan(m.*En .*Em)                  sum2nan(n.*En .*Em)                  sum2nan(En .*Em)                  MNnotnan];
D5=[sum2nan((m.*m.*Em+m.*n.*En) .*(m.*n.*Em+n.*n.*En)) sum2nan(m.*Em .*(m.*n.*Em+n.*n.*En)) sum2nan(n.*Em .*(m.*n.*Em+n.*n.*En)) sum2nan(Em .*(m.*n.*Em+n.*n.*En))...
    sum2nan((m.*n.*Em+n.*n.*En) .*(m.*n.*Em+n.*n.*En)) sum2nan(m.*En .*(m.*n.*Em+n.*n.*En)) sum2nan(n.*En .*(m.*n.*Em+n.*n.*En)) sum2nan(En .*(m.*n.*Em+n.*n.*En)) MNnotnan];
D6=[sum2nan((m.*m.*Em+m.*n.*En) .*m.*En)               sum2nan(m.*Em .*m.*En)               sum2nan(n.*Em .*m.*En)               sum2nan(Em .*m.*En)...
    sum2nan((m.*n.*Em+n.*n.*En) .*m.*En)               sum2nan(m.*En .*m.*En)               sum2nan(n.*En .*m.*En)               sum2nan(En .*m.*En)               MNnotnan];
D7=[sum2nan((m.*m.*Em+m.*n.*En) .*n.*En)               sum2nan(m.*Em .*n.*En)               sum2nan(n.*Em .*n.*En)               sum2nan(Em .*n.*En)...
    sum2nan((m.*n.*Em+n.*n.*En) .*n.*En)               sum2nan(m.*En .*n.*En)               sum2nan(n.*En .*n.*En)               sum2nan(En .*n.*En)               MNnotnan];
D8=[sum2nan((m.*m.*Em+m.*n.*En) .*En)                  sum2nan(m.*Em .*En)                  sum2nan(n.*Em .*En)                  sum2nan(Em .*En)...
    sum2nan((m.*n.*Em+n.*n.*En) .*En)                  sum2nan(m.*En .*En)                  sum2nan(n.*En .*En)                  sum2nan(En .*En)                  MNnotnan];
D9=[sum2nan((m.*m.*Em+m.*n.*En) )                      sum2nan(m.*Em  )                     sum2nan(n.*Em  )                     sum2nan(Em  )...
    sum2nan((m.*n.*Em+n.*n.*En) )                      sum2nan(m.*En  )                     sum2nan(n.*En  )                     sum2nan(En  )                     MNnotnan];

DERIVATIVES = [D1;D2;D3;D4;D5;D6;D7;D8;D9];

Derivatives = ...
 - [sum2nan(El .*(m.*m.*Em+m.*n.*En)); ...     % g
    sum2nan(El .*m.*Em); ...                   % b
    sum2nan(El .*n.*Em); ...                   % c
    sum2nan(El .*Em); ...                      % a
    sum2nan(El .*(m.*n.*Em+n.*n.*En)); ...     % h
    sum2nan(El .*m.*En); ...                   % e
    sum2nan(El .*n.*En); ...                   % f
    sum2nan(El .*En); ...                      % d
    sum2nan(El  ) ...                          % additive gain
   ];                         % vector of derivatives

parameters = DERIVATIVES\Derivatives;

disp('est_pseudo: adding identity and converting from pel to norm''d [0,1) rep.')
parameters(:) = parameters(:) + [0;1;0;0; 0;0;1;0;0];  % convert relative change like u=func(x+deltax) to absolute like u=func(x)
not_sure=1;
parameters = parameters(:) .* [M;1;N/M;1/M; N;M/N;1;1/N; not_sure];

if (nargout == 9) | (nargout == 0)
  % p1 = parameters(4);     % 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(8);
  p5 = parameters(6);
  p6 = parameters(7);
  p7 = parameters(1);
  p8 = parameters(5);
  additive_gain = parameters(9);
end%if
if nargout == 9
  p1 = parameters(4);
end%if

if nargout == 1
  p1 = parameters([4 2 3 8 6 7 1 5 9]);
end%if

if nargout == 0
  disp('est_pseudo: 0 output args given; displaying result to screen:')
  disp(sprintf('%g %g %g %g %g %g %g %g %g',parameters(1),p2,p3,p4,p5,p6,p7,p8,additive_gain))
end%if

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