function [intersection,status] = cv_int_plane_plane(plane1,plane2,dist_eps)
%
% cv_int_plane_plane - determine intersection of 2 planes
% Call: [intersection,status] = cv_int_plane_plane(plane1,plane2);
% On input:
%     plane1: 3 points defining a plane (matrx 3x3)
%     plane2: 3 points defining a plane (matrix 3x3)
%     dist_eps: minimum distance to count as intersection
% On output:
%     intersection: plane    if same plane   (matrix 3x3)
%                   inf      if no intersection
%                   line     if line intersection  (matrix 2x3)
%     status:  -1    if plane intersection
%               0    if no intersection
%               1    if line intersection
% Author:
%     Tom Henderson
%     7 January 2000
% Custom Functions Used:
%     none
% Method:
%     find parameters for points on 2 planes that minimize distance
% Testing:
% > cv_int_plane_plane(p1,p1,0.01)
%
% ans =
%
%     0     0     0
%     1     0     0
%     1     1     0
%
% > cv_int_plane_plane(p1,p3,0.01)
% Warning: Matrix is singular to working precision.
% > In C:\MATLABR11\work\cv_dist_line_plane.m at line 45
%  In C:\MATLABR11\work\cv_int_line_plane.m at line 65
%  In C:\MATLABR11\work\cv_int_plane_plane.m at line 48
%
% ans =
%
%         0         0         0
%         0    1.0000         0
%
% > cv_int_plane_plane(p1,p2,0.01)
%
% ans =
%
%   Inf
%
if ((cv_dist_pt_plane(plane1(1,:),plane2)<dist_eps) ...  % same plane
      & (cv_dist_pt_plane(plane1(2,:),plane2)<dist_eps) ...
      & (cv_dist_pt_plane(plane1(3,:),plane2)<dist_eps))
   intersection = plane1;
   status = -1;
   return
end

pm = mean(plane1);   % get a 4th point and form all pairwise lines with points
l1 = plane1(1:2,:);  % from plane 1
l2 = plane1(2:3,:);
l3 = [plane1(1,:) ; plane1(2,:)];
l4 = [plane1(1,:) ; pm];
l5 = [plane1(2,:) ; pm];
l6 = [plane1(3,:) ; pm];
lines = [l1(1,:) l1(2,:); l2(1,:) l2(2,:); ...
      l3(1,:) l3(2,:); l4(1,:) l4(2,:); ...
      l5(1,:) l5(2,:); l6(1,:) l6(2,:)];
 
for i1 = 1:5     %  find 2 points that are not the same
   [p1,s1] = cv_int_line_plane([lines(i1,1:3);lines(i1,4:6)],plane2,dist_eps);
   for i2 = i1+1:6
      [p2,s2] = cv_int_line_plane([lines(i2,1:3);lines(i2,4:6)],plane2,dist_eps);
      if (s1>0 & s2>0)
         if (cv_dist_pt_pt(p1,p2)>dist_eps)
            status = 1;
            intersection = [p1 ; p2];   % and that gives the line of intersection
            return
         end
      end
   end
end

intersection = inf;    % no intersection
status = 0;
