7.4 The Halton Sampler

The underlying goal of the StratifiedSampler is to generate a well-distributed but non-uniform set of sample points, with no two sample points too close together and no excessively large regions of the sample space that have no samples. As Figure 7.18 showed, a jittered pattern does this much better than a random pattern does, although its quality can suffer when samples in adjacent strata happen to be close to the shared boundary of their two strata.

This section introduces the HaltonSampler, which is based on algorithms that directly generate low-discrepancy point sets. Unlike the points generated by the StratifiedSampler, the HaltonSampler not only generates points that are guaranteed to not clump too closely together, but it also generates points that are simultaneously well distributed over all of the dimensions of the sample vector—not just one or two dimensions at a time, as the StratifiedSampler did.

7.4.1 Hammersley and Halton Sequences

The Halton and Hammersley sequences are two closely related low-discrepancy point sets. Both are based on a construction called the radical inverse, which is based on the fact that a positive integer value a can be expressed in a base  b with a sequence of digits d Subscript m Baseline left-parenthesis a right-parenthesis ellipsis d 2 left-parenthesis a right-parenthesis d 1 left-parenthesis a right-parenthesis uniquely determined by

a equals sigma-summation Underscript i equals 1 Overscript m Endscripts d Subscript i Baseline left-parenthesis a right-parenthesis b Superscript i minus 1 Baseline comma
(7.6)

where all digits d Subscript i Baseline left-parenthesis a right-parenthesis are between 0 and b minus 1 .

The radical inverse function normal upper Phi Subscript b in base b converts a nonnegative integer a to a fractional value in left-bracket 0 comma 1 right-parenthesis by reflecting these digits about the radix point:

normal upper Phi Subscript b Baseline left-parenthesis a right-parenthesis equals 0 period d 1 left-parenthesis a right-parenthesis d 2 left-parenthesis a right-parenthesis ellipsis d Subscript m Baseline left-parenthesis a right-parenthesis period
(7.7)

Thus, the contribution of the digit d Subscript i Baseline left-parenthesis a right-parenthesis to the radical inverse is d Subscript i Baseline left-parenthesis a right-parenthesis slash b Superscript i .

One of the simplest low-discrepancy sequences is the van der Corput sequence, which is a 1D sequence given by the radical inverse function in base 2:

x Subscript a Baseline equals normal upper Phi 2 left-parenthesis a right-parenthesis period

Table 7.3 shows the first few values of the van der Corput sequence. Notice how it recursively splits the intervals of the 1D line in half, generating a sample point at the center of each interval.

Table 7.3: The radical inverse normal upper Phi 2 left-parenthesis a right-parenthesis of the first few non-negative integers, computed in base 2. Notice how successive values of normal upper Phi 2 left-parenthesis a right-parenthesis are not close to any of the previous values of normal upper Phi 2 left-parenthesis a right-parenthesis . As more and more values of the sequence are generated, samples are necessarily closer to previous samples, although with a minimum distance that is guaranteed to be reasonably good.

a Base 2 normal upper Phi 2 left-parenthesis a right-parenthesis
0 0 0
1 1 0.1 equals 1 slash 2
2 10 0.01 equals 1 slash 4
3 11 0.11 equals 3 slash 4
4 100 0.001 equals 1 slash 8
5 101 0.101 equals 5 slash 8
vertical-ellipsis

The discrepancy of this sequence is

upper D Subscript upper N Superscript asterisk Baseline left-parenthesis upper P right-parenthesis equals upper O left-parenthesis StartFraction log upper N Over upper N EndFraction right-parenthesis comma

which matches the best discrepancy that has been attained for infinite sequences in n  dimensions,

upper D Subscript upper N Superscript asterisk Baseline left-parenthesis upper P right-parenthesis equals upper O left-parenthesis StartFraction left-parenthesis log upper N right-parenthesis Superscript n Baseline Over upper N EndFraction right-parenthesis period

To generate points in an n -dimensional Halton sequence, we use the radical inverse base b , with a different base for each dimension of the pattern. The bases used must all be relatively prime to each other, so a natural choice is to use the first n prime numbers left-parenthesis p 1 comma ellipsis comma p Subscript n Baseline ):

x Subscript a Baseline equals left-parenthesis normal upper Phi 2 left-parenthesis a right-parenthesis comma normal upper Phi 3 left-parenthesis a right-parenthesis comma normal upper Phi 5 left-parenthesis a right-parenthesis comma ellipsis comma normal upper Phi Subscript p Sub Subscript n Subscript Baseline left-parenthesis a right-parenthesis right-parenthesis period

One of the most useful characteristics of the Halton sequence is that it can be used even if the total number of samples needed isn’t known in advance; all prefixes of the sequence are well distributed, so as additional samples are added to the sequence low discrepancy will be maintained. (However, its distribution is best when the total number of samples is the product of powers of the bases normal upper Pi left-parenthesis p Subscript i Baseline right-parenthesis Superscript k Super Subscript i for exponents k Subscript i .)

The discrepancy of an n -dimensional Halton sequence is

