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]

Mathematica graphics

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]

Mathematica graphics

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

  • 2
    $\begingroup$ 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). $\endgroup$
    – whuber
    Nov 25, 2012 at 23:07
  • $\begingroup$ 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. $\endgroup$ Nov 26, 2012 at 4:15

3 Answers 3


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]

enter image description here

  • $\begingroup$ 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 $\endgroup$
    – s0rce
    Nov 25, 2012 at 23:05
  • $\begingroup$ for future users do you think you could elaborate your answer into a complete solution? $\endgroup$
    – s0rce
    Nov 25, 2012 at 23:17
  • $\begingroup$ @s0rce Yes, got distracted playing around with it :) $\endgroup$
    – ssch
    Nov 25, 2012 at 23:19
  • $\begingroup$ My solution to the select slowness was more brute force Parallelize@Select[points, rf @@ # &]. You're method is much better. $\endgroup$
    – s0rce
    Nov 25, 2012 at 23:36
  • $\begingroup$ @OleksandrR. Feel free. $\endgroup$
    – ssch
    Nov 26, 2012 at 0:55

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


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

Mathematica graphics


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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.