point_in_triangle

PURPOSE ^

POINT_IN_TRIANGLE tests points for membership in triangles

SYNOPSIS ^

function out = point_in_triangle(P,E,V,epsilon, str)

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Fri 30-Dec-2022 19:44:54 by m2html © 2005