upper D Subscript upper N Superscript asterisk Baseline left-parenthesis x Subscript a Baseline right-parenthesis equals upper O left-parenthesis StartFraction left-parenthesis log upper N right-parenthesis Superscript n Baseline Over upper N EndFraction right-parenthesis comma

which is asymptotically optimal.

If the number of samples upper N is fixed, the Hammersley point set can be used, giving slightly lower discrepancy. Hammersley point sets are defined by

x Subscript a Baseline equals left-parenthesis StartFraction a Over upper N EndFraction comma normal upper Phi Subscript b 1 Baseline left-parenthesis a right-parenthesis comma normal upper Phi Subscript b 2 Baseline left-parenthesis a right-parenthesis comma ellipsis comma normal upper Phi Subscript b Sub Subscript n Subscript Baseline left-parenthesis a right-parenthesis right-parenthesis comma

where upper N is the total number of samples to be taken and as before all of the bases b Subscript i are relatively prime. Figure 7.25(a) shows a plot of the first 216 points of the 2D Halton sequence. Figure 7.25(b) shows the first 256 points of the Hammersley sequence.

Figure 7.25: The First Points of Two Low-Discrepancy Sequences in 2D. (a) Halton (216 points), (b) Hammersley (256 points).

The function RadicalInverse() computes the radical inverse for a given number a using the baseIndexth prime number as the base. The function is implemented using an enormous switch statement, where baseIndex is mapped to the appropriate prime number and then a separate RadicalInverseSpecialized() template function actually computes the radical inverse. (The reason for the curious switch-based structure will be explained shortly.)

<<Low Discrepancy Function Definitions>>= 
Float RadicalInverse(int baseIndex, uint64_t a) { switch (baseIndex) { case 0: <<Compute base-2 radical inverse>> 
return ReverseBits64(a) * 0x1p-64;
case 1: return RadicalInverseSpecialized<3>(a); case 2: return RadicalInverseSpecialized<5>(a); case 3: return RadicalInverseSpecialized<7>(a); <<Remainder of cases for RadicalInverse()>> 
} }

For the base-2 radical inverse, we can take advantage of the fact that numbers in digital computers are already represented in base 2 to compute the radical inverse more efficiently. For a 64-bit value a , we have from Equation (7.6)

a equals sigma-summation Underscript i equals 1 Overscript 64 Endscripts d Subscript i Baseline left-parenthesis a right-parenthesis 2 Superscript i minus 1 Baseline period

First consider the result of reversing the bits of a , still considering it as an integer value, which gives

sigma-summation Underscript i equals 1 Overscript 64 Endscripts d Subscript i Baseline left-parenthesis a right-parenthesis 2 Superscript 64 minus i Baseline period

If we then divide this value by 2 Superscript 64 , we have

sigma-summation Underscript i equals 1 Overscript 64 Endscripts d Subscript i Baseline left-parenthesis a right-parenthesis 2 Superscript negative i Baseline comma

which is normal upper Phi 2 left-parenthesis a right-parenthesis . Thus, the base-2 radical inverse can equivalently be computed with a bit reverse and a power-of-two division.

The bits of an integer quantity can be efficiently reversed with a series of logical bit operations. The first line of the ReverseBits32() function, which reverses the bits of a 32-bit integer, swaps the lower 16 bits with the upper 16 bits of the value. The next line simultaneously swaps the first 8 bits of the result with the second 8 bits and the third 8 bits with the fourth. This process continues until the last line, which swaps adjacent bits. To understand this code, it’s helpful to write out the binary values of the various hexadecimal constants. For example, 0xff00ff00 is 11111111000000001111111100000000 in binary; it’s easy to see that a bitwise OR with this value masks off the first and third 8-bit quantities.

<<Low Discrepancy Inline Functions>>= 
inline uint32_t ReverseBits32(uint32_t n) { n = (n << 16) | (n >> 16); n = ((n & 0x00ff00ff) << 8) | ((n & 0xff00ff00) >> 8); n = ((n & 0x0f0f0f0f) << 4) | ((n & 0xf0f0f0f0) >> 4); n = ((n & 0x33333333) << 2) | ((n & 0xcccccccc) >> 2); n = ((n & 0x55555555) << 1) | ((n & 0xaaaaaaaa) >> 1); return n; }

The bits of a 64-bit value can then be reversed by reversing the two 32-bit components individually and then interchanging them.

<<Low Discrepancy Inline Functions>>+=  
inline uint64_t ReverseBits64(uint64_t n) { uint64_t n0 = ReverseBits32((uint32_t)n); uint64_t n1 = ReverseBits32((uint32_t)(n >> 32)); return (n0 << 32) | n1; }

To compute the base-2 radical inverse, then, we reverse the bits and multiply by 1 slash 2 Superscript 64 , where the hexadecimal floating-point constant 0x1p-64 is used for the value 2 Superscript negative 64 . As explained in Section 3.9.1, implementing a power-of-two division via the corresponding power-of-two multiplication gives the same result with IEEE floating point. (And floating-point multiplication is generally more efficient than floating-point division.)

<<Compute base-2 radical inverse>>= 
return ReverseBits64(a) * 0x1p-64;

