Poisson Disk Sampling

wood logs pile

Poisson disk sampling is a technique for generating random points in the given space such that it ensures that points are distributed more evenly compared to purely random sampling, and it helps to avoid clustering or large gaps between points. The technique balances randomness with uniformity, making it suitable for procedural placement of trees in forest or other such tasks where clumping is undesirable.

Classical Algorithm for Poisson Disk

The old method involves generating random points in space. Each random point is rejected if it is closer to any of the existing points than the given radius, otherwise it is accepted. This approach is often referred to as dart throwing (paper).

This is the most naïve & easiest of the approaches. However, it is exponentially slow as the number of points increase. Also, it does not guarantee that the algorithm will generate maximal poisson disk distribution, since it becomes increasingly difficult to generate a point randomly in the desired location if most of the space is already used up by other points. This also makes it difficult to signal termination of algorithm.

Bridson’s Algorithm

I am going to use this since it is faster & is applicable to n-dimensional space (paper). According to author, it works in linear time.

In Nutshell

In nutshell, Bridson’s algorithm starts with a randomly selected point in space. To ensure that a newly generated point after this one must never be within radius r of this point, this algorithm generates new random point within a range between r and 2r centered on the original point. After new point has been created, it is still not finalized, since it could be within range of some other point across neighboring areas. So we add an extra validation to check all neighboring cells & see if we are not very near to some other point there. If the point passes this validation, it is emitted as a legal sample. This process is iteratively applied until the entire space is adequately sampled.

The Algorithm Steps:

From the original paper, the algorithm steps are as follows & I have simplified each step after quotation:

Initial Step:

“The algorithm takes as input the extent of the sample domain in Rn , the minimum distance r between samples, and a constant k as the limit of samples to choose before rejection in the algorithm (typically k = 30).”

It means we need to create a function and pass extents (width & height in 2D) of our space, the radius & the constant k.

def generate_bridson_sampling_points(width=1.0, height=1.0, radius=0.025, num_neighbors=30):
    pass

Step 0:

“Initialize an n-dimensional background grid for storing samples and accelerating spatial searches. We pick the cell size to be bounded by r/sqrt(n), so that each grid cell will contain at most one sample, and thus the grid can be implemented as a simple n dimensional array of integers: the default −1 value indicates no sample, a non-negative integer gives the index of the sample located in a cell.”

In case of 2D, we can use 2D array of int to create a grid. The size of each cell must be r/sqrt(2) for 2D grid and r/sqrt(n) for nth-dimensional grid. r is radius as mentioned above & x is the cell’s length. The reason why we use grid with this cell size is because we want to make our algorithm fast. If we choose some other value for cell size, say a big value, then our algorithm (which is instructed to sample points r-distance away from existing point) could sample multiple poisson disc points within a same cell, which we are avoiding because we want to iterate faster through the grid. By having a single sample per cell, we can make it easy for computer to iterate over the grid faster.

Now if we kept our cell size to be x = r/sqrt(2), any next sample, which is going to be generated r-distance away from existing point, will surely fall within the next cell.

By default, all grid points will be initialized with -1 value, meaning that there is no sample here.

Step 1:

“Select the initial sample, x0, randomly chosen uniformly from the domain. Insert it into the background grid, and initialize the “active list” (an array of sample indices) with this index (zero).”

We start the process with selecting a random point from our space to begin with & we insert it into the grid as well is in the active list (which is an array of sampled points).

Step 2:

“While the active list is not empty, choose a random index from it (say i). Generate up to k points chosen uniformly from the spherical annulus between radius r and 2r around xi. For each point in turn, check if it is within distance r of existing samples (using the background grid to only test nearby samples). If a point is adequately far from existing samples, emit it as the next sample and add it to the active list. If after k attempts no such point is found, instead remove i from the active list.”

We continue with a loop whose termination condition is that the active list, which is initially filled with our randomly selected point, is not empty. We choose a random index from it (initially, there is only a single index that we created, so select it). Now from this randomly selected point, attempt to generate random points around it within range(r, 2r); and for each those point check if it is at least r distance away from all existing samples. If it is indeed at least r-distance away from existing samples, just add it to the active list. Else if after k attempts, we could not generate a point that is far enough, just discard the current point from the active list around which we are creating samples.

