Cycles: add Sobol-Burley sampling pattern

Based on the paper "Practical Hash-based Owen Scrambling" by Brent Burley,
2020, Journal of Computer Graphics Techniques.

It is distinct from the existing Sobol sampler in two important ways:
* It is Owen scrambled, which gives it a much better convergence rate in many
  situations.
* It uses padding for higher dimensions, rather than using higher Sobol
  dimensions directly. In practice this is advantagous because high-dimensional
  Sobol sequences have holes in their sampling patterns that don't resolve
  until an unreasonable number of samples are taken. (See Burley's paper for
  details.)

The pattern reduces noise in some benchmark scenes, however it is also slower,
particularly on the CPU. So for now Progressive Multi-Jittered sampling remains
the default.

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

4 files changed, 213 insertions, 80 deletions

diff --git a/intern/cycles/kernel/sample/jitter.h b/intern/cycles/kernel/sample/jitter.h
index b5cfa624406..dd170cf2120 100644
--- a/intern/cycles/kernel/sample/jitter.h
+++ b/intern/cycles/kernel/sample/jitter.h
@@ -1,20 +1,12 @@
 /* SPDX-License-Identifier: Apache-2.0
  * Copyright 2011-2022 Blender Foundation */
 
+#include "kernel/sample/util.h"
+#include "util/hash.h"
+
 #pragma once
 CCL_NAMESPACE_BEGIN
 
-ccl_device_inline uint32_t laine_karras_permutation(uint32_t x, uint32_t seed)
-{
-  x += seed;
-  x ^= (x * 0x6c50b47cu);
-  x ^= x * 0xb82f1e52u;
-  x ^= x * 0xc7afe638u;
-  x ^= x * 0x8d22f6e6u;
-
-  return x;
-}
-
 ccl_device_inline uint32_t nested_uniform_scramble(uint32_t x, uint32_t seed)
 {
   x = reverse_integer_bits(x);
@@ -24,46 +16,6 @@ ccl_device_inline uint32_t nested_uniform_scramble(uint32_t x, uint32_t seed)
   return x;
 }
 
