如何高效地确定一个点是否在点云的凸包内？

Question

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

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

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

引发LinAlgError错误：数组必须是方阵。

- AME

Answer 1

这里有一个简单的解决方案，仅需使用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

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

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

print in_hull(tested,cloud)

如果你已经安装了matplotlib，你也可以使用下面这个函数，它会调用第一个函数并绘制结果。对于仅包含2D数据的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')

- Juh_

---

是否可以找到点云凸包的外部点？因为我想从距离计算中删除那些来自外部三角形且通常具有较大距离的点。 - Liwellyen

其实很简单：假设cloud是一个N行K列的数组，其中包含了K维空间中的N个点。使用ConvexHull(cloud).vertices（来自scipy.spatial）可以得到凸包上的点的索引，也就是“外部点”。 - Juh_

你可以放心地假设这是一种可靠的方法，因为它在“Delaunay.find_simplex”的文档中有解释，该方法对于凸壳外的点返回-1。 现在，如果你想要更多的控制，或者想要一个更快的算法，我建议使用@nils提供的解决方案。它更复杂，但只计算所需的内容（我没有测试过，但看起来是这样的）。 - Juh_

是的：ConvexHull没有提供合适的API。在这里，我建议使用一个方法，它可以做更多的事情，但很容易实现。请注意，我几年前就停止使用scipy了，所以它可能已经发展了。 - Juh_

在第hull = Delaunay(hull)行出现了“TypeError：float（）参数必须是字符串或数字”的错误。有什么想法吗？ - Gulzar

Answer 2

我不会使用凸包算法，因为你不需要计算凸包，你只需要检查你的点是否可以表示为定义凸包的一部分点的凸组合。此外，在高维空间中找到凸包是计算上非常昂贵的。

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

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))

举个例子，在10维空间中，我解决了10000点的问题。执行时间在毫秒级别。想必不会想知道使用QHull需要多长时间。

- Nils

---

@Juh_: 假设 {x_1,...,x_n} 是 n 个点的集合，{w_1,...,w_n} 是变量权重，y 是你想通过这些 n 个点的组合来描述的点。那么 \sum_i w_i x_i = y_i，然后你想要... - Nils

@Juh_ 这很棘手。我无法在此处编写数学公式。Scipy假设您有以下问题：min_x {c'w | Aw=b, w>=0}，其中w是变量，c是目标系数，而Aw=b是约束条件（LP中默认为w>=0）。由于c为零，不存在实际优化。求解器只需检查可行性，即是否存在w使得Aw=b被满足。现在，在我们的情况下，b = [y_1,...,y_d,1]，而A = [[x_11 w_1,...,x_n1 w_n],...,[x_1d w_1,...,x_nd w_n],[w_1,...,w_n]]。在上面的代码中，查询点y被称为x，而点集x被称为“points”。 - Nils

@Juh_：“为什么需要添加“缩放”维度（即1）？”这是为了实现凸组合的要求，否则您将检查该点是否位于锥体中，这不是您想要的。 - Nils

@S.A Z 是点云。如果有的话，描述 x 的 Z 中点的权重存储在 lp 的属性中。 - Nils

好的解决方案。但是我发现当我检查 10,000 个点是否属于凸包时，速度非常慢。 - user9562553

Answer 3

你好，我不确定如何使用你的程序库来实现这个功能。但是下面用文字描述了一个简单的算法：

1. 创建一个明确在凸包外的点，将其称为Y。
2. 生成一条连接你要判断的点（X）和新点Y的线段。
3. 循环遍历所有凸包的边缘线段，并检查每个边缘线段是否与XY相交。
4. 如果你计算得到的相交数是偶数（包括0），则X在凸包外；否则X在凸包内。
5. 如果发现XY通过你凸包上的一个顶点，或者直接重叠于你凸包的边界，则将Y稍微移动一下。
6. 以上方法也适用于凹多边形。可以参考下面的插图（绿色点是你要判断的X点。黄色标记表示相交点）。

- firemana

---

+1 很好的方法。对于凸包而言，可能更容易找到一个肯定在凸包内部的点（所有凸包顶点的平均值），然后按照相反的成功条件使用您的方法。 - Jaime

虽然这有点吹毛求疵，但有几种情况会导致算法失败：1）如果你选择的一点与凸壳上一对顶点共线，且测试点也与这些顶点共线，则你会得到无限数量交点。2）如果你的测试点和X、外部点Y与一个奇数个面（三维情况）的交点上的顶点共线，那么你将错误地认为测试点实际上在凸壳内...至少，你需要检查第二种情况。例如，确保XYV不共线。 - wmsmith

还有，注意示例中的某些多边形不是凸包，对于凸包，您最多会找到两个交点。我也不确定如何选择“明显在外”的点。也许更容易找到一个“明显在内”的点（例如重心），然后查看其是否具有一个或零个交点，这也可以消除共线性问题（我假设壳是一个凸多边形）。 - Vincenzooo

这需要先找到凸包（作为多边形）。但是，正如Nils的解决方案所示，这一步骤对于整个任务并不是必要的。 - Keenan Pepper

@Vincenzooo 如果你找到了最小点（按字典顺序），然后在所有维度上减去一些量，那么你肯定在外部。此外，有时你可能对点的范围有额外的了解，这使得任务变得微不足道。 - Andreas Vinter-Hviid

Answer 4

首先，获取您点云的凸包。

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

这个Stack Overflow主题还提供了一种解决方法，可以找到点在线的哪一侧： 确定点位于线的哪一侧

---

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

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

看起来您已经有一种方法可以获取您点云的凸包。但是，如果您发现必须自己实现凸包，维基百科在这里列出了一个很好的凸包算法列表： 凸包算法

- Sildoreth

---

如果有人已经计算出了点的凸包，那么这种方法是最简单的。 - CODError

Answer 5

使用ConvexHull的equations属性：

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

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

演示：

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)

