inpolyhedron - are points inside a triangulated volume?

Version (22.4 KB) by Sven
Test if 3d points are inside a mesh. Or, voxelise a mask from a surface. Mesh can be non-convex too!
Updated 12 Nov 2015

View License

Editor's Note: This file was selected as MATLAB Central Pick of the Week

INPOLYHEDRON Tests if points are inside a 3D triangulated (faces/vertices) surface
User's note:
inpolyhedron adopts the widely used convention that surface face normals point OUT from the object. If your faces point in, simply call inpolyhedron(...,'flipNormals',true).
(see discussion at
IN = INPOLYHEDRON(FV,QPTS) tests if the query points (QPTS) are inside the
patch/surface/polyhedron defined by FV (a structure with fields 'vertices' and
'faces'). QPTS 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,...) takes faces/vertices separately, rather than in
an FV structure.

IN = INPOLYHEDRON(..., X, Y, Z) voxelises a mask of 3D gridded query points
rather than an N-by-3 array of points. X, Y, and Z coordinates of the grid
supplied in XVEC, YVEC, and ZVEC respectively. IN will return as a 3D logical
volume with SIZE(IN) = [LENGTH(YVEC) LENGTH(XVEC) LENGTH(ZVEC)], equivalent to
syntax used by MESHGRID. INPOLYHEDRON handles this input faster and with a lower
memory footprint than using MESHGRID to make full X, Y, Z query points matrices.

INPOLYHEDRON(...,'PropertyName',VALUE,'PropertyName',VALUE,...) tests query
points using the following optional property values:

TOL - Tolerance on the tests for "inside" the surface. You can think of
tol as the distance a point may possibly lie above/below the surface, and still
be perceived as on the surface. Due to numerical rounding nothing can ever be
done exactly here. Defaults to ZERO. Note that in the current implementation TOL
only affects points lying above/below a surface triangle (in the Z-direction).
Points coincident with a vertex in the XY plane are considered INside the surface.
More formal rules can be implemented with input/feedback from users.

GRIDSIZE - Internally, INPOLYHEDRON uses a divide-and-conquer algorithm to
split all faces into a chessboard-like grid of GRIDSIZE-by-GRIDSIZE regions.
Performance will be a tradeoff between a small GRIDSIZE (few iterations, more
data per iteration) and a large GRIDSIZE (many iterations of small data
calculations). The sweet-spot has been experimentally determined (on a win64
system) to be correlated with the number of faces/vertices. You can overwrite
this automatically computed choice by specifying a GRIDSIZE parameter.

FACENORMALS - By default, the normals to the FACE triangles are computed as the
cross-product of the first two triangle edges. You may optionally specify face
normals here if they have been pre-computed.

FLIPNORMALS - (Defaults FALSE). To match a wider convention, triangle
face normals are presumed to point OUT from the object's surface. If
your surface normals are defined pointing IN, then you should set the
FLIPNORMALS option to TRUE to use the reverse of this convention.

tmpvol = zeros(20,20,20); % Empty voxel volume
tmpvol(5:15,8:12,8:12) = 1; % Turn some voxels on
tmpvol(8:12,5:15,8:12) = 1;
tmpvol(8:12,8:12,5:15) = 1;
fv = isosurface(tmpvol, 0.99); % Create the patch object
fv.faces = fliplr(fv.faces); % Ensure normals point OUT
% Test SCATTERED query points
pts = rand(200,3)*12 + 4; % Make some query points
in = inpolyhedron(fv, pts); % Test which are inside the patch
figure, hold on, view(3) % Display the result
plot3(pts(~in,1),pts(~in,2),pts(~in,3),'ro'), axis image
% Test STRUCTURED GRID of query points
gridLocs = 3:2.1:19;
[x,y,z] = meshgrid(gridLocs,gridLocs,gridLocs);
in = inpolyhedron(fv, gridLocs,gridLocs,gridLocs);
figure, hold on, view(3) % Display the result
plot3(x(in), y(in), z(in),'bo','MarkerFaceColor','b')
plot3(x(~in),y(~in),z(~in),'ro'), axis image

Cite As

Sven (2024). inpolyhedron - are points inside a triangulated volume? (, MATLAB Central File Exchange. Retrieved .

MATLAB Release Compatibility
Created with R2012a
Compatible with any release
Platform Compatibility
Windows macOS Linux

Inspired by: Inhull

Inspired: Asteroid Shape Data Explorer, in_polyhedron

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
Version Published Release Notes

Used quoted semicolons ':' inside function handle calls to conform with new 2015b syntax interpreter
Fixed title.

Fixed indeterminate behaviour for when query
points lie *exactly* in line with an "overhanging" vertex.

Put an upper limit on gridsize to prevent memory issues for huge trianulations.

v3.1 - Dropped nested unique call made redundant via v2.1 gridded point handling. Refreshed grid size selection via optimisation.

NEW CONVENTION ADOPTED to expect face normals pointing IN
Vertically oriented faces are now ignored. Results in speed-up and bug-fix.

The icon was lost last upload...

Much improved handling of gridded input. Now no longer unpacks the 3d grid. Saves big memory footprint!

Added logic to consider query points coincident with a vertex as IN (previously considered OUT)

A speedup tweak to distribute facets to grid locations - way to go, accumarray()