我在 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：数组必须是正方形。
【问题讨论】：
标签： python numpy convex-hull
这是一个只需要 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，您还可以使用以下函数调用第一个函数并绘制结果。仅适用于二维数据，由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')
【讨论】：
cloud 是一个 NxK 数组，由 K 维的 N 个点组成，ConvexHull(cloud).vertices（来自 scipy.spatial）给出了凸包上点的索引，即“外点”
Delaunay.find_simplex 的文档中所解释的那样，它为船体外部的点返回-1。现在，如果你想要更多的控制，或者想要更快的算法，我推荐下面@nils 的解决方案。它更复杂，但只计算需要的东西（我没有测试它，但看起来确实如此）
hull = Delaunay(hull).有什么想法吗？
我不会使用凸包算法，因为您不需要计算凸包，您只想检查您的点是否可以表示为一组点的凸组合，其中一个子集定义了一个凸包船体。此外，找到凸包的计算成本很高，尤其是在更高维度上。
事实上，仅仅找出一个点是否可以表示为另一组点的凸组合的问题就可以表述为线性规划问题。
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 需要多长时间。
【讨论】：
您好，我不确定如何使用您的程序库来实现这一点。但是有一个简单的算法可以实现这一点，用文字描述：
【讨论】：
首先，获取点云的凸包。
然后以逆时针顺序循环遍历凸包的所有边缘。对于每条边，检查您的目标点是否位于该边的“左侧”。执行此操作时，将边缘视为围绕凸包逆时针指向的向量。如果目标点在所有向量的“左侧”，则它包含在多边形中；否则，它位于多边形之外。
另一个 Stack Overflow 主题包括一个解决方案来查找点在哪条线的“边”上： Determine Which Side of a Line a Point Lies
请注意，这仅适用于凸多边形。但是你正在处理一个凸包，所以它应该适合你的需要。
看起来您已经有办法为您的点云获取凸包。但是如果你发现你必须自己实现，维基百科在这里有一个很好的凸包算法列表： Convex Hull Algorithms
【讨论】：
使用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)
【讨论】：
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吗？这对我来说不是很清楚。对于单个方程，该点积的物理意义是什么？我猜这意味着“刻面的正常点”，但我不明白为什么会这样
为了完整起见，这里是一个穷人的解决方案：
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 两次。
二维中的两个示例图：
错误：
正确：
【讨论】：
看起来您正在使用 2D 点云，所以我想引导您到 inclusion test 进行凸多边形的多边形内点测试。
Scipy 的凸包算法允许找到 2 维或更多维的凸包，这比 2D 点云所需的要复杂。因此，我建议使用不同的算法，例如this one。这是因为您真正需要的凸包的多边形点测试是按顺时针顺序排列的凸包点列表，以及多边形内部的一个点。
这种方法的时间表现如下：
其中 N 是点云中的点数，h 是点云凸包中的点数。
【讨论】：
背负this answer，一次检查一个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]
【讨论】：
在@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 类。
【讨论】：
如果你想继续使用 scipy，你必须使用凸包（你这样做了）
>>> from scipy.spatial import ConvexHull
>>> points = np.random.rand(30, 2) # 30 random points in 2-D
>>> hull = ConvexHull(points)
然后在船体上建立点列表。这是 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-')
所以从那开始，我建议计算船体上的点列表
pts_hull = [(points[simplex,0], points[simplex,1])
for simplex in hull.simplices]
（虽然我没试过）
您还可以使用自己的代码来计算船体，返回 x,y 点。
如果您想知道原始数据集中的某个点是否在船体上，那么您就完成了。
我想知道一个点是在船体内部还是外部，你必须多做一些工作。你需要做的可能是
对于连接你的船体的两个单纯形的所有边：决定你的点是在上面还是在下面
如果点低于所有线，或高于所有线，则它在船体之外
作为一种加速，只要一个点在一条线上方和另一条线下方，它就在船体内部。
【讨论】：
基于this 帖子，这是我针对具有 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.')
【讨论】：
d = np.roll(pts, 2) - pts