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@inra.fr
0041 % Created: 2011-04-08,    using Matlab 7.9.0.529 (R2009b)
0042 % Copyright 2011 INRA - Cepia Software Platform.
0043 
0044 
0045 %% Default values
0046 
0047 % default return value
0048 point = [NaN NaN NaN];
0049 pos = NaN;
0050 isInside = false;
0051 
0052 tol = 1e-13;
0053 if ~isempty(varargin)
0054     tol = varargin{1};
0055 end
0056 
0057 
0058 %% Process inputs
0059 
0060 % triangle edge vectors
0061 if size(triangle, 2) > 3
0062     % triangle is given as a 1-by-9 row vector
0063     t0  = triangle(1:3);
0064     u   = triangle(4:6) - t0;
0065     v   = triangle(7:9) - t0;
0066 else
0067     % triangle is given as a 3-by-3 array
0068     t0  = triangle(1, 1:3);
0069     u   = triangle(2, 1:3) - t0;
0070     v   = triangle(3, 1:3) - t0;
0071 end
0072 
0073 
0074 %% Compute intersection
0075 
0076 % triangle normal
0077 n   = cross(u, v);
0078 
0079 % test for degenerate case of flat triangle
0080 if vectorNorm3d(n) < tol
0081     return;
0082 end
0083 
0084 % line direction vector
0085 dir = line(4:6);
0086 
0087 % vector between triangle origin and line origin
0088 w0 = line(1:3) - t0;
0089 
0090 % compute projection of each vector on the plane normal
0091 a = -dot(n, w0);
0092 b = dot(n, dir);
0093 
0094 % test case of line parallel to the triangle
0095 if abs(b) < tol
0096     return;    
0097 end
0098 
0099 % compute intersection point of line with supporting plane
0100 % If r < 0: point before ray
0101 % If r > 1: point after edge
0102 pos = a / b;
0103 
0104 % coordinates of intersection point
0105 point = line(1:3) + pos * dir;
0106 
0107 
0108 %% test if intersection point is inside triangle
0109 
0110 % normalize direction vectors of triangle edges
0111 uu  = dot(u, u);
0112 uv  = dot(u, v);
0113 vv  = dot(v, v);
0114 
0115 % coordinates of vector v in triangle basis
0116 w   = point - t0;
0117 wu  = dot(w, u);
0118 wv  = dot(w, v);
0119 
0120 % normalization constant
0121 D = uv^2 - uu * vv;
0122 
0123 % test first coordinate
0124 s = (uv * wv - vv * wu) / D;
0125 if s < 0.0 || s > 1.0
0126     point = [NaN NaN NaN];
0127     pos = NaN;
0128     return;
0129 end
0130 
0131 % test second coordinate, and third triangle edge
0132 t = (uv * wu - uu * wv) / D;
0133 if t < 0.0 || (s + t) > 1.0
0134     point = [NaN NaN NaN];
0135     pos = NaN;
0136     return;
0137 end
0138 
0139 % set the validity flag
0140 isInside = true;
0141