intersectLinePlane

INTERSECTLINEPLANE Intersection point between a 3D line and a plane.

   PT = intersectLinePlane(LINE, PLANE)
   Returns the intersection point of the given line and the given plane.
   LINE:  [x0 y0 z0 dx dy dz]
   PLANE: [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
   PT:    [xi yi zi]
   If LINE and PLANE are parallel, return [NaN NaN NaN].
   If LINE (or PLANE) is a matrix with 6 (or 9) columns and N rows, result
   is an array of points with N rows and 3 columns.
   
   PT = intersectLinePlane(LINE, PLANE, TOL)
   Specifies the tolerance factor to test if a line is parallel to a
   plane. Default is 1e-14.

   Example
     % define horizontal plane through origin
     plane = [0 0 0   1 0 0   0 1 0];
     % intersection with a vertical line
     line = [2 3 4  0 0 1];
     intersectLinePlane(line, plane)
     ans = 
        2   3   0
     % intersection with a line "parallel" to plane
     line = [2 3 4  1 2 0];
     intersectLinePlane(line, plane)
     ans = 
       NaN  NaN  NaN

   See also:
   lines3d, planes3d, points3d, clipLine3d

0001 function point = intersectLinePlane(line, plane, varargin)
0002 %INTERSECTLINEPLANE Intersection point between a 3D line and a plane.
0003 %
0004 %   PT = intersectLinePlane(LINE, PLANE)
0005 %   Returns the intersection point of the given line and the given plane.
0006 %   LINE:  [x0 y0 z0 dx dy dz]
0007 %   PLANE: [x0 y0 z0 dx1 dy1 dz1 dx2 dy2 dz2]
0008 %   PT:    [xi yi zi]
0009 %   If LINE and PLANE are parallel, return [NaN NaN NaN].
0010 %   If LINE (or PLANE) is a matrix with 6 (or 9) columns and N rows, result
0011 %   is an array of points with N rows and 3 columns.
0012 %
0013 %   PT = intersectLinePlane(LINE, PLANE, TOL)
0014 %   Specifies the tolerance factor to test if a line is parallel to a
0015 %   plane. Default is 1e-14.
0016 %
0017 %   Example
0018 %     % define horizontal plane through origin
0019 %     plane = [0 0 0   1 0 0   0 1 0];
0020 %     % intersection with a vertical line
0021 %     line = [2 3 4  0 0 1];
0022 %     intersectLinePlane(line, plane)
0023 %     ans =
0024 %        2   3   0
0025 %     % intersection with a line "parallel" to plane
0026 %     line = [2 3 4  1 2 0];
0027 %     intersectLinePlane(line, plane)
0028 %     ans =
0029 %       NaN  NaN  NaN
0030 %
0031 %   See also:
0032 %   lines3d, planes3d, points3d, clipLine3d
0033 %
0034 
0035 %   ---------
0036 %   author : David Legland
0037 %   INRA - TPV URPOI - BIA IMASTE
0038 %   created the 17/02/2005.
0039 %
0040 
0041 %   HISTORY
0042 %   24/11/2005 add support for multiple input
0043 %   23/06/2006 correction from Songbai Ji allowing different number of
0044 %       lines or plane if other input has one row
0045 %   14/12/2006 correction for parallel lines and plane normals
0046 %   05/01/2007 fixup for parallel lines and plane normals
0047 %   24/04/2007 rename as 'intersectLinePlane'
0048 %   11/19/2010 Added bsxfun functionality for improved speed (Sven Holcombe)
0049 %   01/02/2011 code cleanup, add option for tolerance, update doc
0050 
0051 
0052 % extract tolerance if needed
0053 tol = 1e-14;
0054 if nargin > 2
0055     tol = varargin{1};
0056 end
0057 
0058 % unify sizes of data
0059 nLines  = size(line, 1);
0060 nPlanes = size(plane, 1);
0061 
0062 % N planes and M lines not allowed
0063 if nLines ~= nPlanes && min(nLines, nPlanes) > 1
0064     error('MatGeom:geom3d:intersectLinePlane', ...
0065         'Input must have same number of rows, or one must be 1');
0066 end
0067 
0068 % plane normal
0069 n = crossProduct3d(plane(:,4:6), plane(:,7:9));
0070 
0071 % difference between origins of plane and line
0072 dp = bsxfun(@minus, plane(:, 1:3), line(:, 1:3));
0073 
0074 % dot product of line direction with plane normal
0075 denom = sum(bsxfun(@times, n, line(:,4:6)), 2);
0076 
0077 % relative position of intersection point on line (can be inf in case of a
0078 % line parallel to the plane)
0079 t = sum(bsxfun(@times, n, dp),2) ./ denom;
0080 
0081 % compute coord of intersection point
0082 point = bsxfun(@plus, line(:,1:3),  bsxfun(@times, [t t t], line(:,4:6)));
0083 
0084 % set indices of line and plane which are parallel to NaN
0085 par = abs(denom) < tol;
0086 point(par,:) = NaN;