在Python中查找点是否在3D多边形中

问题描述 投票：0回答：4

我试图找出一个点是否在 3D 多边形中。我使用了在网上找到的另一个脚本来解决许多使用光线投射的 2D 问题。我想知道如何改变它以适用于 3D 多边形。我不会看到带有很多凹面或孔或任何东西的非常奇怪的多边形。这是 python 中的 2D 实现：

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

任何帮助将不胜感激！谢谢。

python 3d polygons point-in-polygon
4个回答
6
投票

这里提出了类似的问题，但重点是效率。

@Brian
和@fatalaccidents建议的scipy.spatial.ConvexHull方法有效，但是如果您需要检查多个点，则非常慢。

嗯，最有效的解决方案，也来自

scipy.spatial
，但利用了
Delaunay
细分：

from scipy.spatial import Delaunay

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

这是可行的，因为如果该点不在任何单纯形中（即在三角剖分之外），则由 

-1
返回 
.find_simplex(point)
。 （注意：它适用于 N 维，而不仅仅是 2/3D。）

性能对比

首先一点：

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

结果：

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)

所以

Delaunay
方法更快。但这取决于多边形的大小！我发现，对于由超过 65 个点组成的多边形，
Delaunay
方法变得越来越慢，而 
ConvexHull
方法的速度几乎保持不变。

对于多个点：

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

结果：

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)

所以 

Delaunay
会带来极大的速度提升；更不用说要等多久才能达到10000点以上。在这种情况下，多边形大小不再有太大的影响。

综上所述，

Delaunay
不仅速度快很多，而且代码也非常简洁。

5
投票

我检查了 QHull 版本（从上面）和线性规划解决方案（例如，参见这个问题）。到目前为止，使用 QHull 似乎是最好的选择，尽管我可能会错过 

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)

哪里

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

产生：

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

和

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

产生

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

感谢所有发表评论的人。对于任何寻找此问题答案的人，我找到了一个适用于某些情况（但不适用于复杂情况）的方法。

我正在做的是像 shongololo 建议的那样使用 scipy.spatial.ConvexHull ，但略有不同。我正在制作点云的 3D 凸包，然后将我要检查的点添加到“新”点云中并制作新的 3D 凸包。如果它们相同，那么我假设它一定在凸包内部。如果有人有更强大的方法来做到这一点，我仍然会很感激，因为我认为这有点黑客。代码如下所示：

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

希望这对将来寻找答案的人有所帮助！谢谢！

0
投票

我简单地尝试使用新点形成的三角形面积之和。如果它是内部的，则总和将等于（有一定公差）主三角形。否则不行。在这里

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

可能效率低下。但确实能完成工作。

以下内容（我不知道如何链接@bottlenick 答案）对我不起作用。这是因为我的三角形位于 3D 平面中。可能 scipy 期望它是一个四面体并且崩溃了。也许作者@bottleNick 可以纠正我错过的内容（因为答案中提到了 2/3D）。我很乐意使用它，因为它看起来如此简单和优雅。

from scipy.spatial import Delaunay
Delaunay(poly).find_simplex(point) >= 0  # True if point lies within poly
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.