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

I tried simply by using sum of area of triangles formed by the new point. If it is interior, the sum will equal (with some tolerance) the master triangle. Else no. Here it is

import numpy as np
def triangleArea3d(pts):
# Find area of a triangle in 3D plane
# triangle is defined by pts=np.array([x1,y1,z1],[x2,y2,z2],[x3,y3,z3])
    a=0
    b=1
    c=2
    dr=np.cross(pts[a]-pts[b],pts[a]-pts[c])
    area=np.linalg.norm(dr)
    return area

def pointInTri3d(poly,pt):
    # use tolerance
    tol=1e-15
    n=poly.shape[0]
    amaster=triangleArea3d(poly)
    tri0=np.concatenate((poly[[1,2]], pt), axis=0)
    a0=triangleArea3d(tri0)
    tri1=np.concatenate((poly[[0,2]], pt), axis=0)
    a1=triangleArea3d(tri1)
    tri2=np.concatenate((poly[[0,1]], pt), axis=0)
    a2=triangleArea3d(tri2)
    asub=a0+a1+a2
    # if asub==amaster:
    if abs(asub - amaster) < tol:
        flag=1
    else:
        flag=0
    return flag,asub,amaster

Probably inefficient. But does the job.

The following (i dont know how to link @bottlenick answer) did not work for me. It is because my triangle was in the 3D plane. Probably scipy expects it to be a tetrahedron and gets crashed. May be the author @bottleNick can correct me what I missed (since 2/3D is mentioned in the answer). I would be happy to use that since it looks so simple and elegant.

from scipy.spatial import Delaunay
Delaunay(poly).find_simplex(point) >= 0  # True if point lies within poly

Answer 2

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.

Answer 3

I checked out the QHull version (from above) and the linear programming solution (e.g. see this question). So far, using QHull seems to be the best bet, although I might be missing some optimizations with the scipy.spatial LP.

import numpy
import numpy.random
from numpy import zeros, ones, arange, asarray, concatenate
from scipy.optimize import linprog

from scipy.spatial import ConvexHull

def pnt_in_cvex_hull_1(hull, pnt):
    '''
    Checks if `pnt` is inside the convex hull.
    `hull` -- a QHull ConvexHull object
    `pnt` -- point array of shape (3,)
    '''
    new_hull = ConvexHull(concatenate((hull.points, [pnt])))
    if numpy.array_equal(new_hull.vertices, hull.vertices): 
        return True
    return False


def pnt_in_cvex_hull_2(hull_points, pnt):
    '''
    Given a set of points that defines a convex hull, uses simplex LP to determine
    whether point lies within hull.
    `hull_points` -- (N, 3) array of points defining the hull
    `pnt` -- point array of shape (3,)
    '''
    N = hull_points.shape[0]
    c = ones(N)
    A_eq = concatenate((hull_points, ones((N,1))), 1).T   # rows are x, y, z, 1
    b_eq = concatenate((pnt, (1,)))
    result = linprog(c, A_eq=A_eq, b_eq=b_eq)
    if result.success and c.dot(result.x) == 1.:
        return True
    return False


points = numpy.random.rand(8, 3)
hull = ConvexHull(points, incremental=True)
hull_points = hull.points[hull.vertices, :]
new_points = 1. * numpy.random.rand(1000, 3)

where

%%time
in_hull_1 = asarray([pnt_in_cvex_hull_1(hull, pnt) for pnt in new_points], dtype=bool)

produces:

CPU times: user 268 ms, sys: 4 ms, total: 272 ms
Wall time: 268 ms

and

%%time
in_hull_2 = asarray([pnt_in_cvex_hull_2(hull_points, pnt) for pnt in new_points], dtype=bool)

produces

CPU times: user 3.83 s, sys: 16 ms, total: 3.85 s
Wall time: 3.85 s

Answer 4

Thanks to all that have commented. For anyone looking for an answer for this, I have found one that works for some cases (but not complicated cases).

What I am doing is using scipy.spatial.ConvexHull like shongololo suggested, but with a slight twist. I am making a 3D convex hull of the point cloud, and then adding the point I am checking into a "new" point cloud and making a new 3D convex hull. If they are identical then I am assuming that it must be inside the convex hull. I would still appreciate if someone has a more robust way of doing this, as I see this as a bit hackish. The code would look something like as follows:

from scipy.spatial import ConvexHull

def pnt_in_pointcloud(points, new_pt):
    hull = ConvexHull(points)
    new_pts = points + new_pt
    new_hull = ConvexHull(new_pts)
    if hull == new_hull: 
        return True
    else:
        return False

Hopefully this helps someone in the future looking for an answer! Thanks!