Python 找到点是否位于点云的凸包中的有效方法是什么？

What's an efficient way to find if a point lies in the convex hull of a point cloud?

提问by AME

I have a point cloud of coordinates in numpy. For a high number of points, I want to find out if the points lie in the convex hull of the point cloud.

我在numpy中有一个坐标点云。对于大量点，我想找出这些点是否位于点云的凸包中。

I tried pyhull but I cant figure out how to check if a point is in the ConvexHull:

我试过 pyhull，但我不知道如何检查一个点是否在ConvexHull：

hull = ConvexHull(np.array([(1, 2), (3, 4), (3, 6)]))
for s in hull.simplices:
    s.in_simplex(np.array([2, 3]))

raises LinAlgError: Array must be square.

引发 LinAlgError: Array must be square。

回答by kiriloff

If you want to keep with scipy, you have to convex hull (you did so)

如果你想保持 scipy，你必须凸包（你这样做了）

>>> from scipy.spatial import ConvexHull
>>> points = np.random.rand(30, 2)   # 30 random points in 2-D
>>> hull = ConvexHull(points)

then build the list of points on the hull. Here is the code from doc to plot the hull

然后在船体上建立点列表。这是 doc 中绘制船体的代码

>>> import matplotlib.pyplot as plt
>>> plt.plot(points[:,0], points[:,1], 'o')
>>> for simplex in hull.simplices:
>>>     plt.plot(points[simplex,0], points[simplex,1], 'k-')

So starting from that, I would propose for computing list of points on the hull

所以从那开始，我会建议计算船体上的点列表

pts_hull = [(points[simplex,0], points[simplex,1]) 
                            for simplex in hull.simplices] 

(although i did not try)

（虽然我没试过）

And you can also come with your own code for computing the hull, returning the x,y points.

您还可以使用自己的代码来计算船体，返回 x,y 点。

If you want to know if a point from your original dataset is on the hull, then you are done.

如果您想知道原始数据集中的一个点是否在 hull 上，那么您就完成了。

I what you want is to know if a any point is inside the hull or outside, you must do a bit of work more. What you will have to do could be

我想要的是知道任何一点是在船体内部还是外部，您必须多做一些工作。你必须做的可能是

• for all edges joining two simplices of your hull: decide whether your point is above or under

• if point is below all lines, or above all lines, it is outside the hull

• 对于连接船体的两个简单的所有边：确定您的点是在上面还是在下面

• 如果点低于所有线，或高于所有线，则它在船体之外

As a speed up, as soon as a point has been above one line and below one another, it is inside the hull.

作为一种加速，只要一个点高于一条线并低于另一条线，它就在船体内部。

回答by Nuclearman

It looks like you are using a 2D point cloud, so I'd like to direct you to the inclusion testfor point-in-polygon testing of convex polygons.

看起来您正在使用 2D 点云，所以我想指导您进行凸多边形的点入多边形测试的包含测试。

Scipy's convex hull algorithm allows for finding convex hulls in 2 or more dimensions which is more complicated than it needs to be for a 2D point cloud. Therefore, I recommend using a different algorithm, such as this one. This is because as you really need for point-in-polygon testing of a convex hull is the list of convex hull points in clockwise order, and a point that is inside of the polygon.

Scipy 的凸包算法允许在 2 维或更多维中查找凸包，这比 2D 点云所需的要复杂。因此，我建议使用不同的算法，例如this one。这是因为您真正需要凸包的多边形内点测试是按顺时针顺序排列的凸包点列表，以及多边形内部的一个点。

The time performance of this approach is as followed:

该方法的时间性能如下：

• O(N log N) to construct the convex hull
• O(h) in preprocessing to calculate (and store) the wedge angles from the interior point
• O(log h) per point-in-polygon query.

• O(N log N) 构造凸包
• O(h) 在预处理中计算（和存储）来自内点的楔角
• 每个多边形点查询的 O(log h)。

Where N is the number of points in the point cloud and h is the number of points in the point clouds convex hull.