-ccl_device_inline 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;
-}
-
-ccl_device_inline uint cmj_hash_simple(uint i, uint p)
-{
-  i = (i ^ 61) ^ p;
-  i += i << 3;
-  i ^= i >> 4;
-  i *= 0x27d4eb2d;
-  return i;
-}
-
-ccl_device_inline float cmj_randfloat(uint i, uint p)
-{
-  return cmj_hash(i, p) * (1.0f / 4294967808.0f);
-}
-
-ccl_device_inline float cmj_randfloat_simple(uint i, uint p)
-{
-  return cmj_hash_simple(i, p) * (1.0f / (float)0xFFFFFFFF);
-}
-
-ccl_device_inline float cmj_randfloat_simple_dist(uint i, uint p, float d)
-{
-  return cmj_hash_simple(i, p) * (d / (float)0xFFFFFFFF);
-}
-
 ccl_device float pmj_sample_1D(KernelGlobals kg, uint sample, uint rng_hash, uint dimension)
 {
   uint hash = rng_hash;
@@ -71,16 +23,12 @@ ccl_device float pmj_sample_1D(KernelGlobals kg, uint sample, uint rng_hash, uin
   if (kernel_data.integrator.scrambling_distance < 1.0f) {
     hash = kernel_data.integrator.seed;
 
-    jitter_x = cmj_randfloat_simple_dist(
-        dimension, rng_hash, kernel_data.integrator.scrambling_distance);
+    jitter_x = hash_wang_seeded_float(dimension, rng_hash) *
+               kernel_data.integrator.scrambling_distance;
   }
 
   /* Perform Owen shuffle of the sample number to reorder the samples. */
-#ifdef _SIMPLE_HASH_
-  const uint rv = cmj_hash_simple(dimension, hash);
-#else /* Use a _REGULAR_HASH_. */
-  const uint rv = cmj_hash(dimension, hash);
-#endif
+  const uint rv = hash_cmj_seeded_uint(dimension, hash);
 #ifdef _XOR_SHUFFLE_
 #  warning "Using XOR shuffle."
   const uint s = sample ^ rv;
@@ -101,11 +49,7 @@ ccl_device float pmj_sample_1D(KernelGlobals kg, uint sample, uint rng_hash, uin
 
 #ifndef _NO_CRANLEY_PATTERSON_ROTATION_
   /* Use Cranley-Patterson rotation to displace the sample pattern. */
-#  ifdef _SIMPLE_HASH_
-  float dx = cmj_randfloat_simple(d, hash);
-#  else
-  float dx = cmj_randfloat(d, hash);
-#  endif
+  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);
@@ -129,18 +73,14 @@ ccl_device void pmj_sample_2D(KernelGlobals kg,
   if (kernel_data.integrator.scrambling_distance < 1.0f) {
     hash = kernel_data.integrator.seed;
 
-    jitter_x = cmj_randfloat_simple_dist(
-        dimension, rng_hash, kernel_data.integrator.scrambling_distance);
-    jitter_y = cmj_randfloat_simple_dist(
-        dimension + 1, rng_hash, kernel_data.integrator.scrambling_distance);
+    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;
   }
 
   /* Perform a shuffle on the sample number to reorder the samples. */
-#ifdef _SIMPLE_HASH_
-  const uint rv = cmj_hash_simple(dimension, hash);
-#else /* Use a _REGULAR_HASH_. */
-  const uint rv = cmj_hash(dimension, hash);
-#endif
+  const uint rv = hash_cmj_seeded_uint(dimension, hash);
 #ifdef _XOR_SHUFFLE_
 #  warning "Using XOR shuffle."
   const uint s = sample ^ rv;
@@ -159,13 +99,8 @@ ccl_device void pmj_sample_2D(KernelGlobals kg,
 
 #ifndef _NO_CRANLEY_PATTERSON_ROTATION_
   /* Use Cranley-Patterson rotation to displace the sample pattern. */
-#  ifdef _SIMPLE_HASH_
-  float dx = cmj_randfloat_simple(d, hash);
-  float dy = cmj_randfloat_simple(d + 1, hash);
-#  else
-  float dx = cmj_randfloat(d, hash);
-  float dy = cmj_randfloat(d + 1, hash);
-#  endif
+  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;
diff --git a/intern/cycles/kernel/sample/pattern.h b/intern/cycles/kernel/sample/pattern.h
index 89500d51872..e8c3acb5cf7 100644
--- a/intern/cycles/kernel/sample/pattern.h
+++ b/intern/cycles/kernel/sample/pattern.h
@@ -4,6 +4,7 @@
 #pragma once
 
 #include "kernel/sample/jitter.h"
+#include "kernel/sample/sobol_burley.h"
 #include "util/hash.h"
 
 CCL_NAMESPACE_BEGIN
@@ -48,6 +49,10 @@ ccl_device_forceinline float path_rng_1D(KernelGlobals kg,
   return (float)drand48();
 #endif
 
+  if (kernel_data.integrator.sampling_pattern == SAMPLING_PATTERN_SOBOL_BURLEY) {
+    return sobol_burley_sample_1D(sample, dimension, rng_hash);
+  }
+
 #ifdef __SOBOL__
   if (kernel_data.integrator.sampling_pattern == SAMPLING_PATTERN_PMJ)
 #endif
@@ -66,7 +71,7 @@ ccl_device_forceinline float path_rng_1D(KernelGlobals kg,
   /* Hash rng with dimension to solve correlation issues.
    * See T38710, T50116.
    */
-  uint tmp_rng = cmj_hash_simple(dimension, rng_hash);
+  uint tmp_rng = hash_wang_seeded_uint(dimension, rng_hash);
   shift = tmp_rng * (kernel_data.integrator.scrambling_distance / (float)0xFFFFFFFF);
 
   return r + shift - floorf(r + shift);
@@ -86,6 +91,11 @@ ccl_device_forceinline void path_rng_2D(KernelGlobals kg,
   return;
 #endif
 
+  if (kernel_data.integrator.sampling_pattern == SAMPLING_PATTERN_SOBOL_BURLEY) {
+    sobol_burley_sample_2D(sample, dimension, rng_hash, fx, fy);
+    return;
+  }
+
 #ifdef __SOBOL__
   if (kernel_data.integrator.sampling_pattern == SAMPLING_PATTERN_PMJ)
 #endif
diff --git a/intern/cycles/kernel/sample/sobol_burley.h b/intern/cycles/kernel/sample/sobol_burley.h
new file mode 100644
index 00000000000..4e041aa075e
--- /dev/null
+++ b/intern/cycles/kernel/sample/sobol_burley.h
@@ -0,0 +1,143 @@
+/* SPDX-License-Identifier: Apache-2.0
+ * Copyright 2011-2022 Blender Foundation */
+
+/*
+ * A shuffled, Owen-scrambled Sobol sampler, implemented with the
+ * techniques from the paper "Practical Hash-based Owen Scrambling"
+ * by Brent Burley, 2020, Journal of Computer Graphics Techniques.
+ *
+ * Note that unlike a standard high-dimensional Sobol sequence, this
+ * Sobol sampler uses padding to achieve higher dimensions, as described
+ * in Burley's paper.
+ */
+
+#pragma once
+
+#include "kernel/sample/util.h"
+#include "util/hash.h"
+#include "util/math.h"
+#include "util/types.h"
+
+CCL_NAMESPACE_BEGIN
+
+/*
+ * Computes a single dimension of a sample from an Owen-scrambled
+ * Sobol sequence.  This is used in the main sampling functions,
+ * sobol_burley_sample_#D(), below.
+ *
+ * - rev_bit_index: the sample index, with reversed order bits.
+ * - dimension:     the sample dimension.
+ * - scramble_seed: the Owen scrambling seed.
+ *
+ * Note that the seed must be well randomized before being
+ * passed to this function.
+ */
+ccl_device_forceinline float sobol_burley(uint rev_bit_index, uint dimension, uint scramble_seed)
+{
+  uint result = 0;
+
+  if (dimension == 0) {
+    // Fast-path for dimension 0, which is just Van der corput.
+    // This makes a notable difference in performance since we reuse
+    // dimensions for padding, and dimension 0 is reused the most.
+    result = reverse_integer_bits(rev_bit_index);
+  }
+  else {
+    uint i = 0;
+    while (rev_bit_index != 0) {
+      uint j = count_leading_zeros(rev_bit_index);
+      result ^= sobol_burley_table[dimension][i + j];
+      i += j + 1;
+
+      // We can't do "<<= j + 1" because that can overflow the shift
+      // operator, which doesn't do what we need on at least x86.
+      rev_bit_index <<= j;
+      rev_bit_index <<= 1;
+    }
+  }
+
+  // Apply Owen scrambling.
+  result = reverse_integer_bits(reversed_bit_owen(result, scramble_seed));
+
+  return uint_to_float_excl(result);
+}
+
+/*
+ * Computes a 1D Owen-scrambled and shuffled Sobol sample.
+ */
+ccl_device float sobol_burley_sample_1D(uint index, uint dimension, uint seed)
+{
+  // Include the dimension in the seed, so we get decorrelated
+  // sequences for different dimensions via shuffling.
+  seed ^= hash_hp_uint(dimension);
+
+  // Shuffle.
+  index = reversed_bit_owen(reverse_integer_bits(index), seed ^ 0xbff95bfe);
+
+  return sobol_burley(index, 0, seed ^ 0x635c77bd);
+}
+
+/*
+ * Computes a 2D Owen-scrambled and shuffled Sobol sample.
+ */
+ccl_device void sobol_burley_sample_2D(
+    uint index, uint dimension_set, uint seed, ccl_private float *x, ccl_private float *y)
+{
+  // Include the dimension set in the seed, so we get decorrelated
+  // sequences for different dimension sets via shuffling.
+  seed ^= hash_hp_uint(dimension_set);
+
+  // Shuffle.
+  index = reversed_bit_owen(reverse_integer_bits(index), seed ^ 0xf8ade99a);
+
+  *x = sobol_burley(index, 0, seed ^ 0xe0aaaf76);
+  *y = sobol_burley(index, 1, seed ^ 0x94964d4e);
+}
+
+/*
+ * Computes a 3D Owen-scrambled and shuffled Sobol sample.
+ */
+ccl_device void sobol_burley_sample_3D(uint index,
+                                       uint dimension_set,
+                                       uint seed,
+                                       ccl_private float *x,
+                                       ccl_private float *y,
+                                       ccl_private float *z)
+{
+  // Include the dimension set in the seed, so we get decorrelated
+  // sequences for different dimension sets via shuffling.
+  seed ^= hash_hp_uint(dimension_set);
+
+  // Shuffle.
+  index = reversed_bit_owen(reverse_integer_bits(index), seed ^ 0xcaa726ac);
+
+  *x = sobol_burley(index, 0, seed ^ 0x9e78e391);
+  *y = sobol_burley(index, 1, seed ^ 0x67c33241);
+  *z = sobol_burley(index, 2, seed ^ 0x78c395c5);
+}
+
+/*
+ * Computes a 4D Owen-scrambled and shuffled Sobol sample.
+ */
+ccl_device void sobol_burley_sample_4D(uint index,
+                                       uint dimension_set,
+                                       uint seed,
+                                       ccl_private float *x,
+                                       ccl_private float *y,
+                                       ccl_private float *z,
+                                       ccl_private float *w)
+{
+  // Include the dimension set in the seed, so we get decorrelated
+  // sequences for different dimension sets via shuffling.
+  seed ^= hash_hp_uint(dimension_set);
+
+  // Shuffle.
+  index = reversed_bit_owen(reverse_integer_bits(index), seed ^ 0xc2c1a055);
+
+  *x = sobol_burley(index, 0, seed ^ 0x39468210);
+  *y = sobol_burley(index, 1, seed ^ 0xe9d8a845);
+  *z = sobol_burley(index, 2, seed ^ 0x5f32b482);
+  *w = sobol_burley(index, 3, seed ^ 0x1524cc56);
+}
+
+CCL_NAMESPACE_END
diff --git a/intern/cycles/kernel/sample/util.h b/intern/cycles/kernel/sample/util.h
new file mode 100644
index 00000000000..33056bb7819
--- /dev/null
+++ b/intern/cycles/kernel/sample/util.h
@@ -0,0 +1,45 @@
+/* SPDX-License-Identifier: Apache-2.0
+ * Copyright 2011-2022 Blender Foundation */
+
+#pragma once
+
+#include "util/types.h"
+
+CCL_NAMESPACE_BEGIN
+
+/*
+ * Performs base-2 Owen scrambling on a reversed-bit 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
+ */
+ccl_device_inline uint reversed_bit_owen(uint n, uint seed)
+{
+  n ^= n * 0x3d20adea;
+  n += seed;
+  n *= (seed >> 16) | 1;
+  n ^= n * 0x05526c56;
+  n ^= n * 0x53a22864;
+
+  return n;
+}
+
+/*
+ * 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
+ */
+ccl_device_inline uint laine_karras_permutation(uint x, uint seed)
+{
+  x += seed;
+  x ^= (x * 0x6c50b47cu);
+  x ^= x * 0xb82f1e52u;
+  x ^= x * 0xc7afe638u;
+  x ^= x * 0x8d22f6e6u;
+
+  return x;
+}
+
+CCL_NAMESPACE_END