For other bases, the RadicalInverseSpecialized() template function computes the radical inverse by computing the digits d Subscript i starting with d 1 and computing a series v Subscript i where v 1 equals d 1 , v 2 equals b d 1 plus d 2 such that

v Subscript n Baseline equals b Superscript n minus 1 Baseline d 1 plus b Superscript n minus 2 Baseline d 2 plus midline-horizontal-ellipsis plus d Subscript n Baseline period

(For example, in base 10, it would convert the value 1234 to 4321.) This value can be found entirely using integer arithmetic, without accumulating any round-off error.

The final value of the radical inverse is then found by converting to floating-point and multiplying by 1 slash b Superscript n , where n is the number of digits in the value, to get the value in Equation (7.7). The term for this multiplication is built up in invBaseN as the digits are processed.

<<Low Discrepancy Static Functions>>= 
template <int base> static Float RadicalInverseSpecialized(uint64_t a) { const Float invBase = (Float)1 / (Float)base; uint64_t reversedDigits = 0; Float invBaseN = 1; while (a) { uint64_t next = a / base; uint64_t digit = a - next * base; reversedDigits = reversedDigits * base + digit; invBaseN *= invBase; a = next; } return std::min(reversedDigits * invBaseN, OneMinusEpsilon); }

A natural question to ask would be why a template function parameterized on the base is used here (rather than, say, a regular function call that took the base as a parameter, which would avoid the generation of a separate code path for each base). The motivation is that integer division is shockingly slow on modern CPUs, and much more efficient approaches are possible for division by a compile-time constant.

For example, integer division of a 32-bit value by 3 can be computed exactly by multiplying this value by 2863311531 to get a 64-bit intermediate and then shifting the result right by 33 bits; these are both fairly efficient operations. (A similar approach can be used for dividing 64-bit values by 3, but the magic constant is much larger; see Warren (2006) for more about these techniques.) Thus, using a template function here allows the compiler to see that the division to compute the value of next in the while loop is actually a division by a constant and gives it a chance to apply this optimization. The code with this optimization runs 5.9 times faster on a 2015-era laptop than an implementation based on integer division instructions.

Another optimization is that we avoid computing a running sum over reversed digits multiplied by the reciprocal base; instead, this multiplication is postponed all the way until the end when the loop terminates. The main issue here is that floating-point and integer units on current processors operate fairly independently from each other. Referencing an integer variable within a floating computation in a tight loop would introduce pipeline bubbles related to the amount of time that is needed to convert and move the values from one unit to the other.

It will be useful to be able to compute the inverse of the radical inverse function; the InverseRadicalInverse() function takes the reversed integer digits in some base, corresponding to reversedDigits in the RadicalInverseSpecialized() template function immediately before being multiplied by the 1 slash b Superscript n factor to convert to a floating-point value in left-bracket 0 comma 1 right-parenthesis . Note that in order to be able to compute the inverse correctly, the total number of digits in the original value must be known: for example, both 1234 and 123400 are converted to 4321 after the integer-only part of the radical inverse algorithm; trailing zeros become leading zeros, which are lost.

<<Low Discrepancy Inline Functions>>+=  
template <int base> inline uint64_t InverseRadicalInverse(uint64_t inverse, int nDigits) { uint64_t index = 0; for (int i = 0; i < nDigits; ++i) { uint64_t digit = inverse % base; inverse /= base; index = index * base + digit; } return index; }

The Hammersley and Halton sequences have the shortcoming that as the base b increases, sample values can exhibit surprisingly regular patterns. This issue can be addressed with scrambled Halton and Hammersley sequences, where a permutation is applied to the digits when computing the radical inverse:

normal upper Psi Subscript b Baseline left-parenthesis a right-parenthesis equals 0 period p left-parenthesis d 1 left-parenthesis a right-parenthesis right-parenthesis p left-parenthesis d 2 left-parenthesis a right-parenthesis right-parenthesis ellipsis p left-parenthesis d Subscript m Baseline left-parenthesis a right-parenthesis right-parenthesis comma
(7.8)

where p is a permutation of the digits left-parenthesis 0 comma 1 comma ellipsis comma b minus 1 right-parenthesis . Note that the same permutation is used for each digit, and the same permutation is used for generating all of the sample points in a given base b . Figure 7.26 shows the effect of scrambling with the Halton sequence.

Figure 7.26: Plot of Halton Sample Values with and without Scrambling. (a) In higher dimensions of the sample vector, projections of sample values start to exhibit regular structure. Here, points from the dimensions left-parenthesis normal upper Phi 29 left-parenthesis a right-parenthesis comma normal upper Phi 31 left-parenthesis a right-parenthesis right-parenthesis are shown. (b) Scrambled sequences, Equation (7.8), break up this structure by randomly permuting the digits of sample indices.

In the following, we will use random permutations, though specific constructions of permutations can give slightly better results; see the “Further Reading” section for more details.

The ComputeRadicalInversePermutations() function computes these random permutation tables. It initializes a single contiguous array for all of the permutations, where the first two values are a permutation of the integers zero and one for b equals 2 , the next three values are a permutation of 0 comma 1 comma 2 for b equals 3 , and so forth for successive prime bases. At entry to the for loop below, p points to the start of the permutation array to initialize for the current prime base.

