function [intersection,status] = cv_int_line_plane(line,plane,dist_eps)
%
% cv_int_line_plane - determine intersection of line and plane
% Call: [intersection,status] = cv_int_line_plane(line,plane);
% On input:
%     line: 2 points defining line (matrix 2x3)
%     plane: 3 points defining plane (matrix 3x3)
%     dist_eps: minimum distance to count as intersection
% On output:
%     intersection: 2 points    if line in plane (matrix 2x3)
%                   inf         if no intersection
%                   1 point     if single point intersection
%     status:  -1    if line in plane
%               0    if no intersection
%               1    if single point intersection
% Author:
%     Tom Henderson
%     7 January 2000
% Custom Functions Used:
%     none
% Method:
%     find parameters for points in line and plane that minimize distance
% Testing:
% > cv_int_line_plane([0 1 0; 1 0 0], [0 0 0; 1 0 0; 1 1 0],0.01)
%
% ans =
%
%     0     1     0
%     1     0     0
%
% > cv_int_line_plane([0 0 1; 1 0 0], [0 0 0; 1 0 0; 1 1 0],0.01)
%
% ans =
%
%     1     0     0
%
% > cv_int_line_plane([0 0 1; 1 0 1], [0 0 0; 1 0 0; 1 1 0],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 33
%
% ans =
%
%   Inf
%
% > cv_int_line_plane([0 0 1; 1 0 0.9], [0 0 0; 1 0 0; 1 1 0],0.01)
%
% ans =
%
%   10.0000         0         0
%
% > cv_int_line_plane([0 0 1; 1 1 0.9], [0 0 0; 1 0 0; 1 1 0],0.01)
%
% ans =
%
%   10.0000   10.0000         0
%
if (cv_dist_pt_plane(line(1,:),plane)<dist_eps) ...   % line in plane
      & (cv_dist_pt_plane(line(2,:),plane)<dist_eps)
   intersection = line;
   status = -1;
   return
end

d = cv_dist_line_plane(line,plane);
if d<dist_eps   % line intersects plane in point
    p0 = plane(1,:);
	p1 = plane(2,:);
	p2 = plane(3,:);
	r0 = line(1,:);
	r1 = line(2,:);
	v = (p1-p0)';
	w = (p2-p0)';
	s = (r1-r0)';
	A = [v w -s];
	b = (p0-r0)';
	u = A\b;
   intersection = r0 - u(3)*s';
   status = 1;
else    %  no intersection
   intersection = inf;
   status = 0;
end
