Dart-Throwing 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:

1. Initialize an empty set \(S\). This set will contain the final sample points.

2. While \(|S| < N\):

   • Generate a random sample \(P\) within \(E\).

   • Check the distance \(d_i\) between \(P\) and each point \(s_i\) in \(S\). 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 might be impossible for any newly generated point to be approved, causing the algorithm to potentially 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>


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 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 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;
    }
                            
public:
    
    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)
    {
        std::vector<std::array<float, D>> samples;
    
        int attempt = 0;
        while((samples.size() < numSamples) && (attempt < maxNumAttempts)){
            std::array<float, D> point = Random(min, max);
            
            bool approved = true;
            for(int i = 0; i < samples.size(); i++){
                if(Length(Subtract(point, samples[i])) < minDistance){
                    approved = false;
                    break;
                }
            }
            if(approved)
                samples.push_back(point);
                
            attempt++;
        }
    
        return samples;
    }
};

The github repository containing the entire code for the Dart-Throwing algorithm can be found here.

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

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