Update from book source

Folded HaltonPixelIndexer back into HaltonSampler, among other small changes.

src/pbrt/film.cpp

@@ -240,7 +240,6 @@ pstd::optional<SquareMatrix<3>> PixelSensor::SolveXYZFromSensorRGB(
     SpectrumHandle sensorIllum, SpectrumHandle outputIllum) const {
     Float rgbCamera[24][3], xyzOutput[24][3];
     // Compute _rgbCamera_ values for training swatches
-    RGB outputWhite = IlluminantToSensorRGB(outputIllum);
     for (size_t i = 0; i < swatchReflectances.size(); ++i) {
         RGB rgb = ProjectReflectance<RGB>(swatchReflectances[i], sensorIllum, &r_bar,
                                           &g_bar, &b_bar);
@@ -264,8 +263,8 @@ pstd::optional<SquareMatrix<3>> PixelSensor::SolveXYZFromSensorRGB(
 }
 
 // Swatch reflectances are taken from Danny Pascale's Macbeth chart measurements
-// BabelColor ColorChecker data: Copyright © 2004-2012 Danny Pascale (www.babelcolor.com);
-// used by permission.
+// BabelColor ColorChecker data: Copyright (c) 2004-2012 Danny Pascale
+// (www.babelcolor.com); used by permission.
 // http://www.babelcolor.com/index_htm_files/ColorChecker_RGB_and_spectra.zip
 std::vector<SpectrumHandle> PixelSensor::swatchReflectances{
     PiecewiseLinearSpectrum::FromInterleaved(

src/pbrt/film.h

@@ -79,8 +79,8 @@ class PixelSensor {
     PBRT_CPU_GPU
     RGB ToSensorRGB(const SampledSpectrum &s, const SampledWavelengths &lambda) const {
         return RGB((r_bar.Sample(lambda) * s).Average(),
-                      (g_bar.Sample(lambda) * s).Average(),
-                      (b_bar.Sample(lambda) * s).Average());
+                   (g_bar.Sample(lambda) * s).Average(),
+                   (b_bar.Sample(lambda) * s).Average());
     }
 
     // PixelSensor Public Members

src/pbrt/samplers.cpp

@@ -31,12 +31,44 @@ std::string SamplerHandle::ToString() const {
 // HaltonSampler Method Definitions
 HaltonSampler::HaltonSampler(int samplesPerPixel, const Point2i &fullResolution,
                              pstd::vector<DigitPermutation> *dp, Allocator alloc)
-    : digitPermutations(dp),
-      samplesPerPixel(samplesPerPixel),
-      haltonPixelIndexer(fullResolution) {
+    : digitPermutations(dp), samplesPerPixel(samplesPerPixel) {
     // Generate random digit permutations for Halton sampler
     if (!dp)
         digitPermutations = ComputeRadicalInversePermutations(0, alloc);
+
+    // Find radical inverse base scales and exponents that cover sampling area
+    for (int i = 0; i < 2; ++i) {
+        int base = (i == 0) ? 2 : 3;
+        int scale = 1, exp = 0;
+        while (scale < std::min(fullResolution[i], MaxHaltonResolution)) {
+            scale *= base;
+            ++exp;
+        }
+        baseScales[i] = scale;
+        baseExponents[i] = exp;
+    }
+
+    // Compute multiplicative inverses for _baseScales_
+    multInverse[0] = multiplicativeInverse(baseScales[1], baseScales[0]);
+    multInverse[1] = multiplicativeInverse(baseScales[0], baseScales[1]);
+}
+
+uint64_t HaltonSampler::multiplicativeInverse(int64_t a, int64_t n) {
+    int64_t x, y;
+    extendedGCD(a, n, &x, &y);
+    return Mod(x, n);
+}
+
+void HaltonSampler::extendedGCD(uint64_t a, uint64_t b, int64_t *x, int64_t *y) {
+    if (b == 0) {
+        *x = 1;
+        *y = 0;
+        return;
+    }
+    int64_t d = a / b, xp, yp;
+    extendedGCD(b, a % b, &xp, &yp);
+    *x = yp;
+    *y = xp - (d * yp);
 }
 
 std::vector<SamplerHandle> HaltonSampler::Clone(int n, Allocator alloc) {
@@ -50,10 +82,9 @@ std::vector<SamplerHandle> HaltonSampler::Clone(int n, Allocator alloc) {
 }
 
 std::string HaltonSampler::ToString() const {
-    return StringPrintf("[ HaltonSampler digitPermutations: %p haltonPixelIndexer: %s "
-                        "pixel: %s sampleIndex: %d dimension: %d samplesPerPixel: %d ]",
-                        digitPermutations, haltonPixelIndexer, pixel, sampleIndex,
-                        dimension, samplesPerPixel);
+    return StringPrintf("[ HaltonSampler digitPermutations: %p "
+                        "sampleIndex: %d dimension: %d samplesPerPixel: %d ]",
+                        digitPermutations, sampleIndex, dimension, samplesPerPixel);
 }
 
 HaltonSampler *HaltonSampler::Create(const ParameterDictionary &parameters,

src/pbrt/samplers.h

@@ -46,12 +46,21 @@ class HaltonSampler {
 
     PBRT_CPU_GPU
     void StartPixelSample(const Point2i &p, int index, int dim) {
-        if (p != pixel)
-            haltonPixelIndexer.SetPixel(p);
-        haltonPixelIndexer.SetPixelSample(index);
-
-        pixel = p;
-        sampleIndex = index;
+        int sampleStride = baseScales[0] * baseScales[1];
+        int64_t pixelSampleForIndex = 0;
+
+        if (sampleStride > 1) {
+            Point2i pm(Mod(p[0], MaxHaltonResolution), Mod(p[1], MaxHaltonResolution));
+            for (int i = 0; i < 2; ++i) {
+                uint64_t dimOffset =
+                    (i == 0) ? InverseRadicalInverse(pm[i], 2, baseExponents[i])
+                             : InverseRadicalInverse(pm[i], 3, baseExponents[i]);
+                pixelSampleForIndex +=
+                    dimOffset * (sampleStride / baseScales[i]) * multInverse[i];
+            }
+            pixelSampleForIndex %= sampleStride;
+        }
+        sampleIndex = pixelSampleForIndex + index * sampleStride;
         dimension = dim;
     }
 
@@ -60,24 +69,23 @@ class HaltonSampler {
         if (dimension >= PrimeTableSize)
             dimension = 2;
         int dim = dimension++;
-        return ScrambledRadicalInverse(dim, haltonPixelIndexer.SampleIndex(),
-                                       (*digitPermutations)[dim]);
+        return ScrambledRadicalInverse(dim, sampleIndex, (*digitPermutations)[dim]);
     }
 
     PBRT_CPU_GPU
     Point2f Get2D() {
         if (dimension == 0) {
             dimension += 2;
-            return haltonPixelIndexer.SampleFirst2D();
+            return {RadicalInverse(0, sampleIndex >> baseExponents[0]),
+                    RadicalInverse(1, sampleIndex / baseScales[1])};
         } else {
             if (dimension + 1 >= PrimeTableSize)
                 dimension = 2;
 
             int dim = dimension;
             dimension += 2;
-            return {ScrambledRadicalInverse(dim, haltonPixelIndexer.SampleIndex(),
-                                            (*digitPermutations)[dim]),
-                    ScrambledRadicalInverse(dim + 1, haltonPixelIndexer.SampleIndex(),
+            return {ScrambledRadicalInverse(dim, sampleIndex, (*digitPermutations)[dim]),
+                    ScrambledRadicalInverse(dim + 1, sampleIndex,
                                             (*digitPermutations)[dim + 1])};
         }
     }
@@ -86,14 +94,18 @@ class HaltonSampler {
     std::string ToString() const;
 
   private:
+    // HaltonSampler Private Methods
+    static uint64_t multiplicativeInverse(int64_t a, int64_t n);
+    static void extendedGCD(uint64_t a, uint64_t b, int64_t *x, int64_t *y);
+
     // HaltonSampler Private Members
     pstd::vector<DigitPermutation> *digitPermutations;
     int samplesPerPixel;
-    HaltonPixelIndexer haltonPixelIndexer;
-    Point2i pixel =
-        Point2i(std::numeric_limits<int>::max(), std::numeric_limits<int>::max());
-    int sampleIndex = 0;
+    static constexpr int MaxHaltonResolution = 128;
+    int64_t sampleIndex = 0;
     int dimension = 0;
+    Point2i baseScales, baseExponents;
+    int multInverse[2];
 };
 
 // PaddedSobolSampler Definition

src/pbrt/util/lowdiscrepancy.cpp

@@ -29,46 +29,6 @@ std::string DigitPermutation::ToString() const {
     return s + " ]";
 }
 
-// HaltonIndexer Local Constants
-constexpr int HaltonPixelIndexer::MaxHaltonResolution;
-
-// HaltonIndexer Method Definitions
-HaltonPixelIndexer::HaltonPixelIndexer(const Point2i &fullResolution) {
-    // Find radical inverse base scales and exponents that cover sampling area
-    for (int i = 0; i < 2; ++i) {
-        int base = (i == 0) ? 2 : 3;
-        int scale = 1, exp = 0;
-        while (scale < std::min(fullResolution[i], MaxHaltonResolution)) {
-            scale *= base;
-            ++exp;
-        }
-        baseScales[i] = scale;
-        baseExponents[i] = exp;
-    }
-
-    // Compute multiplicative inverses for _baseScales_
-    multInverse[0] = multiplicativeInverse(baseScales[1], baseScales[0]);
-    multInverse[1] = multiplicativeInverse(baseScales[0], baseScales[1]);
-}
-
-uint64_t HaltonPixelIndexer::multiplicativeInverse(int64_t a, int64_t n) {
-    int64_t x, y;
-    extendedGCD(a, n, &x, &y);
-    return Mod(x, n);
-}
-
-void HaltonPixelIndexer::extendedGCD(uint64_t a, uint64_t b, int64_t *x, int64_t *y) {
-    if (b == 0) {
-        *x = 1;
-        *y = 0;
-        return;
-    }
-    int64_t d = a / b, xp, yp;
-    extendedGCD(b, a % b, &xp, &yp);
-    *x = yp;
-    *y = xp - (d * yp);
-}
-
 std::string ToString(RandomizeStrategy r) {
     switch (r) {
     case RandomizeStrategy::None:

src/pbrt/util/lowdiscrepancy.h

@@ -30,8 +30,7 @@ class DigitPermutation {
         CHECK_LT(base, 65536);  // uint16_t
         // Compute number of digits needed for _base_
         nDigits = 0;
-        Float invBase = (Float)1 / (Float)base;
-        Float invBaseN = 1;
+        Float invBase = (Float)1 / (Float)base, invBaseN = 1;
         while (1 - invBaseN < 1) {
             ++nDigits;
             invBaseN *= invBase;
@@ -40,10 +39,11 @@ class DigitPermutation {
         permutations = alloc.allocate_object<uint16_t>(nDigits * base);
         for (int digitIndex = 0; digitIndex < nDigits; ++digitIndex) {
             // Compute random permutation for _digitIndex_
-            uint32_t digitSeed = (base * 32 + digitIndex) ^ seed;
-            for (int digitValue = 0; digitValue < base; ++digitValue)
-                Perm(digitIndex, digitValue) =
-                    PermutationElement(digitValue, base, digitSeed);
+            uint32_t digitSeed = MixBits(((base << 8) + digitIndex) ^ seed);
+            for (int digitValue = 0; digitValue < base; ++digitValue) {
+                int index = digitIndex * base + digitValue;
+                permutations[index] = PermutationElement(digitValue, base, digitSeed);
+            }
         }
     }
 
@@ -56,17 +56,9 @@ class DigitPermutation {
 
     std::string ToString() const;
 
-    int base;
-
   private:
-    // DigitPermutation Private Methods
-    PBRT_CPU_GPU
-    uint16_t &Perm(int digitIndex, int digitValue) {
-        return permutations[digitIndex * base + digitValue];
-    }
-
-    int nDigits;
-    // indexed by [digitIndex * base + digitValue]
+    // DigitPermutation Private Members
+    int base, nDigits;
     uint16_t *permutations;
 };
 
@@ -94,9 +86,8 @@ struct NoRandomizer {
 // Low Discrepancy Inline Functions
 PBRT_CPU_GPU inline Float RadicalInverse(int baseIndex, uint64_t a) {
     int base = Primes[baseIndex];
-    Float invBase = (Float)1 / (Float)base;
+    Float invBase = (Float)1 / (Float)base, invBaseN = 1;
     uint64_t reversedDigits = 0;
-    Float invBaseN = 1;
     while (a) {
         // Extract least significant digit from _a_ and update _reversedDigits_
         uint64_t next = a / base;
@@ -122,11 +113,11 @@ inline uint64_t InverseRadicalInverse(uint64_t inverse, int base, int nDigits) {
 PBRT_CPU_GPU inline Float ScrambledRadicalInverse(int baseIndex, uint64_t a,
                                                   const DigitPermutation &perm) {
     int base = Primes[baseIndex];
-    const Float invBase = (Float)1 / (Float)base;
+    Float invBase = (Float)1 / (Float)base, invBaseN = 1;
     uint64_t reversedDigits = 0;
-    Float invBaseN = 1;
     int digitIndex = 0;
     while (1 - invBaseN < 1) {
+        // Permute least significant digit from _a_ and update _reversedDigits_
         uint64_t next = a / base;
         int digitValue = a - next * base;
         reversedDigits = reversedDigits * base + perm.Permute(digitIndex, digitValue);
@@ -230,66 +221,6 @@ struct OwenScrambler {
     uint32_t seed;
 };
 
-// HaltonPixelIndexer Defintion
-class HaltonPixelIndexer {
-  public:
-    // HaltonPixelIndexer Public Methods
-    HaltonPixelIndexer(const Point2i &fullResolution);
-
-    PBRT_CPU_GPU
-    void SetPixel(const Point2i &p) {
-        pixelSampleForIndex = 0;
-
-        int sampleStride = baseScales[0] * baseScales[1];
-        if (sampleStride > 1) {
-            Point2i pm(Mod(p[0], MaxHaltonResolution), Mod(p[1], MaxHaltonResolution));
-            for (int i = 0; i < 2; ++i) {
-                uint64_t dimOffset =
-                    (i == 0) ? InverseRadicalInverse(pm[i], 2, baseExponents[i])
-                             : InverseRadicalInverse(pm[i], 3, baseExponents[i]);
-                pixelSampleForIndex +=
-                    dimOffset * (sampleStride / baseScales[i]) * multInverse[i];
-            }
-            pixelSampleForIndex %= sampleStride;
-        }
-    }
-
-    PBRT_CPU_GPU
-    void SetPixelSample(int pixelSample) {
-        int sampleStride = baseScales[0] * baseScales[1];
-        sampleIndex = pixelSampleForIndex + pixelSample * sampleStride;
-    }
-
-    PBRT_CPU_GPU
-    Point2f SampleFirst2D() const {
-        return {RadicalInverse(0, sampleIndex >> baseExponents[0]),
-                RadicalInverse(1, sampleIndex / baseScales[1])};
-    }
-
-    PBRT_CPU_GPU
-    int64_t SampleIndex() const { return sampleIndex; }
-
-    std::string ToString() const {
-        return StringPrintf("[ HaltonPixelIndexer pixelSampleForIndex: %d "
-                            "sampleIndex: %d baseScales: %s baseExponents: %s "
-                            "multInverse[0]: %d multInverse[1]: %d ]",
-                            pixelSampleForIndex, sampleIndex, baseScales, baseExponents,
-                            multInverse[0], multInverse[1]);
-    }
-
-  private:
-    // HaltonPixelIndexer Private Methods
-    static uint64_t multiplicativeInverse(int64_t a, int64_t n);
-    static void extendedGCD(uint64_t a, uint64_t b, int64_t *x, int64_t *y);
-
-    // HaltonPixelIndexer Private Members
-    static constexpr int MaxHaltonResolution = 128;
-    Point2i baseScales, baseExponents;
-    int multInverse[2];
-    int64_t pixelSampleForIndex;
-    int64_t sampleIndex;
-};
-
 enum class RandomizeStrategy { None, CranleyPatterson, Xor, Owen };
 
 std::string ToString(RandomizeStrategy r);