CODEWITHSUNDEEP
pythonnumpyscipycomputational-geometryconvex-hull
Gulzar
Gulzar

Reputation: 27946

How to get uniformly distributed points in convex hull?

Given a set of points, 

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

I would like to uniformly fill the volume defined by the convex hull in which they lie by a list (np.array of shape nx3) of 3d points.

I can get the convex hull by

hull = scipy.spatial.ConvexHull(points)

What would be the fastest way to get a list of points that uniformly fills this hull's volume?

Upvotes: 6

Views: 3439

Answers (3)

Daniel F
Daniel F

Reputation: 14399

  1. Find delaunay simplices of the hull

  2. randomly sample the simplices based on their area

  3. for each simplex, find uniform distribution of sampled points using dirichelet distribution

  4. 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))

EDIT: functionalized and extended to higher dimensions (Warning: ConvexHull only works up to 9D)

Upvotes: 8

Puco4
Puco4

Reputation: 633

In case it helps someone, here is the code to implement Yves's answer:

  1. Draw points uniformly in the bounding box of the convex hull.
  2. Reject the points falling outside the hull.
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)

enter image description here

Upvotes: 4

user1196549
user1196549

Reputation:

Draw points uniformly in the bounding box and reject those that aren't inside the hull. (As the hull is convex, this can be done in linear time O(F) without preprocessing, and logarithmic time O(log F) after preprocessing, projecting the faces to a plane and considering a planar subdivision).

Upvotes: 2

Related Questions