Checking if a point is in a convex 3D polyhedron

18

7

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

s0rce

Posted 2012-11-25T22:45:52.130

Reputation: 9 292

2Generating 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 – 2012-11-25T23:07:26.397

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 – 2012-11-26T04:15:38.567

Answers

18

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

ssch

Posted 2012-11-25T22:45:52.130

Reputation: 16 150

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. http://i.stack.imgur.com/Xs6dr.png

– s0rce – 2012-11-25T23:05:11.573

for future users do you think you could elaborate your answer into a complete solution? – s0rce – 2012-11-25T23:17:22.753

@s0rce Yes, got distracted playing around with it :) – ssch – 2012-11-25T23:19:47.760

My solution to the select slowness was more brute force Parallelize@Select[points, rf @@ # &]. You're method is much better. – s0rce – 2012-11-25T23:36:03.930

@OleksandrR. Feel free. – ssch – 2012-11-26T00:55:27.657

@OleksandrR, I'm curious. – s0rce – 2012-11-26T03:06:21.430

8

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]

Mathematica graphics

RunnyKine

Posted 2012-11-25T22:45:52.130

Reputation: 32 260

4

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 O'Rourke

Posted 2012-11-25T22:45:52.130

Reputation: 4 389

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 – 2012-11-26T10:51:22.830