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. The solution to this is poisson disc sampling.
Comparison of noise samples with poisson samples:
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).”
RADIUS = 8
K = 30 # Number of samples to choose before rejection in the algorithm
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.
GRID_SIZE = RADIUS / math.sqrt(2)
# Create a grid to store points
cols, rows = int(WIDTH / GRID_SIZE) + 1, int(HEIGHT / GRID_SIZE) + 1
def initialize_grid():
return [[None for _ in range(rows)] for _ in range(cols)]
grid = initialize_grid()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. You can initialize it to be None/null, its your choice.
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).
# List to store active points
active = []
points = []
# restart, called whenever we want to start it over again.
def restart_simulation(start_point):
global grid, active, points
grid = initialize_grid()
points = [start_point]
active = [start_point]
col = int(start_point[0] / GRID_SIZE)
row = int(start_point[1] / GRID_SIZE)
grid[col][row] = start_point
# Initialize with a random point
initial_point = (random.uniform(0, WIDTH), random.uniform(0, HEIGHT))
restart_simulation(initial_point)
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.
# Main loop
running = True
while running:
if active:
rand_index = random.randint(0, len(active) - 1)
point = active[rand_index]
found = False
for _ in range(K):
new_point = generate_point_around(point) # defined below
if in_bounds(new_point) and fits(new_point): # functions defined below
points.append(new_point)
active.append(new_point)
col = int(new_point[0] / GRID_SIZE)
row = int(new_point[1] / GRID_SIZE)
grid[col][row] = new_point
found = True
break
if not found:
active.pop(rand_index)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.
Helper functions that are used above:
# Helper functions
def distance(p1, p2):
return math.hypot(p1[0] - p2[0], p1[1] - p2[1])
def generate_point_around(point):
r = RADIUS * (random.random() + 1)
angle = 2 * math.pi * random.random()
new_x = point[0] + r * math.cos(angle)
new_y = point[1] + r * math.sin(angle)
return new_x, new_y
def in_bounds(point):
return 0 <= point[0] < WIDTH and 0 <= point[1] < HEIGHT
def fits(point):
col = int(point[0] / GRID_SIZE)
row = int(point[1] / GRID_SIZE)
for i in range(max(col - 2, 0), min(col + 3, cols)):
for j in range(max(row - 2, 0), min(row + 3, rows)):
neighbor = grid[i][j]
if neighbor is not None and distance(point, neighbor) < RADIUS:
return False
return TrueComplete Python Poisson Disc Code for PyGame
import pygame
import random
import math
# Pygame setup
pygame.init()
WIDTH, HEIGHT = 800, 800
screen = pygame.display.set_mode((WIDTH, HEIGHT))
pygame.display.set_caption("Poisson Disc Sampling using Bridson's Algorithm")
clock = pygame.time.Clock()
# Parameters for Poisson disc sampling
RADIUS = 8
K = 30 # Number of samples to choose before rejection in the algorithm
GRID_SIZE = RADIUS / math.sqrt(2)
# Create a grid to store points
cols, rows = int(WIDTH / GRID_SIZE) + 1, int(HEIGHT / GRID_SIZE) + 1
def initialize_grid():
return [[None for _ in range(rows)] for _ in range(cols)]
grid = initialize_grid()
# List to store active points
active = []
points = []
# Helper functions
def distance(p1, p2):
return math.hypot(p1[0] - p2[0], p1[1] - p2[1])
def generate_point_around(point):
r = RADIUS * (random.random() + 1)
angle = 2 * math.pi * random.random()
new_x = point[0] + r * math.cos(angle)
new_y = point[1] + r * math.sin(angle)
return new_x, new_y
def in_bounds(point):
return 0 <= point[0] < WIDTH and 0 <= point[1] < HEIGHT
def fits(point):
col = int(point[0] / GRID_SIZE)
row = int(point[1] / GRID_SIZE)
for i in range(max(col - 2, 0), min(col + 3, cols)):
for j in range(max(row - 2, 0), min(row + 3, rows)):
neighbor = grid[i][j]
if neighbor is not None and distance(point, neighbor) < RADIUS:
return False
return True
def restart_simulation(start_point):
global grid, active, points
grid = initialize_grid()
points = [start_point]
active = [start_point]
col = int(start_point[0] / GRID_SIZE)
row = int(start_point[1] / GRID_SIZE)
grid[col][row] = start_point
# Initialize with a random point
initial_point = (random.uniform(0, WIDTH), random.uniform(0, HEIGHT))
restart_simulation(initial_point)
# Main loop
running = True
while running:
screen.fill((30, 30, 30))
for event in pygame.event.get():
if event.type == pygame.QUIT:
running = False
elif event.type == pygame.MOUSEBUTTONDOWN:
mouse_pos = pygame.mouse.get_pos()
restart_simulation(mouse_pos)
if active:
rand_index = random.randint(0, len(active) - 1)
point = active[rand_index]
found = False
for _ in range(K):
new_point = generate_point_around(point)
if in_bounds(new_point) and fits(new_point):
points.append(new_point)
active.append(new_point)
col = int(new_point[0] / GRID_SIZE)
row = int(new_point[1] / GRID_SIZE)
grid[col][row] = new_point
found = True
break
if not found:
active.pop(rand_index)
# Draw points
for p in points:
pygame.draw.circle(screen, (255, 255, 255), (int(p[0]), int(p[1])), 2)
# Draw active points
for a in active:
pygame.draw.circle(screen, (255, 0, 0), (int(a[0]), int(a[1])), 2)
pygame.display.flip()
# clock.tick(1000)
pygame.quit()
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.
- Point cloud from 3D mesh.
More reading
- Godot implementation by someone here: https://github.com/udit/poisson-disc-sampling