% est_dozen                                          Steve Mann, 1992-1994
% (not called est_dozen2 because the 2 is obvious from number of parameters)
%
% estimate the 12 Taylor series of pchirp parameters from a pair of images
%
% The model is:   u = p1 + p2x + p3y +  p4x^2 +  p5y^2 +  p6xy
%                 v = p7 + p8x + p9y + p10x^2 + p11y^2 + p12xy
%
% Use: p = est_dozen(E1,E2);
%  or: [p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12] = bilinear_est(A,B);
%  or : est_dozen(E1,E2);      % nicely formatted screen output
%
% relation to pchirp2: p1=-(ah+bg); p2=(a-cg); p3=(b-ch); p4=c;
%                      p5=-(dh+eg); p6=(d-fg); p7=(e-fh); p8=f;
% sign convention: m (or X) down, n (or Y) to right
% properly handles situation when one or both images contain some NaN points

function [p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12]=f(E1,E2); % 1 or 12 returns

%%%%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
[n m]=meshdom(1:N,M:-1:1);          % variables ``x'' and ``y''

DERIVATIVES = ... % must indent; if you do not indent then the )...sum2nan will not work, you need ) ...sum2nan (e.g. the space)
 [...
  [sum2nan(Em       .*Em) sum2nan(m.*Em       .*Em) sum2nan(n.*Em       .*Em) sum2nan(m.*m.*Em       .*Em) sum2nan(n.*n.*Em       .*Em) sum2nan(m.*n.*Em      .*Em)...
   sum2nan(En ) sum2nan(m.*En ) sum2nan(n.*En ) sum2nan(m.*m.*En ) sum2nan(n.*n.*En ) sum2nan(m.*n.*En )];...
  [sum2nan(Em ) sum2nan(m.*Em ) sum2nan(n.*Em ) sum2nan(m.*m.*Em ) sum2nan(n.*n.*Em ) sum2nan(m.*n.*Em )...
   sum2nan(En ) sum2nan(m.*En ) sum2nan(n.*En ) sum2nan(m.*m.*En ) sum2nan(n.*n.*En ) sum2nan(m.*n.*En )];...
  [sum2nan(Em ) sum2nan(m.*Em ) sum2nan(n.*Em ) sum2nan(m.*m.*Em ) sum2nan(n.*n.*Em ) sum2nan(m.*n.*Em )...
   sum2nan(En ) sum2nan(m.*En ) sum2nan(n.*En ) sum2nan(m.*m.*En ) sum2nan(n.*n.*En ) sum2nan(m.*n.*En )];...
  [sum2nan(Em ) sum2nan(m.*Em ) sum2nan(n.*Em ) sum2nan(m.*m.*Em ) sum2nan(n.*n.*Em ) sum2nan(m.*n.*Em )...
   sum2nan(En ) sum2nan(m.*En ) sum2nan(n.*En ) sum2nan(m.*m.*En ) sum2nan(n.*n.*En ) sum2nan(m.*n.*En )];...
  [sum2nan(Em ) sum2nan(m.*Em ) sum2nan(n.*Em ) sum2nan(m.*m.*Em ) sum2nan(n.*n.*Em ) sum2nan(m.*n.*Em )...
   sum2nan(En ) sum2nan(m.*En ) sum2nan(n.*En ) sum2nan(m.*m.*En ) sum2nan(n.*n.*En ) sum2nan(m.*n.*En )];...
  [sum2nan(Em ) sum2nan(m.*Em ) sum2nan(n.*Em ) sum2nan(m.*m.*Em ) sum2nan(n.*n.*Em ) sum2nan(m.*n.*Em )...
   sum2nan(En ) sum2nan(m.*En ) sum2nan(n.*En ) sum2nan(m.*m.*En ) sum2nan(n.*n.*En ) sum2nan(m.*n.*En )];...
  [sum2nan(Em ) sum2nan(m.*Em ) sum2nan(n.*Em ) sum2nan(m.*m.*Em ) sum2nan(n.*n.*Em ) sum2nan(m.*n.*Em )...
   sum2nan(En ) sum2nan(m.*En ) sum2nan(n.*En ) sum2nan(m.*m.*En ) sum2nan(n.*n.*En ) sum2nan(m.*n.*En )];...
  [sum2nan(Em ) sum2nan(m.*Em ) sum2nan(n.*Em ) sum2nan(m.*m.*Em ) sum2nan(n.*n.*Em ) sum2nan(m.*n.*Em )...
   sum2nan(En ) sum2nan(m.*En ) sum2nan(n.*En ) sum2nan(m.*m.*En ) sum2nan(n.*n.*En ) sum2nan(m.*n.*En )];...
  [sum2nan(Em ) sum2nan(m.*Em ) sum2nan(n.*Em ) sum2nan(m.*m.*Em ) sum2nan(n.*n.*Em ) sum2nan(m.*n.*Em )...
   sum2nan(En ) sum2nan(m.*En ) sum2nan(n.*En ) sum2nan(m.*m.*En ) sum2nan(n.*n.*En ) sum2nan(m.*n.*En )];...
  [sum2nan(Em ) sum2nan(m.*Em ) sum2nan(n.*Em ) sum2nan(m.*m.*Em ) sum2nan(n.*n.*Em ) sum2nan(m.*n.*Em )...
   sum2nan(En ) sum2nan(m.*En ) sum2nan(n.*En ) sum2nan(m.*m.*En ) sum2nan(n.*n.*En ) sum2nan(m.*n.*En )];...
  [sum2nan(Em ) sum2nan(m.*Em ) sum2nan(n.*Em ) sum2nan(m.*m.*Em ) sum2nan(n.*n.*Em ) sum2nan(m.*n.*Em )...
   sum2nan(En ) sum2nan(m.*En ) sum2nan(n.*En ) sum2nan(m.*m.*En ) sum2nan(n.*n.*En ) sum2nan(m.*n.*En )];...
  [sum2nan(Em ) sum2nan(m.*Em ) sum2nan(n.*Em ) sum2nan(m.*m.*Em ) sum2nan(n.*n.*Em ) sum2nan(m.*n.*Em )...
   sum2nan(En ) sum2nan(m.*En ) sum2nan(n.*En ) sum2nan(m.*m.*En ) sum2nan(n.*n.*En ) sum2nan(m.*n.*En )];...
 ];

