7.7 Sobol’ Sampler

The last Sampler in this chapter is based on a series of generator matrices due to Sobol prime . The samples from the sequence that these matrices generate are distinguished by both being very efficient to implement—thanks to being entirely based on base-2 computations—while also being extremely well distributed over all n dimensions of the sample vector. Figure 7.34 shows the first few Sobol prime generator matrices.

Figure 7.34: Generator matrices for the first four dimensions of the Sobol prime sequence. Note their regular structure.

Figure 7.35 compares Sobol prime samples to stratified and Halton points with the depth of field test scene.

Figure 7.35: Comparisons of the Stratified, Halton, and Sobol prime Samplers for Rendering Depth of Field. (1) An image rendered using the StratifiedSampler, (2) an image rendered using the HaltonSampler, and (3) an image using the SobolSampler. Both low-discrepancy samplers are better than the stratified sampler. In spite of the structured grid artifacts visible with this undersampled image with the SobolSampler, the Sobol prime sequence often provides a faster rate of convergence than the Halton sequence.

The weakness of the Sobol prime points is that they are prone to structural grid artifacts before convergence; a sense of this issue can be seen in the image sample points shown in Figure 7.36.

Figure 7.36: A grid of 2 times 2 pixels, sampled with 16 Sobol prime samples each. Note that there is a fair amount of structure as well as many samples close to others. The very good distribution properties of the sequence over all n dimensions of the sample vector generally make up for these shortcomings.

This structure is visible in the images in Figure 7.37. In exchange for this weakness, Sobol prime sequences are extremely well distributed over all n dimensions of the sample vector.

Figure 7.37: Undersampled images rendered with (1) the Halton sampler and (2) the Sobol prime sampler. Both exhibit visible structure, though with different visual characteristics. The Sobol prime sequence in particular exhibits a clearly visible checkerboard structure.

<<SobolSampler Declarations>>= 
class SobolSampler : public GlobalSampler { public: <<SobolSampler Public Methods>> 
std::unique_ptr<Sampler> Clone(int seed); SobolSampler(int64_t samplesPerPixel, const Bounds2i &sampleBounds) : GlobalSampler(RoundUpPow2(samplesPerPixel)), sampleBounds(sampleBounds) { resolution = RoundUpPow2(std::max(sampleBounds.Diagonal().x, sampleBounds.Diagonal().y)); log2Resolution = Log2Int(resolution); } int64_t GetIndexForSample(int64_t sampleNum) const; Float SampleDimension(int64_t index, int dimension) const;
private: <<SobolSampler Private Data>> 
const Bounds2i sampleBounds; int resolution, log2Resolution;
};

The SobolSampler uniformly scales the first two dimensions by the smallest power of 2 that causes the left-bracket 0 comma 1 right-parenthesis squared sample domain to cover the image area to be sampled. As with the HaltonSampler, this specific scaling scheme is chosen in order to make it easier to compute the reverse mapping from pixel coordinates to the sample indices that land in each pixel.

<<SobolSampler Public Methods>>= 
SobolSampler(int64_t samplesPerPixel, const Bounds2i &sampleBounds) : GlobalSampler(RoundUpPow2(samplesPerPixel)), sampleBounds(sampleBounds) { resolution = RoundUpPow2(std::max(sampleBounds.Diagonal().x, sampleBounds.Diagonal().y)); log2Resolution = Log2Int(resolution); }

<<SobolSampler Private Data>>= 
const Bounds2i sampleBounds; int resolution, log2Resolution;

The SobolIntervalToIndex() function returns the index of the sampleNumth sample in the pixel p, if the left-bracket 0 comma 1 right-parenthesis squared sampling domain has been scaled by 2 Superscript normal l normal o normal g Baseline 2 normal upper R normal e normal s normal o normal l normal u normal t normal i normal o normal n to cover the pixel sampling area.

<<Low Discrepancy Declarations>>+= 
inline uint64_t SobolIntervalToIndex(const uint32_t log2Resolution, uint64_t sampleNum, const Point2i &p);

