How to get uniformly distributed points in convex hull?

Question

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?

Answer 1

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)

Answer 2

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)

Answer 3

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