Checking if a point is in a convex 3D polyhedron

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"]

-
Generating convex polyhedron from face planes is closely related because when the planes' normal vectors are consistently oriented, they provide a simple mechanism for determining whether points are inside or outside (as illustrated, for instance, in the code in my answer in that thread: see the argument to the RegionPlot3D example there). –  whuber Nov 25 '12 at 23:07
The RedionFunction approach is quite good. In cases where there are many small faces you might instead consider shooting a ray to the outside and counting intersections. For this to be effective you'd need to bin faces in such a way that most never get tested for intersection, that is, only "reasonable contenders" get tested. –  Daniel Lichtblau Nov 26 '12 at 4:15

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]


-
I just found that in the documentation right as I saw your answer. I had been looking around for a while. I was reading about ray-casting algorithms and hoping there would be an easier way. i.stack.imgur.com/Xs6dr.png –  s0rce Nov 25 '12 at 23:05
for future users do you think you could elaborate your answer into a complete solution? –  s0rce Nov 25 '12 at 23:17
@s0rce Yes, got distracted playing around with it :) –  ssch Nov 25 '12 at 23:19
My solution to the select slowness was more brute force Parallelize@Select[points, rf @@ # &]. You're method is much better. –  s0rce Nov 25 '12 at 23:36
@OleksandrR. Feel free. –  ssch Nov 26 '12 at 0:55

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.

-
Joseph, some Mathematica implementation would enhance this answer no end. Since you´ve done it already in C, this should be a picnic ;-) –  Yves Klett Nov 26 '12 at 10:51

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]


-