Poisson Disk Sampling

Poisson disk sampling has numerous uses in computer graphics and other fields. Here I will describe how to generate Poisson disk sample sets for any sampling problem using the cySampleElimination code release.

What is Poisson Disk Sampling?

In a Poisson disk sample set no two samples are too close together. The closeness is defined by the Poisson disk radius, which is the half distance between the two closest samples. As compared to regular random sampling, Poisson disk sample sets provide a much more uniform sample distribution over the sampling domain.

The figure above shows a comparison of 1000 random samples to 1000 Poisson disk samples. Notice that Poisson disk sample points are more uniformly distributed over the sampling domain and no two sample points are too close together. Random samples, on the other hand, produce clusters and relatively large empty spaces between samples. This is not a deficiency of the particular random number generator used for this example. Uniform random samples provide uniform sample placement probability, but not a uniform coverage of the sampling domain.

Beyond providing a nice distribution of sample points, Poisson disk sample sets also provide superior convergence properties with Monte Carlo sampling. The same number of Poisson disk samples typically produces lower noise, as compared to random sampling.

Generating Poisson Disk Sample Sets

The problem with Poisson disk sample sets is that generating them has been notoriously difficult. This is because samples cannot be generated independently; each sample depends on the positions of the other samples. Simple algorithms (like dart throwing) are computationally expensive. Computationally efficient algorithms are complex and limited to certain sampling domains. High-dimensional sampling was practically impossible for anything beyond 6 or 7 dimensions.

All of these problems of Poisson disk sample set generation has been resolved by the weighted sample elimination algorithm I introduced in 2015. You can find detailed information on this algorithm here: http://www.cemyuksel.com/research/sampleelimination/

In the remainder of this solution I will provide basic information about the weighted sample elimination algorithm and describe how the source code in the cySampleElimination code release can be used for generating Poisson disk sample sets for any sampling domain in any dimension.

Weighted Sample Elimination

Weighted sample elimination is a simple algorithm for generating Poisson disk sample sets. One additional advantage of it is that it does not require specifying a Poisson disk radius property. Its input is a collection of (randomly generated) samples and it eliminates a portion of these samples such that the resulting sample set has Poisson disk property. In other words, the weighted sample elimination algorithm receives a large set of input samples and it returns a smaller set of output samples. It does not generate new samples. It merely selects a subset from a given input set.

The input samples can be generated for any sampling domain in any dimensions. The more input samples you provide, the better quality Poisson disk sample set you get, though beyond a certain percentage you get diminishing returns. I would recommend using an input set size no less than 3 times the output set size. I typically use 5 times the output set size, which provides a good balance between computation speed and sampling quality. That said, it is not a good idea to use a very small input sample set, regardless of the output sample set size. For example, if you would merely generate 10 samples, using 50 input samples may not provide a good input sample distribution. In this case, it would be a good idea to start with a much larger input sample set size.

Using cySampleElimination

The cySampleElimination code release includes a template class that can work with any given sample class type, though it is designed to work with classes in the cyVector code release (Vec2 for 2D, Vec3 for 3D, Vec4 for 4D, and Vec for arbitrary dimensions).

The first thing to do is to generate an input sample set. Let's say that we would like to have 1000 output samples, so we will begin with generating 5000 input samples. How you generate these samples is entirely up to you and your sampling problem. Below I include a simple code for generating random samples in 3D using Vec3f.

std::vector< cy::Vec3f > inputPoints(5000);
for ( size_t i=0; i<inputPoints.size(); i++ ) {
	inputPoints[i].x = (float) rand() / RAND_MAX;
	inputPoints[i].y = (float) rand() / RAND_MAX;
	inputPoints[i].z = (float) rand() / RAND_MAX;
}

Notice that all sample values are between zero and one. This is important. If your sampling domain is different, you will need to specify its dimensions. I will discuss this below. For now, let's assume that we would like all sample values to be between zero and one.

For using weighted sample elimination, I first need to create an instance of the WeightedSampleElimination class.

cy::WeightedSampleElimination< cy::Vec3f, float, 3, int > wse;

The template parameters here define the sample class type (in this case Vec3f), floating point type (float or double), the dimensionality of the point type (in this case 3), and the type for the sample index (in this case int). Unless you are working with more than 2 billion samples, sample index type int will suffice.

The constructor will set the default parameters. You can alter them if you like. I will talk about the parameters below. For now, let's assume that we are happy with the default parameters. For generating a Poisson disk sample set, all I need to do is to call the Eliminate() method of the WeightedSampleElimination class.

std::vector< cy::Vec3f > outputPoints(1000);
wse.Eliminate( inputPoints.data(), inputPoints.size(), 
               outputPoints.data(), outputPoints.size() );

That's it! The output points will have Poisson disk property. Notice that we didn't have to specify a Poisson disk radius. The eliminate method can take a few more parameters. I will discuss them below.

Progressive sampling

