www.gusucode.com > in_polyhedron > in_polyhedron_test.m
%% Tutorial and tests of IN_POLYHEDRON function % *By Jarek Tuszynski* (jaroslaw.w.tuszynski@leidos.com) % % IN_POLYHEDRON tests if points are inside a 3D triangulated surface % (faces/vertices) or volume (tetrahedrals/vertices). There are NO % assumptions about orientation of the face normals. % % IN = INPOLYHEDRON(X,POINTS) tests if the query points (POINTS) are inside % the surface/polyhedron defined by X. X can be a structure with fields % 'vertices' and 'faces' or an object of MATLAB triangulation class. In % case of triangulation class object we will only use the outside % boundary. POINTS is an N-by-3 set of XYZ coordinates. IN is an N-by-1 % logical vector which will be TRUE for each query point inside the surface. % % INPOLYHEDRON(FACES,VERTICES,POINTS) takes faces/vertices separately, % rather than in an FV structure. % %% Algorithm % For each point do: % % # shoot a random ray out of the query point in a random direction % # for each face solve: % $\left[\begin{array}{ccc} -d_{x} & v1_{x}-v0_{x} & v2_{x}-v0_{x} \\ -d_{y} & v1_{y}-v0_{y} & v2_{y}-v0_{y} \\ -d_{z} & v1_{z}-v0_{z} & v2_{z}-v0_{z} \end{array}\right]\*\left[\begin{array}{c} t \\ u \\ v \end{array} \right]=\left[\begin{array}{c} o_{x}-v0_{x} \\ o_{y}-v0_{y} \\ o_{z}-v0_{z} \end{array}\right]$ % for $\left[\begin{array}{c} t \\ u \\ v \end{array} \right]$. % _d_ is the ray direction. Variables _u_ , _v_ are barycentric coordinates % and _t/|d|_ is the distance from the intersection point to the ray origin. % Ray/triangle intersect if all _t_, _u_, _v_ and _w_=1-u-v are positive. % # count ray / surface intersections % # even number means inside and odd mean outside % # in rare case the ray hits one of the surface faces right on the % edge repeat the process with a new ray % %% References % Based on % * "Fast, minimum storage ray-triangle intersection". Tomas M鰈ler and % Ben Trumbore. Journal of Graphics Tools, 2(1):21--28, 1997. % http://www.graphics.cornell.edu/pubs/1997/MT97.pdf (Ray/triangle % intersection) % * http://fileadmin.cs.lth.se/cs/Personal/Tomas_Akenine-Moller/raytri/raytri.c % (Ray/triangle intersection) % * Robert Sedgewick "Algorithms" (point in polygon algorithm) % %% Licence % *The function is distributed under BSD License* format compact; % viewing preference clear variables; close all; rng('shuffle'); type('license.txt') %% Test if random points are inside or outside of the volume % defined by MATLAB test object "tetmesh" load tetmesh; TR = triangulation(tet,X); [S.faces, S.vertices] = freeBoundary(TR); n = 2000; % number of points points = 80*rand(n,3) - repmat([40 40 0], n, 1); tic in1 = in_polyhedron(S, points); fprintf('Number of points inside is %i, outside is %i. Calculation time: %f sec\n', ... nnz(in1), nnz(in1==0), toc); %% Plot results figure, hold on, view(3) % Display the result set(gcf, 'Position', get(gcf, 'Position').*[0 0 1.5 1.5]) patch(S,'FaceColor','g','FaceAlpha',0.2) plot3(points( in1,1),points( in1,2),points( in1,3),'bo','MarkerFaceColor','b') plot3(points(~in1,1),points(~in1,2),points(~in1,3),'r.'), axis image legend({'volume', 'points inside', 'points outside'}, 'Location', 'southoutside') %% Compare the results to the output of similar inpolyhedron function % by Sven Holcombe (http://www.mathworks.com/matlabcentral/fileexchange/37856) % inpolyhedron function is usually faster but requires knowlege about the % face normals. if exist('inpolyhedron.m', 'file') tic in2 = inpolyhedron(S, points); fprintf('Number of points inside is %i, outside is %i. Calculation time: %f sec\n', ... nnz(in1), nnz(in1==0), toc); fprintf('Number of differences is %i\n', sum(in1~=in2)); end %% Flip 50% of face normals and repeat msk = rand(size(S.faces,1),1) > 0.5; S.faces(msk,:) = fliplr(S.faces(msk,:)); in3 = in_polyhedron(S, points); fprintf('Number of differences for in_polyhedron is %i\n', sum(in1~=in3)); if exist('inpolyhedron.m', 'file') in4 = inpolyhedron(S, points); fprintf('Number of differences for inpolyhedron is %i\n', sum(in1~=in4)); end