function h = cvexShowStereoImages(title1,title2,I1,I2,pts1,pts2,f,t1,t2)
%cvexShowStereoImages Display stereo images, corresponding points, and epipolar lines.
%  cvexShowStereoImages(TITLE1,TITLE2,I1,I2,PTS1,PTS2,F) displays stereo
%  images, I1 and I2, corresponding points, PTS1 and PTS2, and epipolar
%  lines computed using the fundamental matrix F.
%
%  cvexShowStereoImages(...,T1,T2) firstly transforms the images using
%  projective transformations, T1 and T2, and the displays the images,
%  points, and lines.
%
%    PTS1 and PTS2 are the coordinates of corresponding points in I1 and
%    I2, respectively. Each column of PTS1 or PTS2 has the format of
%    [x, y] which determines the location of a point. Both PTS1 and
%    PTS2 are M-by-2 matrices, where M is the number of points.
% 
%    F is a 3-by-3 fundamental matrix. If P1, a point in I1, corresponds
%    to P2, a point in I2, then [P2,1] * F * [P1,1]' = 0.
%
%    TITLE1 and TITLE2 are the titles for I1 and I2, respectively.
% 
%  H = cvexShowStereoImages(...) 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 estimateUncalibratedRectification

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

% Check the number of input arguments.
if nargin ~= 7 && nargin ~= 9
  error('vision:cvexShowStereoImages:invalidInputs',...
    ['You must use one of the following syntaxes: ',...
     'cvexShowStereoImages(TITLE1,TITLE2,I1,I2,PTS1,PTS2,F) or ',...
     'cvexShowStereoImages(TITLE1,TITLE2,I1,I2,PTS1,PTS2,F,T1,T2).']);
end

I1 = im2single(I1);
I2 = im2single(I2);
pts1 = single(pts1);
pts2 = single(pts2);
epiLines1 = epipolarLine(f', pts2);
epiLines2 = epipolarLine(f, pts1);

% Transform the image, corresponding points, and epipolar lines.
if nargin > 7
  if isEpipoleInImage(f, size(I1)) || isEpipoleInImage(f', size(I2))
    error('vision:cvexShowStereoImages:epipoleInsideImage', ...
      ['The epipole is inside the image. You cannot apply the '...
       'transformation on the entire image.']);
  end
  
  [I1, pts1, epiLines1] = transformImagePoints(I1, pts1, epiLines1, f, t1);
  [I2, pts2, epiLines2] = transformImagePoints(I2, pts2, epiLines2, f, t2);
end
  
epiLines1 = lineToBorderPoints(epiLines1, size(I1));
epiLines2 = lineToBorderPoints(epiLines2, size(I2));

% Display stereo images, corresponding points, and epipolar lines. 
h = cvexShowImagePair(I1, I2, title1, title2, 'MultipleColors', ...
  pts1, pts2, epiLines1,epiLines2);

%==========================================================================
% Transform the image, corresponding points, and epipolar lines.
%==========================================================================
function [I, pts, epiLines] = transformImagePoints(I, pts, epiLines, f, t)

if isEpipoleInImage(f, size(I))
  error('vision:cvexShowStereoImages:epipoleInsideImage', ...
    ['The epipole is inside the image. You cannot apply the '...
     'transformation on the entire image.']);
end

htrans = vision.GeometricTransformer('TransformMatrixSource', 'Input port');
t = single(t);
I = step(htrans, I, t);

isXY = cvstGetCoordsChoice('fcnCvexShowImagePair');
if isXY
  t = convertGeometricTransform(t);
  pts = flipud(pts');
end

ptsr = bsxfun(@plus,  t(:,[1,2])*pts, t(:,3));
pts = bsxfun(@rdivide, ptsr([1,2],:), ptsr(3,:));

if isXY
  pts = fliplr(pts');
  epiLines = epiLines(:, [2 1 3])';
  epiLines = inv(t)' * epiLines;
  epiLines = epiLines([2 1 3], :)';
else
  epiLines = inv(t)' * epiLines;
end

%==========================================================================
function T = convertGeometricTransform(t)
T = t;
T(1,1) = t(2,2);
%T(1,2) = t(1,2);
T(1,3) = t(3,2);
%T(2,1) = t(2,1);
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);
