% Note: This file was taken from David Lowe's website at
% http://www.cs.ubc.ca/~lowe/keypoints/   Parts have been modified as
% indicated.
%
% num = match(image1, image2)
%
% This function reads two images, finds their SIFT features, and
%   displays lines connecting the matched keypoints.  A match is accepted
%   only if its distance is less than distRatio times the distance to the
%   second closest match.
% It returns the number of matches displayed.
%
% Example: match('scene.pgm','book.pgm');

% MODIFIED BY BRANDON SMITH TO RETURN ALL OF THE MATCHED POINT PAIRS
% point_pairs is in the format point_pairs(<pair index>, <component = 
% 1 (x1), 2 (y1), 3 (x2), 4 (y2). Also image1 and image2 were changed
% to be the images themselves, not filenames.
function [num_matches, point_pairs] = match(image1, image2)

% Find SIFT keypoints for each image
[im1, des1, loc1] = sift(image1);
[im2, des2, loc2] = sift(image2);

% For efficiency in Matlab, it is cheaper to compute dot products between
%  unit vectors rather than Euclidean distances.  Note that the ratio of 
%  angles (acos of dot products of unit vectors) is a close approximation
%  to the ratio of Euclidean distances for small angles.
%
% distRatio: Only keep matches in which the ratio of vector angles from the
%   nearest to second nearest neighbor is less than distRatio.
distRatio = 0.6;   

% For each descriptor in the first image, select its match to second image.
des2t = des2';                          % Precompute matrix transpose

% ADDED BY BRANDON SMITH
match = zeros(1,size(des1,1));
% END ADDED BY B.S.

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

   % Check if nearest neighbor has angle less than distRatio times 2nd.
   if (vals(1) < distRatio * vals(2))
      match(i) = indx(1);
   else
      match(i) = 0;
   end
end

% IMITTED BY BRANDON SMITH
% Create a new image showing the two images side by side.
%im3 = appendimages(im1,im2);
%
% Show a figure with lines joining the accepted matches.
%figure('Position', [100 100 size(im3,2) size(im3,1)]);
%colormap('gray');
%imagesc(im3);
%hold on;
%cols1 = size(im1,2);
% END OMITTED

% ADDED BY BRANDON SMITH
count = 1;
% END ADDED BY B.S.

for i = 1: size(des1,1)
  if (match(i) > 0)
      
    % ADDED BY BRANDON SMITH
    point_pairs(count,1) = loc1(i,2);
    point_pairs(count,2) = loc1(i,1);
    point_pairs(count,3) = loc2(match(i),2);
    point_pairs(count,4) = loc2(match(i),1);
    count = count + 1;
    % END ADDED BY B.S.
    
    % OMITTED BY BRANDON SMITH
%    line([loc1(i,2) loc2(match(i),2)+cols1], ...
%         [loc1(i,1) loc2(match(i),1)], 'Color', 'c');
    % END OMITTED BY B.S.
  end
end
%hold off;
num_matches = sum(match > 0);

fprintf('Found %d matches.\n', num_matches);

