MakeQTMovie('start','cs448a-hw1.mov');
MakeQTMovie('quality',0.75);
MakeQTMovie('size',[320 240]);
MakeQTMovie('framerate',30);

sig = 2;

% compute gaussian density (95% falls in interval abs(x-mean) < 2*sig
x = floor(-3*sig):ceil(3*sig);
GAUSS = exp(-0.5*x.^2/sig^2); GAUSS = GAUSS/sum(GAUSS);
dGAUSS = -x.*GAUSS/sig^2;

% load images
F = double(imread('SRI/1bw1.tif')); 
[rows,cols] = size(F);
[X,Y] = meshgrid(1:cols,1:rows);

% select roi
clf; imagesc(F); colormap(gray); axis image;
mask = double(roipoly); roi = find(mask);
mx = mean(X(roi)); my = mean(Y(roi));
hold on; plot(mx,my,'rx'); hold off;

u = 0; v = 0;

for t = 2:99,
 G = double(imread(sprintf('SRI/1bw%d.tif',t)));

 Fx = conv2(F,GAUSS'*dGAUSS,'same');
 Fy = conv2(F,dGAUSS'*GAUSS,'same');

% iterate
 for i = 1:10,
  % this is over-kill to do it on the entire image, but what the heck..
  Gwarp = interp2(X,Y,G,X+u,Y+v,'bilinear');
  Ft = conv2(Gwarp-F,GAUSS'*GAUSS,'same');
  
  c1 = sum(Fx(roi).^2);
  c2 = sum(Fx(roi).*Fy(roi));
  c3 = sum(Fy(roi).^2);
  C = [c1 c2;c2 c3]; D = [sum(Ft(roi).*Fx(roi)); sum(Ft(roi).*Fy(roi))];
  
  uv_new = inv(C)*D;
  u = u-uv_new(1);
  v = v-uv_new(2);
  
  fprintf('err = %d\n',sum(Ft(roi).^2)); 
 end;

 imagesc(G); axis image; hold on; plot(mx+u,my+v,'rx'); hold off;
 drawnow;
 MakeQTMovie('addframe');
end;
MakeQTMovie('finish');
