Poisson Disk Sampling

Cover image credits: https://github.com/udit/poisson-disc-sampling

Poisson disk sampling is a technique for generating random points in the given space such that the points are distributed more evenly compared to purely random sampling, and it helps to avoid clustering or large gaps between points.

In my old games, I used simple noise for procedural placement of trees in forest. This resulted in many trees clump together as the noise samples are generated too near in some cases.

Table of Contents

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 random
import math
import matplotlib.pyplot as plt

def calculate_squared_distance(point1, point2):
    dx = point1[0] - point2[0]
    dy = point1[1] - point2[1]
    return dx * dx + dy * dy

def generate_random_points_around(center_point, radius, num_points=1):
    new_points = []
    for _ in range(num_points):
        r = random.uniform(radius, 2 * radius)
        a = random.uniform(0, 2 * math.pi)
        x = center_point[0] + r * math.sin(a)
        y = center_point[1] + r * math.cos(a)
        new_points.append((x, y))
    return new_points

def is_point_within_limits(point, width, height):
    if 0 <= point[0] < width and 0 <= point[1] < height:
        return True
    else:
        return False

def get_neighborhood_indices(grid_size, index, n=2):
    row, col = index
    row_start = max(row - n, 0)
    row_end = min(row + n + 1, grid_size[0])
    col_start = max(col - n, 0)
    col_end = min(col + n + 1, grid_size[1])
    indices = []
    for r in range(row_start, row_end):
        for c in range(col_start, col_end):
            if (r, c) != (row, col):
                indices.append((r, c))
    return indices

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

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

def generate_bridson_sampling_points(width=1.0, height=1.0, radius=0.025, num_neighbors=30):
    cell_size = radius / math.sqrt(2)
    num_rows = int(math.ceil(width / cell_size))
    num_cols = int(math.ceil(height / cell_size))
    squared_radius = radius * radius
    points_grid = []
    for _ in range(num_rows):
        row = [(0, 0)] * num_cols
        points_grid.append(row)
    neighborhood_indices = {}
    for i in range(num_rows):
        for j in range(num_cols):
            neighborhood_indices[(i, j)] = get_neighborhood_indices((num_rows, num_cols), (i, j), 2)
    points_list = []
    initial_point = (random.uniform(0, width), random.uniform(0, height))
    add_point_to_grid(initial_point, points_list, points_grid, cell_size)
    while points_list:
        random_index = random.randint(0, len(points_list) - 1)
        current_point = points_list[random_index]
        del points_list[random_index]
        new_points = generate_random_points_around(current_point, radius, num_neighbors)
        for new_point in new_points:
            if is_point_within_limits(new_point, width, height) and not is_point_in_neighborhood(new_point, points_grid, neighborhood_indices, cell_size, squared_radius):
                add_point_to_grid(new_point, points_list, points_grid, cell_size)
    sampled_points = []
    for row in points_grid:
        for point in row:
            if point != (0, 0):
                sampled_points.append(point)
    return sampled_points

if __name__ == '__main__':
    plt.figure()
    plt.subplot(1, 1, 1, aspect=1)
    sampled_points: list = generate_bridson_sampling_points()
    X = []
    Y = []
    for point in sampled_points:
        X.append(point[0])
        Y.append(point[1])
    plt.scatter(X, Y, s=10)
    plt.xlim(0, 1)
    plt.ylim(0, 1)
    plt.show()

Uses

Object placement; so the objects do not overlap. In tree placement, we want trees to be evenly spaced apart while still retaining randomness. If we are making a tool such as level generator, objects placement can be made better this way.
Procedural textures; as a mask for some effects such as soap bubbles, air bubbles in water or juice. Also, textures usually contain repeating patterns that can be made using this.
If we wish to generate texture entirely out of dense dots.
Point cloud from 3D mesh.

gameidea