Poisson-Disk Sampling

Poisson-Disk sampling is an algorithm used to generate random samples within a user-defined domain, such that each two points are seperated by a minimum distance $r$. It produces a well-distributed set of random samples that are tightly-packed, yet not clumped or clustered together.

Dart-Throwing Algorithm

One naive implementation of Poisson-Disk sampling is the dart-throwing algorithm. It is described as follows:

Initialize an empty set $S$. This set will contain the final sample points.
While $|S| < N$:
- Generate a random sample $P$ within $E$.
- Check the distance $d_i$ between each point $s_i$ in $S$ and $P$. If there is at least one point $s_i$ in $S$ where $d_i < r$, then we discard $P$ and try again with a new point. Otherwise, $P$ is added to $S$.

One problem with this algorithm is that at some point when the area is fully saturated, and there are still remaining samples to be generated, it will be difficult for any newly generated point to be approved, causing the algorithm to run forever. Hence, as a precaution, a maximum number of attempts $A$ is enforced, where $A >= N$. If all $A$ attempts are used, then the algorithm will terminate regardless of whether $N$ approved sample points were generated or not.

Below is an C++ implementation of the dart-throwing algorithm:

#include <random>


float Random(float min, float max)
{
    static std::default_random_engine rng(std::random_device{}());
    std::uniform_real_distribution<float> dist(min, max);

    return dist(rng);
}

std::vector<float> Random(std::vector<float> min, std::vector<float> max, int n)
{
    std::vector<float> P(n);
    
    for(int i = 0; i < n; i++)
        P[i] = Random(min[i], max[i]);
    
    return P;
}

float ComputeDistanceSquared(const std::vector<float>& P, const std::vector<float>& Q, int n)
{
    float distance = 0.0f;
    for(int i = 0; i < n; i++)
        distance += ( (P[i] - Q[i]) * (P[i] - Q[i]) );
    
    return distance;
}

std::vector<std::vector<float>> DartThrowingAlgorithm(int N, int n, int A, float r, const std::vector<float>& min, const std::vector<float>& max)
{
    if((min.size() != n) || (max.size() != n))
        return {};
        
    float rSquared = r * r;
    std::vector<std::vector<float>> S;
    
    int a = 0;
    while( (samples.size() < N) && (a < A) ){
        std::vector<float> P = Random(min, max, n);
        bool isApproved = true;
        for(int i = 0; i < samples.size(); i++){
            float dSquared = ComputeDistanceSquared(P, S[i], n);
            if(dSquared < rSquared){
                isApproved = false;
                break;
            }
        }
        if(isApproved)
            S.push_back(P);
        a++;
    }
    
    return S;
}

As discussed before, the Dart-Throwing algorithm is naive and inefficient. The best-case scenario for the algorithm happens when each generated point $P$ is approved and added to the set $S$. That is, the best-case time complexity of the algorithm is

\[T(N) = \sum_{k=1}^N D + Dk\] \[\sum_{k=1}^N D(k + 1) = D \sum_{k=1}^N (k + 1)\] \[= D \frac{N(N+3)}{2} = \Omega(D N^2)\]

The worst-case scenario for the algorithm happens when $N-1$ points are generated and approved, and the $N^{th}$ point keeps getting rejected till the $A-N+1$ attempts left are used and the algorithm terminates. That is, the worst-case time complexity of the algorithm is

\[T(N) = \sum_{k=1}^{N-1} k + (A-N+1)(N-1)\] \[= \frac{(N-1)(N-2)}{2} + (AN - A - N^2 + N + N - 1)\] \[= \frac{N^2}{2} - \frac{3N}{2} + 1 + AN - A - N^2 + 2N - 1\] \[= \frac{N^2}{2} - \frac{3N}{2} + AN - A - N^2 + 2N\] \[= -\frac{N^2}{2} + (\frac{1}{2} + A)N - A\] \[= O(N^2) + O(AN)\] \[= O(AN)\]

Bridson’s Algorithm

Bridson’s algorithm is another technique for implementing Poisson disk sampling that is much more efficient than the dart-throwing algorithm. It is described as follows:

The algorithm takes as input the extent $E$ of the sample domain in $R^n$, the minimum distance $r$ between samples, and a constant $k$.