This way, it starts with a single randomly selected point & progresses all the way by creating more samples around it in range(r, 2r) until the entire space within the given bounds is filled with poisson disc samples.

Python Implementation:

import numpy as np
import matplotlib.pyplot as plt


def generate_bridson_sampling_points(width=1.0, height=1.0, radius=0.025, num_neighbors=30):
    def calculate_squared_distance(point1, point2):
        return (point1[0] - point2[0]) ** 2 + (point1[1] - point2[1]) ** 2

    def generate_random_points_around(center_point, num_points=1):
        # Note: This is not uniformly distributed around the center_point, but it's acceptable for this purpose
        radii = np.random.uniform(radius, 2 * radius, num_points)
        angles = np.random.uniform(0, 2 * np.pi, num_points)
        new_points = np.empty((num_points, 2))
        new_points[:, 0] = center_point[0] + radii * np.sin(angles)
        new_points[:, 1] = center_point[1] + radii * np.cos(angles)
        return new_points

    def is_point_within_limits(point):
        return 0 <= point[0] < width and 0 <= point[1] < height

    def get_neighborhood_indices(shape, index, n=2):
        row, col = index
        row_start, row_end = max(row - n, 0), min(row + n + 1, shape[0])
        col_start, col_end = max(col - n, 0), min(col + n + 1, shape[1])
        indices = np.dstack(np.mgrid[row_start:row_end, col_start:col_end])
        indices = indices.reshape(indices.size // 2, 2).tolist()
        indices.remove([row, col])
        return indices

    def is_point_in_neighborhood(point):
        i, j = int(point[0] / cell_size), int(point[1] / cell_size)
        if grid_mask[i, j]:
            return True
        for (i, j) in neighborhood_indices[(i, j)]:
            if grid_mask[i, j] and calculate_squared_distance(point, points_grid[i, j]) < squared_radius:
                return True
        return False

    def add_point_to_grid(point):
        points_list.append(point)
        i, j = int(point[0] / cell_size), int(point[1] / cell_size)
        points_grid[i, j], grid_mask[i, j] = point, True

    # Calculate cell size based on radius and dimensions
    cell_size = radius / np.sqrt(2)
    num_rows = int(np.ceil(width / cell_size))
    num_cols = int(np.ceil(height / cell_size))

    # Calculate squared radius for squared distance comparison
    squared_radius = radius * radius

    # Initialize arrays for points and grid
    points_grid = np.zeros((num_rows, num_cols, 2), dtype=np.float32)
    grid_mask = np.zeros((num_rows, num_cols), dtype=bool)

    # Generate neighborhood indices for each grid cell
    neighborhood_indices = {}
    for i in range(num_rows):
        for j in range(num_cols):
            neighborhood_indices[(i, j)] = get_neighborhood_indices(grid_mask.shape, (i, j), 2)

    # Initialize list to store points
    points_list = []

    # Add an initial random point to the grid
    add_point_to_grid((np.random.uniform(width), np.random.uniform(height)))

    # Generate points using Bridson's algorithm
    while len(points_list):
        random_index = np.random.randint(len(points_list))
        current_point = points_list[random_index]
        del points_list[random_index]

        new_points = generate_random_points_around(current_point, num_neighbors)
        for new_point in new_points:
            if is_point_within_limits(new_point) and not is_point_in_neighborhood(new_point):
                add_point_to_grid(new_point)

    return points_grid[grid_mask]


if __name__ == '__main__':
    plt.figure()
    plt.subplot(1, 1, 1, aspect=1)

    sampled_points = generate_bridson_sampling_points()
    X = [x for (x, y) in sampled_points]
    Y = [y for (x, y) in sampled_points]
    plt.scatter(X, Y, s=10)
    plt.xlim(0, 1)
    plt.ylim(0, 1)
    plt.show()

More reading

While I tried my best to explain the algorithm, I think there are still some limitations in my explanation that I’ll like to address later. Also, to make it more interesting, I’d like to implement something by making use of poisson disc.

Table of Contents