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