Initialize an $n$-dimensional grid $G$ with cell size $\frac{r}{\sqrt(n)}$. Each cell in the grid $G$ will contain only one sample.
Initialize an empty set $A$. This set will contain the “active samples”.
Initialize another empty set $S$. This set will contain the final sample points.
Generate an initial sample $x \in E$ and add it $A$, $S$ and $G$.
While $A$ is non-empty:
- Randomly choose sample $a_i \in A$.
- Generate $k$ random candidate samples $\{ w_1, w_2, ...., w_k\}$ from the spherical annulus between radius $r$ and $2r$ around $a_i$
- For each candidate $w_j$:
  - Locate the cell in $G$ in which $w_j$ will be stored. Say the cell is $C_j$
  - Around $C_j$, there are $3^n - 1$ distinct neighbouring cells. Each one of these cells is either empty or have one existing sample $q_m$.
  - For each sample $q_m$:
    - Calculate its distance from $w_j$.
    - If the distance is less than $r$, then $w_j$ is rejected.
  - If all the samples $q_m$ were adequately far from $w_j$, then $w_j$ is approved, and is added to $S$, $A$ and $G$.
- If all the candidates $w_j$ were rejected, then remove $a_i$ from $A$.

Each cell in the grid $G$ has $3^n - 1$ distinct neighbors.

Proof

Let $C$ be a cell in $G$ with index $i$, such that $$ i = (i_1, i_2, ......, i_n) $$ A cell with index $j$ is a neighbor to $C$ if and only if $$ \forall_{1 \, \leq \, k \, \leq \, n} \; \; j_k = i_k ± d $$ where $d \in \{ -1, 0, 1\} $. That is,

If $d = -1$, then we get the previous cell in dimension $k$

If $d = 0$, then we get the same cell in dimension $k$

If $d = 1$, then we get the next cell in dimension $k$

This means that for every dimension $1 \, \leq \, k \, \leq \, n$, there are $3$ neighboring cells. Hence, the total number of neighbors to $C$ is $$ 3 \times 3 \times ... \times 3 = 3^n $$ However, this count includes the cell $C$ itself, which is the case when $d = 0$ in all $k$ dimensions. Thus, the total number of distinct neighbors is $$ 3^n - 1 $$ $\blacksquare$

Let $C$ be a cell in $G$ with index $i = (i_1, i_2, ......, i_n)$. We know that $C$ has $3^n$ neighbors $\{N_1, N_2, ..., N_{3^n -1}\}$, including $C$ itself. We can calculate the cell index of $N_j$ using the following function:

std::array<int, n> ComputeNeighborCellIndex(std::array<int, n> C, int n, int j)
{
    std::array<int, n> index = C;
    
    for(int i = 0; i < n; i++){
        int offset = (j % 3) - 1;
        j /= 3;
        index[i] += offset;
    }
    
    return index;
}

The idea of this function is that we represent each index $j$, where $1 \leq j \leq 3^n$ as a base-3 number (Ternary). Then using this ternary representation, we extract the offset $d$ by subtracing $1$ from every digit. This offset will be added to the index of cell $C$ to calculate the index of the neighbouring cell. That is,

Decimal	Ternary	Offset
$0$	$000$	$\{-1, -1, -1\}$
$1$	$001$	$\{-1, -1, 0\}$
$2$	$002$	$\{-1, -1, 1\}$
$3$	$010$	$\{-1, 0, -1\}$
$4$	$011$	$\{-1, 0, 0\}$
$5$	$012$	$\{-1, 0, 1\}$
$6$	$020$	$\{-1, 1, -1\}$
$7$	$021$	$\{-1, 1, 0\}$
$8$	$022$	$\{-1, 1, 1\}$
$9$	$100$	$\{0, -1, -1\}$
$10$	$101$	$\{0, -1, 0\}$

The size of each cell in the grid $G$ must be at most $\frac{r}{\sqrt(n)}$.

Proof

We know that the distance between any two points must be at least $r$. We also know that each cell in the grid $G$ must have only one point. Consider the cell $C$ in $G$. Let $d$ be the length of the diagonal of $C$ and $s$ be its cell size. Suppose for the sake of contradiction that $d > r$. This means that two points can be placed inside $C$ whose distance is greater than or equal to $r$. This violates the requirement that each cell in the grid $G$ must have only one point. Thus, $d \leq r$. Using Pythagoras theorem: $$ d^2 = s_1^2 + s_2^2 + ...... + s_n^2 = n \; s^2 $$ $$ s^2 = \frac{d^2}{n} $$ $$ s = \sqrt(\frac{d^2}{n}) = \frac{d}{\sqrt(n)} $$ Since $d \leq r$, then $$ s \leq \frac{r}{\sqrt(n)} $$ $\blacksquare$

