hull = ConvexHull(np.array([(1, 2), (3, 4), (3, 6)]))
for s in hull.simplices:
s.in_simplex(np.array([2, 3]))
引发LinAlgError：数组必须是正方形 .
9 回答
16
这是一个简单的解决方案，只需要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
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')
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 落入"outside"，那么顶点必须改变 . 这适用于n维，但您必须计算两次 ConvexHull .
2-D中的两个示例图：
假：
真正：
20
使用 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)
9 回答
这是一个简单的解决方案，只需要scipy：
它返回一个布尔数组，其中
True值表示位于给定凸包中的点 . 它可以像这样使用：如果安装了matplotlib，还可以使用以下函数调用第一个函数并绘制结果 . 仅适用于2D数据，由
Nx2数组给出：嗨，我不知道如何使用您的程序库来实现这一目标 . 但是有一个简单的算法来实现这一点：
创造一个绝对在船体之外的点 . 叫它Y.
生成一个线段，将您的问题点（X）连接到新点Y.
环绕凸壳的所有边缘段 . 如果片段与XY相交，则检查每个片段 .
如果您计算的交叉点数是偶数（包括0），则X在船体外 . 否则X在船体内 .
如果发生这种情况，XY会穿过船体上的一个顶点，或者直接与你的一个船体边缘重叠，稍微移动一下 .
以上的工作也适用于凹壳 . 您可以在下图中看到（绿点是您要确定的X点 . 黄色标记交叉点.
首先，获得点 Cloud 的凸包 .
然后以逆时针顺序环绕凸包的所有边缘 . 对于每条边，检查您的目标点是否位于该边的“左侧” . 执行此操作时，将边缘视为逆向指向凸包的向量 . 如果目标点是所有向量的“左”，则它由多边形包含;否则，它位于多边形之外 .
此其他Stack Overflow主题包括一个解决方案，用于查找某个点所在的"side"：Determine Which Side of a Line a Point Lies
这种方法的运行时复杂性（一旦你已经有凸包）就是 O(n) ，其中n是凸包具有的边数 .
请注意，这仅适用于凸多边形 . 但是你正在处理一个凸壳，所以它应该适合你的需要 .
看起来你已经有办法为你的点 Cloud 获得凸包 . 但是如果你发现你必须实现自己的，那么维基百科有一个很好的凸包算法列表：Convex Hull Algorithms
我不会使用凸包算法，因为您不需要计算凸包，您只需要检查您的点是否可以表示为子集定义凸包的点集的凸组合 . 此外，找到凸包在计算上是昂贵的，尤其是在更高的尺寸 .
事实上，找出一个点是否可以表示为另一组点的凸组合的单纯问题可以表述为线性规划问题 .
例如，我解决了10维中10000点的问题 . 执行时间在ms范围内 . 不想知道QHull需要多长时间 .
仅仅为了完整性，这是一个穷人的解决方案：
这个想法很简单：如果你添加一个落在船体上的点
p，那么一组点P的凸包顶点将不会改变;凸包的顶点[P1, P2, ..., Pn]和[P1, P2, ..., Pn, p]是相同的 . 但如果p落入"outside"，那么顶点必须改变 . 这适用于n维，但您必须计算两次ConvexHull.2-D中的两个示例图：
假：
真正：
使用
ConvexHull的equations属性：换言之，当且仅当对于每个方程（描述小平面）时，点与正常向量（
eq[:-1]）之间的点积加上偏移量（eq[-1]）小于或等于零时，点在船体中 . 你可能想要比较一个小的正常数tolerance = 1e-12而不是零因为数值精度的问题（否则，你可能会发现凸包的顶点不在凸面上船体） .示范：
看起来你正在使用2D点 Cloud ，所以我想引导你到inclusion test进行凸多边形的多边形点测试 .
Scipy的凸包算法允许在2维或更多维中找到凸包，这比2D点 Cloud 需要更复杂 . 因此，我建议使用不同的算法，例如this one . 这是因为您真正需要凸包的多边形点测试是顺时针顺序的凸包点列表，以及多边形内部的点 .
这种方法的时间表现如下：
O（N log N）构造凸包
O（h）在预处理中从内点计算（并存储）楔角
每个多边形点查询
O（log h） .
其中N是点 Cloud 中的点数，h是点 Cloud 凸壳中的点数 .
如果你想保持scipy，你必须凸壳（你这样做）
然后在船体上 Build 点列表 . 这是来自doc的用于绘制船体的代码
所以从那开始，我建议计算船体上的点列表
（虽然我没试过）
您还可以使用自己的代码来计算船体，返回x，y点 .
If you want to know if a point from your original dataset is on the hull ，那你就完成了 .
I what you want is to know if a any point is inside the hull or outside ，你必须做更多的工作 . 你需要做的就是
用于连接船体两个单体的所有边缘：确定您的点是高于还是低于
如果点低于所有线，或高于所有线，则它位于船体外
作为加速，一旦一个点超过一条线并且低于一条线，它就在船体内部 .
基于this帖子，这是我的快速而肮脏的解决方案，适用于有4个边的凸面区域（您可以轻松地将其扩展到更多）