<<Low Discrepancy Function Definitions>>+=  
std::vector<uint16_t> ComputeRadicalInversePermutations(RNG &rng) { std::vector<uint16_t> perms; <<Allocate space in perms for radical inverse permutations>> 
int permArraySize = 0; for (int i = 0; i < PrimeTableSize; ++i) permArraySize += Primes[i]; perms.resize(permArraySize);
uint16_t *p = &perms[0]; for (int i = 0; i < PrimeTableSize; ++i) { <<Generate random permutation for i th prime base>> 
for (int j = 0; j < Primes[i]; ++j) p[j] = j; Shuffle(p, Primes[i], 1, rng);
p += Primes[i]; } return perms; }

The total size of the permutation array is given by the sum of the prime numbers up to the end of a precomputed table of prime numbers.

<<Allocate space in perms for radical inverse permutations>>= 
int permArraySize = 0; for (int i = 0; i < PrimeTableSize; ++i) permArraySize += Primes[i]; perms.resize(permArraySize);

<<Low Discrepancy Declarations>>= 
static constexpr int PrimeTableSize = 1000; extern const int Primes[PrimeTableSize];

<<Low Discrepancy Data Definitions>>= 
const int Primes[PrimeTableSize] = { 2, 3, 5, 7, 11, <<Subsequent prime numbers>> 
13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919
};

Generating each permutation is easy: we just initialize p to the identity permutation for the current prime length and then randomly shuffle its values.

<<Generate random permutation for i th prime base>>= 
for (int j = 0; j < Primes[i]; ++j) p[j] = j; Shuffle(p, Primes[i], 1, rng);

The ScrambledRadicalInverse() function is essentially the same as RadicalInverse() except that it puts each digit through the permutation table for the given base. See Exercise 7.9.3 for discussion of a more efficient implementation for the base-2 case, following RadicalInverse().

<<Low Discrepancy Function Definitions>>+= 
Float ScrambledRadicalInverse(int baseIndex, uint64_t a, const uint16_t *perm) { switch (baseIndex) { case 0: return ScrambledRadicalInverseSpecialized<2>(perm, a); case 1: return ScrambledRadicalInverseSpecialized<3>(perm, a); case 2: return ScrambledRadicalInverseSpecialized<5>(perm, a); case 3: return ScrambledRadicalInverseSpecialized<7>(perm, a); <<Remainder of cases for ScrambledRadicalInverse()>> 
} }

The implementation below also accounts for a special case that can arise when perm maps the digit 0 to a nonzero value. In this case, the iteration stops prematurely once a reaches 0, incorrectly missing an infinitely long suffix of digits with value perm[0]. Fortunately, this is a geometric series with a simple analytic solution whose value is added in the last line.

<<Low Discrepancy Static Functions>>+= 
template <int base> static Float ScrambledRadicalInverseSpecialized(const uint16_t *perm, uint64_t a) { const Float invBase = (Float)1 / (Float)base; uint64_t reversedDigits = 0; Float invBaseN = 1; while (a) { uint64_t next = a / base; uint64_t digit = a - next * base; reversedDigits = reversedDigits * base + perm[digit]; invBaseN *= invBase; a = next; } return std::min(invBaseN * (reversedDigits + invBase * perm[0] / (1 - invBase)), OneMinusEpsilon); }

7.4.2 Halton Sampler Implementation

The HaltonSampler generates sample vectors using the Halton sequence. Unlike the StratifiedSampler, it is fully deterministic; it uses no pseudo-random numbers in its operation. However, Halton samples can lead to aliasing if the image isn’t sufficiently well sampled. Figure 7.27 compares the results of sampling a checkerboard texture using a Halton-based sampler to using the stratified sampler from the previous section. Note the unpleasant pattern along edges in the foreground and toward the horizon.

Figure 7.27: Comparison of the Stratified Sampler to a Low-Discrepancy Sampler Based on Halton Points on the Image Plane. (1) The jittered stratified sampler with a single sample per pixel and (2) the HaltonSampler sampler with a single sample per pixel. Note that although the Halton pattern is able to reproduce the checker pattern farther toward the horizon than the stratified pattern, there is a regular structure to the error in the low-discrepancy pattern that is visually distracting; it doesn’t turn aliasing into less objectionable noise as well as the jittered approach.

<<HaltonSampler Declarations>>= 
class HaltonSampler : public GlobalSampler { public: <<HaltonSampler Public Methods>> 
HaltonSampler(int nsamp, const Bounds2i &sampleBounds); int64_t GetIndexForSample(int64_t sampleNum) const; Float SampleDimension(int64_t index, int dimension) const; std::unique_ptr<Sampler> Clone(int seed);
private: <<HaltonSampler Private Data>> 
static std::vector<uint16_t> radicalInversePermutations; Point2i baseScales, baseExponents; int sampleStride; int multInverse[2]; mutable Point2i pixelForOffset = Point2i(std::numeric_limits<int>::max(), std::numeric_limits<int>::max()); mutable int64_t offsetForCurrentPixel;
<<HaltonSampler Private Methods>> 
const uint16_t *PermutationForDimension(int dim) const { if (dim >= PrimeTableSize) Severe("HaltonSampler can only sample %d dimensions.", PrimeTableSize); return &radicalInversePermutations[PrimeSums[dim]]; }
};

