我在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：数组必须是正方形。
问问题
53117 次
11 回答
101
这是一个只需要 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')
于 2013-06-03T14:03:15.297 回答
41
我不会使用凸包算法，因为您不需要计算凸包，您只想检查您的点是否可以表示为一组点的凸组合，其中一个子集定义了一个凸包。此外，找到凸包的计算成本很高，尤其是在更高维度上。
事实上，仅仅找出一个点是否可以表示为另一组点的凸组合的问题就可以表述为线性规划问题。
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 需要多长时间。
于 2017-04-22T21:21:22.653 回答
27
嗨，我不确定如何使用您的程序库来实现这一点。但是有一个简单的算法可以实现这一点，用文字描述：
- 创建一个绝对在您的船体之外的点。叫它Y
- 生成一条线段，将您的问题点 (X) 连接到新点 Y。
- 环绕凸包的所有边缘段。如果线段与 XY 相交，请检查它们中的每一个。
- 如果你统计的交叉点数是偶数（包括 0），那么 X 在船体外。否则 X 在船体内部。
- 如果发生这种情况，XY 穿过船体上的一个顶点，或者直接与船体的一个边缘重叠，请稍微移动 Y。
- 以上也适用于凹形船体。您可以在下图中看到（绿点是您要确定的 X 点。黄色标记交叉点。
于 2013-06-04T04:42:48.033 回答
24
首先，获取点云的凸包。
然后以逆时针顺序循环遍历凸包的所有边缘。对于每条边，检查您的目标点是否位于该边的“左侧”。执行此操作时，将边缘视为围绕凸包逆时针指向的向量。如果目标点在所有向量的“左侧”，则它包含在多边形中；否则，它位于多边形之外。
另一个 Stack Overflow 主题包括一个解决方案来查找点位于直线的 哪一侧：确定点位于直线的哪一侧
这种方法的运行时复杂度（一旦你已经有了凸包）是O(n)，其中 n 是凸包具有的边数。
请注意，这仅适用于凸多边形。但是你正在处理一个凸包，所以它应该适合你的需要。
看起来您已经有办法为您的点云获取凸包。但是如果你发现你必须自己实现，维基百科在这里有一个很好的凸包算法列表： Convex Hull Algorithms
于 2013-06-03T21:39:23.917 回答
15
使用 的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)
换句话说，当且仅当对于每个方程（描述小平面），该点与法向量 ( 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)
于 2017-02-10T17:46:09.563 回答
7
只是为了完整性，这里是一个穷人的解决方案：
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的凸包的顶点是相同的。但是如果落在“外面”，那么顶点必须改变。这适用于 n 维，但您必须计算两次。[P1, P2, ..., Pn][P1, P2, ..., Pn, p]pConvexHull
二维中的两个示例图：
错误的：
真的：
于 2014-04-08T08:34:43.490 回答
2
在@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类来提高效率。
于 2021-05-09T12:40:52.490 回答
2
为了背负这个答案，一次检查一个 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]
于 2021-02-10T14:22:17.687 回答
1
如果你想和 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 点。
如果您想知道原始数据集中的一个点是否在 hull 上，那么您就完成了。
我想知道任何一点是在船体内部还是外部，你必须多做一些工作。你需要做的可能是
对于连接你的船体的两个单纯形的所有边：决定你的点是高于还是低于
如果点低于所有线或高于所有线，则它位于船体外部
作为一种加速，一旦一个点在一条线之上并在另一条线之下，它就在船体内部。
于 2013-05-25T15:40:23.507 回答
0
基于这篇文章，这是我针对具有 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.')
于 2017-03-21T12:51:47.777 回答