On Poisson Disk Sampling

In 2024, a team of nine mathematicians released a monstrous, nearly 1,000 page proof of the geometric Langlands conjecture. It is a crowning achievement in pure mathematics, and I have accepted that I will never understand even the statements that they proved, much less the proof itself.

On the total opposite end of the spectrum, in 2007, Robert Bridson published a one page paper that has nearly 1,000 citations and takes less than 10 minutes to fully understand. It presents a simple solution to a problem that commonly arises in computer graphics and simulations: placing things randomly, but not too close together.

Say you’re trying to procedurally generate a forest and need a way to place the trees. The problem with plain random sampling is obvious:

Some of the trees would be on top of each other! What we need is the ability to set a minimum distance between any two trees. A distribution of trees that obeys this rule is called a Poisson disk distribution. We can try a naive rejection sampling approach where we throw random darts and reject any point that falls within that minimum distance of any other point, but without a more clever data structure, it takes linear time to check collisions for each sample and the rejection rate quickly approaches one. Bridson’s algorithm gives us an efficient way to do this.

Bridson’s Algorithm

Suppose the desired minimum distance between points is $r$ and we are working in a $d$-dimensional space. Bridson’s algorithm goes as follows:

1. Partition the space into a grid of side length $\frac{r}{\sqrt{d}}$. This guarantees that each grid cell can have at most one point inside it.
2. Initialize a list active to have one random point chosen uniformly from the space.
3. While active is non-empty:
   1. Select an element $p$ from active uniformly at random.
   2. Uniformly sample the annulus centered at $p$ of inner radius $r$ and outer radius $2r$ at most $k$ times. If a valid Poisson disk sample is found, using the grid for efficient collision detection, add it to active and pick a new $p$. If no valid point is found within $k$ attempts, remove $p$ from active. Bridson recommends setting $k = 30$.

The easiest way to uniformly sample the annulus is to generate a random unit vector $\vec{v} \in \mathbb{R}^d$ and a number $x$ chosen uniformly from the interval $\left[\frac{1}{2^d}, 1\right)$, and then your final sample is $2rx^{1/d}\cdot\vec{v}$. A few years ago, I made a video that explains why this works. In two dimensions, picking a unit vector is equivalent to picking an angle $\theta \in [0, 2\pi)$. In higher dimensions, you can normalize a vector where each component is sampled from a normal distribution.

Improvements

There are two simple improvements to Bridson’s algorithm that I have found drastically reduce the number of iterations required to generate the same number of points. The first works in two dimensions, but the second works in higher dimensions as well.

Let’s start with the two dimensional improvement. Consider when the algorithm places a point $p$ and then samples its annulus to get a new point $q$. We call $p$ the parent of $q$. There is valuable information stored in the relation between these points. When we inevitably sample the annulus centered at $q$, there is an entire range of angles that we need not consider because the points within them would be too close to $p$. This range is represented by the dotted lines in the following figure.

While the visual intuition is easy to grasp, translating it into a formula is a tedious trigonometry exercise. I’ll spare you the details and claim without proof that the cone formed by the dotted lines is centered at angle $\alpha$ and its width is $2\beta$ where

The only interesting part of this formula is the minimum that appears in the equation for $\beta$. This accounts for the fact that either the inner or outer circle of the annulus can bound the cone, depending on the distance between $p$ and $q$. We have to pick the minimum to guarantee that the intersection of the cone and the annulus is entirely contained in the circle. Note how the boundary points of the cone jump from the outer circle to the inner when the distance between the points crosses $\sqrt{3} \cdot r$. The first term in the minimum is the angle of intersection with the outer circle and the second term is that with the inner circle.

Implementing this just requires storing the parent of each point. Then you can calculate the angles of the cone and generate the angle $\theta$ for the next sample in the range outside the cone.

The following graph shows how many points were generated by Bridson’s algorithm with and without the parental optimization. Each datapoint is the average of 100 trials on a grid of side length $\ell = 100$ with $r = 1$.

This improvement could likely be generalized to higher dimensions, but it would require more space to store the contact vectors of the annuli, and the returns would likely diminish because the volume of the intersection of an annulus with a sphere becomes proportionally insignificant in higher dimensions. Similarly, there is nothing stopping us from storing the children of a point in addition to its parent to eliminate even more sections of the annulus, but this would also require more storage and it would make selecting $\theta$ far slower.

Instead of changing how we pick the angle to the next sample, the second improvement changes how we pick the distance to the next sample. Consider the distribution of the distances from each point in the annulus to its center. Its cumulative distribution function (CDF) is proportional to $x^d$ on the interval $[r, 2r]$. For an explanation of this, I again defer to my video. But what happens if we change the exponent to be some constant $c$ other than $d$? Then we can move the points closer to or further away from the center. For $c \neq 0$, the exact CDF is