<<HaltonSampler Method Definitions>>= 
HaltonSampler::HaltonSampler(int samplesPerPixel, const Bounds2i &sampleBounds) : GlobalSampler(samplesPerPixel) { <<Generate random digit permutations for Halton sampler>>  <<Find radical inverse base scales and exponents that cover sampling area>> 
Vector2i res = sampleBounds.pMax - sampleBounds.pMin; for (int i = 0; i < 2; ++i) { int base = (i == 0) ? 2 : 3; int scale = 1, exp = 0; while (scale < std::min(res[i], kMaxResolution)) { scale *= base; ++exp; } baseScales[i] = scale; baseExponents[i] = exp; }
<<Compute stride in samples for visiting each pixel area>>  <<Compute multiplicative inverses for baseScales>> 
multInverse[0] = multiplicativeInverse(baseScales[1], baseScales[0]); multInverse[1] = multiplicativeInverse(baseScales[0], baseScales[1]);
}

The permutation tables for the scrambled radical inverses are shared across all HaltonSampler instances and are computed the first time the constructor runs. For pbrt’s requirements, this approach is fine: the current implementation only uses different sampler instances for different tiles of the image, where we’d like to always use the same permutations anyway. For other uses, it could be worthwhile to have more control over when different permutations are used.

<<Generate random digit permutations for Halton sampler>>= 

<<HaltonSampler Private Data>>= 
static std::vector<uint16_t> radicalInversePermutations;

The utility method PermutationForDimension() returns a pointer to the start of the permutation array for the given dimension.

<<HaltonSampler Private Methods>>= 
const uint16_t *PermutationForDimension(int dim) const { if (dim >= PrimeTableSize) Severe("HaltonSampler can only sample %d dimensions.", PrimeTableSize); return &radicalInversePermutations[PrimeSums[dim]]; }

To be able to quickly find the offset for a given dimension, it’s helpful to have the sums of the prime numbers preceding each prime.

