function h = cvexShowMatches(I1, I2, matched_points1, matched_points2, ...
    title1, title2, coordinate, t1, t2)
%cvexShowMatches Display corresponding feature points.
%  cvexShowMatches(I1,I2) displays an overlay of a pair of images.
%
%  cvexShowMatches(I1,I2,MATCHED_POINTS1,MATCHED_POINTS2) or
%  cvexShowMatches(I1,I2,MATCHED_POINTS1,MATCHED_POINTS2,TITLE1,TTILE2)
%  displays an overlay of a pair of images with a color-coded plot of
%  the corresponding points connected by a line.
%
%    I1 and I2 are the images to display.
%
%    MATCHED_POINTS1 and MATCHED_POINTS2 are the coordinates of
%    corresponding points in I1 and I2, respectively. The format of
%    MATCHED_POINTS1 and MATCHED_POINTS2 is [r1,r2,...,rN1;
%    c1,c2,...,cN1], where r and c are the row and column coordinates,
%    respectively, and N is the number of points. The coordinate of the
%    top-left pixel is assumed to be [0; 0].
% 
%    TITLE1 and TITLE2 are the legend titles for I1 and I2, respectively.
%
%  cvexShowMatches(I1,I2,MATCHED_POINTS1,MATCHED_POINTS2,TITLE1,TTILE2,'RC')
%  is the same as cvexShowMatches(I1,I2,MATCHED_POINTS1,MATCHED_POINTS2,
%  TITLE1,TTILE2)
%
%  cvexShowMatches(I1,I2,MATCHED_POINTS1,MATCHED_POINTS2,TITLE1,TTILE2,'XY')
%  The format of MATCHED_POINTS1 and MATCHED_POINTS2 is [x1,x2,...,xN1;
%  y1,y2,...,yN1], where x and y are the column and row coordinates,
%  respectively, and N is the number of points. The coordinate of the
%  top-left pixel is assumed to be [1; 1].
%
%  cvexShowMatches(...,T1,T2) transforms the images and matched points using
%  projective transformations, T1 and T2, and displays them.
%
%   H = cvexShowMatches(...) returns the handle to the figure
%   object. 
%
%   Notes
%   -----
%   This is a helper function in support of examples and is subject to
%   change or removal in a future release.
%
% See also vision.CornerDetector, detectSURFFeatures, matchFeatures,
%     vision.GeometricTransformEstimator

% Copyright 2010 The MathWorks, Inc.
% $Revision: 1.1.6.5 $  $Date: 2011/05/09 03:43:35 $

assert (nargin == 2 || nargin == 4 || nargin == 6 || ...
  ((nargin == 7 || nargin == 9) && ...
  (strcmp(coordinate,'RC') || strcmp(coordinate,'XY'))), ...
  'vision:cvexShowMatches:invalidInputs',...
  ['You must use one of the following syntaxes: ', ...
  'CVEXSHOWMATCHES(I1,I2) ',...
  'CVEXSHOWMATCHES(I1,I2,MATCHED_POINTS1,MATCHED_POINTS2) ',...
  'CVEXSHOWMATCHES(I1,I2,MATCHED_POINTS1,MATCHED_POINTS2,TITLE1,TTILE2) ',...
  'CVEXSHOWMATCHES(I1,I2,MATCHED_POINTS1,MATCHED_POINTS2,TITLE1,TTILE2,''RC'') ',...
  'CVEXSHOWMATCHES(I1,I2,MATCHED_POINTS1,MATCHED_POINTS2,TITLE1,TTILE2,''XY'') ',...
  'CVEXSHOWMATCHES(I1,I2,MATCHED_POINTS1,MATCHED_POINTS2,TITLE1,TTILE2,''RC'',T1,T2) ',...
  'CVEXSHOWMATCHES(I1,I2,MATCHED_POINTS1,MATCHED_POINTS2,TITLE1,TTILE2,''XY'',T1,T2)']);

I1 = im2double(I1);
I2 = im2double(I2);