Bridson’s algorithm is much faster than the dart-throwing algorithm. The number of iterations of the loop in Step $5$ is exactly $2N-1$, where $N$ is the number of samples. Each iteration takes $O(k)$ time. This means that the time complexity $T(N)$ of the algorithm is

\[T(N) = k(2N-1) = 2kN - k = O(2kN) - O(k) = O(kN)\]

Since $k$ is constant, then $T(N) = O(N)$ time. That is, the algorithm is linear.

Step $5$ is executed exaclty $2N-1$ times.

Proof

The loop in step $5$ of the algorithm terminates only if the active set $A$ is empty. For each iteration of the loop, exactly one of the following two events occur:

A new sample is added to $A$

An existing active sample is removed from $A$

Each sample point $s_i \in S$ must once be added to $A$ and once be removed from $A$. After it is removal from $A$, it can never be readded back into $A$. This means that the total number of insertions is $N$ and the total number of removals is $N$. Since only one event happens in each iteration of the loop, then the number of iterations of the loop is $2N$. However, the insertion of the first sample $x$ into $A$ is done in step $4$ before the loop. Thus, the total number of iterations of the loop in step $5$ is exactly $2N-1$.

$\blacksquare$

Below is a full C++ implementation of Bridson’s algorithm:

#include <random>
#include <limits>



template<int D>
class PoissonDisk{
    
private:
    
    static float Length(const std::array<float, D>& v)
    {
        float length = 0.0f;
        
        for(int i = 0; i < D; i++)
            length += (v[i]*v[i]);
        
        return std::sqrt(length);
    }
    
    static float Random(float min, float max)
    {
        static std::default_random_engine rng(std::random_device{}());
        std::uniform_real_distribution<float> dist(min, max);
        
        return dist(rng);
    }
    
    static std::array<float, D> Random(const std::array<float, D>& min, const std::array<float, D>& max)
    {
        std::array<float, D> p;
        
        for(int i = 0; i < D; i++)
            p[i] = Random(min[i], max[i]);
        
        return p;
    }
    
    static std::array<float, D> Normalize(const std::array<float, D>& v)
    {
        std::array<float, D> u = v;
        
        float length = Length(u);
        
        if(length == 0.0f)
            return u;
        
        for(int i = 0; i < D; i++)
            u[i] /= length;
        
        return u;
    }
    
    static std::array<float, D> RandomWithinSphericalAnnulus(const std::array<float, D>& sample, float minDistance)
    {
        std::array<float, D> minDir;
        for(int i = 0; i < D; i++)
            minDir[i] = -1.0f;
        
        std::array<float, D> maxDir;
        for(int i = 0; i < D; i++)
            maxDir[i] = 1.0f;
        
        std::array<float, D> randomDir = Normalize(Random(minDir, maxDir));
        
        float t = Random(minDistance, 2.0f * minDistance);
        
        std::array<float, D> newSample;
        for(int i = 0; i < D; i++)
            newSample[i] = sample[i] + (randomDir[i] * t);
        
        return newSample;
    }
        
    static std::array<float, D> Subtract(const std::array<float, D>& a, const std::array<float, D>& b)
    {
        std::array<float, D> c;
        
        for(int i = 0; i < D; i++)
            c[i] = a[i] - b[i];
        
        return c;
    }
    
    static std::array<int, D> ComputeGridIndex(const std::array<float, D>& sample, const std::array<float, D>& min, float cellSize)
    {
        std::array<int, D> gridIndex;
        
        for(int i = 0; i < D; i++)
            gridIndex[i] = (sample[i] - min[i]) / cellSize;
        
        return gridIndex;
    }
    
    static std::array<int, D> ComputeNeighborGridIndex(const std::array<int, D>& index, int rank)
    {
        int d, offset, tmp = rank;
        std::array<int, D> neighborIndex = index;
        
        for(d = 0; d < D; d++){
            offset = (tmp % 3) - 1;
            tmp /= 3;
            neighborIndex[d] += offset;
        }
        
        return neighborIndex;
    }
    
    static int ComputeFlattenedGridIndex(const std::array<int, D>& gridIndex, const std::array<int, D>& gridResolutions)
    {
        int flattenedGridIndex = gridIndex[0];
        int stride = gridResolutions[0];
        for(int i = 1; i < D; i++){
            flattenedGridIndex += (gridIndex[i] * stride);
            stride *= gridResolutions[i];
        }
        
        return flattenedGridIndex;
    }
    