<<Low Discrepancy Data Definitions>>+= 
const int PrimeSums[PrimeTableSize] = { 0, 2, 5, 10, 17, <<Subsequent prime sums>> 
28, 41, 58, 77, 100, 129, 160, 197, 238, 281, 328, 381, 440, 501, 568, 639, 712, 791, 874, 963, 1060, 1161, 1264, 1371, 1480, 1593, 1720, 1851, 1988, 2127, 2276, 2427, 2584, 2747, 2914, 3087, 3266, 3447, 3638, 3831, 4028, 4227, 4438, 4661, 4888, 5117, 5350, 5589, 5830, 6081, 6338, 6601, 6870, 7141, 7418, 7699, 7982, 8275, 8582, 8893, 9206, 9523, 9854, 10191, 10538, 10887, 11240, 11599, 11966, 12339, 12718, 13101, 13490, 13887, 14288, 14697, 15116, 15537, 15968, 16401, 16840, 17283, 17732, 18189, 18650, 19113, 19580, 20059, 20546, 21037, 21536, 22039, 22548, 23069, 23592, 24133, 24680, 25237, 25800, 26369, 26940, 27517, 28104, 28697, 29296, 29897, 30504, 31117, 31734, 32353, 32984, 33625, 34268, 34915, 35568, 36227, 36888, 37561, 38238, 38921, 39612, 40313, 41022, 41741, 42468, 43201, 43940, 44683, 45434, 46191, 46952, 47721, 48494, 49281, 50078, 50887, 51698, 52519, 53342, 54169, 54998, 55837, 56690, 57547, 58406, 59269, 60146, 61027, 61910, 62797, 63704, 64615, 65534, 66463, 67400, 68341, 69288, 70241, 71208, 72179, 73156, 74139, 75130, 76127, 77136, 78149, 79168, 80189, 81220, 82253, 83292, 84341, 85392, 86453, 87516, 88585, 89672, 90763, 91856, 92953, 94056, 95165, 96282, 97405, 98534, 99685, 100838, 102001, 103172, 104353, 105540, 106733, 107934, 109147, 110364, 111587, 112816, 114047, 115284, 116533, 117792, 119069, 120348, 121631, 122920, 124211, 125508, 126809, 128112, 129419, 130738, 132059, 133386, 134747, 136114, 137487, 138868, 140267, 141676, 143099, 144526, 145955, 147388, 148827, 150274, 151725, 153178, 154637, 156108, 157589, 159072, 160559, 162048, 163541, 165040, 166551, 168074, 169605, 171148, 172697, 174250, 175809, 177376, 178947, 180526, 182109, 183706, 185307, 186914, 188523, 190136, 191755, 193376, 195003, 196640, 198297, 199960, 201627, 203296, 204989, 206686, 208385, 210094, 211815, 213538, 215271, 217012, 218759, 220512, 222271, 224048, 225831, 227618, 229407, 231208, 233019, 234842, 236673, 238520, 240381, 242248, 244119, 245992, 247869, 249748, 251637, 253538, 255445, 257358, 259289, 261222, 263171, 265122, 267095, 269074, 271061, 273054, 275051, 277050, 279053, 281064, 283081, 285108, 287137, 289176, 291229, 293292, 295361, 297442, 299525, 301612, 303701, 305800, 307911, 310024, 312153, 314284, 316421, 318562, 320705, 322858, 325019, 327198, 329401, 331608, 333821, 336042, 338279, 340518, 342761, 345012, 347279, 349548, 351821, 354102, 356389, 358682, 360979, 363288, 365599, 367932, 370271, 372612, 374959, 377310, 379667, 382038, 384415, 386796, 389179, 391568, 393961, 396360, 398771, 401188, 403611, 406048, 408489, 410936, 413395, 415862, 418335, 420812, 423315, 425836, 428367, 430906, 433449, 435998, 438549, 441106, 443685, 446276, 448869, 451478, 454095, 456716, 459349, 461996, 464653, 467312, 469975, 472646, 475323, 478006, 480693, 483382, 486075, 488774, 491481, 494192, 496905, 499624, 502353, 505084, 507825, 510574, 513327, 516094, 518871, 521660, 524451, 527248, 530049, 532852, 535671, 538504, 541341, 544184, 547035, 549892, 552753, 555632, 558519, 561416, 564319, 567228, 570145, 573072, 576011, 578964, 581921, 584884, 587853, 590824, 593823, 596824, 599835, 602854, 605877, 608914, 611955, 615004, 618065, 621132, 624211, 627294, 630383, 633492, 636611, 639732, 642869, 646032, 649199, 652368, 655549, 658736, 661927, 665130, 668339, 671556, 674777, 678006, 681257, 684510, 687767, 691026, 694297, 697596, 700897, 704204, 707517, 710836, 714159, 717488, 720819, 724162, 727509, 730868, 734229, 737600, 740973, 744362, 747753, 751160, 754573, 758006, 761455, 764912, 768373, 771836, 775303, 778772, 782263, 785762, 789273, 792790, 796317, 799846, 803379, 806918, 810459, 814006, 817563, 821122, 824693, 828274, 831857, 835450, 839057, 842670, 846287, 849910, 853541, 857178, 860821, 864480, 868151, 871824, 875501, 879192, 882889, 886590, 890299, 894018, 897745, 901478, 905217, 908978, 912745, 916514, 920293, 924086, 927883, 931686, 935507, 939330, 943163, 947010, 950861, 954714, 958577, 962454, 966335, 970224, 974131, 978042, 981959, 985878, 989801, 993730, 997661, 1001604, 1005551, 1009518, 1013507, 1017508, 1021511, 1025518, 1029531, 1033550, 1037571, 1041598, 1045647, 1049698, 1053755, 1057828, 1061907, 1065998, 1070091, 1074190, 1078301, 1082428, 1086557, 1090690, 1094829, 1098982, 1103139, 1107298, 1111475, 1115676, 1119887, 1124104, 1128323, 1132552, 1136783, 1141024, 1145267, 1149520, 1153779, 1158040, 1162311, 1166584, 1170867, 1175156, 1179453, 1183780, 1188117, 1192456, 1196805, 1201162, 1205525, 1209898, 1214289, 1218686, 1223095, 1227516, 1231939, 1236380, 1240827, 1245278, 1249735, 1254198, 1258679, 1263162, 1267655, 1272162, 1276675, 1281192, 1285711, 1290234, 1294781, 1299330, 1303891, 1308458, 1313041, 1317632, 1322229, 1326832, 1331453, 1336090, 1340729, 1345372, 1350021, 1354672, 1359329, 1363992, 1368665, 1373344, 1378035, 1382738, 1387459, 1392182, 1396911, 1401644, 1406395, 1411154, 1415937, 1420724, 1425513, 1430306, 1435105, 1439906, 1444719, 1449536, 1454367, 1459228, 1464099, 1468976, 1473865, 1478768, 1483677, 1488596, 1493527, 1498460, 1503397, 1508340, 1513291, 1518248, 1523215, 1528184, 1533157, 1538144, 1543137, 1548136, 1553139, 1558148, 1563159, 1568180, 1573203, 1578242, 1583293, 1588352, 1593429, 1598510, 1603597, 1608696, 1613797, 1618904, 1624017, 1629136, 1634283, 1639436, 1644603, 1649774, 1654953, 1660142, 1665339, 1670548, 1675775, 1681006, 1686239, 1691476, 1696737, 1702010, 1707289, 1712570, 1717867, 1723170, 1728479, 1733802, 1739135, 1744482, 1749833, 1755214, 1760601, 1765994, 1771393, 1776800, 1782213, 1787630, 1793049, 1798480, 1803917, 1809358, 1814801, 1820250, 1825721, 1831198, 1836677, 1842160, 1847661, 1853164, 1858671, 1864190, 1869711, 1875238, 1880769, 1886326, 1891889, 1897458, 1903031, 1908612, 1914203, 1919826, 1925465, 1931106, 1936753, 1942404, 1948057, 1953714, 1959373, 1965042, 1970725, 1976414, 1982107, 1987808, 1993519, 1999236, 2004973, 2010714, 2016457, 2022206, 2027985, 2033768, 2039559, 2045360, 2051167, 2056980, 2062801, 2068628, 2074467, 2080310, 2086159, 2092010, 2097867, 2103728, 2109595, 2115464, 2121343, 2127224, 2133121, 2139024, 2144947, 2150874, 2156813, 2162766, 2168747, 2174734, 2180741, 2186752, 2192781, 2198818, 2204861, 2210908, 2216961, 2223028, 2229101, 2235180, 2241269, 2247360, 2253461, 2259574, 2265695, 2271826, 2277959, 2284102, 2290253, 2296416, 2302589, 2308786, 2314985, 2321188, 2327399, 2333616, 2339837, 2346066, 2352313, 2358570, 2364833, 2371102, 2377373, 2383650, 2389937, 2396236, 2402537, 2408848, 2415165, 2421488, 2427817, 2434154, 2440497, 2446850, 2453209, 2459570, 2465937, 2472310, 2478689, 2485078, 2491475, 2497896, 2504323, 2510772, 2517223, 2523692, 2530165, 2536646, 2543137, 2549658, 2556187, 2562734, 2569285, 2575838, 2582401, 2588970, 2595541, 2602118, 2608699, 2615298, 2621905, 2628524, 2635161, 2641814, 2648473, 2655134, 2661807, 2668486, 2675175, 2681866, 2688567, 2695270, 2701979, 2708698, 2715431, 2722168, 2728929, 2735692, 2742471, 2749252, 2756043, 2762836, 2769639, 2776462, 2783289, 2790118, 2796951, 2803792, 2810649, 2817512, 2824381, 2831252, 2838135, 2845034, 2851941, 2858852, 2865769, 2872716, 2879665, 2886624, 2893585, 2900552, 2907523, 2914500, 2921483, 2928474, 2935471, 2942472, 2949485, 2956504, 2963531, 2970570, 2977613, 2984670, 2991739, 2998818, 3005921, 3013030, 3020151, 3027278, 3034407, 3041558, 3048717, 3055894, 3063081, 3070274, 3077481, 3084692, 3091905, 3099124, 3106353, 3113590, 3120833, 3128080, 3135333, 3142616, 3149913, 3157220, 3164529, 3171850, 3179181, 3186514, 3193863, 3201214, 3208583, 3215976, 3223387, 3230804, 3238237, 3245688, 3253145, 3260604, 3268081, 3275562, 3283049, 3290538, 3298037, 3305544, 3313061, 3320584, 3328113, 3335650, 3343191, 3350738, 3358287, 3365846, 3373407, 3380980, 3388557, 3396140, 3403729, 3411320, 3418923, 3426530, 3434151, 3441790, 3449433, 3457082, 3464751, 3472424, 3480105, 3487792, 3495483, 3503182, 3510885, 3518602, 3526325, 3534052, 3541793, 3549546, 3557303, 3565062, 3572851, 3580644, 3588461, 3596284, 3604113, 3611954, 3619807, 3627674, 3635547, 3643424, 3651303, 3659186, 3667087, 3674994,
};

