% est_affine1                                          Steve Mann, 1991-1993
%
% 1D estimate the 2 affine parameters from a pair of 1d ``images''
%
% The model is:   u-x = ax + b
%                 u   = ax + b - x
%                     = (a-1)x + b
%
% Use: p = est_affine1(g1,g2);
%  or: [p1,p2] = est_affine1(g1,g2);
%  or : est_affine1(g1,g2);      % nicely formatted screen output
%
%     : est_affine1(g1,g2,0);    % user-specified Ex exponent for `affine fit'
%                                % default is 1 which is for `affine flow'
%
% properly handles situation when one or both images contain some NaN points

function [p1,p2]=f(E1,E2,w); % 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

if nargin==2
  w=1;
  disp('est_affine1: weighting (exponent in Em^2) defaulting to 1 for `flow''')
end%if

if (w~=0)&(w~=1)
  disp('est_affine1: not sure how to handle weights other than 0 or 1')
end%if

[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((Em.*Em).^w .*m.*m) sum2nan((Em.*Em).^w .*m)];...
  [sum2nan((Em.*Em).^w    .*m) sum2nan((Em.*Em).^w    )]...
 ];

Derivatives = ...
 - [sum2nan(El   .*m    .*Em.^(2*w-1)); ...
    sum2nan(El          .*Em.^(2*w-1)) ...
   ];                         % vector of derivatives

parameters = DERIVATIVES\Derivatives;

disp('est_affine: adding identity')
parameters(:) = parameters(:) + [1;0];  % convert relative change like u=func(x+deltax) to absolute like u=func(x)

if (nargout == 2) | (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
                                  % no semicolon
end%if
if nargout == 2
  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',parameters(1),p2))
end%if

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