Finding if point is in 3D poly in python

Question

I am trying to find out whether a point is in a 3D poly. I had used another script I found online to take care of a lot of the 2D problems using ray casting. I was wondering how this could be changed to work for 3D polygons. I'm not going to be looking at really strange polygons with a lot of concavity or holes or anything. Here is the 2D implementation in python:

def point_inside_polygon(x,y,poly):

    n = len(poly)
    inside =False

    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xinters:
                        inside = not inside
        p1x,p1y = p2x,p2y

    return inside

Any help would be greatly appreciated! Thank you.

Answer 1

A similar question was posed here, but with focus on efficiency.

The scipy.spatial.ConvexHull approach suggested here by @Brian and @fatalaccidents works, but get's very slow if you need to check more than one point.

Well, the most efficient solution, also comes from scipy.spatial, but makes use of Delaunay tesselation:

from scipy.spatial import Delaunay

Delaunay(poly).find_simplex(point) >= 0  # True if point lies within poly

This works, because -1 is returned by .find_simplex(point) if the point is not in any of the simplices (i.e. outside of the triangulation). (Note: it works in N dimensions, not only 2/3D.)

---

Performance comparison

First for one point:

import numpy
from scipy.spatial import ConvexHull, Delaunay

def in_poly_hull_single(poly, point):
    hull = ConvexHull(poly)
    new_hull = ConvexHull(np.concatenate((poly, [point])))
    return np.array_equal(new_hull.vertices, hull.vertices)

poly = np.random.rand(65, 3)
point = np.random.rand(3)

%timeit in_poly_hull_single(poly, point)
%timeit Delaunay(poly).find_simplex(point) >= 0

Result:

2.63 ms ± 280 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
1.49 ms ± 153 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

So the Delaunay approach is faster. But this depends on the polygon size! I found that for a polygon consisting of more than ~65 points, the Delaunay approach becomes increasingly slower, while the ConvexHull approach remains almost constant in speed.

For multiple points:

def in_poly_hull_multi(poly, points):
    hull = ConvexHull(poly)
    res = []
    for p in points:
        new_hull = ConvexHull(np.concatenate((poly, [p])))
        res.append(np.array_equal(new_hull.vertices, hull.vertices))
    return res

points = np.random.rand(10000, 3)

%timeit in_poly_hull_multi(poly, points)
%timeit Delaunay(poly).find_simplex(points) >= 0

Result:

155 ms ± 9.42 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
1.81 ms ± 106 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

So Delaunay gives an extreme speed increase; not to mention how long we have to wait for 10'000 points or more. In such case, the polygon size doesn't have a too large influence anymore.

---

In summary, Delaunay is not only much faster, but also very concise in the code.