Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
2  
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

2 Answers 2

up vote 12 down vote accepted

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

share|improve this answer
    
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.

share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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