其中 N 是点云中的点数，h 是点云凸包中的点数。

回答by Juh_

Here is an easy solution that requires only scipy:

这是一个简单的解决方案，只需要 scipy：

def in_hull(p, hull):
    """
    Test if points in `p` are in `hull`

    `p` should be a `NxK` coordinates of `N` points in `K` dimensions
    `hull` is either a scipy.spatial.Delaunay object or the `MxK` array of the 
    coordinates of `M` points in `K`dimensions for which Delaunay triangulation
    will be computed
    """
    from scipy.spatial import Delaunay
    if not isinstance(hull,Delaunay):
        hull = Delaunay(hull)

    return hull.find_simplex(p)>=0

It returns a boolean array where Truevalues indicate points that lie in the given convex hull. It can be used like this:

它返回一个布尔数组，其中True值表示位于给定凸包中的点。它可以像这样使用：

tested = np.random.rand(20,3)
cloud  = np.random.rand(50,3)

print in_hull(tested,cloud)

If you have matplotlib installed, you can also use the following function that calls the first one and plots the results. For 2D data only, given by Nx2arrays:

如果您安装了 matplotlib，您还可以使用以下函数调用第一个函数并绘制结果。仅对于二维数据，由Nx2数组给出：

def plot_in_hull(p, hull):
    """
    plot relative to `in_hull` for 2d data
    """
    import matplotlib.pyplot as plt
    from matplotlib.collections import PolyCollection, LineCollection

    from scipy.spatial import Delaunay
    if not isinstance(hull,Delaunay):
        hull = Delaunay(hull)

    # plot triangulation
    poly = PolyCollection(hull.points[hull.vertices], facecolors='w', edgecolors='b')
    plt.clf()
    plt.title('in hull')
    plt.gca().add_collection(poly)
    plt.plot(hull.points[:,0], hull.points[:,1], 'o', hold=1)


    # plot the convex hull
    edges = set()
    edge_points = []

    def add_edge(i, j):
        """Add a line between the i-th and j-th points, if not in the list already"""
        if (i, j) in edges or (j, i) in edges:
            # already added
            return
        edges.add( (i, j) )
        edge_points.append(hull.points[ [i, j] ])

    for ia, ib in hull.convex_hull:
        add_edge(ia, ib)

    lines = LineCollection(edge_points, color='g')
    plt.gca().add_collection(lines)
    plt.show()    

    # plot tested points `p` - black are inside hull, red outside
    inside = in_hull(p,hull)
    plt.plot(p[ inside,0],p[ inside,1],'.k')
    plt.plot(p[-inside,0],p[-inside,1],'.r')

回答by Sildoreth

First, obtain the convex hull for your point cloud.

首先，获得点云的凸包。

Then loop over all of the edges of the convex hull in counter-clockwise order. For each of the edges, check whether your target point lies to the "left" of that edge. When doing this, treat the edges as vectors pointing counter-clockwise around the convex hull. If the target point is to the "left" of all of the vectors, then it is contained by the polygon; otherwise, it lies outside the polygon.

然后按逆时针顺序遍历凸包的所有边缘。对于每条边，检查您的目标点是否位于该边的“左侧”。执行此操作时，将边缘视为围绕凸包逆时针指向的向量。如果目标点在所有向量的“左侧”，则它包含在多边形中；否则，它位于多边形之外。

This other Stack Overflow topic includes a solution to finding which "side" of a line a point is on: Determine Which Side of a Line a Point Lies

另一个 Stack Overflow 主题包括一个解决方案，用于查找点位于线的哪一侧：确定点位于线的 哪一侧

---

---

这种方法的运行时复杂度（一旦你已经有了凸包）是 O(n)O(n)其中 n 是凸包的边数。

Note that this will work only for convex polygons. But you're dealing with a convex hull, so it should suit your needs.

请注意，这仅适用于凸多边形。但是您正在处理凸包，因此它应该适合您的需求。

