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 % E-mail: david.legland@inrae.fr
0038 % Created: 2005-02-17
0039 % Copyright 2005-2024 INRA - TPV URPOI - BIA IMASTE
0040 
0041 % extract tolerance if needed
0042 tol = 1e-14;
0043 if nargin > 2
0044     tol = varargin{1};
0045 end
0046 
0047 % unify sizes of data
0048 nLines  = size(line, 1);
0049 nPlanes = size(plane, 1);
0050 
0051 % N planes and M lines not allowed
0052 if nLines ~= nPlanes && min(nLines, nPlanes) > 1
0053     error('MatGeom:geom3d:intersectLinePlane', ...
0054         'Input must have same number of rows, or one must be 1');
0055 end
0056 
0057 % plane normal
0058 n = crossProduct3d(plane(:,4:6), plane(:,7:9));
0059 
0060 % difference between origins of plane and line
0061 dp = bsxfun(@minus, plane(:, 1:3), line(:, 1:3));
0062 
0063 % dot product of line direction with plane normal
0064 denom = sum(bsxfun(@times, n, line(:,4:6)), 2);
0065 
0066 % relative position of intersection point on line (can be inf in case of a
0067 % line parallel to the plane)
0068 t = sum(bsxfun(@times, n, dp),2) ./ denom;
0069 
0070 % compute coord of intersection point
0071 point = bsxfun(@plus, line(:,1:3),  bsxfun(@times, [t t t], line(:,4:6)));
0072 
0073 % set indices of line and plane which are parallel to NaN
0074 par = abs(denom) < tol;
0075 point(par,:) = NaN;