Cycles: add Progressive Multi-Jitter sampling pattern

This sampling pattern is particularly suited to adaptive sampling, and will
be used for that upcoming feature.

Based on "Progressive Multi-Jittered Sample Sequences" by Per Christensen,
Andrew Kensler and Charlie Kilpatrick.

Ref D4686

intern/cycles/blender/addon/properties.py

@@ -112,6 +112,7 @@ enum_use_layer_samples = (
 enum_sampling_pattern = (
     ('SOBOL', "Sobol", "Use Sobol random sampling pattern"),
     ('CORRELATED_MUTI_JITTER', "Correlated Multi-Jitter", "Use Correlated Multi-Jitter random sampling pattern"),
+    ('PROGRESSIVE_MUTI_JITTER', "Progressive Multi-Jitter", "Use Progressive Multi-Jitter random sampling pattern"),
 )
 
 enum_integrator = (

intern/cycles/kernel/kernel_jitter.h

@@ -195,4 +195,35 @@ ccl_device void cmj_sample_2D(int s, int N, int p, float *fx, float *fy)
 }
 #endif
 
+ccl_device float pmj_sample_1D(KernelGlobals *kg, int sample, int rng_hash, int dimension)
+{
+  /* Fallback to random */
+  if (sample > NUM_PMJ_SAMPLES) {
+    int p = rng_hash + dimension;
+    return cmj_randfloat(sample, p);
+  }
+  uint tmp_rng = cmj_hash_simple(dimension, rng_hash);
+  int index = ((dimension % NUM_PMJ_PATTERNS) * NUM_PMJ_SAMPLES + sample) * 2;
+  return __uint_as_float(kernel_tex_fetch(__sample_pattern_lut, index) ^ (tmp_rng & 0x007fffff)) -
+         1.0f;
+}
+
+ccl_device void pmj_sample_2D(
+    KernelGlobals *kg, int sample, int rng_hash, int dimension, float *fx, float *fy)
+{
+  if (sample > NUM_PMJ_SAMPLES) {
+    int p = rng_hash + dimension;
+    *fx = cmj_randfloat(sample, p);
+    *fy = cmj_randfloat(sample, p + 1);
+  }
+  uint tmp_rng = cmj_hash_simple(dimension, rng_hash);
+  int index = ((dimension % NUM_PMJ_PATTERNS) * NUM_PMJ_SAMPLES + sample) * 2;
+  *fx = __uint_as_float(kernel_tex_fetch(__sample_pattern_lut, index) ^ (tmp_rng & 0x007fffff)) -
+        1.0f;
+  tmp_rng = cmj_hash_simple(dimension + 1, rng_hash);
+  *fy = __uint_as_float(kernel_tex_fetch(__sample_pattern_lut, index + 1) ^
+                        (tmp_rng & 0x007fffff)) -
+        1.0f;
+}
+
 CCL_NAMESPACE_END

intern/cycles/kernel/kernel_random.h

@@ -43,7 +43,7 @@ ccl_device uint sobol_dimension(KernelGlobals *kg, int index, int dimension)
   uint i = index + SOBOL_SKIP;
   for (int j = 0, x; (x = find_first_set(i)); i >>= x) {
     j += x;
-    result ^= kernel_tex_fetch(__sobol_directions, 32 * dimension + j - 1);
+    result ^= kernel_tex_fetch(__sample_pattern_lut, 32 * dimension + j - 1);
   }
   return result;
 }
