当前位置： 文江博客话题详情

查找一个点是否位于一组点的凸包内，而不计算包本身

发布于 2024-10-16 02:25:11 字数 99 浏览 7 评论 0原文

测试点 P 是否位于由一组点 X 形成的凸包内的最简单方法是什么？

我想要一种在高维空间（例如，最多 40 维）中工作的算法，该算法不会显式计算凸包本身。有什么想法吗？

原文

What is the simplest way to test if a point P is inside a convex hull formed by a set of points X?

I'd like an algorithm that works in a high-dimensional space (say, up to 40 dimensions) that doesn't explicitly compute the convex hull itself. Any ideas?

如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论， 你可以免费 注册 一个本站的账号。

评论（7）

凉薄对峙 2024-10-23 02:25:11

该问题可以通过找到线性规划的可行点来解决。如果您对完整细节感兴趣，而不是仅仅将 LP 插入现有求解器，我建议您阅读 Boyd 和 Vandenberghe 关于凸优化的优秀书籍。

设A = (X[1] X[2] ... X[n])，即第一列为v1，第二列为v2 等。

解决以下 LP 问题，

minimize (over x): 1
s.t.     Ax = P
         x^T * [1] = 1
         x[i] >= 0  \forall i

其中

  1. x^T 是 x 的转置
  2. [1] 是 all- 1 个向量。

问题有解当且仅当该点位于凸包内。

The problem can be solved by finding a feasible point of a Linear Program. If you're interested in the full details, as opposed to just plugging an LP into an existing solver, I'd recommend reading Chapter 11.4 in Boyd and Vandenberghe's excellent book on convex optimization.

Set A = (X[1] X[2] ... X[n]), that is, the first column is v1, the second v2, etc.

Solve the following LP problem,

minimize (over x): 1
s.t.     Ax = P
         x^T * [1] = 1
         x[i] >= 0  \forall i

where

  1. x^T is the transpose of x
  2. [1] is the all-1 vector.

The problem has a solution iff the point is in the convex hull.

回复 原文
谁的年少不轻狂 2024-10-23 02:25:11

当且仅当从该点到其他点的所有向量的方向小于其周围圆/球/超球面的一半时，该点位于其他点的凸包之外。

下面是两点情况的草图，蓝色的点在凸包（绿色）内部，红色的点在凸包外部：

点的草图

对于红色的，存在圆的二等分，使得向量从该点到凸包上的点只相交半个圆。
对于蓝点，不可能找到这样的二等分。

The point lies outside of the convex hull of the other points if and only if the direction of all the vectors from it to those other points are on less than one half of a circle/sphere/hypersphere around it.

Here is a sketch for the situation of two points, a blue one inside the convex hull (green) and the red one outside:

Sketch of the points

For the red one, there exist bisections of the circle, such that the vectors from the point to the points on the convex hull intersect only one half of the circle.
For the blue point, it is not possible to find such a bisection.

回复 原文
乖乖 2024-10-23 02:25:11

您不必计算凸包本身，因为它在多维空间中看起来相当麻烦。凸包有

一个众所周知的属性：任意向量（点）v 点 [v1, v2, .., vn] 的凸包内部可以表示为 sum(ki*vi)，其中 0 0 ;= ki <= 1 且 sum(ki) = 1。相应地，凸包之外的任何点都不会具有这样的表示。

在 m 维空间中，这将为我们提供具有 n 个未知数的 m 个线性方程组。

编辑
我不确定一般情况下这个新问题的复杂性，但对于 m = 2 来说，它似乎是线性的。也许，在这方面有更多经验的人会纠正我。

You don't have to compute convex hull itself, as it seems quite troublesome in multidimensional spaces. There's a well-known property of convex hulls:

Any vector (point) v inside convex hull of points [v1, v2, .., vn] can be presented as sum(ki*vi), where 0 <= ki <= 1 and sum(ki) = 1. Correspondingly, no point outside of convex hull will have such representation.

In m-dimensional space, this will give us the set of m linear equations with n unknowns.

edit
I'm not sure about complexity of this new problem in general case, but for m = 2 it seems linear. Perhaps, somebody with more experience in this area will correct me.

回复 原文
囍笑 2024-10-23 02:25:11

我在 16 维上也遇到了同样的问题。由于即使 qhull 也无法正常工作，因为必须生成太多面，所以我通过测试开发了自己的方法，是否可以在新点和参考数据之间找到分离超平面（我称之为“HyperHull”；）） 。

寻找分离超平面的问题可以转化为凸二次规划问题（参见：SVM） 。我使用 cvxopt 在 python 中完成了此操作，代码不到 170 行（包括 I/O）。即使存在问题，该算法在任何维度上都无需修改即可工作，即维度越高，船体上的点数越高（请参阅：关于多面体中随机点的凸包）。由于船体不是明确构造的，而是仅检查点是否在内部，因此与快速船体等相比，该算法在更高维度上具有很大的优势。

该算法可以“自然地”并行化，并且加速应该等于处理器的数量。

I had the same problem with 16 dimensions. Since even qhull didn't work properly as too much faces had to be generated, I developed my own approach by testing, whether a separating hyperplane can be found between the new point and the reference data (I call this "HyperHull" ;) ).

