inpolyhedron - are points inside a triangulated volume?

Test if 3d points are inside a mesh. Or, voxelise a mask from a surface. Mesh can be non-convex too!

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 http://blogs.mathworks.com/pick/2013/09/06/inpolyhedron/)
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.

Example:
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
patch(fv,'FaceColor','g','FaceAlpha',0.2)
plot3(pts(in,1),pts(in,2),pts(in,3),'bo','MarkerFaceColor','b')
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
patch(fv,'FaceColor','g','FaceAlpha',0.2)
plot3(x(in), y(in), z(in),'bo','MarkerFaceColor','b')
plot3(x(~in),y(~in),z(~in),'ro'), axis image