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.
解决方案
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