The following figure lets you see what the CDF and 500 random samples in the annulus look like as $c$ changes. Remember that $c = 2$ gives a uniform distribution.

You may have noticed that the slider allows you to set $c = 0$ even though $F_0(x)$ is undefined due to a division by zero. To rememdy this, we define $F_0(x) = \lim_{c \to 0} F_c(x)$, which is

In order to sample the radius using an arbitrary value of $c$, we can apply inverse transform sampling. When $c = 0$, the radius should be $r \cdot 2^x$ where $x$ is a uniform random variable on the interval $[0, 1)$. Otherwise, we use $2ry^{1/c}$ where $y$ is a uniform random variable between 1 and $\frac{1}{2^c}$. The bounds of the interval swap depending on whether $c$ is positive or negative.

The following heatmap shows the impact $c$ has on the number of points generated by Bridson’s algorithm. We use the same experimental setup as before with 100 trials on a grid of side length $\ell = 100$ with $r = 1$.

This seems to suggest that we should set $c$ to be some very negative number, or even take the limit as $c$ approaches negative infinity, forcing each point to be a distance of exactly $r$ from its parent. While it is true that this would maximize the number of points generated and create more tightly packed configurations, it would do so at the expense of the distribution “feeling” random. In the extreme case of $c = -\infty$, you get many artifacts like long strings of points and gaps where the restricted distance cannot reach. It is also not difficult to reconstruct the tree of how the points were generated after the fact.

So we need to strike a balance between maximizing the density of points and preserving randomness. If you fix some $15 \leq k \leq 40$, I have found it best to set $c = -1.4 - \frac{17}{\sqrt{k}}$ when using the parental optimization. This equation was derived empirically to make the number of points generated roughly match the expected output from a uniform and maximal Poisson disk sampler (more on this later). For every $k$ between 15 and 40, I did a binary search to find the value of $c$ that would make it so if you drew a circle of radius $r/2$ around every point, those circles would take up 54.7% of the entire area. This percentage is the saurated coverage of circular disks under the random sequential adsorption model. Of course, this only applies in two dimensions and $c$ will need to be tuned differently in higher dimensions.

Stippling

So far, we have kept $r$ constant, but there is no reason for this. We can dynamically set the minimum distance between points according to a function $r\colon \mathbb{R}^d \to \mathbb{R}$. So after placing a point $p$, we sample an annulus of inner radius $r(p)$. A fun application of this is to define $r$ as the brightness of each pixel in an image to produce a stippling effect. The figure on the left does this in black and white, and the figure on the right combines three sets of Poisson disk samples, one for each color channel.

Birdson’s algorithm is an inherently sequential one, but there are others that are designed to be executed in parallel, resulting in massive performance boosts. My favorite of these is PixelPie which runs entirely on the GPU. I used this algorithm to create Poisson Cam, a realtime video stippler using Poisson disk sampling. This was one of my favorite projects to work on because it taught me shader programming, Rust, stream compaction algorithms, and of course the PixelPie algorithm itself.

Maximality and Uniformity

In 2022, Scott A. Mitchell published a really cool paper introducing a beautiful new way of generating Poisson disk samples in two dimensions that, to my knowledge, has received no attention since its publication. This last section is dedicated to it.

Before I describe the algorithm, I want to explain the three things that make it better than Bridson’s:

1. Maximality. After the algorithm finishes, it is guaranteed to be impossible to fit another point without violating the Poisson disk property.
2. Uniformity. The algorithm samples a uniform distribution over all maximal sets of Poisson disk samples.
3. Determinism. The algorithm does not rely on rejection sampling, meaning there are no failed attempts to place points.

Maximality and uniformity have been achieved by many algorithms, the most popular of which is hierarchical dart throwing, but Mitchell’s is the first algorithm to do this without some kind of rejection sampling. It is also very performant. While I have not done any rigorous benchmarking, I have found that Mitchell’s implementation of the algorithm runs about as quickly and generates roughly the same number of points as my implementation of Bridson’s algorithm (i.e., with the parental optimization) with $k = 20$ and $c = -5.2$. As far as I can tell, the only downside to this algorithm is that it is substantially more complex. The following is a greatly simplified summary:

1. Partition the space into a grid of side length $\frac{r}{\sqrt{d}}$.
2. While there is room to place another point:
   1. Randomly select a cell $c$ from the grid weighted by the remaining areas.
   2. Decompose $c$ into disjoint triangles and chocks.¹
   3. Select a random triangle or chock $t$ weighted by area.
   4. Uniformly sample a point $p$ from $t$. Add it to the final set.
   5. Carve out a circle of radius $r$ centered at $p$ from the grid.

The original paper does an excellent job of explaining the algorithm in detail, so the most helpful contribution I can make is the following visualization: