Cycles: improve Progressive Multi-Jittered sampling

Fix two issues in the previous implementation:
* Only power-of-two prefixes were progressively stratified, not suffixes.
  This resulted in unnecessarily increased noise when using non-power-of-two
  sample counts.
* In order to try to get away with just a single sample pattern, the code
  used a combination of sample index shuffling and Cranley-Patterson rotation.
  Index shuffling is normally fine, but due to the sample patterns themselves
  not being quite right (as described above) this actually resulted in
  additional increased noise. Cranley-Patterson, on the other hand, always
  increases noise with randomized (t,s) nets like PMJ02, and should be avoided
  with these kinds of sequences.

Addressed with the following changes:
* Replace the sample pattern generation code with a much simpler algorithm
  recently published in the paper "Stochastic Generation of (t, s) Sample
  Sequences". This new implementation is easier to verify, produces fully
  progressively stratified PMJ02, and is *far* faster than the previous code,
  being O(N) in the number of samples generated.
* It keeps the sample index shuffling, which works correctly now due to the
  improved sample patterns. But it now uses a newer high-quality hash instead
  of the original Laine-Karras hash.
* The scrambling distance feature cannot (to my knowledge) be implemented with
  any decorrelation strategy other than Cranley-Patterson, so Cranley-Patterson
  is still used when that feature is enabled. But it is now disabled otherwise,
  since it increases noise.
* In place of Cranley-Patterson, multiple independent patterns are generated
  and randomly chosen for different pixels and dimensions as described in the
  original PMJ paper. In this patch, the pattern selection is done via
  hash-based shuffling to ensure there are no repeats within a single pixel
  until all patterns have been used.

The combination of these fixes brings the quality of Cycles' PMJ sampler in
line with the previously submitted Sobol-Burley sampler in D15679. They are
essentially indistinguishable in terms of quality/noise, which is expected
since they are both randomized (0,2) sequences.

Differential Revision: https://developer.blender.org/D15746

7 files changed, 166 insertions, 388 deletions

