intersectLineMesh3d

 Intersection points of a 3D line with a mesh.

   INTERS = intersectLineMesh3d(LINE, VERTICES, FACES)
   Compute the intersection points between a 3D line and a 3D mesh defined
   by vertices and faces. The mesh data is provided as a pair of arrays,
   with VERTICES being a NV-by-3 array of vertex coordinates, and FACES
   being a NF-by-3 or NF-by-4 array of face vertex indices.
   The LINE data correspond to a 1-by-6 row vector concatenating the line
   origin and direction. LINE can also be a NL-by-6 array representing a
   collection of lines with various origins and directions.

   INTERS = intersectLineMesh3d(LINE, MESH)
   Provides the mesh data as a struct with the fields 'vertices' and
   'faces'.

   [INTERS, POS, IFACES] = intersectLineMesh3d(...)
   Also returns the position of each intersection point on the input line,
   and the index of the intersected faces.
   If POS > 0, the point is also on the ray corresponding to the line. 
   
   [INTERS, POS, IFACES, ILINES] = intersectLineMesh3d(...)
   Also returns the index of the line each intersection point belong to.

   Example
     [V, F] = createCube;
     line = [.2 .3 .4 1 0 0];
     pts = intersectLineMesh3d(line, V, F)
     pts =
         1.0000    0.3000    0.4000
              0    0.3000    0.4000

   See also
   meshes3d, triangulateFaces, intersectLineTriangle3d

0001 function [points, pos, faceInds, lineInds] = intersectLineMesh3d(line, vertices, varargin)
0002 % Intersection points of a 3D line with a mesh.
0003 %
0004 %   INTERS = intersectLineMesh3d(LINE, VERTICES, FACES)
0005 %   Compute the intersection points between a 3D line and a 3D mesh defined
0006 %   by vertices and faces. The mesh data is provided as a pair of arrays,
0007 %   with VERTICES being a NV-by-3 array of vertex coordinates, and FACES
0008 %   being a NF-by-3 or NF-by-4 array of face vertex indices.
0009 %   The LINE data correspond to a 1-by-6 row vector concatenating the line
0010 %   origin and direction. LINE can also be a NL-by-6 array representing a
0011 %   collection of lines with various origins and directions.
0012 %
0013 %   INTERS = intersectLineMesh3d(LINE, MESH)
0014 %   Provides the mesh data as a struct with the fields 'vertices' and
0015 %   'faces'.
0016 %
0017 %   [INTERS, POS, IFACES] = intersectLineMesh3d(...)
0018 %   Also returns the position of each intersection point on the input line,
0019 %   and the index of the intersected faces.
0020 %   If POS > 0, the point is also on the ray corresponding to the line.
0021 %
0022 %   [INTERS, POS, IFACES, ILINES] = intersectLineMesh3d(...)
0023 %   Also returns the index of the line each intersection point belong to.
0024 %
0025 %   Example
0026 %     [V, F] = createCube;
0027 %     line = [.2 .3 .4 1 0 0];
0028 %     pts = intersectLineMesh3d(line, V, F)
0029 %     pts =
0030 %         1.0000    0.3000    0.4000
0031 %              0    0.3000    0.4000
0032 %
0033 %   See also
0034 %   meshes3d, triangulateFaces, intersectLineTriangle3d
0035 %
0036 
0037 % ------
0038 % Author: David Legland
0039 % e-mail: david.legland@inrae.fr
0040 % Created: 2011-12-20,    using Matlab 7.9.0.529 (R2009b)
0041 % Copyright 2011 INRA - Cepia Software Platform.
0042 
0043 % tolerance for detecting if a point is on line or within edge bounds
0044 tol = 1e-12;
0045 
0046 % parsing
0047 if ~isempty(varargin)
0048     if isscalar(varargin{1})
0049         tol = varargin{1};
0050         [vertices, faces] = parseMeshData(vertices);
0051     else
0052         faces = varargin{1};
0053         varargin(1) = [];
0054         if ~isempty(varargin)
0055             tol = varargin{1};
0056         end
0057     end
0058 else
0059     [vertices, faces] = parseMeshData(vertices);
0060 end
0061 
0062 % ensure the mesh has triangular faces
0063 tri2Face = [];
0064 if iscell(faces) || size(faces, 2) ~= 3
0065     [faces, tri2Face] = triangulateFaces(faces);
0066 end
0067 
0068 % find triangle edge vectors
0069 t0  = vertices(faces(:,1), :);
0070 u   = vertices(faces(:,2), :) - t0;
0071 v   = vertices(faces(:,3), :) - t0;
0072 
0073 % triangle normal
0074 n   = crossProduct3d(u, v);
0075 
0076 % direction vectors of lines and origins of lines
0077 dv = permute(line(:,4:6), [3 2 1]);
0078 d0 = permute(line(:,1:3), [3 2 1]);
0079 
0080 % vector between triangle origin and line origin
0081 w0 = d0 - t0;
0082 
0083 a = -sum(n .* w0, 2); % negative dot product
0084 b = sum(n .* dv, 2);  % dot product
0085 
0086 valid = abs(b) > tol & vectorNorm3d(n) > tol;
0087 
0088 % compute intersection point of line with supporting plane
0089 % If pos < 0: point before ray
0090 % IF pos > |dir|: point after edge
0091 pos = a ./ b;
0092 
0093 % coordinates of intersection point
0094 points = d0 + (pos .* dv);
0095 
0096 
0097 %% test if intersection point is inside triangle
0098 
0099 % normalize direction vectors of triangle edges
0100 uu  = dot(u, u, 2);
0101 uv  = dot(u, v, 2);
0102 vv  = dot(v, v, 2);
0103 
0104 % coordinates of vector v in triangle basis
0105 w   = points - t0;
0106 wu  = sum(w .* u, 2);
0107 wv  = sum(w .* v, 2);
0108 
0109 % normalization constant
0110 D = uv.^2 - uu .* vv;
0111 
0112 % test first coordinate
0113 s = (uv .* wv - vv .* wu) ./ D;
0114 ind1 = s < -tol | s > (1.0 + tol);
0115 
0116 % test second coordinate, and third triangle edge
0117 t = (uv .* wu - uu .* wv) ./ D;
0118 ind2 = t < -tol | (s + t) > (1.0 + tol);
0119 
0120 % keep only interesting points
0121 inds = ~ind1 & ~ind2 & valid;
0122 [faceInds, lineInds] = find(permute(inds, [1 3 2]));
0123 
0124 % Bit of an indexing trick to get points in appropriate order
0125 points = points(sub2ind(size(points), ...
0126     faceInds+[0 0 0], faceInds*0+(1:3), lineInds+[0 0 0]) );
0127 
0128 if nargout > 1
0129     pos = pos(inds);
0130     
0131     % convert to face indices of original mesh
0132     if ~isempty(tri2Face)
0133         faceInds = tri2Face(faceInds);
0134     end
0135 end