    static std::array<int, D> ComputeGridResolutions(float cellSize, const std::array<float, D>& min, const std::array<float, D>& max)
    {
        std::array<int, D> gridResolutions;
        for(int i = 0; i < D; i++)
            gridResolutions[i] = std::ceil((max[i] - min[i]) / cellSize);
        
        return gridResolutions;
    }
    
    static int ComputeTotalNumGridCells(const std::array<int, D>& gridResolutions)
    {
        int totalNumGridCells = 1;
        for(int i = 0; i < D; i++)
            totalNumGridCells *= gridResolutions[i];

        return totalNumGridCells;
    }
        
public:
    
    PoissonDisk(){}
    
    static std::vector<std::array<float, D>> Generate(int numSamples, int maxNumAttempts, float minDistance, const std::array<float, D>& min, const std::array<float, D>& max)
    {
        int numNeighbors = pow(3, D);
        float cellSize = minDistance / sqrt(D);
        
        std::array<int, D> gridResolutions = ComputeGridResolutions(cellSize, min, max);
        
        int totalNumGridCells = ComputeTotalNumGridCells(gridResolutions);
        
        std::array<float, D> initialSample;
        for(int i = 0; i < D; i++)
            initialSample[i] = (max[i] + min[i]) / 2.0f;                
        std::array<int, D> initialSampleGridIndex = ComputeGridIndex(initialSample, min, cellSize);
        int initialSampleFlattenedGridIndex = ComputeFlattenedGridIndex(initialSampleGridIndex, gridResolutions);
        
        std::vector<std::array<float, D>> samples;
        samples.push_back(initialSample);
        
        std::vector<int> activeSamples;
        activeSamples.push_back(0);
        
        std::vector<int> grid(totalNumGridCells, -1);
        grid[initialSampleFlattenedGridIndex] = 0;
        
        while(!activeSamples.empty() && samples.size() < numSamples){
            std::array<float, D> sample = samples[activeSamples[activeSamples.size()-1]];
            
            bool retire = true;
            
            for(int i = 0; i < maxNumAttempts; i++){
                std::array<float, D> candidate = RandomWithinSphericalAnnulus(sample, minDistance);
                                
                std::array<int, D> candidateGridIndex = ComputeGridIndex(candidate, min, cellSize);
                int candidateFlattenedGridIndex = ComputeFlattenedGridIndex(candidateGridIndex, gridResolutions);
                
                bool candidateApproved = true;
                
                for(int m = 0; m < numNeighbors; m++){
                    
                    std::array<int, D> neighborGridIndex = ComputeNeighborGridIndex(candidateGridIndex, m);
                    
                    if(neighborGridIndex == candidateGridIndex)
                        continue;
                    
                    int neighborFlattenedGridIndex = ComputeFlattenedGridIndex(neighborGridIndex, gridResolutions);
                                        
                    if(grid[neighborFlattenedGridIndex] == -1)
                        continue;
                    
                    std::array<float, D> neighbor = samples[grid[neighborFlattenedGridIndex]];
                    
                    if(Length(Subtract(neighbor, candidate)) < minDistance){
                        candidateApproved = false;
                        break;
                    }
                }
                
                if(candidateApproved){
                    samples.push_back(candidate);
                    activeSamples.push_back(samples.size() - 1);
                    grid[candidateFlattenedGridIndex] = samples.size() - 1;
                    retire = false;
                }
            }
            
            if(retire)
                activeSamples.pop_back();
        }
        
        return samples;
    }
    
    ~PoissonDisk(){}
};

References

[1] Fast Poisson Disk Sampling in Arbitrary Dimensions by Robert Bridson

Decimal	Ternary	Offset
\(0\)	\(000\)	\(\{-1, -1, -1\}\)
\(1\)	\(001\)	\(\{-1, -1, 0\}\)
\(2\)	\(002\)	\(\{-1, -1, 1\}\)
\(3\)	\(010\)	\(\{-1, 0, -1\}\)
\(4\)	\(011\)	\(\{-1, 0, 0\}\)
\(5\)	\(012\)	\(\{-1, 0, 1\}\)
\(6\)	\(020\)	\(\{-1, 1, -1\}\)
\(7\)	\(021\)	\(\{-1, 1, 0\}\)
\(8\)	\(022\)	\(\{-1, 1, 1\}\)
\(9\)	\(100\)	\(\{0, -1, -1\}\)
\(10\)	\(101\)	\(\{0, -1, 0\}\)