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

Question

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

4

我正在尝试确定一个点是否在三维多边形内。我使用了另一个我在网上找到的脚本来处理很多二维问题，使用光线投射。我想知道如何将其更改为适用于三维多边形。我不会查看具有许多凹度或孔洞等奇怪多边形。这是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

非常感谢您的帮助！谢谢。

- fatalaccidents

2

根据你所说的3D多边形是什么意思，问题可能在于定义什么是3D多边形，以及如何通过算法来定义它。你的多边形是否在平面上，因此可以投影到2D平面上？如果不是，那么你是否指的是3D网格（像盒子或金字塔）？如果不是这样的话，如果它是一个扭曲的煎饼形状之类的东西，那么我无法想象如何准确地定义一个点是否在多边形内部。 - songololo

我希望poly是一个(x,y,z)点的列表。就像上面的代码一样，它只处理基本形状，因为连通性是假定的，所以任何凹度或类似的东西可能会破坏算法。例如，我可能有一个从球体方程或某种圆柱体、锥体或平行六面体生成的点列表。希望这能澄清问题。 - fatalaccidents

听起来你想知道一个点是否在3D网格内。（不确定如何实现，但一定是可能的。）也许可以看看 qhull，它可以做凸包，它可能有在 Python 中暴露的功能，让你检查一个点是否在凸包内。（虽然这对于凸缺陷是行不通的。） - songololo

为了完整性...可以在http://geospatialpython.com/2011/01/point-in-polygon.html找到发布的代码和相关信息。 - sancho.s ReinstateMonicaCellio

3个回答

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

- Brian

1

感谢所有评论过的人。对于任何寻找答案的人，我已经找到了一个适用于某些情况（但不是复杂情况）的解决方案。

我的做法是使用类似shongololo建议的scipy.spatial.ConvexHull，但稍微有所改动。我正在制作一个点云的三维凸包，然后将我要检查的点添加到“新”点云中，并制作一个新的三维凸包。如果它们是相同的，则我认为它必须在凸包内。如果有更稳健的方法，请告诉我，因为我认为这有点不正规。代码如下：

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

希望这能帮助到未来有类似问题的人！谢谢！

- fatalaccidents

好的想法。对于大多数情况应该都适用，但不适用于具有凸性的对象。我看到scipy.spatial.ConvexHull有一个add_points方法，这可能让您消除new_pts行，并改为执行以下操作：new_hull = hull.add_points(new_pt)...（也许比从头开始创建凸包更有效？不过请检查它是否直接编辑原始数据，否则会使您的真实比较出现问题...） - songololo

2

这不起作用是因为 == 运算符尚未实现，以便检查所有点是否对应。演示它不起作用的简单方法是运行："from copy import deepcopy; ConvexHull(vertices) == deepcopy(ConvexHull(vertices))"，这将返回 False。 - Christian O'Reilly

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- BottleNick · Accepted Answer

有类似的问题被提出这里，但是重点在于效率。

scipy.spatial.ConvexHull的方法建议由@Brian和@fatalaccidents提出，虽然可以工作，但是如果需要检查多个点，则会变得非常慢。

嗯，最有效的解决方案同样来自scipy.spatial，但利用了Delaunay剖分：

from scipy.spatial import Delaunay

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

这有效，因为如果点不在任何单纯形内（即在三角剖分外部），则.find_simplex(point)将返回-1。（注意：它适用于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 可以极大地提高速度；更不用说我们需要等待 10,000 个或更多点的时间有多长了。在这种情况下，多边形的大小对结果没有太大的影响。

总之，Delaunay 不仅速度更快，而且代码非常简洁。