- Charlie Brummitt

---

你能解释一下为什么“一个点在凸壳内，当且仅当对于每个方程（描述面）点积与法向量（eq[:-1]）加上偏移量（eq[-1]）小于或等于零”吗？我不太清楚。单个方程中的点积有什么物理意义？我猜它意味着“面的法线指向该点”，但我不明白为什么是这样。 - Gulzar

这个语句是从定义凸包的一种方式中得出的。根据Qhull文档（由scipy使用的代码）：“点集P的凸包是包含P的最小凸集。如果P是有限的，则凸包定义了一个矩阵A和一个向量b，使得对于P中的所有x，Ax+b <= [0,...]”_A_的行是单位法线；_b_的元素是偏移量。 - Charlie Brummitt

这是一个不错的解决方案。但对于10,000个二维点的凸包成员测试来说，速度有点慢。 - user9562553

我已经分享了函数的向量化版本。 - Philippe Miron

Answer 6

为了完整起见，这里提供一个简单的解决方案：

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()

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

下面是两个二维示例图：

错误： 正确：

- feqwix

---

我很喜欢它！但是我要说一下：维度灾难。超过8个维度，核函数就会分裂。 - Ulf Aslak

Answer 7

看起来您正在使用2D点云，因此我想引导您前往凸多边形点内测试 链接。

Scipy的凸包算法可以用于在2个或更多维度中找到凸包，这比2D点云需要更复杂。因此，我建议使用其他算法，例如 此算法。这是因为您真正需要的是按顺时针顺序列出的凸包点列表和一个在多边形内部的点，以进行凸多边形内点测试。

这种方法的时间性能如下：

- O（N log N）用于构建凸包 - 预处理中O（h）计算（和存储）来自内部点的楔角 - 每个点对多边形查询为O（log h）。

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

- Nuclearman

Answer 8

在@Charlie Brummitt的工作基础上，我实现了一个更高效的版本，使其能够同时检查多个点是否在凸包内，并用更快的线性代数替换任何循环。

import numpy as np
from scipy.spatial.qhull import _Qhull

def in_hull(points, queries):
    hull = _Qhull(b"i", points,
                  options=b"",
                  furthest_site=False,
                  incremental=False, 
                  interior_point=None)
    equations = hull.get_simplex_facet_array()[2].T
    return np.all(queries @ equations[:-1] < - equations[-1], axis=1)

# ============== Demonstration ================

points = np.random.rand(8, 2)
queries = np.random.rand(3, 2)
print(in_hull(points, queries))

请注意，我正在使用更低级别的_Qhull类以提高效率。

- milembar

---

我遇到了这个错误：无法从'scipy.spatial.qhull'导入名称'_Qhull'。 - zxdawn

是的，这是因为你的 scipy 版本非常新。 _Qhull 不再被公开。所以你别无选择，只能使用旧版本的 scipy 或者稍微慢一些的 python 接口。 - milembar

Answer 9

参考此答案，若想同时检查numpy数组中的所有点，则以下方式适用于我：

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))

# get array of boolean values indicating in hull if True
in_hull = np.all(np.add(np.dot(random_points, hull.equations[:,:-1].T),
                        hull.equations[:,-1]) <= tolerance, axis=1)

random_points_in_hull = random_points[in_hull]

- Viyaleta

Answer 10

对于那些感兴趣的人，我制作了一个向量化版本的 @charlie-brummit answer：

def points_in_hull(p, hull, tol=1e-12):
    return np.all(hull.equations[:,:-1] @ p.T + np.repeat(hull.equations[:,-1][None,:], len(p), axis=0).T <= tol, 0)

其中p现在是一个[N,2]数组。它比推荐的解决方案(@Sildoreth answer)快约6倍，比原始方法快约10倍。

没有for循环的改编演示:(从下面粘贴以避免在线程中搜索)

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')

in_hull = points_in_hull(random_points, hull)
plt.scatter(random_points[in_hull,0], random_points[in_hull,1], marker='x', color='g')
plt.scatter(random_points[~in_hull,0], random_points[~in_hull,1], marker='d', color='m')

- Philippe Miron