The problem of finding a separating hyperplane can be transformed to a convex quadratic programming problem (see: SVM). I did this in python using cvxopt with less then 170 lines of code (including I/O). The algorithm works without modification in any dimension even if there exists the problem, that as higher the dimension as higher the number of points on the hull (see: On the convex hull of random points in a polytope). Since the hull isn't explicitely constructed but only checked, whether a point is inside or not, the algorithm has very big advantages in higher dimensions compared to e.g. quick hull.

This algorithm can 'naturally' be parallelized and speed up should be equal to number of processors.

回复 原文
家住魔仙堡 2024-10-23 02:25:11

虽然最初的帖子是三年前的，但也许这个答案仍然会有帮助。 Gilbert-Johnson-Keerthi (GJK) 算法找到两个凸多胞体之间的最短距离，每个凸多胞体被定义为一组生成器的凸包——值得注意的是，凸包本身不需要计算。在一种特殊情况下，即所询问的情况，其中一个多面体只是一个点。为什么不尝试使用GJK算法来计算P与点X的凸包之间的距离？如果该距离为 0，则 P 在 X 内部（或至少在其边界上）。 Octave/Matlab 中的 GJK 实现（称为 ClosestPointInConvexPolytopeGJK.m）以及支持代码可在 http://www.99main.com/~centore/MunsellAndKubelkaMunkToolbox/MunsellAndKubelkaMunkToolbox.html 。 GJK 算法的简单描述可以在第 4 节中找到。 2 篇论文，位于 http://www.99main.com/~centore/ColourSciencePapers /GJKinConstrainedLeastSquares.pdf 。我已经对 31 维空间中的一些非常小的集合 X 使用了 GJK 算法，并且得到了很好的结果。 GJK 的性能与其他人推荐的线性规划方法相比如何还不确定（尽管任何比较都会很有趣）。 GJK 方法确实避免了计算凸包或用线性不等式表示包，这两者都可能非常耗时。希望这个答案有帮助。

Though the original post was three years ago, perhaps this answer will still be of help. The Gilbert-Johnson-Keerthi (GJK) algorithm finds the shortest distance between two convex polytopes, each of which is defined as the convex hull of a set of generators---notably, the convex hull itself does not have to be calculated. In a special case, which is the case being asked about, one of the polytopes is just a point. Why not try using the GJK algorithm to calculate the distance between P and the convex hull of the points X? If that distance is 0, then P is inside X (or at least on its boundary). A GJK implementation in Octave/Matlab, called ClosestPointInConvexPolytopeGJK.m, along with supporting code, is available at http://www.99main.com/~centore/MunsellAndKubelkaMunkToolbox/MunsellAndKubelkaMunkToolbox.html . A simple description of the GJK algorithm is available in Sect. 2 of a paper, at http://www.99main.com/~centore/ColourSciencePapers/GJKinConstrainedLeastSquares.pdf . I've used the GJK algorithm for some very small sets X in 31-dimensional space, and had good results. How the performance of GJK compares to the linear programming methods that others are recommending is uncertain (although any comparisons would be interesting). The GJK method does avoid computing the convex hull, or expressing the hull in terms of linear inequalities, both of which might be time-consuming. Hope this answer helps.

回复 原文
神爱温柔 2024-10-23 02:25:11

您是否愿意接受通常应该有效但不能保证有效的启发式答案？如果你是，那么你可以尝试这个随机的想法。

令 f(x) 为到 P 的距离的立方乘以 X 中的事物数量，减去到 X 中所有点的距离的立方之和。从随机的某个位置开始，并使用爬山算法来最大化f(x) 对于距离 P 很远的球体中的 x。除了退化情况外，如果 P 不在凸包中，则应该很有可能找到 P 位于其一侧的超平面的法线，X 中的所有内容都在 的另一边。

Are you willing to accept a heuristic answer that should usually work but is not guaranteed to? If you are then you could try this random idea.

Let f(x) be the cube of the distance to P times the number of things in X, minus the sum of the cubes of the distance to all of the points in X. Start somewhere random, and use a hill climbing algorithm to maximize f(x) for x in a sphere that is very far away from P. Excepting degenerate cases, if P is not in the convex hull this should have a very good probability of finding the normal to a hyperplane which P is on one side of, and everything in X is on the other side of.

回复 原文
笑忘罢 2024-10-23 02:25:11

使用 scipy.optimize.minimize 测试点是否位于外壳空间中的文章。

基于用户1071136的回答。

如果计算凸包，它确实会快得多，所以我为想要这样做的人添加了几行。我从 graham 扫描（仅限 2D）切换到 scipy qhull 算法。

scipy.optimize.minimize 文档：
https://docs.scipy.org/doc/scipy/reference/ Optimize.nonlin.html

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