Derivatives = ...
 - [sum2nan(El       .*Em); ...
    sum2nan(El    .*m.*Em); ...
    sum2nan(El    .*n.*Em); ...
    sum2nan(El .*m.*m.*Em); ...
    sum2nan(El .*n.*n.*Em); ...
    sum2nan(El .*m.*n.*Em); ...
    sum2nan(El       .*En); ...
    sum2nan(El    .*m.*En); ...
    sum2nan(El    .*n.*En); ...
    sum2nan(El .*m.*m.*En); ...
    sum2nan(El .*n.*n.*En); ...
    sum2nan(El .*m.*n.*En) ...
   ];                         % vector of derivatives

parameters = DERIVATIVES\Derivatives;

if nargout == 12
  p1 = parameters(1);
  p2 = parameters(2);
  p3 = parameters(3);
  p4 = parameters(4);
  p5 = parameters(5);
  p6 = parameters(6);
  p7 = parameters(7);
  p8 = parameters(8);
  p9 = parameters(9);
  p10= parameters(10);
  p11= parameters(11);
  p12= parameters(12);
end%if

if nargout == 1
  p1 = parameters;
end%if

if nargout == 0
  disp('est_bilinear: 0 output arguments were given so I am displaying result to screen:')
  disp(sprintf('%g %g %g %g %g %g %g %g %g %g %g %g',p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12))
end%if