POINT_IN_TRIANGLE tests points for membership in triangles
POINT_IN_TRIANGLE(P, E, V)
POINT_IN_TRIANGLE(P, E, V, epsilon)
Inputs:
P - [p x 2] or [p x 3] list of point coordinates
V - [v x 2] or [v x 3] list of triangle vertices
E - [e x 3] matrix specifying the vertices V consituting each of the
e triangles
epsilon - threshold to counteract numerical instability:
epsilon > 0 makes the triangle bigger so points on edges
and vertices are correctly classified as inside. May
contain false positives.
epsilon = 0 may miss some points on edges and vertices due
to numerical issues.
epsilon < 0 makes the trianle smaller so points on edges
and vertices are classified as outside. May miss more
points than intended.
Default: epsilon = eps
POINT_IN_TRIANGLE(..., 'match') specifies that each point is to be tested
only against its corresponding triangle (number of points must equal
number of triangles). Option has no effect in 2D.0001 function out = point_in_triangle(P,E,V,epsilon, str) 0002 %POINT_IN_TRIANGLE tests points for membership in triangles 0003 % POINT_IN_TRIANGLE(P, E, V) 0004 % POINT_IN_TRIANGLE(P, E, V, epsilon) 0005 % Inputs: 0006 % P - [p x 2] or [p x 3] list of point coordinates 0007 % V - [v x 2] or [v x 3] list of triangle vertices 0008 % E - [e x 3] matrix specifying the vertices V consituting each of the 0009 % e triangles 0010 % epsilon - threshold to counteract numerical instability: 0011 % epsilon > 0 makes the triangle bigger so points on edges 0012 % and vertices are correctly classified as inside. May 0013 % contain false positives. 0014 % epsilon = 0 may miss some points on edges and vertices due 0015 % to numerical issues. 0016 % epsilon < 0 makes the trianle smaller so points on edges 0017 % and vertices are classified as outside. May miss more 0018 % points than intended. 0019 % Default: epsilon = eps 0020 % 0021 % POINT_IN_TRIANGLE(..., 'match') specifies that each point is to be tested 0022 % only against its corresponding triangle (number of points must equal 0023 % number of triangles). Option has no effect in 2D. 0024 0025 0026 % (C) 2012-2015 Bartlomiej Grychtol 0027 % License: GPL version 2 or 3 0028 % $Id: point_in_triangle.m 5034 2015-05-26 15:09:20Z bgrychtol $ 0029 0030 switch nargin 0031 case {1 2 3} 0032 epsilon = eps; 0033 str = []; 0034 case 4 0035 if ischar(epsilon) 0036 str = epsilon; 0037 epsilon = eps; 0038 else 0039 str = []; 0040 end 0041 end 0042 0043 0044 match = strcmp(str,'match'); 0045 0046 switch size(P,2) 0047 case 2 0048 out = point_in_triangle_2d(P,E,V,epsilon); 0049 0050 case 3 0051 [u, v] = point_in_triangle_3d_wrapper(P,E,V,match); 0052 out= test(u,v,epsilon); 0053 0054 otherwise 0055 error('EIDORS:WrongInput','Points must be 2D or 3D'); 0056 end 0057 0058 0059 function out = test(u,v,epsilon) 0060 out = u >= -epsilon & v >= -epsilon & (u+v-epsilon) <= 1; 0061 0062 function [u, v] = point_in_triangle_3d_wrapper(P,E,V, match) 0063 nPts = size(P,1); 0064 nTri = size(E,1); 0065 0066 if nPts == 1 || match 0067 [u, v] = point_in_triangle_3d(P,E,V); 0068 else 0069 u = zeros(nPts, nTri); 0070 v = zeros(nPts, nTri); 0071 for i = 1:nPts 0072 [u(i,:), v(i,:)] = point_in_triangle_3d(P(i,:),E,V); 0073 end 0074 0075 end 0076 0077 %-------------------------------------------------------------------------% 0078 % Decide if point P is in triangles E indexing nodes V in 3D 0079 function [u, v] = point_in_triangle_3d(p, E, V) 0080 %http://www.blackpawn.com/texts/pointinpoly/default.html 0081 % vectors 0082 v0 = V(E(:,3),:) - V(E(:,1),:); 0083 v1 = V(E(:,2),:) - V(E(:,1),:); 0084 v2 = bsxfun(@minus,p, V(E(:,1),:)); 0085 0086 % dot products 0087 dot00 = dot(v0, v0, 2); 0088 dot01 = dot(v0, v1, 2); 0089 dot02 = dot(v0, v2, 2); 0090 dot11 = dot(v1, v1, 2); 0091 dot12 = dot(v1, v2, 2); 0092 0093 % barycentric coordinates 0094 invDenom = 1 ./ (dot00 .* dot11 - dot01 .* dot01); 0095 u = (dot11 .* dot02 - dot01 .* dot12) .* invDenom; 0096 v = (dot00 .* dot12 - dot01 .* dot02) .* invDenom; 0097 0098 0099 %-------------------------------------------------------------------------% 0100 % Decide if points P are in triangles E indexing nodes V in 2D 0101 function out = point_in_triangle_2d(P,E,V,epsilon) 0102 X = reshape(V(E,1),size(E)); 0103 Y = reshape(V(E,2),size(E)); 0104 T = [ bsxfun(@minus, X(:,1:2), X(:,3)) bsxfun(@minus, Y(:,1:2), Y(:,3))]; 0105 invdetT = 1./(T(:,1).*T(:,4) - T(:,2).*T(:,3)); 0106 y23 = Y(:,2) - Y(:,3); 0107 y31 = Y(:,3) - Y(:,1); 0108 x32 = X(:,3) - X(:,2); 0109 x13 = X(:,1) - X(:,3); 0110 nPts = size(P,1); 0111 out = sparse(size(E,1),nPts); 0112 for i = 1:nPts 0113 p1X = P(i,1) - X(:,3); 0114 p2Y = P(i,2) - Y(:,3); 0115 u = (y23 .* p1X + x32 .* p2Y) .* invdetT; 0116 v = (y31 .* p1X + x13 .* p2Y) .* invdetT; 0117 out(:,i) = test(u,v,epsilon); 0118 end 0119 out = out'; 0120 0121 % u = ( repmat(Y(:,2)-Y(:,3),1,nPts).*bsxfun(@minus,P(:,1)',X(:,3)) ... 0122 % + repmat(X(:,3)-X(:,2),1,nPts).*bsxfun(@minus,P(:,2)',Y(:,3)) )... 0123 % .* repmat(invdetT,1,nPts); 0124 % v = ( repmat(Y(:,3)-Y(:,1),1,nPts).*bsxfun(@minus,P(:,1)',X(:,3)) ... 0125 % + repmat(X(:,1)-X(:,3),1,nPts).*bsxfun(@minus,P(:,2)',Y(:,3)) )... 0126 % .* repmat(invdetT,1,nPts); 0127 % u = u'; 0128 % v = v'; 0129 % 0130