Update from book source.

Removed support for double-precision Sobol' samples.

src/pbrt/cmd/pspec.cpp

@@ -195,14 +195,9 @@ GenerateSamples(std::string samplerName, int nPoints, int iter) {
         RNG rng(Options->seed, iter);
         uint32_t r[2] = {rng.Uniform<uint32_t>(), rng.Uniform<uint32_t>()};
 
-        for (int i = 0; i < nPoints; ++i) {
-            Float u2 = SobolSample(i, 0, NoRandomizer());
-            uint32_t uu = OwenScrambleBinaryFull(uint32_t(u2 * 0x1p32), r[0]);
-            u2 = uu * 0x1p-32;
-
-            points.push_back(
-                Point2f(u2, OwenScrambledRadicalInverse(1, i, r[1])));
-        }
+        for (int i = 0; i < nPoints; ++i)
+            points.push_back(Point2f(SobolSample(i, 0, OwenScrambler(r[0])),
+                                     OwenScrambledRadicalInverse(1, i, r[1])));
     } else {
         Sampler sampler = [&]() -> Sampler {
             if (samplerName == "random")

src/pbrt/samplers.h

@@ -76,12 +76,6 @@ class HaltonSampler {
         return SampleDimension(dimension++);
     }
 
-    PBRT_CPU_GPU
-    Point2f GetPixel2D() {
-        return {RadicalInverse(0, haltonIndex >> baseExponents[0]),
-                RadicalInverse(1, haltonIndex / baseScales[1])};
-    }
-
     PBRT_CPU_GPU
     Point2f Get2D() {
         if (dimension + 1 >= PrimeTableSize)
@@ -91,6 +85,12 @@ class HaltonSampler {
         return {SampleDimension(dim), SampleDimension(dim + 1)};
     }
 
+    PBRT_CPU_GPU
+    Point2f GetPixel2D() {
+        return {RadicalInverse(0, haltonIndex >> baseExponents[0]),
+                RadicalInverse(1, haltonIndex / baseScales[1])};
+    }
+
     std::vector<Sampler> Clone(int n, Allocator alloc);
     std::string ToString() const;
 
@@ -116,18 +116,15 @@ class HaltonSampler {
 
     PBRT_CPU_GPU
     Float SampleDimension(int dimension) const {
-        switch (randomize) {
-        case RandomizeStrategy::None:
+        if (randomize == RandomizeStrategy::None)
             return RadicalInverse(dimension, haltonIndex);
-        case RandomizeStrategy::PermuteDigits:
+        else if (randomize == RandomizeStrategy::PermuteDigits)
             return ScrambledRadicalInverse(dimension, haltonIndex,
                                            (*digitPermutations)[dimension]);
-        case RandomizeStrategy::Owen:
+        else {
+            DCHECK_EQ(randomize, RandomizeStrategy::Owen);
             return OwenScrambledRadicalInverse(dimension, haltonIndex,
-                                               MixBits(1 + (uint64_t(dimension) << 32)));
-        default:
-            LOG_FATAL("Unhandled randomization strategy");
-            return {};
+                                               MixBits(1 + (dimension << 4)));
         }
     }
 
@@ -205,19 +202,14 @@ class PaddedSobolSampler {
     // PaddedSobolSampler Private Methods
     PBRT_CPU_GPU
     Float SampleDimension(int dimension, uint32_t a, uint32_t hash) const {
-        switch (randomize) {
-        case RandomizeStrategy::None:
+        if (randomize == RandomizeStrategy::None)
             return SobolSample(a, dimension, NoRandomizer());
-        case RandomizeStrategy::PermuteDigits:
+        else if (randomize == RandomizeStrategy::PermuteDigits)
             return SobolSample(a, dimension, BinaryPermuteScrambler(hash));
-        case RandomizeStrategy::FastOwen:
+        else if (randomize == RandomizeStrategy::FastOwen)
             return SobolSample(a, dimension, FastOwenScrambler(hash));
-        case RandomizeStrategy::Owen:
+        else
             return SobolSample(a, dimension, OwenScrambler(hash));
-        default:
-            LOG_FATAL("Unhandled randomization strategy");
-            return {};
-        }
     }
 
     // PaddedSobolSampler Private Members
@@ -230,6 +222,7 @@ class PaddedSobolSampler {
 // ZSobolSampler Definition
 class ZSobolSampler {
   public:
+    // ZSobolSampler Public Methods
     ZSobolSampler(int samplesPerPixel, Point2i fullResolution,
                   RandomizeStrategy randomize = RandomizeStrategy::PermuteDigits,
                   int seed = 0)
@@ -315,36 +308,62 @@ class ZSobolSampler {
     PBRT_CPU_GPU
     uint64_t GetSampleIndex() const {
         static const uint8_t permutations[24][4] = {
-            {0, 1, 2, 3}, {0, 1, 3, 2}, {0, 2, 1, 3}, {0, 2, 3, 1}, {0, 3, 2, 1},
-            {0, 3, 1, 2}, {1, 0, 2, 3}, {1, 0, 3, 2}, {1, 2, 0, 3}, {1, 2, 3, 0},
-            {1, 3, 2, 0}, {1, 3, 0, 2}, {2, 1, 0, 3}, {2, 1, 3, 0}, {2, 0, 1, 3},
-            {2, 0, 3, 1}, {2, 3, 0, 1}, {2, 3, 1, 0}, {3, 1, 2, 0}, {3, 1, 0, 2},
-            {3, 2, 1, 0}, {3, 2, 0, 1}, {3, 0, 2, 1}, {3, 0, 1, 2}};
-
+            {0, 1, 2, 3},
+            {0, 1, 3, 2},
+            {0, 2, 1, 3},
+            {0, 2, 3, 1},
+            // Define remaining 20 4-way permutations
+            {0, 3, 2, 1},
+            {0, 3, 1, 2},
+            {1, 0, 2, 3},
+            {1, 0, 3, 2},
+            {1, 2, 0, 3},
+            {1, 2, 3, 0},
+            {1, 3, 2, 0},
+            {1, 3, 0, 2},
+            {2, 1, 0, 3},
+            {2, 1, 3, 0},
+            {2, 0, 1, 3},
+            {2, 0, 3, 1},
+            {2, 3, 0, 1},
+            {2, 3, 1, 0},
+            {3, 1, 2, 0},
+            {3, 1, 0, 2},
+            {3, 2, 1, 0},
+            {3, 2, 0, 1},
+            {3, 0, 2, 1},
+            {3, 0, 1, 2}
+
+        };
         uint64_t sampleIndex = 0;
         bool pow2Samples = log2SamplesPerPixel & 1;
         int lastDigit = pow2Samples ? 1 : 0;
         for (int i = nBase4Digits - 1; i >= lastDigit; --i) {
+            // Randomly permute $i$th base 4 digit in _mortonIndex_
             int digitShift = 2 * i;
             int digit = (mortonIndex >> digitShift) & 3;
             int p = HashPerm(mortonIndex >> (digitShift + 2));
             digit = permutations[p][digit];
             sampleIndex |= uint64_t(digit) << digitShift;
         }
+        // Handle power-of-2 (but not 4) sample count
         if (pow2Samples) {
             sampleIndex |= (mortonIndex & 3);
             sampleIndex >>= 1;
             sampleIndex ^= MixBits((mortonIndex >> 2) ^ (0x55555555 * dimension)) & 1;
         }
+
         return sampleIndex;
     }
 
   private:
+    // ZSobolSampler Private Methods
     PBRT_CPU_GPU
     int HashPerm(uint64_t index) const {
         return uint32_t(MixBits(index ^ (0x55555555 * dimension)) >> 24) % 24;
     }
 
+    // ZSobolSampler Private Members
     RandomizeStrategy randomize;
     int log2SamplesPerPixel, seed, nBase4Digits;
     uint64_t mortonIndex;
@@ -503,28 +522,29 @@ class SobolSampler {
         return SampleDimension(dimension++);
     }
 
+    PBRT_CPU_GPU
+    Point2f Get2D() {
+        if (dimension + 1 >= NSobolDimensions)
+            dimension = 2;
+        Point2f u(SampleDimension(dimension), SampleDimension(dimension + 1));
+        dimension += 2;
+        return u;
+    }
+
     PBRT_CPU_GPU
     Point2f GetPixel2D() {
-        Point2f u(SampleDimension(0), SampleDimension(1));
+        Point2f u(SobolSample(sobolIndex, 0, NoRandomizer()),
+                  SobolSample(sobolIndex, 1, NoRandomizer()));
         // Remap Sobol$'$ dimensions used for pixel samples
         for (int dim = 0; dim < 2; ++dim) {
-            CHECK_RARE(1e-7, u[dim] * scale - pixel[dim] < 0);
-            CHECK_RARE(1e-7, u[dim] * scale - pixel[dim] > 1);
+            DCHECK_RARE(1e-7, u[dim] * scale - pixel[dim] < 0);
+            DCHECK_RARE(1e-7, u[dim] * scale - pixel[dim] > 1);
             u[dim] = Clamp(u[dim] * scale - pixel[dim], 0, OneMinusEpsilon);
         }
 
         return u;
     }
 
-    PBRT_CPU_GPU
-    Point2f Get2D() {
-        if (dimension + 1 >= NSobolDimensions)
-            dimension = 2;
-        Point2f u(SampleDimension(dimension), SampleDimension(dimension + 1));
-        dimension += 2;
-        return u;
-    }
-
     std::vector<Sampler> Clone(int n, Allocator alloc);
     std::string ToString() const;
 
@@ -532,11 +552,11 @@ class SobolSampler {
     // SobolSampler Private Methods
     PBRT_CPU_GPU
     Float SampleDimension(int dimension) const {
-        // Return un-randomized Sobol sample if appropriate
-        if (dimension < 2 || randomize == RandomizeStrategy::None)
+        // Return un-randomized Sobol$'$ sample if appropriate
+        if (randomize == RandomizeStrategy::None)
             return SobolSample(sobolIndex, dimension, NoRandomizer());
 
-        // Return randomized Sobol sample using _randomize_
+        // Return randomized Sobol$'$ sample using _randomize_
         uint32_t hash = MixBits((uint64_t(dimension) << 32) ^ GetOptions().seed);
         if (randomize == RandomizeStrategy::PermuteDigits)
             return SobolSample(sobolIndex, dimension, BinaryPermuteScrambler(hash));

src/pbrt/util/lowdiscrepancy.h

@@ -63,8 +63,8 @@ class DigitPermutation {
 };
 
 // Low Discrepancy Declarations
-inline PBRT_CPU_GPU uint64_t SobolIntervalToIndex(const uint32_t log2Scale,
-                                                  uint64_t sampleIndex, const Point2i &p);
+inline PBRT_CPU_GPU uint64_t SobolIntervalToIndex(uint32_t log2Scale,
+                                                  uint64_t sampleIndex, Point2i p);
 
 PBRT_CPU_GPU inline Float BlueNoiseSample(Point2i p, int instance);
 
@@ -154,41 +154,8 @@ PBRT_CPU_GPU inline uint32_t MultiplyGenerator(pstd::span<const uint32_t> C, uin
     return v;
 }
 
-// Laine et al., Stratified Sampling for Stochastic Transparency, Sec 3.1...
-PBRT_CPU_GPU inline uint32_t FastOwenBinaryScramble(uint32_t v, uint32_t hash) {
-    v = ReverseBits32(v);
-    v += hash;
-    v ^= v * 0x6c50b47cu;
-    v ^= v * 0xb82f1e52u;
-    v ^= v * 0xc7afe638u;
-    v ^= v * 0x8d22f6e6u;
-    return ReverseBits32(v);
-}
-
-PBRT_CPU_GPU inline uint32_t OwenScrambleBinaryFull(uint32_t v, uint32_t hash) {
-    if (hash & 1)
-        v ^= 1u << 31;
-
-    for (int b = 1; b < 32; ++b) {
-        uint32_t mask = (~0u) << (32 - b);
-        if (MixBits((v & mask) ^ hash) & (1u << b))
-            v ^= 1u << (31 - b);
-    }
-
-    return v;
-}
-
-template <typename R>
-PBRT_CPU_GPU inline Float SobolSample(int64_t index, int dimension, R randomizer) {
-#ifdef PBRT_FLOAT_AS_DOUBLE
-    return SobolSampleDouble(index, dimension, randomizer);
-#else
-    return SobolSampleFloat(index, dimension, randomizer);
-#endif
-}
-
 template <typename R>
-PBRT_CPU_GPU inline float SobolSampleFloat(int64_t a, int dimension, R randomizer) {
+PBRT_CPU_GPU inline Float SobolSample(int64_t a, int dimension, R randomizer) {
     DCHECK_LT(dimension, NSobolDimensions);
     DCHECK(a >= 0 && a < (1ull << SobolMatrixSize));
     // Compute initial Sobol sample _v_ using generator matrices
@@ -197,6 +164,7 @@ PBRT_CPU_GPU inline float SobolSampleFloat(int64_t a, int dimension, R randomize
         if (a & 1)
             v ^= SobolMatrices32[i];
 
+    // Randomize Sobol sample and return floating-point value
     v = randomizer(v);
     return std::min(v * 0x1p-32f, FloatOneMinusEpsilon);
 }
@@ -233,7 +201,7 @@ PBRT_CPU_GPU inline Float BlueNoiseSample(Point2i p, int instance) {
 // BinaryPermuteScrambler Definition
 struct BinaryPermuteScrambler {
     PBRT_CPU_GPU
-    BinaryPermuteScrambler(uint32_t permutation) : permutation(permutation) {}
+    BinaryPermuteScrambler(uint32_t perm) : permutation(perm) {}
     PBRT_CPU_GPU
     uint32_t operator()(uint32_t v) const { return permutation ^ v; }
     uint32_t permutation;
@@ -243,8 +211,19 @@ struct BinaryPermuteScrambler {
 struct FastOwenScrambler {
     PBRT_CPU_GPU
     FastOwenScrambler(uint32_t seed) : seed(seed) {}
+    // FastOwenScrambler Public Methods
+    // Laine et al., Stratified Sampling for Stochastic Transparency, Sec 3.1...
     PBRT_CPU_GPU
-    uint32_t operator()(uint32_t v) const { return FastOwenBinaryScramble(v, seed); }
+    uint32_t operator()(uint32_t v) const {
+        v = ReverseBits32(v);
+        v += seed;
+        v ^= v * 0x6c50b47cu;
+        v ^= v * 0xb82f1e52u;
+        v ^= v * 0xc7afe638u;
+        v ^= v * 0x8d22f6e6u;
+        return ReverseBits32(v);
+    }
+
     uint32_t seed;
 };
 
@@ -252,8 +231,20 @@ struct FastOwenScrambler {
 struct OwenScrambler {
     PBRT_CPU_GPU
     OwenScrambler(uint32_t seed) : seed(seed) {}
+    // OwenScrambler Public Methods
     PBRT_CPU_GPU
-    uint32_t operator()(uint32_t v) const { return OwenScrambleBinaryFull(v, seed); }
+    uint32_t operator()(uint32_t v) const {
+        if (seed & 1)
+            v ^= 1u << 31;
+        for (int b = 1; b < 32; ++b) {
+            // Apply Owen scrambling to binary digit _b_ in _v_
+            uint32_t mask = (~0u) << (32 - b);
+            if (MixBits((v & mask) ^ seed) & (1u << b))
+                v ^= 1u << (31 - b);
+        }
+        return v;
+    }
+
     uint32_t seed;
 };
 
@@ -263,7 +254,7 @@ enum class RandomizeStrategy { None, PermuteDigits, FastOwen, Owen };
 std::string ToString(RandomizeStrategy r);
 
 PBRT_CPU_GPU
-inline uint64_t SobolIntervalToIndex(uint32_t m, uint64_t frame, const Point2i &p) {
+inline uint64_t SobolIntervalToIndex(uint32_t m, uint64_t frame, Point2i p) {
     if (m == 0)
         return frame;
 
@@ -285,26 +276,6 @@ inline uint64_t SobolIntervalToIndex(uint32_t m, uint64_t frame, const Point2i &
     return index;
 }
 
-PBRT_CPU_GPU
-inline uint64_t SobolSampleBits64(int64_t a, int dimension) {
-    CHECK_LT(dimension, NSobolDimensions);
-    DCHECK(a >= 0 && a < (1ull << SobolMatrixSize));
-    uint64_t v = 0;
-    for (int i = dimension * SobolMatrixSize; a != 0; a >>= 1, i++)
-        if (a & 1)
-            v ^= SobolMatrices64[i];
-    return v;
-}
-
-template <typename R>
-PBRT_CPU_GPU inline double SobolSampleDouble(int64_t a, int dimension, R randomizer) {
-    uint64_t v = SobolSampleBits64(a, dimension);
-    // FIXME? We just scramble the high bits here...
-    uint32_t vs = randomizer(v >> 32);
-    v = (uint64_t(vs) << 32) | (v & 0xffffffff);
-    return std::min(v * (1.0 / (1ULL << SobolMatrixSize)), DoubleOneMinusEpsilon);
-}
-
 }  // namespace pbrt
 
 #endif  // PBRT_UTIL_LOWDISCREPANCY_H

src/pbrt/util/sampling_test.cpp

@@ -144,18 +144,10 @@ TEST(LowDiscrepancy, RadicalInverse) {
     }
 }
 
-TEST(LowDiscrepancy, Sobol) {
-    // Check that float and double variants match (as float values).
-    for (int i = 0; i < 256; ++i) {
-        for (int dim = 0; dim < 100; ++dim) {
-            EXPECT_EQ(SobolSampleFloat(i, dim, NoRandomizer()),
-                      (float)SobolSampleDouble(i, dim, NoRandomizer()));
-        }
-    }
-
+TEST(LowDiscrepancy, SobolFirstDimension) {
     // Make sure first dimension is the regular base 2 radical inverse
     for (int i = 0; i < 8192; ++i) {
-        EXPECT_EQ(SobolSampleFloat(i, 0, NoRandomizer()),
+        EXPECT_EQ(SobolSample(i, 0, NoRandomizer()),
                   ReverseBits32(i) * 2.3283064365386963e-10f);
     }
 }

src/pbrt/util/sobolmatrices.h

@@ -15,7 +15,6 @@ namespace pbrt {
 static constexpr int NSobolDimensions = 1024;
 static constexpr int SobolMatrixSize = 52;
 extern PBRT_CONST uint32_t SobolMatrices32[NSobolDimensions * SobolMatrixSize];
-extern PBRT_CONST uint64_t SobolMatrices64[NSobolDimensions * SobolMatrixSize];
 
 extern PBRT_CONST uint64_t VdCSobolMatrices[][SobolMatrixSize];
 extern PBRT_CONST uint64_t VdCSobolMatricesInv[][SobolMatrixSize];