%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  A function for video stabilization     
%  
%  Hao Jiang, Oct, 2007
%    
%  Notes: The code uses Dr. David Lowe's SIFT feaure
%         and the matching function in the demo code.
%         Execute compfeature(.) first to 
%         extract the features in each video frame. 
%         The function assumes that the video is in ./video
%         and the features are in ./data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function stabled(start_frame_num, end_frame_num)

for vframe = start_frame_num : 1 : end_frame_num
 T = zeros(3,3);
 count = 0;
 vframe
 % use a sliding window of length 31
 for tframe = vframe-15 : 1 : vframe+15 
  count = count + 1;
  filename1 = sprintf('./video/%d.jpg', vframe);
  filename2 = sprintf('./video/%d.jpg', tframe);

  % the sift features are pre-computed
  cmd = sprintf('load ./data/des_%d.mat', vframe);
  eval(cmd);
  des1 = des;
  cmd = sprintf('load ./data/des_%d.mat', tframe);
  eval(cmd);
  des2 = des;
  cmd = sprintf('load ./data/loc_%d.mat', vframe);
  eval(cmd);
  loc1 = loc;
  cmd = sprintf('load ./data/loc_%d.mat', tframe);
  eval(cmd);
  loc2 = loc;

  % match the feature points
  distRatio = 0.6;   

  des2t = des2';                          % Precompute matrix transpose
  for i = 1 : size(des1,1)
     dotprods = des1(i,:) * des2t;        % Computes vector of dot products
     [vals,indx] = sort(acos(dotprods));  % Take inverse cosine and sort results

     if (vals(1) < distRatio * vals(2))
       match(i) = indx(1);
     else
       match(i) = 0;
     end
  end

  inp = [];
  outp = [];
  for i = 1: size(des1,1)
    if (match(i) > 0)
      inp = [inp; loc1(i,2), loc1(i,1)];
      outp = [outp; loc2(match(i),2), loc2(match(i),1)];
    end
  end

  % compute the similarity transformation
  tt = cp2tform(inp, outp, 'linear conformal');
  L = tt.tdata.T;
  scale = sqrt(L(1,1)^2 + L(1,2)^2);
  L(:,1:2) = L(:,1:2)/scale;
  % sum the transforms up 
  T = T + L;
 end

  % get the average transformation
  T = T/count;

  % warping by backprojection.
  % inv(T) is the transform from target back to the original image 
  T = inv(T); 

  im = imread(filename1);
  im = im2double(im);
  width = size(im, 2);
  height = size(im, 1);
  % x and y are coordinates on the warped image
 [x,y] = meshgrid(1:width, 1:height);
  nxy = [x(:), y(:), ones(length(x(:)), 1)] * T;
  % newx and newx are corresponding point coordinates in the original image
  newx = reshape(nxy(:,1), height, width);
  newy = reshape(nxy(:,2), height, width);

  im1 = imread(filename1);
  im1 = im2double(im1);
  % get the value of each pixel in the warp image by bi-linear interpolation
  imt1 = interp2(im(:,:,1), newx, newy, 'linear');
  imt2 = interp2(im(:,:,2), newx, newy, 'linear');
  imt3 = interp2(im(:,:,3), newx, newy, 'linear');
  imt = cat(3, imt1, imt2, imt3);
  imshow(imt);
  outname = sprintf('%d.bmp', vframe);
  imwrite(imt, outname);
end