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
produces
CPU times: user 3.83 s, sys: 16 ms, total: 3.85 s Wall time: 3.85 s
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
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) >= 0Result:
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)
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 FalseIn summary, Delaunay is not only much anycodings_3d faster, but also very concise in the anycodings_3d code.,Thanks to all that have commented. For anycodings_3d anyone looking for an answer for this, I anycodings_3d have found one that works for some cases anycodings_3d (but not complicated cases). ,This works, because -1 is returned by anycodings_3d .find_simplex(point) if the point is not anycodings_3d in any of the simplices (i.e. outside of anycodings_3d the triangulation). (Note: it works in N anycodings_3d dimensions, not only 2/3D.),I am trying to find out whether a point is anycodings_polygons in a 3D poly. I had used another script I anycodings_polygons found online to take care of a lot of the 2D anycodings_polygons problems using ray casting. I was wondering anycodings_polygons how this could be changed to work for 3D anycodings_polygons polygons. I'm not going to be looking at anycodings_polygons really strange polygons with a lot of anycodings_polygons concavity or holes or anything. Here is the anycodings_polygons 2D implementation in python:
I am trying to find out whether a point is anycodings_polygons in a 3D poly. I had used another script I anycodings_polygons found online to take care of a lot of the 2D anycodings_polygons problems using ray casting. I was wondering anycodings_polygons how this could be changed to work for 3D anycodings_polygons polygons. I'm not going to be looking at anycodings_polygons really strange polygons with a lot of anycodings_polygons concavity or holes or anything. Here is the anycodings_polygons 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 insideI checked out the QHull version (from anycodings_3d above) and the linear programming anycodings_3d solution (e.g. see this question). So anycodings_3d far, using QHull seems to be the best anycodings_3d bet, although I might be missing some anycodings_3d 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
produces
CPU times: user 3.83 s, sys: 16 ms, total: 3.85 s Wall time: 3.85 s
Well, the most efficient solution, also anycodings_3d comes from scipy.spatial, but makes use anycodings_3d of Delaunay tesselation:
from scipy.spatial import Delaunay Delaunay(poly).find_simplex(point) >= 0 # True if point lies within poly
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) >= 0Result:
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)
What I am doing is using anycodings_3d scipy.spatial.ConvexHull like shongololo anycodings_3d suggested, but with a slight twist. I anycodings_3d am making a 3D convex hull of the point anycodings_3d cloud, and then adding the point I am anycodings_3d checking into a "new" point cloud and anycodings_3d making a new 3D convex hull. If they anycodings_3d are identical then I am assuming that it anycodings_3d must be inside the convex hull. I would anycodings_3d still appreciate if someone has a more anycodings_3d robust way of doing this, as I see this anycodings_3d as a bit hackish. The code would look anycodings_3d 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 False20th January, 2016: Initial version
gp = [p1, p2, p3, p4, p5, p6, p7, p8]; gpInst = GeoPolygon(gp); procInst = GeoPolygonProc(gpInst);
if (x > procInst.x0 and x < procInst.x1 and y > procInst.y0 and y < procInst.y1 and z > procInst.z0 and z < procInst.z1):
if (procInst.PointInside3DPolygon(x, y, z)):class GeoPoint(object):
def __init__(self, x, y, z):
self.x = x
self.y = y
self.z = z
def __add__(self, p):
return GeoPoint(self.x + p.x, self.y + p.y, self.z + p.z)from GeoVector import * class GeoPlane(object): def __init__(self, a, b, c, d): self.a = a self.b = b self.c = c self.d = d @staticmethod def Create(p0, p1, p2): # p0, p1, p2: GeoPoint v = GeoVector(p0, p1) u = GeoVector(p0, p2) # normal vector n = u * v; a = n.x b = n.y c = n.z d = -(a * p0.x + b * p0.y + c * p0.z) return GeoPlane(a, b, c, d) def __neg__(self): return GeoPlane(-self.a, -self.b, -self.c, -self.d) def __mul__(self, pt): # pt: GeoPoint return (self.a * pt.x + self.b * pt.y + self.c * pt.z + self.d)
class GeoFace(object):
# ptList: vertice GeoPoint Array
# idxList: vertice index Integer Array
def __init__(self, ptList, idxList):
#alloc instance array memory
self.__v = [] # vertice point array
self.__idx = [] # vertice index array
self.__n = len(ptList) # number of vertices
for i in range(0, self.__n):
self.__v.append(ptList[i])
self.__idx.append(idxList[i])
@property
def v(self):
return self.__v
@property
def idx(self):
return self.__idx
@property
def n(self):
return self.__nLast Updated : 26 Jul, 2022
No Yes Yes Yes No No
How to check if point is inside a polygon?,It is also possible to do PIP other way around, i.e. to check if polygon contains a point:,Let’s check if those points are within the polygon,using a function called .within() that checks if a point is within a polygon
from shapely.geometry import Point, Polygon # Create Point objects p1 = Point(24.952242, 60.1696017) p2 = Point(24.976567, 60.1612500) # Create a Polygon coords = [(24.950899, 60.169158), (24.953492, 60.169158), (24.953510, 60.170104), (24.950958, 60.169990)] poly = Polygon(coords)
# Let 's check what we have In[1]: print(p1) POINT(24.952242 60.1696017) In[2]: print(p2) POINT(24.976567 60.16125) In[3]: print(poly) POLYGON((24.950899 60.169158, 24.953492 60.169158, 24.95351 60.170104, 24.950958 60.16999, 24.950899 60.169158))
# Check
if p1 is within the polygon using the within
function
In[4]: p1.within(poly)
Out[4]: True
# Check
if p2 is within the polygon
In[5]: p2.within(poly)
Out[5]: False# Our point
In[6]: print(p1)
POINT(24.952242 60.1696017)
# The centroid
In[7]: print(poly.centroid)
POINT(24.95224242849236 60.16960179038188)# Does polygon contain p1 ?
In[8] : poly.contains(p1)
Out[8]: True
# Does polygon contain p2 ?
In[9] : poly.contains(p2)
Out[9]: Falsefrom shapely.geometry import LineString, MultiLineString # Create two lines line_a = LineString([(0, 0), (1, 1)]) line_b = LineString([(1, 1), (0, 2)])