For some sampling problems it is desirable to introduce samples one by one, and stop introducing new samples once a desired property is reached. This is called progressive sampling (or adaptive sampling). The Eliminate() method can easily produce a progressive sample set by setting its progressive argument as true.

wse.Eliminate( inputPoints.data(), inputPoints.size(), 
               outputPoints.data(), outputPoints.size(),
               true );

In this case the samples in the output set is ordered, such that when the samples are introduced one by one in this order, each subset in the sequence becomes a Poisson disk sample set. The animation below shows an example sequence.

Sampling Arbitrary Volumes

As indicated above, our example sample points had values between zero and one. If you would like to sample a different domain, you must specify the bounds of your domain. For example, let's assume that our sampling domain is a box in 3D, where the minimum bound is (-2,-1,0) and the maximum bound is (2,1,1). We can set these bounds, as shown below.

cy::Vec3f minimumBounds(-2,-1, 0);
cy::Vec3f maximumBounds( 2, 1, 1);
wse.SetBoundsMin( minimumBounds );
wse.SetBoundsMax( maximumBounds );

Note that these bounds specify a box-shaped domain. If you would like to sample an arbitrary volume in 3D, setting these bounds won't help. Instead, we must directly specify the d_max argument of the Eliminate() method. The d_max argument determines the maximum possible distance between two samples in the sampling domain for the desired number of output samples. This is a well-defined value in 2D and 3D, and it can be approximated in higher dimensions. The GetMaxPoissonDiskRadius() method helps with this parameter. The d_max argument should be twice the radius value that this method returns.

float d_max = 2 * wse.GetMaxPoissonDiskRadius( 3, outputPoints.size(), volume );
wse.Eliminate( inputPoints.data(), inputPoints.size(), 
               outputPoints.data(), outputPoints.size(),
               isProgressive,
               d_max );

The parameters of the GetMaxPoissonDiskRadius() method are the dimensions, the output sample size, and the volume of the sampling domain. If the volume is given as zero (the default argument value), the volume is computed using the bounds. Otherwise, the bounds are completely ignored.

Sampling Surfaces in 3D

The last argument of the Eliminate() method can be used for sampling surfaces in 3D (or low-dimensional manifolds in a high-dimensional space). First of all, the input points must be generated on the desired surface. Then, while calling the Eliminate() method, we must specify both the correct d_max argument and the dimensionality of the sampling domain. In our example, since we are sampling a surface in 3D, the sampling domain should be considered 2D.

float d_max = 2 * wse.GetMaxPoissonDiskRadius( 2, outputPoints.size(), area );
wse.Eliminate( inputPoints.data(), inputPoints.size(), 
               outputPoints.data(), outputPoints.size(),
               isProgressive,
               d_max, 2 );

In this case, the area variable holds the total surface area of the sampling domain.

Notice that the WeightedSampleElimination object was created as a 3D sampler (working on samples in 3D), but we specify the dimensionality of the sampling domain as 2D, since we would like to sample a surface. That said, if d_max is specified and progressive sampling is not used, the last argument that specifies the dimensionality argument is not used.

Tiling

It is sometimes desirable to generate a sample set and then tile it over a domain. Also, when using statistical analysis using power spectrum, the generated sample set must be tileable. Tiling involves some extra computation; therefore, it is off by default. You can turn tiling on using the SetTiling() method.

Note that if tiling is off, weighted sample elimination would produce more output sample sets near the boundaries of the sampling domain. Therefore, it is a good idea to turn on tiling for box-shaped sampling domains, even if the sampling problem does not involve tiling.

Adaptive Sampling

Sometimes it is desirable to have a different sampling density at different parts of a domain. We can achieve this by adjusting the weight function accordingly. I will not provide more details on this beyond saying that you can use a custom weight function. This custom weight function should return a larger value for the parts of the domain that should be sampled more sparsely and should return a smaller value for the parts of the domain with a denser desired sample distribution. You will also need to provide a custom d_max argument, which should be the value for the most sparsely sampled area.

Additional Parameters

The default weight function uses three parameters: alpha, beta, and gamma. The alpha parameter is the exponent used for the weight function. Its default value is 8. You probably don't need to modify this parameter at all.

The beta parameter specifies the amount of weight limiting. Weight limiting is a method that can provide higher quality sample sets, but aggressive weight limiting can sometimes generate a few samples that are closely placed, which is often undesirable. By default, weight limiting is on and the beta parameter is set to 0.65. The value of the beta parameter should be between zero and one. Larger values provide more aggressive weight limiting and smaller values diminish the impact of weight limiting. Setting beta to zero is the same as turning off weight limiting.

Finally, the gamma parameter is used for determining the limit for weight limiting based on the relative ratio of input and output set sizes.

You can find more details on these parameters in the original publication. For any additional information on weighted sample elimination, please refer to the related project page: http://www.cemyuksel.com/research/sampleelimination/

www.cemyuksel.com