intersectLineTriangle3d

INTERSECTLINETRIANGLE3D Intersection point of a 3D line and a 3D triangle.

   POINT = intersectLineTriangle3d(LINE, TRI)
   Compute coordinates of the intersection point between the line LINE and
   the triangle TRI.
   LINE is a 1-by-6 row vector given as: [X0 Y0 Z0 DX DY DZ]
   TRI is given either as a row vector [X1 Y1 Z1 X2 Y2 Z2 X3 Y3 Z3], or as
   a 3-by-3 array, each row containing coordinates of one of the triangle
   vertices.
   The result is a 1-by-3 array containing coordinates of the intesection
   point, or [NaN NaN NaN] if the line and the triangle do not intersect.

   [POINT POS] = intersectLineTriangle3d(LINE, TRI)
   Also returns the position of the intersection point on the line, or NaN
   if the line and the supporting plane of the triangle are parallel.

   [POINT POS ISINSIDE] = intersectLineTriangle3d(LINE, TRI)
   Also returns a boolean value, set to true if the line and the triangle
   intersect each other. Can be used for testing validity of result.

   Example
     line = [1 1 0 0 0 1];
     tri = [0 0 5;5 0 0;0 5 0];
     intersectLineTriangle3d(line, tri)
     ans = 
         1   1   3

   See also
   points3d, lines3d, polygons3d, intersectRayPolygon3d,
   distancePointTriangle3d

   References
   Algorithm adapted from SoftSurfer Ray/Segment-Triangle intersection
   http://softsurfer.com/Archive/algorithm_0105/algorithm_0105.htm

0001 function [point, pos, isInside] = intersectLineTriangle3d(line, triangle, varargin)
0002 %INTERSECTLINETRIANGLE3D Intersection point of a 3D line and a 3D triangle.
0003 %
0004 %   POINT = intersectLineTriangle3d(LINE, TRI)
0005 %   Compute coordinates of the intersection point between the line LINE and
0006 %   the triangle TRI.
0007 %   LINE is a 1-by-6 row vector given as: [X0 Y0 Z0 DX DY DZ]
0008 %   TRI is given either as a row vector [X1 Y1 Z1 X2 Y2 Z2 X3 Y3 Z3], or as
0009 %   a 3-by-3 array, each row containing coordinates of one of the triangle
0010 %   vertices.
0011 %   The result is a 1-by-3 array containing coordinates of the intesection
0012 %   point, or [NaN NaN NaN] if the line and the triangle do not intersect.
0013 %
0014 %   [POINT POS] = intersectLineTriangle3d(LINE, TRI)
0015 %   Also returns the position of the intersection point on the line, or NaN
0016 %   if the line and the supporting plane of the triangle are parallel.
0017 %
0018 %   [POINT POS ISINSIDE] = intersectLineTriangle3d(LINE, TRI)
0019 %   Also returns a boolean value, set to true if the line and the triangle
0020 %   intersect each other. Can be used for testing validity of result.
0021 %
0022 %   Example
0023 %     line = [1 1 0 0 0 1];
0024 %     tri = [0 0 5;5 0 0;0 5 0];
0025 %     intersectLineTriangle3d(line, tri)
0026 %     ans =
0027 %         1   1   3
0028 %
0029 %   See also
0030 %   points3d, lines3d, polygons3d, intersectRayPolygon3d,
0031 %   distancePointTriangle3d
0032 %
0033 %   References
0034 %   Algorithm adapted from SoftSurfer Ray/Segment-Triangle intersection
0035 %   http://softsurfer.com/Archive/algorithm_0105/algorithm_0105.htm
0036 %
0037 
0038 % ------
0039 % Author: David Legland
0040 % E-mail: david.legland@inrae.fr
0041 % Created: 2011-04-08, using Matlab 7.9.0.529 (R2009b)
0042 % Copyright 2011-2024 INRA - Cepia Software Platform
0043 
0044 %% Default values
0045 
0046 % default return value
0047 point = [NaN NaN NaN];
0048 pos = NaN;
0049 isInside = false;
0050 
0051 tol = 1e-13;
0052 if ~isempty(varargin)
0053     tol = varargin{1};
0054 end
0055 
0056 
0057 %% Process inputs
0058 
0059 % triangle edge vectors
0060 if size(triangle, 2) > 3
0061     % triangle is given as a 1-by-9 row vector
0062     t0  = triangle(1:3);
0063     u   = triangle(4:6) - t0;
0064     v   = triangle(7:9) - t0;
0065 else
0066     % triangle is given as a 3-by-3 array
0067     t0  = triangle(1, 1:3);
0068     u   = triangle(2, 1:3) - t0;
0069     v   = triangle(3, 1:3) - t0;
0070 end
0071 
0072 
0073 %% Compute intersection
0074 
0075 % triangle normal
0076 n   = cross(u, v);
0077 
0078 % test for degenerate case of flat triangle
0079 if vectorNorm3d(n) < tol
0080     return;
0081 end
0082 
0083 % line direction vector
0084 dir = line(4:6);
0085 
0086 % vector between triangle origin and line origin
0087 w0 = line(1:3) - t0;
0088 
0089 % compute projection of each vector on the plane normal
0090 a = -dot(n, w0);
0091 b = dot(n, dir);
0092 
0093 % test case of line parallel to the triangle
0094 if abs(b) < tol
0095     return;    
0096 end
0097 
0098 % compute intersection point of line with supporting plane
0099 % If r < 0: point before ray
0100 % If r > 1: point after edge
0101 pos = a / b;
0102 
0103 % coordinates of intersection point
0104 point = line(1:3) + pos * dir;
0105 
0106 
0107 %% test if intersection point is inside triangle
0108 
0109 % normalize direction vectors of triangle edges
0110 uu  = dot(u, u);
0111 uv  = dot(u, v);
0112 vv  = dot(v, v);
0113 
0114 % coordinates of vector v in triangle basis
0115 w   = point - t0;
0116 wu  = dot(w, u);
0117 wv  = dot(w, v);
0118 
0119 % normalization constant
0120 D = uv^2 - uu * vv;
0121 
0122 % test first coordinate
0123 s = (uv * wv - vv * wu) / D;
0124 if s < 0.0 || s > 1.0
0125     point = [NaN NaN NaN];
0126     pos = NaN;
0127     return;
0128 end
0129 
0130 % test second coordinate, and third triangle edge
0131 t = (uv * wu - uu * wv) / D;
0132 if t < 0.0 || (s + t) > 1.0
0133     point = [NaN NaN NaN];
0134     pos = NaN;
0135     return;
0136 end
0137 
0138 % set the validity flag
0139 isInside = true;
0140