【问题标题】:Checking if a point is in ConvexHull?检查一个点是否在 ConvexHull 中?
【发布时间】:2019-01-17 04:02:53
【问题描述】:

我无法理解如何计算一个 n 维点是否在一个 n 维 ConvexHull 内。

这里问了一个非常相似的问题(相同): What's an efficient way to find if a point lies in the convex hull of a point cloud?

但是,答案让我感到困惑或似乎对我不起作用,我不知道为什么。

def in_hull(p, hull):
    """ Copied and from the Top Original answer """
    from scipy.spatial import Delaunay
    if not isinstance(hull,Delaunay):
        hull = Delaunay(hull)

    return hull.find_simplex(p)>=0

这个函数给了我很多错误或不想要的结果,我正在使用它。但是,在调试时,我编写了一个简单的脚本来测试一些明显的期望:

如果我用一组点构造一个 ConvexHull, 当我检查那组积分的“会员资格”时,它们都应该是 “成员”。

results_all = []
for _ in range(5000):
    cloud = np.random.rand(5000, 2)
    result = in_hull(cloud, cloud)
    results_all.append(np.all(result))

arr = np.array(results_all)
print(np.sum(np.logical_not(arr)))

虽然这种情况很少见,但在随机生成的数据(5000 个中有 3 个)上似乎会失败,但在实际数据上问题更大。我所说的失败是指我实际上遇到了一些情况,并不是所有的点都被视为成员。

我做错了什么吗?或者也许完全是误解?在这一点上我很困惑,所以很想解释一下发生了什么。

最后,我想要;给定一个在前一阶段计算的 ConvexHull;能够确定点是否位于船体内。

【问题讨论】:

    标签: python numpy scipy


    【解决方案1】:

    对于几乎平坦的单纯形(三角形),Delaunay 对象的 find_simplex 方法似乎是一个边缘情况问题。

    下面是一个代码,用于查找和绘制只有 3 个点的错误案例:

    import matplotlib.pylab as plt
    from scipy.spatial import Delaunay
    from scipy.spatial import delaunay_plot_2d
    
    for _ in range(5000):
        cloud = np.random.rand(3, 2)
    
        tri = Delaunay(cloud)
    
        if np.any( tri.find_simplex(cloud)<0 ):
            print('break at', _)
    
            delaunay_plot_2d(tri);
            id_break = np.where(tri.find_simplex(cloud)<0)
            plt.plot( *cloud[id_break].ravel(), 'or' );
            break
    

    here 提出的另一种方法似乎效果很好:

    hull = ConvexHull(cloud)
    
    def point_in_hull(point, hull, tolerance=1e-12):
        return all(
            (np.dot(eq[:-1], point) + eq[-1] <= tolerance)
            for eq in hull.equations)
    
    [ point_in_hull(point, hull) for point in cloud ]
    # [True, True, True]
    

    【讨论】:

    • 非常感谢!这解决了我的问题,是的,这正是边境案件。
    猜你喜欢
    • 2015-10-07
    • 2014-12-12
    • 1970-01-01
    • 2011-05-11
    • 2013-07-20
    • 2011-04-18
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多