如何在凸包内获取均匀分布的点？

Question

如何在凸包内获取均匀分布的点？

pythonnumpyscipycomputational-geometryconvex-hull

4

给定一组点，

points = np.random.randn(...) # n 3d points

我希望能将一个由一组三维点列表（形如nx3的np.array）定义的凸包体积均匀填充。

我可以通过以下方式获取凸包：

hull = scipy.spatial.ConvexHull(points)

如何以最快的方式获取一个点列表，使这些点均匀地填充该凸壳的体积？

- Gulzar

泊松盘采样怎么样？ - Felipe Gutierrez

3个回答

4

如果有人需要帮助，这里是实现Yves答案的代码：

在凸包的边界框中均匀绘制点。
拒绝掉落在凸包外部的点。

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull
from matplotlib.patches import Rectangle
from matplotlib.path import Path

#Initial positions
pos = np.random.uniform(-30, 10, (20, 2))

#Convex hull of initial positions
hull = ConvexHull( pos )

#Bounding box
bbox = [hull.min_bound, hull.max_bound]

#Hull path
hull_path = Path( hull.points[hull.vertices] )

#Draw n random points inside the convex hull
n = 10**3
rand_points = np.empty((n, 2))
for i in range(n):
    #Draw a random point in the bounding box of the convex hull
    rand_points[i] = np.array([np.random.uniform(bbox[0][0], bbox[1][0]), np.random.uniform(bbox[0][1], bbox[1][1])])

    #We check if the random point is inside the convex hull, otherwise we draw it again            
    while hull_path.contains_point(rand_points[i]) == False:
        rand_points[i] = np.array([np.random.uniform(bbox[0][0], bbox[1][0]), np.random.uniform(bbox[0][1], bbox[1][1])])

#Plot
plt.scatter(pos[:, 0],pos[:, 1], marker='o',  c='blue', alpha = 1, label ='Initial points')
for simplex in hull.simplices:
        plt.plot(hull.points[simplex, 0], hull.points[simplex, 1], '-k')
plt.gca().add_patch(Rectangle((bbox[0][0], bbox[0][1]), bbox[1][0] - bbox[0][0], bbox[1][1] - bbox[0][1],facecolor = 'None', edgecolor = 'cyan'))
plt.scatter(rand_points[:, 0],rand_points[:, 1], marker='o',  c='red', alpha = 0.31, label ='Random points inside hull')
plt.legend()
plt.savefig("fig.png", dpi = 300)

- Puco4

这是未矢量化的，速度非常慢。您能制作一个矢量化版本吗？当前状态下，无法使用。 - Gulzar

你说的“向量化”是什么意思，@Gulzar？我不是程序员，也不深入了解Python，但我需要运行这个算法，最终成功了并决定分享出来。如果你能指出如何改进它，我会进行修改。 - Puco4

我不能教授所有相关的内容，但如果你感兴趣，可以从这里开始 https://realpython.com/numpy-array-programming/ - Gulzar

重点是Python循环速度非常慢。Numpy向量操作是用C实现的，并具有许多底层优化。在生产级别代码中使用它们是必须的。 - Gulzar

2

@Gulzar 我尝试使用列表推导式将for循环向量化，但由于while循环的原因，我没有成功。如果有人能想到更好的方法，我会很乐意改进答案。无论如何，这个算法在我的特定问题中起作用，我相信它仍然可以帮助其他人。 - Puco4

2

在边界框中均匀绘制点，并拒绝那些不在凸壳内的点。（由于凸壳是凸的，因此可以在线性时间O（F）内完成而无需预处理，在预处理后可以在对数时间O（log F）内完成，将面投影到平面上并考虑平面细分）。

- user1196549

抱歉，我不明白该怎么做。F是什么？预处理是什么？我应该将哪个平面投影到哪些点上，为什么这样就足够了？ - Gulzar

@Gulzar：F是面数。您可以在任何平面上进行投影。请参阅https://en.wikipedia.org/wiki/Point_location。 - user1196549

哦，是的，我已经实现了那个。我认为这可以更快地完成，因为实际上并不是所有点都必须被分类为内部或外部。你如何实现log(f)？ - Gulzar

1

这有点危险，因为您可能需要生成大量的点。除非您旋转坐标轴以匹配外壳的主轴，否则与其边界框相比，您的外壳可以任意小（想象一下由绕线x = y = z旋转的瘦椭圆表示的外壳），并且您最终可能会拒绝大多数生成的点。当然，通过轴变换和3D，您的最坏情况拒绝率为pi / 6 = 52％，但随着维度的增加，您的最坏情况拒绝率将逐渐恶化。4d = pi ^ 2/32 = 30％，5d：16％，6d：8％。 - Daniel F

据我所知，数字3不是一个增长的数字。当然，你可以找到一些极端情况，即使在二维空间中效率也为零。我的回答是针对实际目的的，尽管高效的解决方案并不容易实现，只有在面数较多的情况下才有意义。 - user1196549

1

@DanielF：时间稳定性？请详细说明。 - user1196549

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Daniel F · Accepted Answer

Find delaunay simplices of the hull
randomly sample the simplices based on their area
for each simplex, find uniform distribution of sampled points using dirichelet distribution

multiply the distributions with the simplices to find final points.

import numpy as np
from numpy.linalg import det
from scipy.stats import dirichlet


def dist_in_hull(points, n):
    dims = points.shape[-1]
    hull = points[ConvexHull(points).vertices]
    deln = hull[Delaunay(hull).simplices]

    vols = np.abs(det(deln[:, :dims, :] - deln[:, dims:, :])) / np.math.factorial(dims)    
    sample = np.random.choice(len(vols), size = n, p = vols / vols.sum())

    return np.einsum('ijk, ij -> ik', deln[sample], dirichlet.rvs([1]*(dims + 1), size = n))

编辑：功能化并扩展到更高的维度（警告：ConvexHull仅适用于9D以下）