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

Question

我试图找出一个点是否在 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

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

Answer 1

我检查了 QHull 版本（从上面）和线性规划解决方案（例如，参见 this question）。到目前为止，使用 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

Answer 2

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

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

嗯，most efficient solution，也来自scipy.spatial，但是使用了Delaunay tesselation：

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不仅速度快很多，而且代码也非常简洁。

Answer 3

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

我正在做的是像 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

希望这可以帮助将来寻找答案的人！谢谢！