function q = dummy(a,b,exponent)                         Steve Mann, 1991

% motion from fft2       2-D version
% fast conv; see how much image is shifted over
% motion_estimation by fft methods: cross spectrum, CROSS SPECTRUM
%
% Q = est_trans_fft2(A,B); % estimates motion
%
% Q = est_trans_fft2(A,B,.5);  % uses alternate exponent in the division by
%                              % the magnitude
%
% exponent should be chosen from 0 (autocorrelation) to 1 (phase only)
%
% experiments on using [trash,m,n]=max2(Q) as a motion estimator show that:
% exponents from 0 to .3 give the same wrong answer
% exponents from .5 to 1 generally give the same right answer
%                .4 gave one wrong answer; .45 another wrong answer
%                .49 gave right answer
%
% Example application: p_trans=est_trans_fft2(A,B);
%                      pchirp2(pinverse(p_trans),B); % will look alot like A
%                      p=est_pchirp2(A,B,pinverse(p_trans)); % starting guess
%                      B=pchirp2(p,B);  % B almost exactly like A by now
%
% Bugs: RUNS REALLY SLOW IF THERE IS NaN IN DATA (due to bug in Matlab's FFT)
%       workaround by using: est_trans_fft2(A,chnan(B,0)), or the like

if nargin < 3
  disp('est_trans_fft2: exponent defaulting to 1, so we have phase only')
  exponent = 1;
end%if

%disp('est_trans_fft2: about to do fft2 of each input')
A = fft2(a);
B = fft2(b);
% conj is like reversal in time domain
% Q = conj(B).*A./sqrt(A.*conj(A).*B.*conj(B));

if exponent == 1  % Girod method
  disp('est_trans_fft2: about to mult FFTs and div by magnitudes (phase only)')
  Q = conj(B).*A./(abs(A).*abs(B));
end%if

if exponent == 0  % Autocorrelation method
  disp('est_trans_fft2: about to mult FFTs (equiv to compute autocorrelation function)')
  Q = conj(B).*A;
end%if

if exponent == 2 % this is too much; what is the meaning???
		 % exponent should normally be 0 to 1 to be meaningful
  disp('exponent should normally be between 0 and 1 but i will let you try 2')
  disp('est_trans_fft2: about to mult conj and div by power spectra (too much div.)')
  Q = conj(B).*A./(A.*conj(A).*B.*conj(B));
end%if

if (exponent ~= 2) & (exponent ~= 1) & (exponent ~= 0)
  disp(sprintf('est_trans_fft2: about to mult by conj and div by magnitudes^%g',exponent))
  Q = conj(B).*A./(abs(A).*abs(B)).^exponent;
end%if

if nargout > 0
  %disp('est_trans_fft2: about to do inv fft')
  q = ifft2(Q);     % q should be real (e.g. imag part should be nearly 0)
  q = real(q);
end%if

if nargout == 0
  disp('est_motion_fft2: no output arg. given so displaying images to screen')
  %dispmatb(log(real(q+max2(q)/1000000)))  % display as log image

  disp('est_trans_fft2: about to take real part (imag part should be small)')
  q = ifft2(Q);     % q should be real (e.g. imag part should be nearly 0)
  q = real(q);

  clg
  subplot(221)
  tvs(a)
  title('Reference frame')
  subplot(222)
  tvs(b)
  title('Comparison Frame')
  subplot(223)
  tvs(q)
  title('est_trans_fft2')
  subplot(224)
  disp('est_trans_fft2: about to display fftshifted version')
  % log scale doesn't help         tvs(log(q+max2(q)/1000000))
  tvs(fftshift(q))
  title('est_trans_fft2 fftshifted')
end%if