function q = dummy(a,b,exponent) % motion from fft2 2-D version Steve Mann (steve@media.mit.edu) % fast conv; see how much image is shifted over % motion_estimation by fft methods: cross spectrum, CROSS SPECTRUM % % Example use: q=motion_cross_fft2(a,b); % should have a peak where motion is... % q=motion_cross_fft2(a,b,.5); % uses alternate exponent % default exponent is 1 (Girod method) % % 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 % % Bugs: RUNS REALLY SLOW IF THERE IS NaN IN DATA (due to bug in Matlab's FFT) % workaround by using: motion_cross_fft2(A,chnan(B,0)), or the like % % used by est_trans_fft2 % (e.g. if you just want 2 numbers you want to use est_trans_fft2 and not this) if nargin < 3 disp('motion_cross_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); disp('motion_cross_fft2.m: using add-hoc window (found far better than rect)') [M,N]=size(a); window=hanning(M)*hanning(N).'; A = fft2(a.*window); B = fft2(b.*window); %%%%B = fft2(rev(b)); % conj doesnt work quite right in 2D % the phase shift seems to have effect??? % but this should not be a problem % conj is like reversal in time domain % Q = conj(B).*A./sqrt(A.*conj(A).*B.*conj(B)); if exponent == 1 % Girod method disp('motion_cross_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('motion_cross_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('motion_cross_fft: exponent normally be between 0 and 1 but will let you try 2') disp('motion_cross_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('motion_cross_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('motion_cross_fft2: no output arg. given so displaying images to screen') %dispmatb(log(real(q+max2(q)/1000000))) % display as log image disp('motion_cross_fft2: about to take real part (imag part should be small)') q_noans = ifft2(Q); % q should be real (e.g. imag part should be nearly 0) q_noans = real(q_noans); clg subplot(221) tvs(a) title('Reference frame') subplot(222) tvs(b) title('Comparison Frame') subplot(223) tvs(-q_noans) title('motion_cross_fft2 negated') subplot(224) disp('motion_cross_fft2: about to display fftshifted version') % log scale doesn't help tvs(log(q+max2(q)/1000000)) tvs(-fftshift(q_noans)) title('motion_cross_fft2 fftshifted negated') end%if if nargout > 0 q = ifft2(Q); % q should be real (e.g. imag part should be nearly 0) q = real(q); end%if %dispmatb(log(real(q+max2(q)/1000000))) % display as log image