@@ -56,7 +56,9 @@ ccl_device_forceinline float path_rng_1D(
 #ifdef __DEBUG_CORRELATION__
   return (float)drand48();
 #endif
-
+  if (kernel_data.integrator.sampling_pattern == SAMPLING_PATTERN_PMJ) {
+    return pmj_sample_1D(kg, sample, rng_hash, dimension);
+  }
 #ifdef __CMJ__
 #  ifdef __SOBOL__
   if (kernel_data.integrator.sampling_pattern == SAMPLING_PATTERN_CMJ)
@@ -99,7 +101,10 @@ ccl_device_forceinline void path_rng_2D(KernelGlobals *kg,
   *fy = (float)drand48();
   return;
 #endif
-
+  if (kernel_data.integrator.sampling_pattern == SAMPLING_PATTERN_PMJ) {
+    pmj_sample_2D(kg, sample, rng_hash, dimension, fx, fy);
+    return;
+  }
 #ifdef __CMJ__
 #  ifdef __SOBOL__
   if (kernel_data.integrator.sampling_pattern == SAMPLING_PATTERN_CMJ)
@@ -284,4 +289,28 @@ ccl_device float lcg_step_float_addrspace(ccl_addr_space uint *rng)
   return (float)*rng * (1.0f / (float)0xFFFFFFFF);
 }
 
+ccl_device_inline bool sample_is_even(int pattern, int sample)
+{
+  if (pattern == SAMPLING_PATTERN_PMJ) {
+    /* See Section 10.2.1, "Progressive Multi-Jittered Sample Sequences", Christensen et al.
+     * We can use this to get divide sample sequence into two classes for easier variance
+     * estimation. There must be a more elegant way of writing this? */
+#if defined(__GNUC__) && !defined(__KERNEL_GPU__)
+    return __builtin_popcount(sample & 0xaaaaaaaa) & 1;
+#elif defined(__NVCC__)
+    return __popc(sample & 0xaaaaaaaa) & 1;
+#else
+    int i = sample & 0xaaaaaaaa;
+    i = i - ((i >> 1) & 0x55555555);
+    i = (i & 0x33333333) + ((i >> 2) & 0x33333333);
+    i = (((i + (i >> 4)) & 0xF0F0F0F) * 0x1010101) >> 24;
+    return i & 1;
+#endif
+  }
+  else {
+    /* TODO(Stefan): Are there reliable ways of dividing CMJ and Sobol into two classes? */
+    return sample & 0x1;
+  }
+}
+
 CCL_NAMESPACE_END

intern/cycles/kernel/kernel_textures.h

@@ -77,7 +77,7 @@ KERNEL_TEX(KernelShader, __shaders)
 KERNEL_TEX(float, __lookup_table)
 
 /* sobol */
-KERNEL_TEX(uint, __sobol_directions)
+KERNEL_TEX(uint, __sample_pattern_lut)
 
 /* image textures */
 KERNEL_TEX(TextureInfo, __texture_info)

intern/cycles/kernel/kernel_types.h

@@ -267,6 +267,7 @@ enum PathTraceDimension {
 enum SamplingPattern {
   SAMPLING_PATTERN_SOBOL = 0,
   SAMPLING_PATTERN_CMJ = 1,
+  SAMPLING_PATTERN_PMJ = 2,
 
   SAMPLING_NUM_PATTERNS,
 };
@@ -1667,6 +1668,10 @@ typedef struct WorkTile {
   ccl_global float *buffer;
 } WorkTile;
 
+/* Precoumputed sample table sizes for PMJ02 sampler. */
+#define NUM_PMJ_SAMPLES 64 * 64
+#define NUM_PMJ_PATTERNS 48
+
 CCL_NAMESPACE_END
 
 #endif /*  __KERNEL_TYPES_H__ */

intern/cycles/render/CMakeLists.txt

@@ -24,6 +24,7 @@ set(SRC
   hair.cpp
   image.cpp
   integrator.cpp
+  jitter.cpp
   light.cpp
   merge.cpp
   mesh.cpp
@@ -62,6 +63,7 @@ set(SRC_HEADERS
   image.h
   integrator.h
   light.h
+  jitter.h
   merge.h
   mesh.h
   nodes.h

intern/cycles/render/integrator.cpp

@@ -18,11 +18,14 @@
 #include "render/background.h"
 #include "render/integrator.h"
 #include "render/film.h"
+#include "render/jitter.h"
 #include "render/light.h"
 #include "render/scene.h"
 #include "render/shader.h"
 #include "render/sobol.h"
 
+#include "kernel/kernel_types.h"
+
 #include "util/util_foreach.h"
 #include "util/util_hash.h"
 
@@ -78,6 +81,7 @@ NODE_DEFINE(Integrator)
   static NodeEnum sampling_pattern_enum;
   sampling_pattern_enum.insert("sobol", SAMPLING_PATTERN_SOBOL);
   sampling_pattern_enum.insert("cmj", SAMPLING_PATTERN_CMJ);
+  sampling_pattern_enum.insert("pmj", SAMPLING_PATTERN_PMJ);
   SOCKET_ENUM(sampling_pattern, "Sampling Pattern", sampling_pattern_enum, SAMPLING_PATTERN_SOBOL);
 
   return type;
@@ -203,18 +207,34 @@ void Integrator::device_update(Device *device, DeviceScene *dscene, Scene *scene
   int dimensions = PRNG_BASE_NUM + max_samples * PRNG_BOUNCE_NUM;
   dimensions = min(dimensions, SOBOL_MAX_DIMENSIONS);
 
-  uint *directions = dscene->sobol_directions.alloc(SOBOL_BITS * dimensions);
+  if (sampling_pattern == SAMPLING_PATTERN_SOBOL) {
+    uint *directions = dscene->sample_pattern_lut.alloc(SOBOL_BITS * dimensions);
 
-  sobol_generate_direction_vectors((uint(*)[SOBOL_BITS])directions, dimensions);
+    sobol_generate_direction_vectors((uint(*)[SOBOL_BITS])directions, dimensions);
 
-  dscene->sobol_directions.copy_to_device();
+    dscene->sample_pattern_lut.copy_to_device();
+  }
+  else {
+    constexpr int sequence_size = NUM_PMJ_SAMPLES;
+    constexpr int num_sequences = NUM_PMJ_PATTERNS;
+    float2 *directions = (float2 *)dscene->sample_pattern_lut.alloc(sequence_size * num_sequences *
+                                                                    2);
+    TaskPool pool;
+    for (int j = 0; j < num_sequences; ++j) {
+      float2 *sequence = directions + j * sequence_size;
+      pool.push(
+          function_bind(&progressive_multi_jitter_02_generate_2D, sequence, sequence_size, j));
+    }
+    pool.wait_work();
+    dscene->sample_pattern_lut.copy_to_device();
+  }
 
   need_update = false;
 }
 
 void Integrator::device_free(Device *, DeviceScene *dscene)
 {
-  dscene->sobol_directions.free();
+  dscene->sample_pattern_lut.free();
 }
 
 bool Integrator::modified(const Integrator &integrator)

intern/cycles/render/jitter.cpp

@@ -0,0 +1,287 @@
+/*
+ * Copyright 2019 Blender Foundation
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/* 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"
+ */
+
+#include "render/jitter.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)
+{
+  /* Offset samples by 1.0 for faster scrambling in kernel_random.h */
+  for (int i = 0; i < size; ++i) {
+    points[i].x += 1.0f;
+    points[i].y += 1.0f;
+  }
+
+  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;
+    }
+  }
+}
+
+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

intern/cycles/render/jitter.h

@@ -0,0 +1,29 @@
+/*
+ * Copyright 2019 Blender Foundation
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef __JITTER_H__
+#define __JITTER_H__
+
+#include "util/util_types.h"
+
+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
+
+#endif /* __JITTER_H__ */

intern/cycles/render/scene.cpp

@@ -77,7 +77,7 @@ DeviceScene::DeviceScene(Device *device)
       svm_nodes(device, "__svm_nodes", MEM_TEXTURE),
       shaders(device, "__shaders", MEM_TEXTURE),
       lookup_table(device, "__lookup_table", MEM_TEXTURE),
-      sobol_directions(device, "__sobol_directions", MEM_TEXTURE),
+      sample_pattern_lut(device, "__sample_pattern_lut", MEM_TEXTURE),
       ies_lights(device, "__ies", MEM_TEXTURE)
 {
   memset((void *)&data, 0, sizeof(data));

intern/cycles/render/scene.h

@@ -119,7 +119,7 @@ class DeviceScene {
   device_vector<float> lookup_table;
 
   /* integrator */
-  device_vector<uint> sobol_directions;
+  device_vector<uint> sample_pattern_lut;
 
   /* ies lights */
   device_vector<float> ies_lights;

intern/cycles/render/session.cpp

@@ -976,7 +976,7 @@ bool Session::update_scene()
   Integrator *integrator = scene->integrator;
   BakeManager *bake_manager = scene->bake_manager;
 
-  if (integrator->sampling_pattern == SAMPLING_PATTERN_CMJ || bake_manager->get_baking()) {
+  if (integrator->sampling_pattern != SAMPLING_PATTERN_SOBOL || bake_manager->get_baking()) {
     int aa_samples = tile_manager.num_samples;
 
     if (aa_samples != integrator->aa_samples) {