Description of distancePointMesh

0001 function [dist, proj] = distancePointMesh(point, vertices, faces, varargin)
0002 %DISTANCEPOINTMESH  Shortest distance between a (3D) point and a triangle mesh.
0003 %
0004 %   DIST = distancePointMesh(POINT, VERTICES, FACES)
0005 %   Returns the shortest distance between the query point POINT and the
0006 %   triangular mesh defined by the set of vertex coordinates VERTICES and
0007 %   the set of faces FACES. VERTICES is a NV-by-3 array, and FACES is a
0008 %   NF-by-3 array of vertex indices.
0009 %   If FACES is NF-by-4 array, it is converted to a (NF*2)-by-3 array.
0010 %
0011 %   [DIST, PROJ] = distancePointMesh(...)
0012 %   Also returns the projection of the query point on the triangular mesh.
0013 %
0014 %   ... = distancePointMesh(..., 'algorithm', ALGO)
0015 %   Allows to choose the type of algorithm. Options are:
0016 %   * sequential:   process each face sequentially, using the function
0017 %               distancePointTriangle3d
0018 %   * vectorized:   vectorized algorithm, usually faster for large number
0019 %               of faces
0020 %   * auto:         (default) automatically choose the most appropriate
0021 %               between sequential and vectorized.
0022 %
0023 %   Example
0024 %     [V, F] = torusMesh();
0025 %     F2 = triangulateFaces(F);
0026 %     P = [10 20 30];
0027 %     [D, PROJ] = distancePointMesh(P, V, F2);
0028 %     figure; drawMesh(V, F)
0029 %     view(3); axis equal; lighting gouraud; light;
0030 %     drawPoint3d(P);
0031 %     drawPoint3d(PROJ, 'm*');
0032 %     drawEdge3d([P PROJ], 'linewidth', 2, 'color', 'b');
0033 %
0034 %   See also
0035 %     distancePointTriangle3d
0036 %
0037 %   References
0038 %   * "Distance Between Point and Triangle in 3D", David Eberly (1999)
0039 %   https://www.geometrictools.com/Documentation/DistancePoint3Triangle3.pdf
0040 %   * <a href="matlab:
0041 %     web('https://fr.mathworks.com/matlabcentral/fileexchange/22857-distance-between-a-point-and-a-triangle-in-3d')
0042 %   ">Distance between a point and a triangle in 3d</a>, by Gwendolyn Fischer.
0043 %   * <a href="matlab:
0044 %     web('https://fr.mathworks.com/matlabcentral/fileexchange/52882-point2trimesh------distance%C2%A0between-point-and-triangulated-surface')
0045 %   ">Distance Between Point and Triangulated Surface</a>, by Daniel Frisch.
0046  
0047 % ------
0048 % Author: David Legland
0049 % e-mail: david.legland@inra.fr
0050 % Created: 2018-03-08,    using Matlab 9.3.0.713579 (R2017b)
0051 % Copyright 2018 INRA - Cepia Software Platform.
0052 
0053 
0054 %% Parses input arguments
0055 
0056 % check the case of mesh given as structure
0057 if isstruct(vertices)
0058     faces = vertices.faces;
0059     vertices = vertices.vertices;
0060 end
0061 
0062 % default option
0063 algo = 'auto';
0064 
0065 % check optional arguments
0066 while length(varargin) > 1
0067     varName = varargin{1};
0068     if ~ischar(varName)
0069         error('Require options as parameter name-value pairs');
0070     end
0071     
0072     if strcmpi(varName, 'algorithm')
0073         algo = varargin{2};
0074     else
0075         error(['Unknown option name: ' varName]);
0076     end
0077     
0078     varargin(1:2) = [];
0079 end
0080 
0081 % number of faces
0082 nFaces = size(faces, 1);
0083 
0084 if size(faces, 2) > 3 || iscell(faces)
0085     faces = triangulateFaces(faces);
0086 end
0087 
0088 % If algorithm is chosen automatically, choose depending on face number
0089 if strcmpi(algo, 'auto')
0090     if size(faces, 1) > 30
0091         algo = 'vectorized';
0092     else
0093         algo = 'sequential';
0094     end
0095 end
0096 
0097 % switch to vectorized algorithm if necessary
0098 if strcmpi(algo, 'vectorized')
0099     if nargout > 1
0100         [dist, proj] = distancePointTrimesh_vectorized(point, vertices, faces);
0101     else
0102         dist = distancePointTrimesh_vectorized(point, vertices, faces);
0103     end
0104     return;
0105 end
0106 
0107 
0108 %% Sequential algorithm
0109 % For each point, iterates over the triangular faces
0110 
0111 % allocate memory for result
0112 nPoints = size(point, 1);
0113 dist = zeros(nPoints, 1);
0114 
0115 if nargout > 1
0116     projPoints = zeros(nPoints, 3);
0117 end
0118 
0119 % iterate over points
0120 for i = 1:nPoints
0121     % % min distance and projection for current point
0122     minDist = inf;
0123     proj = [0 0 0];
0124     
0125     % iterate over faces
0126     for iFace = 1:nFaces
0127         % create triange for current face
0128         face = faces(iFace, :);
0129         triangle = vertices(face, :);
0130         
0131         [distf, projf] = distancePointTriangle3d(point(i,:), triangle);
0132         
0133         if distf < minDist
0134             minDist = distf;
0135             proj = projf;
0136         end
0137     end
0138     
0139     dist(i) = minDist;
0140     if nargout > 1
0141         projPoints(i,:) = proj;
0142     end
0143 end
0144 end
0145 
0146 function [dist, proj] = distancePointTrimesh_vectorized(point, vertices, faces)
0147 %DISTANCEPOINTTRIMESH Vectorized version of the distancePointTrimesh function
0148 %
0149 %   output = distancePointTrimesh_vectorized(input)
0150 %
0151 %   This version is  vectorized over faces: for each query point, the
0152 %   minimum distance to each triangular face is computed in parallel.
0153 %   Then the minimum distance over faces is kept.
0154 %
0155 %   Example
0156 %   distancePointTrimesh
0157 %
0158 %   See also
0159 %
0160 
0161 % ------
0162 % Author: David Legland
0163 % e-mail: david.legland@inra.fr
0164 % Created: 2018-03-08,    using Matlab 9.3.0.713579 (R2017b)
0165 % Copyright 2018 INRA - Cepia Software Platform.
0166 
0167 % Regions are not numbered as in the original paper of D. Eberly to allow
0168 % automated computation of regions from the 3 conditions on lines.
0169 % Region indices are computed as follow:
0170 %   IND = b2 * 2^2 + b1 * 2 + b0
0171 % with:
0172 %   b0 = 1 if s < 0, 0 otherwise
0173 %   b1 = 1 if t < 0, 0 otherwise
0174 %   b2 = 1 if s+t > 1, 0 otherwise
0175 % resulting ion the following region indices:
0176 %        /\ t
0177 %        |
0178 %   \ R5 |
0179 %    \   |
0180 %     \  |
0181 %      \ |
0182 %       \| P3
0183 %        *
0184 %        |\
0185 %        | \
0186 %   R1   |  \   R4
0187 %        |   \
0188 %        | R0 \
0189 %        |     \
0190 %        | P1   \ P2
0191 %  ------*-------*------> s
0192 %        |        \
0193 %   R3   |   R2    \   R6
0194 
0195 % allocate memory for result
0196 nPoints = size(point, 1);
0197 dist = zeros(nPoints, 1);
0198 proj = zeros(nPoints, 3);
0199 
0200 % triangle origins and direction vectors
0201 p1  = vertices(faces(:,1),:);
0202 v12 = vertices(faces(:,2),:) - p1;
0203 v13 = vertices(faces(:,3),:) - p1;
0204 
0205 % identify coefficients of second order equation that do not depend on
0206 % query point
0207 a = dot(v12, v12, 2);
0208 b = dot(v12, v13, 2);
0209 c = dot(v13, v13, 2);
0210 
0211 % iterate on query points
0212 for i = 1:nPoints
0213     % coefficients of second order equation that depend on query point
0214     diffP = bsxfun(@minus, p1, point(i, :));
0215     d = dot(v12, diffP, 2);
0216     e = dot(v13, diffP, 2);
0217     
0218     % compute position of projected point in the plane of the triangle
0219     det = a .* c - b .* b;
0220     s   = b .* e - c .* d;
0221     t   = b .* d - a .* e;
0222     
0223     % compute region index (one for each face)
0224     regIndex = (s < 0) + 2 * (t < 0) + 4 * (s + t > det);
0225     
0226     % for each region, process all faces whose projection fall within it
0227     
0228     % region 0
0229     % the minimum distance occurs inside the triangle
0230     inds = regIndex == 0;
0231     s(inds) = s(inds) ./ det(inds);
0232     t(inds) = t(inds) ./ det(inds);
0233     
0234 
0235     % region 1 (formerly region 3)
0236     % The minimum distance must occur on the line s = 0
0237     inds = find(regIndex == 1);
0238     s(inds) = 0;
0239     t(inds(e(inds) >= 0)) = 0;
0240     
0241     inds2 = inds(e(inds) < 0);
0242     bool3 = c(inds2) <= -e(inds2);
0243     t(inds2(bool3)) = 1;
0244     inds3 = inds2(~bool3);
0245     t(inds3) = -e(inds3) ./ c(inds3);
0246     
0247     
0248     % region 2 (formerly region 5)
0249     % The minimum distance must occur on the line t = 0
0250     inds = find(regIndex == 2);
0251     t(inds) = 0;
0252     s(inds(d(inds) >= 0)) = 0;
0253 
0254     inds2 = inds(d(inds) < 0);
0255     bool3 = a(inds2) <= -d(inds2);
0256     s(inds2(bool3)) = 1;
0257     inds3 = inds2(~bool3);
0258     s(inds3) = -d(inds3) ./ a(inds3);
0259 
0260 
0261     % region 3 (formerly region 4)
0262     % The minimum distance must occur
0263     % * on the line t = 0
0264     % * on the line s = 0 with t >= 0
0265     % * at the intersection of the two lines
0266     inds = find(regIndex == 3);
0267     
0268     inds2 = inds(d(inds) < 0);
0269     % minimum on edge t = 0 with s > 0.
0270     t(inds2) = 0;
0271     bool3 = a(inds2) <= -d(inds2);
0272     s(inds2(bool3)) = 1;
0273     inds3 = inds2(~bool3);
0274     s(inds3) = -d(inds3) ./ a(inds3);
0275       
0276     inds2 = inds(d(inds) >= 0);
0277     % minimum on edge s = 0
0278     s(inds2) = 0;
0279     bool3 = e(inds2) >= 0;
0280     t(inds2(bool3)) = 0;
0281     bool3 = e(inds2) < 0 & c(inds2) <= e(inds2);
0282     t(inds2(bool3)) = 1;
0283     bool3 = e(inds2) < 0 & c(inds2) > e(inds2);
0284     inds3 = inds2(bool3);
0285     t(inds3) = -e(inds3) ./ c(inds3);
0286     
0287 
0288     % region 4 (formerly region 1)
0289     % The minimum distance must occur on the line s + t = 1
0290     inds = find(regIndex == 4);
0291     numer = (c(inds) + e(inds)) - (b(inds) + d(inds));
0292     s(inds(numer <= 0)) = 0;
0293     inds2 = inds(numer > 0);
0294     numer = numer(numer > 0);
0295     denom = a(inds2) - 2 * b(inds2) + c(inds2);
0296     s(inds2(numer > denom)) = 1;
0297     bool3 = numer <= denom;
0298     s(inds2(bool3)) = numer(bool3) ./ denom(bool3);
0299     t(inds) = 1 - s(inds);
0300 
0301 
0302     % Region 5 (formerly region 2)
0303     % The minimum distance must occur:
0304     % * on the line s + t = 1
0305     % * on the line s = 0 with t <= 1
0306     % * or at the intersection of the two (s=0; t=1)
0307     inds = find(regIndex == 5);
0308     tmp0 = b(inds) + d(inds);
0309     tmp1 = c(inds) + e(inds);
0310     
0311     % minimum on edge s+t = 1, with s > 0
0312     bool2 = tmp1 > tmp0;
0313     inds2 = inds(bool2);
0314     numer = tmp1(bool2) - tmp0(bool2);
0315     denom = a(inds2) - 2 * b(inds2) + c(inds2);
0316     bool3 = numer < denom;
0317     s(inds2(~bool3)) = 1;
0318     inds3 = inds2(bool3);
0319     s(inds3) = numer(bool3) ./ denom(bool3);
0320     t(inds2) = 1 - s(inds2);
0321     
0322     % minimum on edge s = 0, with t <= 1
0323     inds2 = inds(~bool2);
0324     s(inds2) = 0;
0325     t(inds2(tmp1(~bool2) <= 0)) = 1;
0326     t(inds2(tmp1(~bool2) > 0 & e(inds2) >= 0)) = 0;
0327     inds3 = inds2(tmp1(~bool2) > 0 & e(inds2) < 0);
0328     t(inds3) = -e(inds3) ./ c(inds3);
0329 
0330         
0331     % region 6 (formerly region 6)
0332     % The minimum distance must occur
0333     % * on the line s + t = 1
0334     % * on the line t = 0, with s <= 1
0335     % * at the intersection of the two lines (s=1; t=0)
0336     inds = find(regIndex == 6);
0337     tmp0 = b(inds) + e(inds);
0338     tmp1 = a(inds) + d(inds);
0339     
0340     % minimum on edge s+t=1, with t > 0
0341     bool2 = tmp1 > tmp0;
0342     inds2 = inds(bool2);
0343     numer = tmp1(bool2) - tmp0(bool2);
0344     denom = a(inds2) - 2 * b(inds2) + c(inds2);
0345     bool3 = numer <= denom;
0346     t(inds2(~bool3)) = 1;
0347     inds3 = inds2(bool3);
0348     t(inds3) = numer(bool3) ./ denom(bool3);
0349     s(inds2) = 1 - t(inds2);
0350 
0351     % minimum on edge t = 0 with s <= 1
0352     inds2 = inds(~bool2);
0353     t(inds2) = 0;
0354     s(inds2(tmp1(~bool2) <= 0)) = 1;
0355     s(inds2(tmp1(~bool2) > 0 & d(inds2) >= 0)) = 0;
0356     inds3 = inds2(tmp1(~bool2) > 0 & d(inds2) < 0);
0357     s(inds3) = -d(inds3) ./ a(inds3);
0358     
0359 
0360     % compute coordinates of closest point on plane
0361     projList = p1 + bsxfun(@times, s, v12) + bsxfun(@times, t, v13);
0362     
0363     % squared distance between point and closest point on plane
0364     [dist(i), ind] = min(sum((bsxfun(@minus, point(i,:), projList)).^2, 2));
0365 
0366     % keep the valid projection
0367     proj(i, :) = projList(ind,:);
0368 end
0369 
0370 % convert squared distance to distance
0371 dist = sqrt(dist);
0372 
0373 end
distancePointMesh

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
