Checking if a point is in a convex 3D polyhedron

Question

Extending from these questions How to check if a 3D point is in a planar polygon? and How to check if a 2D point is in a polygon?.

I'm trying to do this to render specific shapes made up of spheres.

For a sphere it is easy:

(* generate a grid of points*)
d = 20;
points = Table[{x, y, z}, {x, -d, d}, {y, -d, d}, {z, -d, d}]~Flatten~2;

(*check if they are inside a spherical shell*)
points2 = Select[points, 8.5 < Norm[#] < 10 &];

(*render the spheres that are inside the spherical shell*)
Graphics3D[{Sphere[#, 0.75]} & /@ points2, Boxed -> False]

However, I want to try and do the same rendering for other shapes, for example a pentagonal dipyramid

Graphics3D[{Lighter@Lighter@Blue, Opacity[.8], EdgeForm[Thickness[ 0.005]], PolyhedronData["PentagonalDipyramid", "Faces"]}, Boxed -> False]

However, I'm not sure how to check if the points on the grid are with in the polyhedron.

I can access the faces... PolyhedronData["PentagonalDipyramid", "Faces"]

Answer 1

From the Mathematica documentation for PolyhedronData (see Coordinate-related properties under "More information")

RegionFunction – pure function giving True in the interior of the polyhedron.

PolyhedronData["PentagonalDipyramid", "RegionFunction"]

As I was playing around with this I noticed the Select part was very slow when running on lots of points, but Compile sorted that right out

d = 0.04;
points = Table[{x, y, z}, {x, -1, 1, d}, {y, -1, 1, d}, {z, -1, 1, d}]~Flatten~2;

inDipyramidQ = With[{
    rf = PolyhedronData["PentagonalDipyramid", "RegionFunction"][x, y, z]
   },
   Compile[{{pt, _Real, 1}},
    Block[{x = pt[[1]], y = pt[[2]], z = pt[[3]]}, If[#, 1, 0] &[rf]],
    Parallelization -> True, RuntimeAttributes -> Listable
   ]
  ];

(* 0.0624s versus uncompiled inPyramidQ that I aborted after >1min *)
points2 = Pick[points, inDipyramidQ[points], 1]; // AbsoluteTiming

Graphics3D[{Sphere[points2, d/2]}, Boxed -> False]

Answer 2

Version 10 approach:

d = 0.04;
points = Table[{x, y, z}, {x, -1, 1, d}, {y, -1, 1, d}, {z, -1, 1, d}] ~Flatten~ 2;
region = BoundaryDiscretizeGraphics @ PolyhedronData["PentagonalDipyramid"];
rm = RegionMember[region];

Select points in the region:

pin = Pick[points, rm @ points, True];

Visualize:

Graphics3D[{Sphere[pin, d/2]}, Boxed -> False]

Answer 3

Pardon me in advance if I am not directly addressing your question. A point $p$ is in a convex polyhedron if it is "left-of" each of its faces $F$, where "left-of" is defined by the signed volume of $p$. If the polyhedron is triangulated, then $F$ is a triangle, and the key computation is the signed volume of a tetrahedron formed by $F$ and the point $p$. This is all over the Internet, and in many books, including my own, Computational Geometry in C; that link will lead you to explicit code for this computation. If $F$ is not a triangle, then it is easily triangulated, and you can proceed as above.