The general approach used to derive the algorithm it implements is similar to that used by the Halton sampler in its GetIndexForSample() method. Here, scaling by a power of two means that the base-2 logarithm of the scale gives the number of digits of the bold upper C left-bracket d Subscript i Baseline left-parenthesis a right-parenthesis right-bracket Superscript upper T product that form the scaled sample’s integer component. To find the values of a that give a particular integer value after scaling, we can compute the inverse of bold upper C : given

v equals bold upper C left-bracket d Subscript i Baseline left-parenthesis a right-parenthesis right-bracket Superscript upper T Baseline comma

then equivalently

bold upper C Superscript negative 1 Baseline v equals left-bracket d Subscript i Baseline left-parenthesis a right-parenthesis right-bracket Superscript upper T Baseline period

We won’t include the implementation of this method here.

<<SobolSampler Method Definitions>>= 
int64_t SobolSampler::GetIndexForSample(int64_t sampleNum) const { return SobolIntervalToIndex(log2Resolution, sampleNum, Point2i(currentPixel - sampleBounds.pMin)); }

Computing the sample value for a given sample index and dimension is straightforward given the SobolSample() function.

<<SobolSampler Method Definitions>>+= 
Float SobolSampler::SampleDimension(int64_t index, int dim) const { Float s = SobolSample(index, dim); <<Remap Sobol prime dimensions used for pixel samples>> 
if (dim == 0 || dim == 1) { s = s * resolution + sampleBounds.pMin[dim]; s = Clamp(s - currentPixel[dim], (Float)0, OneMinusEpsilon); }
return s; }

The code for computing Sobol prime sample values takes different paths for 32- and 64-bit floating-point values. Different generator matrices are used for these two cases, giving more bits of precision for 64-bit doubles.

<<Low Discrepancy Inline Functions>>+=  
inline Float SobolSample(int64_t index, int dimension, uint64_t scramble = 0) { #ifdef PBRT_FLOAT_AS_DOUBLE return SobolSampleDouble(index, dimension, scramble); #else return SobolSampleFloat(index, dimension, scramble); #endif }

The implementation of the SobolSampleFloat() function is quite similar to that of MultiplyGenerator(), with the differences that it takes a 64-bit index and that the matrices it uses have size 32 times 52 . These larger matrices allow it to generate distinct sample values up to a equals 2 Superscript 52 Baseline minus 1 , rather than 2 Superscript 32 Baseline minus 1 , as with the 32 times 32 matrices used previously.

<<Low Discrepancy Inline Functions>>+= 
inline float SobolSampleFloat(int64_t a, int dimension, uint32_t scramble) { uint32_t v = scramble; for (int i = dimension * SobolMatrixSize; a != 0; a >>= 1, ++i) if (a & 1) v ^= SobolMatrices32[i]; return v * 0x1p-32f; /* 1/2^32 */ }

<<Sobol Matrix Declarations>>= 
static constexpr int NumSobolDimensions = 1024; static constexpr int SobolMatrixSize = 52; extern const uint32_t SobolMatrices32[NumSobolDimensions * SobolMatrixSize];

The SobolSampleDouble() function is similar, except that it uses 64-bit Sobol prime matrices. It is not included in the text here.

Because the SobolSampler is a GlobalSampler, the values returned for the first two dimensions need to be adjusted so that they are offsets from the current pixel. Here, the sample value is scaled up by the power-of-two scale computed in the constructor and then offset by the lower corner of the sample bounds to find the corresponding raster sample location. The current integer pixel coordinate is subtracted to get a result in left-bracket 0 comma 1 right-parenthesis .

<<Remap Sobol prime dimensions used for pixel samples>>= 
if (dim == 0 || dim == 1) { s = s * resolution + sampleBounds.pMin[dim]; s = Clamp(s - currentPixel[dim], (Float)0, OneMinusEpsilon); }