function [index_pairs, match_metric] = matchFeatures(varargin)
%matchFeatures Find matching features
%   INDEX_PAIRS = matchFeatures(FEATURES1, FEATURES2) returns a P-by-2
%   matrix, INDEX_PAIRS, containing indices to the features most likely to
%   correspond between the two input feature matrices. The function takes
%   two inputs, FEATURES1, an M1-by-N matrix, and FEATURES2, an M2-by-N
%   matrix.
%
%   [INDEX_PAIRS, MATCH_METRIC] = matchFeatures(FEATURES1, FEATURES2, ...)
%   also returns the metric values that correspond to the associated
%   features indexed by INDEX_PAIRS in a P-by-1 matrix MATCH_METRIC.
%
%   [INDEX_PAIRS, MATCH_METRIC] =
%   matchFeatures(...,PARAM1,VAL1,PARAM2,VAL2,...) Additional control over
%   the matching algorithm can be accomplished by specifying parameters and
%   corresponding values.
% 
%   Parameters include:
%   -------------------       
%   'Metric'           A string used to specify the feature matching
%                      metric. Possible values are:
%                        'SAD'         : Sum of absolute differences
%                        'SSD'         : Sum of squared differences 
%                        'normxcorr'   : Normalized cross-correlation
%
%                      Default: 'SSD'
%
%   'MatchThreshold'   A scalar T, 0 < T <= 100, specifying a threshold
%                      for selecting the strongest matches. Matches having
%                      metrics more than T percent from a perfect match 
%                      are rejected. Increase T to return more matches.
% 
%                      Default: 1.0
% 
%   Notes
%   ----- 
%   The range of values of MATCH_METRIC varies as a function of the feature
%   matching metric being used. Prior to computation of SAD and SSD
%   metrics, the feature vectors are normalized to unit vectors. The table
%   below summarizes the metric ranges and perfect match values:
%
%      Metric      Range                            Perfect Match Value
%      ----------  -------------------------------  ------------------- 
%      SAD         [0, 2*sqrt(size(FEATURES1, 1))]          0
%      SSD         [0, 4]                                   0
%      normxcorr   [-1, 1]                                  1
%
%   Class Support
%   -------------
%   FEATURES1 and FEATURES2 can be logical, int8, uint8, int16, uint16,
%   int32, uint32, single, or double, and they must be real and nonsparse.
%
%   The output class of INDEX_PAIRS is uint32. MATCH_METRIC is double when
%   FEATURES1 and FEATURES2 are double. Otherwise, it is single.
%
%   Example 1
%   ---------
%   % Find corresponding interest points between a pair of images using
%   % local neighborhoods.
% 
%   I1 = rgb2gray(imread('viprectification_deskLeft.png'));
%   I2 = rgb2gray(imread('viprectification_deskRight.png'));
% 
%   % Find corners  
%   hcornerdet = vision.CornerDetector;
%   points1 = step(hcornerdet, I1);
%   points2 = step(hcornerdet, I2);
% 
%   % Extract neighborhood features
%   [features1, valid_points1] = extractFeatures(I1, points1);
%   [features2, valid_points2] = extractFeatures(I2, points2);
% 
%   % Match features
%   index_pairs = matchFeatures(features1, features2);
% 
%   % Retrieve locations of corresponding points for each image
%   matched_points1 = valid_points1(index_pairs(:, 1), :);
%   matched_points2 = valid_points2(index_pairs(:, 2), :);
% 
%   % Visualize corresponding points
%   cvexShowMatches(I1, I2, matched_points1, matched_points2);
%
%   Example 2
%   ---------
%   % Use SURF features to find corresponding points between two images
%   % rotated and scaled with respect to each other
%   I1 = imread('cameraman.tif');
%   I2 = imresize(imrotate(I1,-20), 1.2);
% 
%   points1 = detectSURFFeatures(I1);
%   points2 = detectSURFFeatures(I2);
% 
%   [f1, vpts1] = extractFeatures(I1, points1);
%   [f2, vpts2] = extractFeatures(I2, points2);
%         
%   index_pairs = matchFeatures(f1, f2) ;
%   matched_pts1 = vpts1(index_pairs(:, 1));
%   matched_pts2 = vpts2(index_pairs(:, 2));
%   
%   % Note that there are still several outliers present in the data, but
%   % otherwise you can clearly see the effects of rotation and scaling on 
%   % the display of matched features.
%   cvexShowMatches(I1,I2,matched_pts1,matched_pts2); % show matching points
%
% See also vision.CornerDetector, detectSURFFeatures, extractFeatures,
%          estimateFundamentalMatrix, vision.GeometricTransformEstimator.