def hull_test(P, X, use_hull=True, verbose=True, hull_tolerance=1e-5, return_hull=True):
    if use_hull:
        hull = ConvexHull(X)
        X = X[hull.vertices]

    n_points = len(X)

    def F(x, X, P):
        return np.linalg.norm( np.dot( x.T, X ) - P )

    bnds = [[0, None]]*n_points # coefficients for each point must be > 0
    cons = ( {'type': 'eq', 'fun': lambda x: np.sum(x)-1} ) # Sum of coefficients must equal 1
    x0 = np.ones((n_points,1))/n_points # starting coefficients

    result = scipy.optimize.minimize(F, x0, args=(X, P), bounds=bnds, constraints=cons)

    if result.fun < hull_tolerance:
        hull_result = True
    else:
        hull_result = False

    if verbose:
        print( '# boundary points:', n_points)
        print( 'x.T * X - P:', F(result.x,X,P) )
        if hull_result: 
            print( 'Point P is in the hull space of X')
        else: 
            print( 'Point P is NOT in the hull space of X')

    if return_hull:
        return hull_result, X
    else:
        return hull_result

对一些示例数据进行测试：

n_dim = 3
n_points = 20
np.random.seed(0)

P = np.random.random(size=(1,n_dim))
X = np.random.random(size=(n_points,n_dim))

_, X_hull = hull_test(P, X, use_hull=True, hull_tolerance=1e-5, return_hull=True)

输出：

# boundary points: 14
x.T * X - P: 2.13984259782e-06
Point P is in the hull space of X

可视化它：

rows = max(1,n_dim-1)
cols = rows
plt.figure(figsize=(rows*3,cols*3))
for row in range(rows):
    for col in range(row, cols):
        col += 1
        plt.subplot(cols,rows,row*rows+col)
        plt.scatter(P[:,row],P[:,col],label='P',s=300)
        plt.scatter(X[:,row],X[:,col],label='X',alpha=0.5)
        plt.scatter(X_hull[:,row],X_hull[:,col],label='X_hull')
        plt.xlabel('x{}'.format(row))
        plt.ylabel('x{}'.format(col))
plt.tight_layout()

船体测试可视化

A write-up to test if a point is in a hull space, using scipy.optimize.minimize.

Based on user1071136's answer.

It does go a lot faster if you compute the convex hull, so I added a couple of lines for people who want to do that. I switched from graham scan (2D only) to the scipy qhull algorithm.

scipy.optimize.minimize documentation:
https://docs.scipy.org/doc/scipy/reference/optimize.nonlin.html

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


def hull_test(P, X, use_hull=True, verbose=True, hull_tolerance=1e-5, return_hull=True):
    if use_hull:
        hull = ConvexHull(X)
        X = X[hull.vertices]

    n_points = len(X)

    def F(x, X, P):
        return np.linalg.norm( np.dot( x.T, X ) - P )

    bnds = [[0, None]]*n_points # coefficients for each point must be > 0
    cons = ( {'type': 'eq', 'fun': lambda x: np.sum(x)-1} ) # Sum of coefficients must equal 1
    x0 = np.ones((n_points,1))/n_points # starting coefficients

    result = scipy.optimize.minimize(F, x0, args=(X, P), bounds=bnds, constraints=cons)

    if result.fun < hull_tolerance:
        hull_result = True
    else:
        hull_result = False

    if verbose:
        print( '# boundary points:', n_points)
        print( 'x.T * X - P:', F(result.x,X,P) )
        if hull_result: 
            print( 'Point P is in the hull space of X')
        else: 
            print( 'Point P is NOT in the hull space of X')

    if return_hull:
        return hull_result, X
    else:
        return hull_result

Test on some sample data:

n_dim = 3
n_points = 20
np.random.seed(0)

P = np.random.random(size=(1,n_dim))
X = np.random.random(size=(n_points,n_dim))

_, X_hull = hull_test(P, X, use_hull=True, hull_tolerance=1e-5, return_hull=True)

Output:

# boundary points: 14
x.T * X - P: 2.13984259782e-06
Point P is in the hull space of X

Visualize it:

rows = max(1,n_dim-1)
cols = rows
plt.figure(figsize=(rows*3,cols*3))
for row in range(rows):
    for col in range(row, cols):
        col += 1
        plt.subplot(cols,rows,row*rows+col)
        plt.scatter(P[:,row],P[:,col],label='P',s=300)
        plt.scatter(X[:,row],X[:,col],label='X',alpha=0.5)
        plt.scatter(X_hull[:,row],X_hull[:,col],label='X_hull')
        plt.xlabel('x{}'.format(row))
        plt.ylabel('x{}'.format(col))
plt.tight_layout()

Visualization of hull test

回复 原文
~没有更多了~

关于作者

鱼窥荷

暂无简介

0 文章
0 评论
22 人气
发私信

相关话题

更多

热门标签

操作系统 程序设计 IT运维 Linux系统管理 JavaScript 服务器应用 solaris C/C++ PHP Shell BSD Vue.js aix Oracle Python HTML 系统管理 HTML5 CSS 前端
更多

推荐作者

daid

文章 0 评论 0

我心依旧

文章 0 评论 0

晒暮凉

文章 0 评论 0

微信用户

文章 0 评论 0

DS

文章 0 评论 0

〆凄凉。

文章 0 评论 0

更多

友情链接

文江博客
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站，即表示您同意使用 Cookies 和您的相关数据。
    原文