Updated 3 Nov 2022
A fast function to check which of a set of 3D-points on a grid are inside and which are outside of one or more closed surfaces defined by a polyhedron. Written in C++ using the Matlab mex interface.
Unlike other point-in-polyhedron functions currently on the Matlab file exchange, this function requires that the points to be checked are arranged on a 3D-grid. This makes the function a lot faster. See the end of this description for benchmark results.
To compile, either open the project in Visual Studio and compile using the "Release mex" configuration, or compile directly from Matlab:
mex -R2018a FindInsideOfPolyhedron.cpp insidepoly_mexfunction.cpp -output InsidePolyhedron
I = InsidePolyhedron(V, F, x, y, z) returns a 3-dimensional boolean array where I(i, j, k) is true if and only if the point (y(i), x(j), z(k)) is inside the polyhedron surface. The dimensions of I are (length(y), length(x), length(z)), which corresponds to the size of the arrays returned by meshgrid(x, y, z).
I = InsidePolyhedron(V, F, x, y, z, useDithering) adds a tiny bit of noise to the vertex coordinates if useDithering is True. This can usually solve issues where one or more of the triangular faces are parallel to the traced ray, which in other circumstances would make the ray-tracing algorithm fail. This problem may occur when the vertices of the surface are themselves aligned on a grid. Points that lie exactly on (or numerically indistingushable from) the surface will end up being randomly assigned as inside or outside the surface.
V is an n by 3 matrix of vertices on the surface. Each row consists of the x, y and z coordinates of a vertex. F is an m by 3 matrix of indices into the rows of V. Each row of F defines a triangular face defined by three vertices. This is the same format as returned by surf = isosurface(...), which returns a struct with the fields surf.vertices and surf.faces.
The arguments x, y, and z are vectors of x, y and z coordinates respectively. This function exploits the fact that the points to be checked are on a grid in order to speed up processing significantly. The points checked by InsidePolyhedron are those that are on the grid defined by [X, Y, Z] = meshgrid(x, y, z). Note, however, that the input arguments are the vectors x, y and z, not the 3-dimensional arrays X, Y and Z.
dx = 0.1;
x = -7:dx:7;y = -7:dx:7; z = -7:dx:7;
[X, Y, Z] = meshgrid(x, y, z);
%Create a logical array where points inside a radius of 5 are assigned true.
I0 = sqrt(X.^2 + Y.^2 + Z.^2) < 5;
%Turn it into a smooth function for use in isosurface
I0smooth = smooth3(I0, 'box', 5);
%Find the isosurface as a polyhedron
surface = isosurface(X, Y, Z, I0smooth, 0.5);
%Use InsidePolyhedron to find points inside and outside the surface.
%Since the grid we check is the same that we used to generate the
%surface, we will need to use dithering to get a good result.
I1 = InsidePolyhedron(surface.vertices, surface.faces, x, y, z, true);
%Show a slice of the result
imagesc(x, y, I1(:,:,50));
The function relies on ray-tracing. For every point, a ray extending from this point in any direction will cross an odd number of faces if the point is inside the surface, and an even number of faces if it is outside. This technique will work even if the structure defines multiple surfaces, but all surfaces must be closed.
The key to the speedup of checking points on a grid compared to checking points individually comes from reducing the number of faces we need to check for each point. For every plane (x, y, z0) , we need check only those faces that cross z0, and for every line (x, y0, z0) in that plane we need only to check those faces that also cross y0. Finally, we only need to check where each of these faces crosses the line once in order to check every point on the line.
Benchmarking with a large grid (~300 M points) showed a speedup of a factor of 250 over "Pointsinpolyhedron" (https://se.mathworks.com/matlabcentral/fileexchange/47909-pointsinpolyhedron-test-if-points-are-in-polyhedron) and a factor 30 over "inpolyhedron" (https://se.mathworks.com/matlabcentral/fileexchange/37856-inpolyhedron-are-points-inside-a-triangulated-volume). On a moderately sized grid (88 by 71 by 61 = 381128 points), the speedup was a factor 30 and 25, respectively.
If you are having trouble compiling or using the function, please do not hesitate to leave a comment.
Øyvind Rørtveit (2023). InsidePolyhedron (https://github.com/BergenParticleTherapy/InsidePolyhedron), GitHub. Retrieved .
MATLAB Release Compatibility
Platform CompatibilityWindows macOS Linux
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!
Versions that use the GitHub default branch cannot be downloaded
Changed illustration image.