• 当前位置：首页
  • IT科技
  • 有什么有效的方法可以确定某个点是否位于点云的凸包中？

有什么有效的方法可以确定某个点是否位于点云的凸包中？

2025-02-27 09:07:00
admin
原创
25
摘要：问题描述：我在 numpy 中有一个坐标点云。对于大量的点，我想找出这些点是否位于点云的凸包中。我尝试了 pyhull，但我不知道如何检查某个点是否在ConvexHull：hull = ConvexHull(np.array([(1, 2), (3, 4), (3, 6)])) for s in hull.s...

问题描述：

我在 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：数组必须是正方形。

解决方案 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，您还可以使用以下函数调用第一个函数并绘制结果。仅适用于二维数据，由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')

解决方案 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 需要多长时间。

解决方案 3：

你好，我不确定如何使用你的程序库来实现这一点。但有一个简单的算法可以用文字描述来实现这一点：

  1. 创建一个绝对位于船体外部的点。称之为 Y

  2. 生成一条连接你所讨论的点（X）和新点 Y 的线段。

  3. 环绕凸包的所有边缘段。检查每个边段是否与 XY 相交。

  4. 如果计算出的交点数为偶数（包括 0），则 X 位于船体外部。否则，X 位于船体内部。

  5. 如果发生这种情况，XY 会穿过船体上的一个顶点，或直接与船体的某个边缘重叠，请稍微移动 Y。

  6. 上述方法也适用于凹壳。您可以在下图中看到（绿点是您要确定的 X 点。黄色标记交叉点。
    插图

解决方案 4：

首先，获取点云的凸包。

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

循环并检查点是否始终在“左边”

这个 Stack Overflow 上的另一个主题包含一个解决方案，用于查找某个点位于一条线的哪“侧”：
确定某个点位于一条线的哪一侧

这种方法的运行时间复杂度（一旦您已经有凸包）为O（n），其中 n 是凸包具有的边数。
请注意，这仅适用于凸多边形。但您正在处理凸包，因此它应该适合您的需求。

看起来你已经有办法获取点云的凸包了。但是如果你发现你必须自己实现，维基百科上有一个很好的凸包算法列表：
凸包算法

解决方案 5：

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

演示输出

解决方案 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两次。

两个二维示例图：

错误的：

新点（红色）落在凸包之外

真的：

新点（红色）落入凸包内

解决方案 7：

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

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

该方法的时间性能如下：

  • 构建凸包的时间为 O(N log N)

  • 预处理中的 O(h) 用于计算（并存储）内部点的楔角

  • 每个多边形内点的查询耗时 O(log h)。

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

解决方案 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请注意，为了提高效率，我使用了较低级别的类。

解决方案 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]

解决方案 10：

对于那些感兴趣的人，我制作了@charlie-brummit答案的矢量化版本：

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答案) 快约 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')

解决方案 11：

如果你想继续使用 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 点。

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

如果你想知道任何一点是在船体内部还是外部，你必须做更多的工作。你需要做的可能是

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

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

为了加速，一旦一个点位于一条线之上和另一条线之下，它就在船体内部。

解决方案 12：

矢量化，带有绘图功能（基于此处的其他精彩答案：1、2）：

import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import ConvexHull


def points_in_hull(points: np.ndarray, hull: ConvexHull, tolerance: float = 1e-12):
    return np.all(np.add(points @ hull.equations[:, :-1].T, hull.equations[:, -1]) <= tolerance, axis=1)


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, (1000, 2))

in_hull = points_in_hull(random_points, hull)
random_points_in_hull = random_points[in_hull]

# plot
fig, ax = plt.subplots(figsize=(20, 10))
plt.scatter(random_points[:, 0], random_points[:, 1], color='blue')
plt.scatter(points[:, 0], points[:, 1], color='green')
plt.scatter(random_points_in_hull[:, 0], random_points_in_hull[:, 1], color='red')
for simplex in hull.simplices:
    plt.plot(points[simplex, 0], points[simplex, 1], 'k-')
plt.show()

解决方案 13：

根据这篇文章，这是我针对具有 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.')

在此处输入图片描述

解决方案 14：

