% create a one dimensional ``pan'' movie and compute how much it shifted over
% each time, and put it back together
% motion from fft (fast conv)

L = 1000;  % length whole big image
%e = rand(1,L); % original image
%E = [9 -8 7 -6 5 -4 3 -2 1];
%E(length(E)+1:L) = zeros(1,L-length(E));
E = rand(1,L);  % white to start with
E = E-mean(E);
fractal_exponent = 1;
not_zero = 1;
f = (0+not_zero:L/2-1+not_zero);  %  1/0 or 1/eps first too big
f = [f(L/2:-1:1) f];                % fftshifted order; symmetric
E = E./f.^fractal_exponent;
E = fftshift(E);  % bring back to unfftshifted order
e = real(ifft(E));
plot(e);

% make the movie, which will be 2d because image is 1d
M = 20;
N = 7;
movie = NaN*ones(M,N);
panvector = [1 1 3 4 4 6 7 9 12 12 13 15 15 16 17 18 20 20 20 21 23 23 24]; % hypothetical pans
for m = 1:M
  movie(m,:) = e(panvector(m):panvector(m)+N-1);
end%for

%-------------- BEGIN UNWRAPPING --------------------------------
qmax = NaN*ones(1,M-1);  % motion for each pair of frames
for m = 1:M-1  % there are M-1 pairs of movies
  a = movie(m,:);
  b = movie(m+1,:);
  A = fft(a);
  B = fft(b);
  BA = conj(B).*A;       % conj is like reversal in time domain
  Q = BA(:)./sqrt(A(:).*conj(A(:)).*B(:).*conj(B(:)));
  q = ifft(Q);   % ``histogram'' of shifts present
  [y qmax(m)] = max(q);  % strongest shift present; index 1 = shift of 0
end%for

e = NaN*ones(1,1000);  % re initilize e so we dont know what it was
e(1:N) = movie(1,:);   % starting at index 1
for m = 1:M-1
  e(qmax(m):qmax(m+1)) = movie(m+1,:);  % overwrites, takes most recent not median
end%for
plot(e)