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

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:数组必须是正方形.

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

Answer 2

嗨,我不知道如何使用您的程序库来实现这一目标.但是有一个简单的算法来实现这一点:

1. 创造一个绝对在船体外面的点.叫它Y.
2. 生成一个线段,将您的问题点(X)连接到新点Y.
3. 环绕凸壳的所有边缘段.如果片段与XY相交,则检查每个片段.
4. 如果您计算的交叉点数是偶数(包括0),则X在船体外.否则X在船体内.
5. 如果发生这种情况,XY会穿过船体上的一个顶点,或者直接与船体的一个边缘重叠,稍微移动一下.
6. 以上工作也适用于凹形船体.您可以在下图中看到(绿点是您要确定的X点.黄色标记交叉点.

Answer 3

首先,获得点云的凸包.

然后以逆时针顺序环绕凸包的所有边缘.对于每条边,检查您的目标点是否位于该边的"左侧".执行此操作时,将边缘视为逆向指向凸包的向量.如果目标点是所有向量的"左",则它由多边形包含; 否则,它位于多边形之外.

这个其他Stack Overflow主题包括一个解决方案,用于查找点所在线的"侧面":确定点的 哪一侧是一个点

这种方法的运行时复杂性(一旦你已经有凸包)就是O(n),其中n是凸包的边数.

请注意,这仅适用于凸多边形.但是你正在处理一个凸壳,所以它应该适合你的需要.

看起来你已经有办法为你的点云获得凸包.但是如果你发现你必须实现自己的,那么维基百科有一个很好的凸包算法列表: 凸壳算法

Answer 4

我不会使用凸包算法,因为您不需要计算凸包,您只需要检查您的点是否可以表示为子集定义凸包的点集的凸组合.此外,找到凸包在计算上是昂贵的,尤其是在更高的尺寸.

事实上,找出一个点是否可以表示为另一组点的凸组合的单纯问题可以表述为线性规划问题.

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点的问题.执行时间在ms范围内.不想知道QHull需要多长时间.

Answer 5

仅仅为了完整性,这是一个穷人的解决方案:

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两次.

2-D中的两个示例图:

假:

真正:

Answer 6

看起来您正在使用 2D 点云，因此我想引导您进行凸多边形的多边形内点测试的包含测试。

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

该方法的时间性能如下：

• O(N log N) 构造凸包
• 预处理中计算（并存储）与内点的楔角的 O(h)
• 每个多边形内点查询的时间复杂度为 O(log h)。

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

Answer 7

使用以下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)

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类是为了提高效率。

Answer 9

如果你想与 scipy 保持一致，你必须使用凸包（你这样做了）

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

然后在船体上建立点列表。这是来自文档的绘制船体的代码

>>> 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 点。

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

我你想要的是知道某个点是在船体内部还是外部，你必须做更多的工作。你需要做的可能是

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

• 如果点位于所有线下方或所有线上方，则它位于船体外部

作为加速，一旦一个点位于一条线上方且另一条线下方，它就位于船体内部。