% est_pchirp1direct                                      Steve Mann, 1991-1993
%
% 1D estimate the 3 pchirp parameters from a pair of 1d ``images''
%
% The model is:   u = (ax+b)/(cx+1)-x
%
% Use: p = est_affine1(g1,g2);
%  or: [p1,p2] = est_affine1(g1,g2);
%  or : est_1(g1,g2);            % nicely formatted screen output
%
% properly handles situation when one or both images contain some NaN points

function [p1,p2,p3]=f(E1,E2); % 1 or 2 return args

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

[M,N]= size(E1);
if min(M,N) ~= 1
  error('est_affine1: first ``image'' must be 1-D (function of 1 variable)')
end%if
[M2,N2]= size(E2);
if min(M2,N2) ~= 1
  error('est_affine1: second ``image'' must be 1-D (function of 1 variable)')
end%if
if (M~=M2)&(N~=N2)
  error('est_affine1: ``images'' must be the same length')
end%if

MN=max(M,N);  % image length

% below both need to have same length, hence first is cropped (last element)
El = E2(1:MN-1) - E1(1:MN-1);   % temporal derivative, cropped 1 off
Em = (diff(E1) + diff(E2))/2;   % (average) spatial direction derivative

% define variable ``x'' and make sure El and Em are also row vectors like ``x''
m=1:MN-1;                       % spatial variable (``x''), same length as El,Em
El=El(:).';
Em=Em(:).';

DERIVATIVES = ... % must indent; if you do not indent...
 [...
  [sum2nan(m.*Em .*            m.*Em) sum2nan(Em .*            m.*Em) sum2nan((m.*El-m.*m.*Em) .*            m.*Em)];...
  [sum2nan(m.*Em .*               Em) sum2nan(Em .*               Em) sum2nan((m.*El-m.*m.*Em) .*               Em)];...
  [sum2nan(m.*Em .* (m.*El-m.*m.*Em)) sum2nan(Em .* (m.*El-m.*m.*Em)) sum2nan((m.*El-m.*m.*Em) .* (m.*El-m.*m.*Em))]...
 ];

Derivatives = ...
   [sum2nan((m.*Em-El) .*            m.*Em); ...       %%%  (pork)   [bacon]    [[bacon]   [deletem]e]  [deleteme]
    sum2nan((m.*Em-El) .*               Em); ...
    sum2nan((m.*Em-El) .* (m.*El-m.*m.*Em)) ...
   ];                         % vector of derivatives

parameters = DERIVATIVES\Derivatives;

if (nargout == 3) | (nargout == 0)
  % p1 = parameters(1);           % if nargin==0 must never make ref to p1
  p2 = parameters(2);             % or you end up with ans=[], when called with
  p3 = parameters(3);             % no semicolon
end%if
if nargout == 3
  p1 = parameters(1);
end%if

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

if nargout == 0   % no output argument
  disp('est_affine1: 0 output arguments given so I am displaying to screen:')
  disp(sprintf('%g %g  %g',parameters(1),p2,p3))
end%if

if (nargout ~= 0) & (nargout ~= 1) & (nargout ~= 3)
  error('est_affine1: you must have 0, 1, or 3 output arguments')
end%if