scipy.spatial对于那些可能感兴趣的人，我编写了自己的实现，它完全不依赖于（尽管如果需要，它用于计算切比雪夫中心），它既更快，又提供了额外的功能。然而，它只支持二维空间。该模块在Githubscipy.optimize上公开可用。

相关推荐
置顶 2024年20款好用的项目管理软件推荐，项目管理提效的20个工具和技巧
2024年20款好用的项目管理软件推荐，项目管理提效的20个工具和技巧
  为什么项目管理通常仍然耗时且低效？您是否还在反复更新电子表格、淹没在便利贴中并参加每周更新会议？这确实是耗费时间和精力。借助软件工具的帮助，您可以一目了然地全面了解您的项目。如今，国内外有足够多优秀的项目管理软件可以帮助您掌控每个项目。什么是项目管理软件？项目管理软件是广泛行业用于项目规划、资源分配和调度的软件。它使项...
项目管理软件   1289  
IPD流程中TR的创新思维与应用实践
IPD流程中TR的创新思维与应用实践
  IPD（Integrated Product Development）流程作为一种先进的产品开发管理模式，在众多企业中得到了广泛应用。其中，技术评审（TR，Technical Review）环节至关重要，它不仅是对技术方案的评估，更是激发创新思维、推动产品创新的关键节点。深入理解TR在IPD流程中的创新思维及其应用实践...
IPD流程中TR   9  
IPD产品开发流程如何提升市场竞争力
IPD产品开发流程如何提升市场竞争力
  IPD（Integrated Product Development）产品开发流程作为一种先进的产品开发管理模式，在众多企业中得到了广泛应用。它打破了传统产品开发过程中部门之间的壁垒，将市场、研发、生产、销售等各个环节有机整合在一起，形成一个高效协同的整体。通过这种方式，企业能够更快速、更精准地开发出满足市场需求的产品...
IPD管理流程   11  
IPD流程中TR与项目成功的关系探讨
IPD流程中TR与项目成功的关系探讨
  IPD（Integrated Product Development）流程作为一种先进的产品开发管理模式，在众多企业中得到了广泛应用。其中，技术评审（TR，Technical Review）环节在整个IPD流程里占据着关键位置，对项目的成功有着深远影响。深入探讨TR与项目成功的关系，有助于企业更好地运用IPD流程，提升...
IPD项目管理   8  
IPD研发管理体系下的产品创新三步骤
IPD研发管理体系下的产品创新三步骤
  IPD研发管理体系旨在打破部门墙，实现跨部门协同，确保产品开发以市场和客户需求为导向，高效、高质量地推出满足市场需求的产品。在这一体系下，产品创新可拆解为三个关键步骤，它们环环相扣，共同推动企业的产品不断迭代升级，在激烈的市场竞争中占据优势。这三个步骤分别聚焦于洞察市场机会、规划产品战略以及执行开发与验证，每一步都蕴含...
IPD框架   10  
热门文章
2024年项目管理软件排行榜
2024年20款好用的项目管理软件推荐，项目管理提效的20个工具和技巧
img020
使用 Python 'Requests' 模块的代理
img044
10个项目管理常见问题解答：包括难点和解决方法
项目管理软件
项目延期的原因有哪些？如何有效预防和管理？
项目里程碑计划模板
项目里程碑计划模板怎么写?设定项目里程碑的5个方法
项目管理软件
项目管理指南：制定项目实施计划的6大步骤有哪些？
img081
工程项目管理中的团队建设与激励策略
项目管理软件有哪些？
热门标签
更多>
IPD流程 项目管理 IT项目管理 敏捷开发 项目管理软件 Scrum 敏捷项目管理 IPD流程管理 项目风险管理 项目管理工具 项目经理 极限编程  产品管理  项目经理职责 项目路线图 项目进度管理 测试管理工具 项目预算管理 PMP认证 瀑布式管理 项目范围管理 IT科技 项目管理平台 工程项目管理 项目资源管理 pmp证书考试
云禅道AD
禅道项目管理软件

云端的项目管理软件

尊享禅道项目软件收费版功能

无需维护，随时随地协同办公

内置subversion和git源码管理

每天备份，随时转为私有部署

免费试用