To map the first two dimensions of samples from left-bracket 0 comma 1 right-parenthesis squared to pixel coordinates, the HaltonSampler finds the smallest scale factor left-parenthesis 2 Superscript j Baseline comma 3 Superscript k Baseline right-parenthesis that is larger than the lower of either the image resolution or kMaxResolution in each dimension. (We will see shortly how this specific choice of scales makes it easy to see which pixel a sample lands in.) After scaling, any samples outside the image extent will be simply ignored.

For images with resolution greater than kMaxResolution in one or both dimensions, a tile of Halton points is repeated across the image. This resolution limit helps maintain sufficient floating-point precision in the computed sample values.

<<Find radical inverse base scales and exponents that cover sampling area>>= 
Vector2i res = sampleBounds.pMax - sampleBounds.pMin; for (int i = 0; i < 2; ++i) { int base = (i == 0) ? 2 : 3; int scale = 1, exp = 0; while (scale < std::min(res[i], kMaxResolution)) { scale *= base; ++exp; } baseScales[i] = scale; baseExponents[i] = exp; }

For each dimension, baseScales holds the scale factor, 2 Superscript j or 3 Superscript k , and baseExponents holds the exponents j and k .

<<HaltonSampler Private Data>>+=  
Point2i baseScales, baseExponents;

<<HaltonSampler Local Constants>>= 
static constexpr int kMaxResolution = 128;

To see why the HaltonSampler uses this scheme to map samples to pixel coordinates, consider the effect of scaling a value computed with the radical inverse base b by a factor b Superscript n . If the digits of a expressed in base b are d Subscript i Baseline left-parenthesis a right-parenthesis , then recall that the radical inverse is the value 0 period d 1 left-parenthesis a right-parenthesis d 2 left-parenthesis a right-parenthesis ellipsis comma base b . If we multiply this value by b squared , for example, we have d 1 left-parenthesis a right-parenthesis d 2 left-parenthesis a right-parenthesis period d 3 left-parenthesis a right-parenthesis ellipsis semicolon the first two digits have moved to the left of the radix point, and the fractional component of the value starts with d 3 left-parenthesis a right-parenthesis .