if nargin > 2
  isXY = cvstGetCoordsChoice('fcnShowMatches');
  if isa(matched_points1, 'SURFPoints')
      matched_points1 = flipud(matched_points1.Location'-1);
  elseif isXY
      matched_points1 = flipud(matched_points1' - 1);
  end

  if isa(matched_points2, 'SURFPoints')
      matched_points2 = flipud(matched_points2.Location'-1);
  elseif isXY
      matched_points2 = flipud(matched_points2' - 1);
  end
    
  matched_points1 = double(matched_points1);
  matched_points2 = double(matched_points2);
  
else
  matched_points1 = [];
  matched_points2 = [];
end

if nargin == 9
    htrans = vision.GeometricTransformer(...
        'TransformMatrixSource', 'Input port');
    
    if isXY
      matched_points1 = transformPoints( ...
        convertGeometricTransform(t1), matched_points1);
      matched_points2 = transformPoints( ...
        convertGeometricTransform(t2), matched_points2);
    else
      matched_points1 = transformPoints(t1, matched_points1);
      matched_points2 = transformPoints(t2, matched_points2);
    end
    I1 = step(htrans, I1, t1);
    I2 = step(htrans, I2, t2);
end

% Create a centered overlay image
imgOverlay = zeros(max([size(I1); size(I2)]));
I1_rowLocation = round((size(imgOverlay, 1) - size(I1, 1))/2);
I1_colLocation = round((size(imgOverlay, 2) - size(I1, 2))/2);
I2_rowLocation = round((size(imgOverlay, 1) - size(I2, 1))/2);
I2_colLocation = round((size(imgOverlay, 2) - size(I2, 2))/2);

[nrows, ncols] = size(I1);
imgOverlay(I1_rowLocation+(1:nrows), I1_colLocation+(1:ncols), 1) = I1;
[nrows, ncols] = size(I2);
imgOverlay(I2_rowLocation+(1:nrows), I2_colLocation+(1:ncols), 2) = I2;
imgOverlay(I2_rowLocation+(1:nrows), I2_colLocation+(1:ncols), 3) = I2;

% Show overlay image
h = figure;
imshow(imgOverlay, [])
axis on
hold on

% If the points are using XY coordinates, convert them to use RC
% coordinates
if nargin == 7 && strcmp(coordinate,'XY')
    matched_points1 = flipud(matched_points1) - 1;
    matched_points2 = flipud(matched_points2) - 1;
end

% Plot color-coded correspondence points connected by lines 
if (~isempty(matched_points1) && ~isempty(matched_points2))
    % Add one to coordinates to adjust zero-based point locations to
    % a one-based image plot
    plot(I1_colLocation + (matched_points1(2,:) + 1), ...
        I1_rowLocation + (matched_points1(1,:) + 1), 'ro');
    
    plot(I2_colLocation + (matched_points2(2,:) + 1), ...
        I2_rowLocation + (matched_points2(1,:) + 1), 'g+');
    
    plot([I1_colLocation + (matched_points1(2,:) + 1); ...
        I2_colLocation + (matched_points2(2,:) + 1)], ...
        [I1_rowLocation + (matched_points1(1,:) + 1); ...
        I2_rowLocation + (matched_points2(1,:) + 1)], 'y-');
    
    if nargin > 4
      legend(title1, title2)
    else
      % turn interpreter off so that characters are displayed
      % without an additional layer of interpretation by the
      % legend command
      legend({inputname(3), inputname(4)}, 'Interpreter', 'none');
    end
end

%==========================================================================
function t_points = transformPoints(t, points)
pts = bsxfun(@plus,  t(:,[1,2])*points, t(:,3));
t_points = bsxfun(@rdivide, pts([1,2],:), pts(3,:));

%==========================================================================
function T = convertGeometricTransform(t)
T = t;
T(1,1) = t(2,2);
% T(1,2) = t(1,2); no change
T(1,3) = t(3,2);
% T(2,1) = t(2,1); no change
T(2,2) = t(1,1);
T(2,3) = t(3,1);
T(3,1) = t(2,3);
T(3,2) = t(1,3);
% T(3,3) = t(3,3); no change