It looks like you already have a way to get the convex hull for your point cloud. But if you find that you have to implement your own, Wikipedia has a nice list of convex hull algorithms here: Convex Hull Algorithms

看起来您已经有办法获得点云的凸包。但是如果你发现你必须实现你自己的，维基百科有一个很好的凸包算法列表： 凸包算法

回答by firemana

Hi I am not sure about how to use your program library to achieve this. But there is a simple algorithm to achieve this described in words:

嗨，我不确定如何使用您的程序库来实现这一点。但是有一个简单的算法可以用文字描述来实现这一点：

1. create a point that is definitely outside your hull. Call it Y
2. produce a line segment connecting your point in question (X) to the new point Y.
3. loop around all edge segments of your convex hull. check for each of them if the segment intersects with XY.
4. If the number of intersection you counted is even (including 0), X is outside the hull. Otherwise X is inside the hull.
5. if so occurs XY pass through one of your vertexes on the hull, or directly overlap with one of your hull's edge, move Y a little bit.
6. the above works for concave hull as well. You can see in below illustration (Green dot is the X point you are trying to determine. Yellow marks the intersection points.

1. 创建一个绝对在您的船体之外的点。叫它Y
2. 生成一条线段，将您的问题点 (X) 连接到新点 Y。
3. 环绕凸包的所有边缘段。如果线段与 XY 相交，请检查它们中的每一个。
4. 如果你统计的交点数为偶数（包括0），则X在船体外。否则 X 在船体内部。
5. 如果发生这种情况，XY 通过船体上的一个顶点，或直接与船体的一个边缘重叠，请稍微移动 Y。
6. 以上也适用于凹壳。您可以在下图中看到（绿点是您要确定的 X 点。黄色标记交叉点。

回答by feqwix

Just for completness, here is a poor's man solution:

为了完整起见，这是一个穷人的解决方案：

import pylab
import numpy
from scipy.spatial import ConvexHull

def is_p_inside_points_hull(points, p):
    global hull, new_points # Remove this line! Just for plotting!
    hull = ConvexHull(points)
    new_points = numpy.append(points, p, axis=0)
    new_hull = ConvexHull(new_points)
    if list(hull.vertices) == list(new_hull.vertices):
        return True
    else:
        return False

# Test:
points = numpy.random.rand(10, 2)   # 30 random points in 2-D
# Note: the number of points must be greater than the dimention.
p = numpy.random.rand(1, 2) # 1 random point in 2-D
print is_p_inside_points_hull(points, p)

# Plot:
pylab.plot(points[:,0], points[:,1], 'o')
for simplex in hull.simplices:
    pylab.plot(points[simplex,0], points[simplex,1], 'k-')
pylab.plot(p[:,0], p[:,1], '^r')
pylab.show()

The idea is simple: the vertices of the convex hull of a set of points Pwon't change if you add a point pthat falls "inside" the hull; the vertices of the convex hull for [P1, P2, ..., Pn]and [P1, P2, ..., Pn, p]are the same. But if pfalls "outside", then the vertices must change. This works for n-dimensions, but you have to compute the ConvexHulltwice.

这个想法很简单：P如果你添加一个p落在凸包“内部”的点，一组点的凸包的顶点不会改变；凸包为的顶点[P1, P2, ..., Pn]和[P1, P2, ..., Pn, p]是相同的。但是如果p落在“外面”，那么顶点必须改变。这适用于 n 维，但您必须计算ConvexHull两次。

Two example plots in 2-D:

二维中的两个示例图：

False:

错误的：

True:

真的：

回答by Charlie Brummitt

Use the equationsattribute of ConvexHull:

使用 的equations属性ConvexHull：

def point_in_hull(point, hull, tolerance=1e-12):
    return all(
        (np.dot(eq[:-1], point) + eq[-1] <= tolerance)
        for eq in hull.equations)

In words, a point is in the hull if and only if for every equation (describing the facets) the dot product between the point and the normal vector (eq[:-1]) plus the offset (eq[-1]) is less than or equal to zero. You may want to compare to a small, positive constant tolerance = 1e-12rather than to zero because of issues of numerical precision (otherwise, you may find that a vertex of the convex hull is not in the convex hull).

换句话说，当且仅当对于每个方程（描述面），点与法向量 ( eq[:-1]) 加上偏移量 ( eq[-1])之间的点积小于或等于零时，点在外壳中。tolerance = 1e-12由于数值精度问题，您可能希望与一个小的正常数而不是零进行比较（否则，您可能会发现凸包的顶点不在凸包中）。

Demonstration:

示范：

import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import ConvexHull

points = np.array([(1, 2), (3, 4), (3, 6), (2, 4.5), (2.5, 5)])
hull = ConvexHull(points)

np.random.seed(1)
random_points = np.random.uniform(0, 6, (100, 2))

for simplex in hull.simplices:
    plt.plot(points[simplex, 0], points[simplex, 1])

plt.scatter(*points.T, alpha=.5, color='k', s=200, marker='v')

for p in random_points:
    point_is_in_hull = point_in_hull(p, hull)
    marker = 'x' if point_is_in_hull else 'd'
    color = 'g' if point_is_in_hull else 'm'
    plt.scatter(p[0], p[1], marker=marker, color=color)

回答by NetSmoothMF

Based on thispost, here is my quick-and-dirty solution for a convex regions with 4 sides (you can easily extend it to more)

基于这篇文章，这是我针对 4 边凸区域的快速解决方案（您可以轻松地将其扩展到更多）

def same_sign(arr): return np.all(arr > 0) if arr[0] > 0 else np.all(arr < 0)

def inside_quad(pts, pt):
    a =  pts - pt
    d = np.zeros((4,2))
    d[0,:] = pts[1,:]-pts[0,:]
    d[1,:] = pts[2,:]-pts[1,:]
    d[2,:] = pts[3,:]-pts[2,:]
    d[3,:] = pts[0,:]-pts[3,:]
    res = np.cross(a,d)
    return same_sign(res), res

points = np.array([(1, 2), (3, 4), (3, 6), (2.5, 5)])

np.random.seed(1)
random_points = np.random.uniform(0, 6, (1000, 2))

print wlk1.inside_quad(points, random_points[0])
res = np.array([inside_quad(points, p)[0] for p in random_points])
print res[:4]
plt.plot(random_points[:,0], random_points[:,1], 'b.')
plt.plot(random_points[res][:,0], random_points[res][:,1], 'r.')

回答by Nils

I would not use a convex hull algorithm, because you do not need to compute the convex hull, you just want to check whether your point can be expressed as a convex combination of the set of points of whom a subset defines a convex hull. Moreover, finding the convex hull is computationally expensive, especially in higher dimensions.

我不会使用凸包算法，因为您不需要计算凸包，您只想检查您的点是否可以表示为一组点的凸组合，其中一个子集定义了凸包。此外，找到凸包在计算上是昂贵的，尤其是在更高维度上。

In fact, the mere problem of finding out whether a point can be expressed as a convex combination of another set of points can be formulated as a linear programming problem.

事实上，仅仅找出一个点是否可以表示为另一组点的凸组合的问题就可以表述为线性规划问题。

import numpy as np
from scipy.optimize import linprog

def in_hull(points, x):
    n_points = len(points)
    n_dim = len(x)
    c = np.zeros(n_points)
    A = np.r_[points.T,np.ones((1,n_points))]
    b = np.r_[x, np.ones(1)]
    lp = linprog(c, A_eq=A, b_eq=b)
    return lp.success

n_points = 10000
n_dim = 10
Z = np.random.rand(n_points,n_dim)
x = np.random.rand(n_dim)
print(in_hull(Z, x))

For the example, I solved the problem for 10000 points in 10 dimensions. The executions time is in the ms range. Would not want to know how long this would take with QHull.

例如，我解决了 10 个维度中 10000 个点的问题。执行时间在ms范围内。不想知道 QHull 需要多长时间。