This operation—scaling by b Superscript n —forms the core of being able to determine which sample indices land in which pixels. Considering the first two digits in the above example, we can see that the integer component of the scaled value ranges from 0 to b squared minus 1 and that as a increases, its last two digits in base b take on any particular value once in every b squared values in this range.

Given a value x , 0 less-than-or-equal-to x less-than-or-equal-to b squared minus 1 , we can find the first value a that gives the value x in the integer components. By definition, the digits of x in base b are d 2 left-parenthesis x right-parenthesis d 1 left-parenthesis x right-parenthesis . Thus, if d 1 left-parenthesis a right-parenthesis equals d 2 left-parenthesis x right-parenthesis and d 2 left-parenthesis a right-parenthesis equals d 1 left-parenthesis x right-parenthesis , then the scaled value of a ’s radical inverse will have an integer component equal to x .

Because the bases b equals 2 and b equals 3 used in the HaltonSampler for pixel samples are relatively prime, it follows that if the sample values are scaled by some left-parenthesis 2 Superscript j Baseline comma 3 Superscript k Baseline right-parenthesis , then any particular pixel in the range left-parenthesis 0 comma 0 right-parenthesis right-arrow left-parenthesis 2 Superscript j Baseline minus 1 comma 3 Superscript k Baseline minus 1 right-parenthesis will be visited once every 2 Superscript j Baseline 3 Superscript k samples. This product is stored in sampleStride.

<<Compute stride in samples for visiting each pixel area>>= 

<<HaltonSampler Private Data>>+=  
int sampleStride;

<<HaltonSampler Private Data>>+=  
int multInverse[2];

The sample index for the first Halton sample that lands in currentPixel is stored in offsetForCurrentPixel. After this offset has first been computed for the first sample in the current pixel, subsequent samples in the pixel are found at increments of sampleStride samples in the Halton sequence.

<<HaltonSampler Method Definitions>>+=  
int64_t HaltonSampler::GetIndexForSample(int64_t sampleNum) const { if (currentPixel != pixelForOffset) { <<Compute Halton sample offset for currentPixel>> 
offsetForCurrentPixel = 0; if (sampleStride > 1) { Point2i pm(Mod(currentPixel[0], kMaxResolution), Mod(currentPixel[1], kMaxResolution)); for (int i = 0; i < 2; ++i) { uint64_t dimOffset = (i == 0) ? InverseRadicalInverse<2>(pm[i], baseExponents[i]) : InverseRadicalInverse<3>(pm[i], baseExponents[i]); offsetForCurrentPixel += dimOffset * (sampleStride / baseScales[i]) * multInverse[i]; } offsetForCurrentPixel %= sampleStride; }
pixelForOffset = currentPixel; } return offsetForCurrentPixel + sampleNum * sampleStride; }

<<HaltonSampler Private Data>>+= 
mutable Point2i pixelForOffset = Point2i(std::numeric_limits<int>::max(), std::numeric_limits<int>::max()); mutable int64_t offsetForCurrentPixel;

Computing the index of the first sample in a given pixel left-parenthesis x comma y right-parenthesis where the samples have been scaled by left-parenthesis 2 Superscript j Baseline comma 3 Superscript k Baseline right-parenthesis involves computing the inverse radical inverse of the last j digits of x in base 2, which we’ll denote by x Subscript r , and of the last k digits of y in base 3, y Subscript r . This gives us a system of equations

StartLayout 1st Row 1st Column x Subscript r 2nd Column identical-to left-parenthesis i normal m normal o normal d 2 Superscript j Baseline right-parenthesis 2nd Row 1st Column y Subscript r 2nd Column identical-to left-parenthesis i normal m normal o normal d 3 Superscript k Baseline right-parenthesis comma EndLayout

where the index i that satisfies these equations is the index of a sample that lies within the given pixel, after scaling. We don’t include the code that solves for i <<Compute Halton sample offset for currentPixel>> here in the book; see Grünschloß et al. (2012) for details of the algorithm used to find i .

The computation of sample offsets doesn’t account for random digit permutations, so those aren’t included in the sample values computed here. Also, because the low baseExponents[i] digits of the first two dimensions are used to select which pixel is sampled, these digits must be discarded before computing the radical inverse for the first two dimensions of the sample vector, since the SampleDimension() method is supposed to return the fractional offset within the pixel being sampled. Higher dimensions are just sampled directly, including the random permutations.

<<HaltonSampler Method Definitions>>+= 
Float HaltonSampler::SampleDimension(int64_t index, int dim) const { if (dim == 0) return RadicalInverse(dim, index >> baseExponents[0]); else if (dim == 1) return RadicalInverse(dim, index / baseScales[1]); else return ScrambledRadicalInverse(dim, index, PermutationForDimension(dim)); }