%   We use the normalized cross-correlation defined in this paper:
%   
%      "Fast Normalized Cross-Correlation", by J. P. Lewis, Industrial
%      Light & Magic.
%      http://www.idiom.com/~zilla/Papers/nvisionInterface/nip.html
%
%  Copyright 2010 The MathWorks, Inc.

isXY = cvstGetCoordsChoice('fcnMatchFeatures');

% Parse and check inputs
[features1, features2, metric, match_thresh] = parseInputs(isXY, varargin{:});

% Match features
if isXY
    features1 = features1';
    features2 = features2';
end

[index_pairs, match_metric] = cvalgMatchFeatures(features1, features2, ...
    metric, match_thresh);
  
if isXY
    index_pairs = index_pairs';
    match_metric = match_metric';
end

%==========================================================================
% Parse and check inputs
%==========================================================================
function [features1, features2, metric, match_thresh] =...
    parseInputs(isXY, varargin)

assert(nargin ~= 4, 'vision:matchFeatures:NotEnoughArgs', ...
    'Not enough input arguments.');

% Value of string enums. The first string is the default.
metricOptions = {'SSD', 'SAD', 'normxcorr'};

% Setup parser
parser = inputParser;
parser.CaseSensitive = true;
parser.addRequired('features1', @checkFeatures);
parser.addRequired('features2', @checkFeatures);
parser.addParamValue('Metric', metricOptions{1}, ...
  @(x)(checkMetric(metricOptions, x)));
parser.addParamValue('MatchThreshold', 1.0, @checkMatchThreshold);

% Parse input
parser.parse(varargin{:});

% Assign outputs
features1 = parser.Results.features1;
features2 = parser.Results.features2;
metric = lower(parser.Results.Metric);
match_thresh = parser.Results.MatchThreshold;

assert(isequal(class(features1), class(features2)), ...
    'vision:matchFeatures:featuresNotSameClass', ...
    'Expected FEATURES1 and FEATURES2 to be the same class.');
  
if isXY
    assert(size(features1, 2) == size(features2, 2), ...
        'vision:matchFeatures:featuresNotSameColumns', ...
        ['Expected FEATURES1 and FEATURES2 to have the same number ', ...
        'of columns.']);
else
    assert(size(features1, 1) == size(features2, 1), ...
        'vision:matchFeatures:featuresNotSameRows', ...
        ['Expected FEATURES1 and FEATURES2 to have the same number ', ...
        'of rows.']);
end

%==========================================================================
function tf = checkFeatures(features)
validateattributes(features, {'logical', 'int8', 'uint8', 'int16', ...
    'uint16', 'int32', 'uint32', 'single', 'double'}, {'2d', ...
    'nonsparse', 'real'}, mfilename, 'FEATURES');
tf = true;

%========================================================================== 
function tf = checkMetric(list, value)
message = 'Expected ''Metric'' to match one of these strings: ';
numStr = length(list);
for idx = 1:(numStr-1)
    message = [message sprintf('%s, ', list{idx})]; %#ok<AGROW>
end
message = [message list{end}];

assert(ischar(value) && any(strcmpi(value, list)), ...
    'vision:matchFeatures:invalidMetric', [message '.']);
tf = true;

%==========================================================================
function tf = checkMatchThreshold(threshold)
validateattributes(threshold, {'numeric'}, {'nonempty', 'nonnan', ...
    'finite', 'nonsparse', 'real', 'positive', 'scalar', '<=', 100}, ...
    mfilename);
tf = true;