diff --git a/intern/cycles/kernel/integrator/subsurface_random_walk.h b/intern/cycles/kernel/integrator/subsurface_random_walk.h
index baca0d745e8..e0c69c96fc6 100644
--- a/intern/cycles/kernel/integrator/subsurface_random_walk.h
+++ b/intern/cycles/kernel/integrator/subsurface_random_walk.h
@@ -229,7 +229,7 @@ ccl_device_inline bool subsurface_random_walk(KernelGlobals kg,
   const float phase_log = logf((diffusion_length + 1.0f) / (diffusion_length - 1.0f));
 
   /* Modify state for RNGs, decorrelated from other paths. */
-  rng_state.rng_hash = hash_cmj_seeded_uint(rng_state.rng_hash + rng_state.rng_offset, 0xdeadbeef);
+  rng_state.rng_hash = hash_hp_seeded_uint(rng_state.rng_hash + rng_state.rng_offset, 0xdeadbeef);
 
   /* Random walk until we hit the surface again. */
   bool hit = false;
diff --git a/intern/cycles/kernel/sample/jitter.h b/intern/cycles/kernel/sample/jitter.h
index dd170cf2120..6a9ff1beec5 100644
--- a/intern/cycles/kernel/sample/jitter.h
+++ b/intern/cycles/kernel/sample/jitter.h
@@ -7,57 +7,40 @@
 #pragma once
 CCL_NAMESPACE_BEGIN
 
-ccl_device_inline uint32_t nested_uniform_scramble(uint32_t x, uint32_t seed)
-{
-  x = reverse_integer_bits(x);
-  x = laine_karras_permutation(x, seed);
-  x = reverse_integer_bits(x);
-
-  return x;
-}
-
 ccl_device float pmj_sample_1D(KernelGlobals kg, uint sample, uint rng_hash, uint dimension)
 {
-  uint hash = rng_hash;
-  float jitter_x = 0.0f;
-  if (kernel_data.integrator.scrambling_distance < 1.0f) {
-    hash = kernel_data.integrator.seed;
+  uint seed = rng_hash;
 
-    jitter_x = hash_wang_seeded_float(dimension, rng_hash) *
-               kernel_data.integrator.scrambling_distance;
+  /* Use the same sample sequence seed for all pixels when using
+   * scrambling distance. */
+  if (kernel_data.integrator.scrambling_distance < 1.0f) {
+    seed = kernel_data.integrator.seed;
   }
 
-  /* Perform Owen shuffle of the sample number to reorder the samples. */
-  const uint rv = hash_cmj_seeded_uint(dimension, hash);
-#ifdef _XOR_SHUFFLE_
-#  warning "Using XOR shuffle."
-  const uint s = sample ^ rv;
-#else /* Use _OWEN_SHUFFLE_ for reordering. */
-  const uint s = nested_uniform_scramble(sample, rv);
-#endif
-
-  /* Based on the sample number a sample pattern is selected and offset by the dimension. */
-  const uint sample_set = s / NUM_PMJ_SAMPLES;
-  const uint d = (dimension + sample_set);
-  const uint dim = d % NUM_PMJ_PATTERNS;
-
-  /* The PMJ sample sets contain a sample with (x,y) with NUM_PMJ_SAMPLES so for 1D
-   *  the x part is used for even dims and the y for odd. */
-  int index = 2 * ((dim >> 1) * NUM_PMJ_SAMPLES + (s % NUM_PMJ_SAMPLES)) + (dim & 1);
-
-  float fx = kernel_data_fetch(sample_pattern_lut, index);
-
-#ifndef _NO_CRANLEY_PATTERSON_ROTATION_
-  /* Use Cranley-Patterson rotation to displace the sample pattern. */
-  float dx = hash_cmj_seeded_float(d, hash);
-  /* Jitter sample locations and map back into [0 1]. */
-  fx = fx + dx + jitter_x;
-  fx = fx - floorf(fx);
-#else
-#  warning "Not using Cranley-Patterson Rotation."
-#endif
+  /* Shuffle the pattern order and sample index to better decorrelate
+   * dimensions and make the most of the finite patterns we have.
+   * The funky sample mask stuff is to ensure that we only shuffle
+   * *within* the current sample pattern, which is necessary to avoid
+   * early repeat pattern use. */
+  uint pattern_i = hash_shuffle_uint(dimension, NUM_PMJ_PATTERNS, seed);
+  /* NUM_PMJ_SAMPLES should be a power of two, so this results in a mask. */
+  uint sample_mask = NUM_PMJ_SAMPLES - 1;
+  uint sample_shuffled = nested_uniform_scramble(sample, hash_wang_seeded_uint(dimension, seed));
+  sample = (sample & ~sample_mask) | (sample_shuffled & sample_mask);
+
+  /* Fetch the sample. */
+  uint index = ((pattern_i * NUM_PMJ_SAMPLES) + sample) % (NUM_PMJ_SAMPLES * NUM_PMJ_PATTERNS);
+  float x = kernel_data_fetch(sample_pattern_lut, index * 2);
+
+  /* Do limited Cranley-Patterson rotation when using scrambling distance. */
+  if (kernel_data.integrator.scrambling_distance < 1.0f) {
+    float jitter_x = hash_wang_seeded_float(dimension, rng_hash) *
+                     kernel_data.integrator.scrambling_distance;
+    x += jitter_x;
+    x -= floorf(x);
+  }
 
-  return fx;
+  return x;
 }
 
 ccl_device void pmj_sample_2D(KernelGlobals kg,
@@ -67,51 +50,41 @@ ccl_device void pmj_sample_2D(KernelGlobals kg,
                               ccl_private float *x,
                               ccl_private float *y)
 {
-  uint hash = rng_hash;
-  float jitter_x = 0.0f;
-  float jitter_y = 0.0f;
-  if (kernel_data.integrator.scrambling_distance < 1.0f) {
-    hash = kernel_data.integrator.seed;
+  uint seed = rng_hash;
 
-    jitter_x = hash_wang_seeded_float(dimension, rng_hash) *
-               kernel_data.integrator.scrambling_distance;
-    jitter_y = hash_wang_seeded_float(dimension + 1, rng_hash) *
-               kernel_data.integrator.scrambling_distance;
+  /* Use the same sample sequence seed for all pixels when using
+   * scrambling distance. */
+  if (kernel_data.integrator.scrambling_distance < 1.0f) {
+    seed = kernel_data.integrator.seed;
   }
 
-  /* Perform a shuffle on the sample number to reorder the samples. */
-  const uint rv = hash_cmj_seeded_uint(dimension, hash);
-#ifdef _XOR_SHUFFLE_
-#  warning "Using XOR shuffle."
-  const uint s = sample ^ rv;
-#else /* Use _OWEN_SHUFFLE_ for reordering. */
-  const uint s = nested_uniform_scramble(sample, rv);
-#endif
-
-  /* Based on the sample number a sample pattern is selected and offset by the dimension. */
-  const uint sample_set = s / NUM_PMJ_SAMPLES;
-  const uint d = dimension + sample_set;
-  uint dim = d % NUM_PMJ_PATTERNS;
-  int index = 2 * (dim * NUM_PMJ_SAMPLES + (s % NUM_PMJ_SAMPLES));
-
-  float fx = kernel_data_fetch(sample_pattern_lut, index);
-  float fy = kernel_data_fetch(sample_pattern_lut, index + 1);
-
-#ifndef _NO_CRANLEY_PATTERSON_ROTATION_
-  /* Use Cranley-Patterson rotation to displace the sample pattern. */
-  float dx = hash_cmj_seeded_float(d, hash);
-  float dy = hash_cmj_seeded_float(d + 1, hash);
-  /* Jitter sample locations and map back to the unit square [0 1]x[0 1]. */
-  float sx = fx + dx + jitter_x;
-  float sy = fy + dy + jitter_y;
-  sx = sx - floorf(sx);
-  sy = sy - floorf(sy);
-#else
-#  warning "Not using Cranley Patterson Rotation."
-#endif
-
-  (*x) = sx;
-  (*y) = sy;
+  /* Shuffle the pattern order and sample index to better decorrelate
+   * dimensions and make the most of the finite patterns we have.
+   * The funky sample mask stuff is to ensure that we only shuffle
+   * *within* the current sample pattern, which is necessary to avoid
+   * early repeat pattern use. */
+  uint pattern_i = hash_shuffle_uint(dimension, NUM_PMJ_PATTERNS, seed);
+  /* NUM_PMJ_SAMPLES should be a power of two, so this results in a mask. */
+  uint sample_mask = NUM_PMJ_SAMPLES - 1;
+  uint sample_shuffled = nested_uniform_scramble(sample, hash_wang_seeded_uint(dimension, seed));
+  sample = (sample & ~sample_mask) | (sample_shuffled & sample_mask);
+
+  /* Fetch the sample. */
+  uint index = ((pattern_i * NUM_PMJ_SAMPLES) + sample) % (NUM_PMJ_SAMPLES * NUM_PMJ_PATTERNS);
+  (*x) = kernel_data_fetch(sample_pattern_lut, index * 2);
+  (*y) = kernel_data_fetch(sample_pattern_lut, index * 2 + 1);
+
+  /* Do limited Cranley-Patterson rotation when using scrambling distance. */
+  if (kernel_data.integrator.scrambling_distance < 1.0f) {
+    float jitter_x = hash_wang_seeded_float(dimension, rng_hash) *
+                     kernel_data.integrator.scrambling_distance;
+    float jitter_y = hash_wang_seeded_float(dimension, rng_hash ^ 0xca0e1151) *
+                     kernel_data.integrator.scrambling_distance;
+    (*x) += jitter_x;
+    (*y) += jitter_y;
+    (*x) -= floorf(*x);
+    (*y) -= floorf(*y);
+  }
 }
 
 CCL_NAMESPACE_END
diff --git a/intern/cycles/kernel/sample/util.h b/intern/cycles/kernel/sample/util.h
index 33056bb7819..29cda179aa2 100644
--- a/intern/cycles/kernel/sample/util.h
+++ b/intern/cycles/kernel/sample/util.h
@@ -8,7 +8,7 @@
 CCL_NAMESPACE_BEGIN
 
 /*
- * Performs base-2 Owen scrambling on a reversed-bit integer.
+ * Performs base-2 Owen scrambling on a reversed-bit unsigned integer.
  *
  * This is equivalent to the Laine-Karras permutation, but much higher
  * quality.  See https://psychopath.io/post/2021_01_30_building_a_better_lk_hash
@@ -25,21 +25,11 @@ ccl_device_inline uint reversed_bit_owen(uint n, uint seed)
 }
 
 /*
- * Performs base-2 Owen scrambling on a reversed-bit integer.
- *
- * This is here for backwards-compatibility, and can be replaced
- * with reversed_bit_owen() above at some point.
- * See https://developer.blender.org/D15679#426304
+ * Performs base-2 Owen scrambling on an unsigned integer.
  */
-ccl_device_inline uint laine_karras_permutation(uint x, uint seed)
+ccl_device_inline uint nested_uniform_scramble(uint i, uint seed)
 {
-  x += seed;
-  x ^= (x * 0x6c50b47cu);
-  x ^= x * 0xb82f1e52u;
-  x ^= x * 0xc7afe638u;
-  x ^= x * 0x8d22f6e6u;
-
-  return x;
+  return reverse_integer_bits(reversed_bit_owen(reverse_integer_bits(i), seed));
 }
 
 CCL_NAMESPACE_END
diff --git a/intern/cycles/kernel/types.h b/intern/cycles/kernel/types.h
index f55ace1a227..655c9c5503b 100644
--- a/intern/cycles/kernel/types.h
+++ b/intern/cycles/kernel/types.h
@@ -1364,10 +1364,14 @@ typedef struct KernelShaderEvalInput {
 } KernelShaderEvalInput;
 static_assert_align(KernelShaderEvalInput, 16);
 
-/* Pre-computed sample table sizes for PMJ02 sampler. */
+/* Pre-computed sample table sizes for PMJ02 sampler.
+ *
+ * Note: divisions *must* be a power of two, and patterns
+ * ideally should be as well.
+ */
 #define NUM_PMJ_DIVISIONS 32
 #define NUM_PMJ_SAMPLES ((NUM_PMJ_DIVISIONS) * (NUM_PMJ_DIVISIONS))
-#define NUM_PMJ_PATTERNS 1
+#define NUM_PMJ_PATTERNS 64
 
 /* Device kernels.
  *
diff --git a/intern/cycles/scene/jitter.cpp b/intern/cycles/scene/jitter.cpp
index 4f9c4fb9418..8287a029752 100644
--- a/intern/cycles/scene/jitter.cpp
+++ b/intern/cycles/scene/jitter.cpp
@@ -2,267 +2,56 @@
  * Copyright 2019-2022 Blender Foundation */
 
 /* This file is based on "Progressive Multi-Jittered Sample Sequences"
- * by Per Christensen, Andrew Kensler and Charlie Kilpatrick.
- * http://graphics.pixar.com/library/ProgressiveMultiJitteredSampling/paper.pdf
- *
- * Performance can be improved in the future by implementing the new
- * algorithm from Matt Pharr in  http://jcgt.org/published/0008/01/04/
- * "Efficient Generation of Points that Satisfy Two-Dimensional Elementary Intervals"
+ * by Christensen, Kensler, and Kilpatrick, but with a much simpler and
+ * faster implementation based on "Stochastic Generation of (t, s)
+ * Sample Sequences" by Helmer, Christensen, and Kensler.
  */
 
 #include "scene/jitter.h"
+#include "util/hash.h"
 
 #include <math.h>
 #include <vector>
 
 CCL_NAMESPACE_BEGIN
 
-static uint cmj_hash(uint i, uint p)
-{
-  i ^= p;
-  i ^= i >> 17;
-  i ^= i >> 10;
-  i *= 0xb36534e5;
-  i ^= i >> 12;
-  i ^= i >> 21;
-  i *= 0x93fc4795;
-  i ^= 0xdf6e307f;
-  i ^= i >> 17;
-  i *= 1 | p >> 18;
-
-  return i;
-}
-
-static float cmj_randfloat(uint i, uint p)
-{
-  return cmj_hash(i, p) * (1.0f / 4294967808.0f);
-}
-
-class PMJ_Generator {
- public:
-  static void generate_2D(float2 points[], int size, int rng_seed_in)
-  {
-    PMJ_Generator g(rng_seed_in);
-    points[0].x = g.rnd();
-    points[0].y = g.rnd();
-    int N = 1;
-    while (N < size) {
-      g.extend_sequence_even(points, N);
-      g.extend_sequence_odd(points, 2 * N);
-      N = 4 * N;
-    }
-  }
-
- protected:
-  PMJ_Generator(int rnd_seed_in) : num_samples(1), rnd_index(2), rnd_seed(rnd_seed_in)
-  {
-  }
-
-  float rnd()
-  {
-    return cmj_randfloat(++rnd_index, rnd_seed);
-  }
-
-  virtual void mark_occupied_strata(float2 points[], int N)
-  {
-    int NN = 2 * N;
-    for (int s = 0; s < NN; ++s) {
-      occupied1Dx[s] = occupied1Dy[s] = false;
-    }
-    for (int s = 0; s < N; ++s) {
-      int xstratum = (int)(NN * points[s].x);
-      int ystratum = (int)(NN * points[s].y);
-      occupied1Dx[xstratum] = true;
-      occupied1Dy[ystratum] = true;
-    }
-  }
-
-  virtual void generate_sample_point(
-      float2 points[], float i, float j, float xhalf, float yhalf, int n, int N)
-  {
-    int NN = 2 * N;
-    float2 pt;
-    int xstratum, ystratum;
-    do {
-      pt.x = (i + 0.5f * (xhalf + rnd())) / n;
-      xstratum = (int)(NN * pt.x);
-    } while (occupied1Dx[xstratum]);
-    do {
-      pt.y = (j + 0.5f * (yhalf + rnd())) / n;
-      ystratum = (int)(NN * pt.y);
-    } while (occupied1Dy[ystratum]);
-    occupied1Dx[xstratum] = true;
-    occupied1Dy[ystratum] = true;
-    points[num_samples] = pt;
-    ++num_samples;
-  }
-
-  void extend_sequence_even(float2 points[], int N)
-  {
-    int n = (int)sqrtf(N);
-    occupied1Dx.resize(2 * N);
-    occupied1Dy.resize(2 * N);
-    mark_occupied_strata(points, N);
-    for (int s = 0; s < N; ++s) {
-      float2 oldpt = points[s];
-      float i = floorf(n * oldpt.x);
-      float j = floorf(n * oldpt.y);
-      float xhalf = floorf(2.0f * (n * oldpt.x - i));
-      float yhalf = floorf(2.0f * (n * oldpt.y - j));
-      xhalf = 1.0f - xhalf;
-      yhalf = 1.0f - yhalf;
-      generate_sample_point(points, i, j, xhalf, yhalf, n, N);
-    }
-  }
-
-  void extend_sequence_odd(float2 points[], int N)
-  {
-    int n = (int)sqrtf(N / 2);
-    occupied1Dx.resize(2 * N);
-    occupied1Dy.resize(2 * N);
-    mark_occupied_strata(points, N);
-    std::vector<float> xhalves(N / 2);
-    std::vector<float> yhalves(N / 2);
-    for (int s = 0; s < N / 2; ++s) {
-      float2 oldpt = points[s];
-      float i = floorf(n * oldpt.x);
-      float j = floorf(n * oldpt.y);
-      float xhalf = floorf(2.0f * (n * oldpt.x - i));
-      float yhalf = floorf(2.0f * (n * oldpt.y - j));
-      if (rnd() > 0.5f) {
-        xhalf = 1.0f - xhalf;
-      }
-      else {
-        yhalf = 1.0f - yhalf;
-      }
-      xhalves[s] = xhalf;
-      yhalves[s] = yhalf;
-      generate_sample_point(points, i, j, xhalf, yhalf, n, N);
-    }
-    for (int s = 0; s < N / 2; ++s) {
-      float2 oldpt = points[s];
-      float i = floorf(n * oldpt.x);
-      float j = floorf(n * oldpt.y);
-      float xhalf = 1.0f - xhalves[s];
-      float yhalf = 1.0f - yhalves[s];
-      generate_sample_point(points, i, j, xhalf, yhalf, n, N);
-    }
-  }
-
-  std::vector<bool> occupied1Dx, occupied1Dy;
-  int num_samples;
-  int rnd_index, rnd_seed;
-};
-
-class PMJ02_Generator : public PMJ_Generator {
- protected:
-  void generate_sample_point(
-      float2 points[], float i, float j, float xhalf, float yhalf, int n, int N) override
-  {
-    int NN = 2 * N;
-    float2 pt;
-    do {
-      pt.x = (i + 0.5f * (xhalf + rnd())) / n;
-      pt.y = (j + 0.5f * (yhalf + rnd())) / n;
-    } while (is_occupied(pt, NN));
-    mark_occupied_strata1(pt, NN);
-    points[num_samples] = pt;
-    ++num_samples;
-  }
-
-  void mark_occupied_strata(float2 points[], int N) override
-  {
-    int NN = 2 * N;
-    int num_shapes = (int)log2f(NN) + 1;
-    occupiedStrata.resize(num_shapes);
-    for (int shape = 0; shape < num_shapes; ++shape) {
-      occupiedStrata[shape].resize(NN);
-      for (int n = 0; n < NN; ++n) {
-        occupiedStrata[shape][n] = false;
-      }
-    }
-    for (int s = 0; s < N; ++s) {
-      mark_occupied_strata1(points[s], NN);
-    }
-  }
-
-  void mark_occupied_strata1(float2 pt, int NN)
-  {
-    int shape = 0;
-    int xdivs = NN;
-    int ydivs = 1;
-    do {
-      int xstratum = (int)(xdivs * pt.x);
-      int ystratum = (int)(ydivs * pt.y);
-      size_t index = ystratum * xdivs + xstratum;
-      assert(index < NN);
-      occupiedStrata[shape][index] = true;
-      shape = shape + 1;
-      xdivs = xdivs / 2;
-      ydivs = ydivs * 2;
-    } while (xdivs > 0);
-  }
-
-  bool is_occupied(float2 pt, int NN)
-  {
-    int shape = 0;
-    int xdivs = NN;
-    int ydivs = 1;
-    do {
-      int xstratum = (int)(xdivs * pt.x);
-      int ystratum = (int)(ydivs * pt.y);
-      size_t index = ystratum * xdivs + xstratum;
-      assert(index < NN);
-      if (occupiedStrata[shape][index]) {
-        return true;
-      }
-      shape = shape + 1;
-      xdivs = xdivs / 2;
-      ydivs = ydivs * 2;
-    } while (xdivs > 0);
-    return false;
-  }
-
- private:
-  std::vector<std::vector<bool>> occupiedStrata;
-};
-
-static void shuffle(float2 points[], int size, int rng_seed)
+void progressive_multi_jitter_02_generate_2D(float2 points[], int size, int rng_seed)
 {
-  if (rng_seed == 0) {
-    return;
-  }
-
-  constexpr int odd[8] = {0, 1, 4, 5, 10, 11, 14, 15};
-  constexpr int even[8] = {2, 3, 6, 7, 8, 9, 12, 13};
-
-  int rng_index = 0;
-  for (int yy = 0; yy < size / 16; ++yy) {
-    for (int xx = 0; xx < 8; ++xx) {
-      int other = (int)(cmj_randfloat(++rng_index, rng_seed) * (8.0f - xx) + xx);
-      float2 tmp = points[odd[other] + yy * 16];
-      points[odd[other] + yy * 16] = points[odd[xx] + yy * 16];
-      points[odd[xx] + yy * 16] = tmp;
-    }
-    for (int xx = 0; xx < 8; ++xx) {
-      int other = (int)(cmj_randfloat(++rng_index, rng_seed) * (8.0f - xx) + xx);
-      float2 tmp = points[even[other] + yy * 16];
-      points[even[other] + yy * 16] = points[even[xx] + yy * 16];
-      points[even[xx] + yy * 16] = tmp;
+  /* Xor values for generating the PMJ02 sequence.  These permute the
+   * order we visit the strata in, which is what makes the code below
+   * produce the PMJ02 sequence.  Other choices are also possible, but
+   * result in different sequences. */
+  static uint xors[2][32] = {
+      {0x00000000, 0x00000000, 0x00000002, 0x00000006, 0x00000006, 0x0000000e, 0x00000036,
+       0x0000004e, 0x00000016, 0x0000002e, 0x00000276, 0x000006ce, 0x00000716, 0x00000c2e,
+       0x00003076, 0x000040ce, 0x00000116, 0x0000022e, 0x00020676, 0x00060ece, 0x00061716,
+       0x000e2c2e, 0x00367076, 0x004ec0ce, 0x00170116, 0x002c022e, 0x02700676, 0x06c00ece,
+       0x07001716, 0x0c002c2e, 0x30007076, 0x4000c0ce},
+      {0x00000000, 0x00000001, 0x00000003, 0x00000003, 0x00000007, 0x0000001b, 0x00000027,
+       0x0000000b, 0x00000017, 0x0000013b, 0x00000367, 0x0000038b, 0x00000617, 0x0000183b,
+       0x00002067, 0x0000008b, 0x00000117, 0x0001033b, 0x00030767, 0x00030b8b, 0x00071617,
+       0x001b383b, 0x00276067, 0x000b808b, 0x00160117, 0x0138033b, 0x03600767, 0x03800b8b,
+       0x06001617, 0x1800383b, 0x20006067, 0x0000808b}};
+
+  uint rng_i = rng_seed;
+
+  points[0].x = hash_hp_float(rng_i++);
+  points[0].y = hash_hp_float(rng_i++);
+
+  /* Subdivide the domain into smaller and smaller strata, filling in new
+   * points as we go. */
+  for (int log_N = 0, N = 1; N < size; log_N++, N *= 2) {
+    float strata_count = (float)(N * 2);
+    for (int i = 0; i < N && (N + i) < size; i++) {
+      /* Find the strata that are already occupied in this cell. */
+      uint occupied_x_stratum = (uint)(points[i ^ xors[0][log_N]].x * strata_count);
+      uint occupied_y_stratum = (uint)(points[i ^ xors[1][log_N]].y * strata_count);
+
+      /* Generate a new point in the unoccupied strata. */
+      points[N + i].x = ((float)(occupied_x_stratum ^ 1) + hash_hp_float(rng_i++)) / strata_count;
+      points[N + i].y = ((float)(occupied_y_stratum ^ 1) + hash_hp_float(rng_i++)) / strata_count;
     }
   }
 }
 
-void progressive_multi_jitter_generate_2D(float2 points[], int size, int rng_seed)
-{
-  PMJ_Generator::generate_2D(points, size, rng_seed);
-  shuffle(points, size, rng_seed);
-}
-
-void progressive_multi_jitter_02_generate_2D(float2 points[], int size, int rng_seed)
-{
-  PMJ02_Generator::generate_2D(points, size, rng_seed);
-  shuffle(points, size, rng_seed);
-}
-
 CCL_NAMESPACE_END
diff --git a/intern/cycles/scene/jitter.h b/intern/cycles/scene/jitter.h
index 2fda0e556c0..8d497d5e9f5 100644
--- a/intern/cycles/scene/jitter.h
+++ b/intern/cycles/scene/jitter.h
@@ -8,7 +8,6 @@
 
 CCL_NAMESPACE_BEGIN
 
-void progressive_multi_jitter_generate_2D(float2 points[], int size, int rng_seed);
 void progressive_multi_jitter_02_generate_2D(float2 points[], int size, int rng_seed);
 
 CCL_NAMESPACE_END
diff --git a/intern/cycles/util/hash.h b/intern/cycles/util/hash.h
index 351b8796be7..4f83f331229 100644
--- a/intern/cycles/util/hash.h
+++ b/intern/cycles/util/hash.h
@@ -4,6 +4,7 @@
 #ifndef __UTIL_HASH_H__
 #define __UTIL_HASH_H__
 
+#include "util/math.h"
 #include "util/types.h"
 
 CCL_NAMESPACE_BEGIN
@@ -407,42 +408,16 @@ ccl_device_inline uint hash_hp_seeded_uint(uint i, uint seed)
   return hash_hp_uint(i ^ seed);
 }
 
-/* Outputs [0.0, 1.0]. */
-ccl_device_inline float hash_hp_seeded_float(uint i, uint seed)
+/* Outputs [0.0, 1.0). */
+ccl_device_inline float hash_hp_float(uint i)
 {
-  return uint_to_float_incl(hash_hp_seeded_uint(i, seed));
+  return uint_to_float_excl(hash_hp_uint(i));
 }
 
-/* ***** CMJ Hash Functions *****
- *
- * These are based on one of the hash functions in the paper
- * "Correlated Multi-Jittered Sampling" by Andrew Kensler, 2013.
- *
- * These are here for backwards-compatibility, and can be replaced
- * by the Hash Prospector hashes above at some point.
- * See https://developer.blender.org/D15679#426304
- */
-
-ccl_device_inline uint hash_cmj_seeded_uint(uint i, uint seed)
-{
-  i ^= seed;
-  i ^= i >> 17;
-  i ^= i >> 10;
-  i *= 0xb36534e5;
-  i ^= i >> 12;
-  i ^= i >> 21;
-  i *= 0x93fc4795;
-  i ^= 0xdf6e307f;
-  i ^= i >> 17;
-  i *= 1 | seed >> 18;
-
-  return i;
-}
-
-/* Outputs [0.0, 1.0]. */
-ccl_device_inline float hash_cmj_seeded_float(uint i, uint seed)
+/* Outputs [0.0, 1.0). */
+ccl_device_inline float hash_hp_seeded_float(uint i, uint seed)
 {
-  return uint_to_float_excl(hash_cmj_seeded_uint(i, seed));
+  return uint_to_float_excl(hash_hp_seeded_uint(i, seed));
 }
 
 /* ***** Modified Wang Hash Functions *****
@@ -463,10 +438,58 @@ ccl_device_inline uint hash_wang_seeded_uint(uint i, uint seed)
   return i;
 }
 
-/* Outputs [0.0, 1.0]. */
+/* Outputs [0.0, 1.0). */
 ccl_device_inline float hash_wang_seeded_float(uint i, uint seed)
 {
-  return uint_to_float_incl(hash_wang_seeded_uint(i, seed));
+  return uint_to_float_excl(hash_wang_seeded_uint(i, seed));
+}
+
+/* ***** Index Shuffling Hash Function *****
+ *
+ * This function takes an index, the length of the thing the index points
+ * into, and returns a shuffled index.  For example, if you pass indices
+ * 0 through 19 to this function with a length parameter of 20, it will
+ * return the indices in a shuffled order with no repeats.  Indices
+ * larger than the length parameter will simply repeat the same shuffled
+ * pattern over and over.
+ *
+ * This is useful for iterating over an array in random shuffled order
+ * without having to shuffle the array itself.
+ *
+ * Passing different seeds results in different random shuffles.
+ *
+ * This function runs in average O(1) time.
+ *
+ * See https://andrew-helmer.github.io/permute/ for details on how this
+ * works.
+ */
+ccl_device_inline uint hash_shuffle_uint(uint i, uint length, uint seed)
+{
+  i = i % length;
+  uint mask = (1 << (32 - count_leading_zeros(length - 1))) - 1;
+
+  do {
+    i ^= seed;
+    i *= 0xe170893d;
+    i ^= seed >> 16;
+    i ^= (i & mask) >> 4;
+    i ^= seed >> 8;
+    i *= 0x0929eb3f;
+    i ^= seed >> 23;
+    i ^= (i & mask) >> 1;
+    i *= 1 | seed >> 27;
+    i *= 0x6935fa69;
+    i ^= (i & mask) >> 11;
+    i *= 0x74dcb303;
+    i ^= (i & mask) >> 2;
+    i *= 0x9e501cc3;
+    i ^= (i & mask) >> 2;
+    i *= 0xc860a3df;
+    i &= mask;
+    i ^= i >> 5;
+  } while (i >= length);
+
+  return i;
 }
 
 /* ********** */