Cycles: Shade volume with null scattering

The distance sampling is mostly based on weighted delta tracking from [Monte Carlo Methods for Volumetric Light Transport Simulation] (http://iliyan.com/publications/VolumeSTAR/VolumeSTAR_EG2018.pdf). The recursive Monte Carlo estimation of the Radiative Transfer Equation is \[\langle L \rangle=\frac{\bar T(x\rightarrow y)}{\bar p(x\rightarrow y)}(L_e+\sigma_s L_s + \sigma_n L).\] where \(\bar T(x\rightarrow y) = e^{-\bar\sigma\Vert x-y\Vert}\) is the majorant transmittance between points \(x\) and \(y\), \(p(x\rightarrow y) = \bar\sigma e^{-\bar\sigma\Vert x-y\Vert}\) is the probability of sampling point \(y\) from point \(x\) following exponential distribution. At each recursive step, we randomly pick one of the two events proportional to their weights: * If \(\xi < \frac{\sigma_s}{\sigma_s+\vert\sigma_n\vert}\), we sample scatter event and evaluate \(L_s\). * Otherwise, no real collision happens and we continue the recursive process. The emission \(L_e\) is evaluated at each step. This also removes some unused volume settings from the UI: * "Max Steps" is removed, because the step size is automatically specified by the volume octree. There is a hard-coded threshold `VOLUME_MAX_STEPS` to prevent numerical issues. * "Homogeneous" is automatically detected during density evaluation An option "Unbiased" is added to the UI. When enabled, densities above the majorant are clamped.
2025-02-11 19:29:51 +01:00
parent 8c36f9ce49
commit a7283fc1d5
26 changed files with 413 additions and 572 deletions
--- a/intern/cycles/blender/addon/properties.py
+++ b/intern/cycles/blender/addon/properties.py
@@ -791,28 +791,11 @@ class CyclesRenderSettings(bpy.types.PropertyGroup):
        default=8,
    )

-    volume_step_rate: FloatProperty(
-        name="Step Rate",
-        description="Globally adjust detail for volume rendering, on top of automatically estimated step size. "
-                    "Higher values reduce render time, lower values render with more detail",
-        default=1.0,
-        min=0.01, max=100.0, soft_min=0.1, soft_max=10.0, precision=2
-    )
-
-    volume_preview_step_rate: FloatProperty(
-        name="Step Rate",
-        description="Globally adjust detail for volume rendering, on top of automatically estimated step size. "
-                    "Higher values reduce render time, lower values render with more detail",
-        default=1.0,
-        min=0.01, max=100.0, soft_min=0.1, soft_max=10.0, precision=2
-    )
-
-    volume_max_steps: IntProperty(
-        name="Max Steps",
-        description="Maximum number of steps through the volume before giving up, "
-        "to avoid extremely long render times with big objects or small step sizes",
-        default=1024,
-        min=2, max=65536
+    volume_unbiased: BoolProperty(
+        name="Unbiased",
+        description="If enabled, volume rendering converges to the correct result with sufficiently large numbers"
+        "of samples, but might appear noisier in the process",
+        default=False,
    )

    dicing_rate: FloatProperty(
@@ -1150,12 +1133,6 @@ class CyclesMaterialSettings(bpy.types.PropertyGroup):
        description="Apply corrections to solve shadow terminator artifacts caused by bump mapping",
        default=True,
    )
-    homogeneous_volume: BoolProperty(
-        name="Homogeneous Volume",
-        description="When using volume rendering, assume volume has the same density everywhere "
-        "(not using any textures), for faster rendering",
-        default=False,
-    )
    volume_sampling: EnumProperty(
        name="Volume Sampling",
        description="Sampling method to use for volumes",
@@ -1170,14 +1147,6 @@ class CyclesMaterialSettings(bpy.types.PropertyGroup):
        default='LINEAR',
    )

-    volume_step_rate: FloatProperty(
-        name="Step Rate",
-        description="Scale the distance between volume shader samples when rendering the volume "
-                    "(lower values give more accurate and detailed results, but also increased render time)",
-        default=1.0,
-        min=0.001, max=1000.0, soft_min=0.1, soft_max=10.0, precision=4
-    )
-
    @classmethod
    def register(cls):
        bpy.types.Material.cycles = PointerProperty(
@@ -1260,12 +1229,6 @@ class CyclesWorldSettings(bpy.types.PropertyGroup):
        min=0, max=1024,
        default=1024,
    )
-    homogeneous_volume: BoolProperty(
-        name="Homogeneous Volume",
-        description="When using volume rendering, assume volume has the same density everywhere "
-        "(not using any textures), for faster rendering",
-        default=False,
-    )
    volume_sampling: EnumProperty(
        name="Volume Sampling",
        description="Sampling method to use for volumes",
@@ -1278,13 +1241,6 @@ class CyclesWorldSettings(bpy.types.PropertyGroup):
        items=enum_volume_interpolation,
        default='LINEAR',
    )
-    volume_step_size: FloatProperty(
-        name="Step Size",
-        description="Distance between volume shader samples when rendering the volume "
-                    "(lower values give more accurate and detailed results, but also increased render time)",
-        default=1.0,
-        min=0.0000001, max=100000.0, soft_min=0.1, soft_max=100.0, precision=4
-    )

    @classmethod
    def register(cls):
--- a/intern/cycles/blender/addon/ui.py
+++ b/intern/cycles/blender/addon/ui.py
@@ -569,11 +569,8 @@ class CYCLES_RENDER_PT_volumes(CyclesButtonsPanel, Panel):
        scene = context.scene
        cscene = scene.cycles

-        col = layout.column(align=True)
-        col.prop(cscene, "volume_step_rate", text="Step Rate Render")
-        col.prop(cscene, "volume_preview_step_rate", text="Viewport")
-
-        layout.prop(cscene, "volume_max_steps", text="Max Steps")
+        col = layout.column()
+        col.prop(cscene, "volume_unbiased", text="Unbiased")


 class CYCLES_RENDER_PT_light_paths(CyclesButtonsPanel, Panel):
@@ -1817,10 +1814,6 @@ class CYCLES_WORLD_PT_settings_volume(CyclesButtonsPanel, Panel):
        sub = col.column()
        col.prop(cworld, "volume_sampling", text="Sampling")
        col.prop(cworld, "volume_interpolation", text="Interpolation")
-        col.prop(cworld, "homogeneous_volume", text="Homogeneous")
-        sub = col.column()
-        sub.active = not cworld.homogeneous_volume
-        sub.prop(cworld, "volume_step_size")


 class CYCLES_WORLD_PT_settings_light_group(CyclesButtonsPanel, Panel):
@@ -1993,10 +1986,6 @@ class CYCLES_MATERIAL_PT_settings_volume(CyclesButtonsPanel, Panel):
        sub = col.column()
        col.prop(cmat, "volume_sampling", text="Sampling")
        col.prop(cmat, "volume_interpolation", text="Interpolation")
-        col.prop(cmat, "homogeneous_volume", text="Homogeneous")
-        sub = col.column()
-        sub.active = not cmat.homogeneous_volume
-        sub.prop(cmat, "volume_step_rate")

    def draw(self, context):
        self.draw_shared(self, context, context.material)
--- a/intern/cycles/blender/shader.cpp
+++ b/intern/cycles/blender/shader.cpp
@@ -1578,10 +1578,8 @@ void BlenderSync::sync_materials(BL::Depsgraph &b_depsgraph, bool update_all)
      shader->set_emission_sampling_method(get_emission_sampling(cmat));
      shader->set_use_transparent_shadow(b_mat.use_transparent_shadow());
      shader->set_use_bump_map_correction(get_boolean(cmat, "use_bump_map_correction"));
-      shader->set_heterogeneous_volume(!get_boolean(cmat, "homogeneous_volume"));
      shader->set_volume_sampling_method(get_volume_sampling(cmat));
      shader->set_volume_interpolation_method(get_volume_interpolation(cmat));
-      shader->set_volume_step_rate(get_float(cmat, "volume_step_rate"));
      shader->set_displacement_method(get_displacement_method(b_mat));

      shader->set_graph(std::move(graph));
@@ -1646,10 +1644,8 @@ void BlenderSync::sync_world(BL::Depsgraph &b_depsgraph, BL::SpaceView3D &b_v3d,

      /* volume */
      PointerRNA cworld = RNA_pointer_get(&b_world.ptr, "cycles");
-      shader->set_heterogeneous_volume(!get_boolean(cworld, "homogeneous_volume"));
      shader->set_volume_sampling_method(get_volume_sampling(cworld));
      shader->set_volume_interpolation_method(get_volume_interpolation(cworld));
-      shader->set_volume_step_rate(get_float(cworld, "volume_step_size"));
    }
    else if (new_viewport_parameters.use_scene_world && b_world) {
      BackgroundNode *background = graph->create_node<BackgroundNode>();
--- a/intern/cycles/blender/sync.cpp
+++ b/intern/cycles/blender/sync.cpp
@@ -349,15 +349,10 @@ void BlenderSync::sync_integrator(BL::ViewLayer &b_view_layer,
  integrator->set_max_glossy_bounce(get_int(cscene, "glossy_bounces"));
  integrator->set_max_transmission_bounce(get_int(cscene, "transmission_bounces"));
  integrator->set_max_volume_bounce(get_int(cscene, "volume_bounces"));
-
+  integrator->set_volume_unbiased(get_boolean(cscene, "volume_unbiased"));
  integrator->set_transparent_min_bounce(get_int(cscene, "min_transparent_bounces"));
  integrator->set_transparent_max_bounce(get_int(cscene, "transparent_max_bounces"));

-  integrator->set_volume_max_steps(get_int(cscene, "volume_max_steps"));
-  const float volume_step_rate = (preview) ? get_float(cscene, "volume_preview_step_rate") :
-                                             get_float(cscene, "volume_step_rate");
-  integrator->set_volume_step_rate(volume_step_rate);
-
  integrator->set_caustics_reflective(get_boolean(cscene, "caustics_reflective"));
  integrator->set_caustics_refractive(get_boolean(cscene, "caustics_refractive"));
  integrator->set_filter_glossy(get_float(cscene, "blur_glossy"));
--- a/intern/cycles/blender/volume.cpp
+++ b/intern/cycles/blender/volume.cpp
@@ -284,7 +284,6 @@ static void sync_volume_object(BL::BlendData &b_data,

  BL::VolumeRender b_render(b_volume.render());

-  volume->set_step_size(b_render.step_size());
  volume->set_object_space((b_render.space() == BL::VolumeRender::space_OBJECT));

  float velocity_scale = b_volume.velocity_scale();
--- a/intern/cycles/kernel/closure/volume.h
+++ b/intern/cycles/kernel/closure/volume.h
@@ -125,8 +125,10 @@ ccl_device float volume_phase_get_g(const ccl_private ShaderVolumeClosure *svc)

 /* Volume sampling utilities. */

-/* todo: this value could be tweaked or turned into a probability to avoid
- * unnecessary work in volumes and subsurface scattering. */
+/* Ignore paths that have volume throughput below this value, to avoid unnecessary work
+ * and precision issues.
+ * TODO: this value could be tweaked or turned into a probability to avoid unnecessary work in
+ * volumes and subsurface scattering. */
 #define VOLUME_THROUGHPUT_EPSILON 1e-6f

 ccl_device Spectrum volume_color_transmittance(Spectrum sigma, const float t)
--- a/intern/cycles/kernel/data_arrays.h
+++ b/intern/cycles/kernel/data_arrays.h
@@ -23,7 +23,6 @@ KERNEL_DATA_ARRAY(KernelObject, objects)
 KERNEL_DATA_ARRAY(Transform, object_motion_pass)
 KERNEL_DATA_ARRAY(DecomposedTransform, object_motion)
 KERNEL_DATA_ARRAY(uint, object_flag)
-KERNEL_DATA_ARRAY(float, object_volume_step)
 KERNEL_DATA_ARRAY(uint, object_prim_offset)

 /* cameras */
--- a/intern/cycles/kernel/data_template.h
+++ b/intern/cycles/kernel/data_template.h
@@ -25,7 +25,6 @@ KERNEL_STRUCT_MEMBER(background, int, use_sun_guiding)
 /* Only shader index. */
 KERNEL_STRUCT_MEMBER(background, int, surface_shader)
 KERNEL_STRUCT_MEMBER(background, int, volume_shader)
-KERNEL_STRUCT_MEMBER(background, float, volume_step_size)
 KERNEL_STRUCT_MEMBER(background, int, transparent)
 KERNEL_STRUCT_MEMBER(background, float, transparent_roughness_squared_threshold)
 /* Sun sampling. */
@@ -45,6 +44,7 @@ KERNEL_STRUCT_MEMBER(background, int, light_index)
 KERNEL_STRUCT_MEMBER(background, int, object_index)
 /* Padding. */
 KERNEL_STRUCT_MEMBER(background, int, pad1)
+KERNEL_STRUCT_MEMBER(background, int, pad2)
 KERNEL_STRUCT_END(KernelBackground)

 /* BVH: own BVH2 if no native device acceleration struct used. */
@@ -201,8 +201,7 @@ KERNEL_STRUCT_MEMBER_DONT_SPECIALIZE
 KERNEL_STRUCT_MEMBER(integrator, int, blue_noise_sequence_length)
 /* Volume render. */
 KERNEL_STRUCT_MEMBER(integrator, int, use_volumes)
-KERNEL_STRUCT_MEMBER(integrator, int, volume_max_steps)
-KERNEL_STRUCT_MEMBER(integrator, float, volume_step_rate)
+KERNEL_STRUCT_MEMBER(integrator, int, volume_unbiased)
 /* Shadow catcher. */
 KERNEL_STRUCT_MEMBER(integrator, int, has_shadow_catcher)
 /* Closure filter. */
--- a/intern/cycles/kernel/geom/object.h
+++ b/intern/cycles/kernel/geom/object.h
@@ -352,7 +352,7 @@ ccl_device_inline float3 object_dupli_uv(KernelGlobals kg, const int object)
  return make_float3(kobject->dupli_uv[0], kobject->dupli_uv[1], 0.0f);
 }

-/* Volume step size */
+/* Volume density */

 ccl_device_inline float object_volume_density(KernelGlobals kg, const int object)
 {
@@ -363,15 +363,6 @@ ccl_device_inline float object_volume_density(KernelGlobals kg, const int object
  return kernel_data_fetch(objects, object).volume_density;
 }

-ccl_device_inline float object_volume_step_size(KernelGlobals kg, const int object)
-{
-  if (object == OBJECT_NONE) {
-    return kernel_data.background.volume_step_size;
-  }
-
-  return kernel_data_fetch(object_volume_step, object);
-}
-
 /* Pass ID for shader */

 ccl_device int shader_pass_id(KernelGlobals kg, const ccl_private ShaderData *sd)
--- a/intern/cycles/kernel/integrator/shade_volume.h
+++ b/intern/cycles/kernel/integrator/shade_volume.h
@@ -52,13 +52,13 @@ struct VolumeIntegrateResult {
  ShaderVolumePhases indirect_phases;
 };

-/* Ignore paths that have volume throughput below this value, to avoid unnecessary work
- * and precision issues.
- * todo: this value could be tweaked or turned into a probability to avoid unnecessary
- * work in volumes and subsurface scattering. */
-#  define VOLUME_THROUGHPUT_EPSILON 1e-6f
+/* We use both volume octree and volume stack, sometimes they disagree on whether a point is inside
+ * a volume or not. We accept small numerical precision issues, above this threshold the volume
+ * stack shall prevail. */
 /* TODO(weizhen): tweak this value. */
 #  define OVERLAP_EXP 5e-4f
+/* Restrict the number of steps in case of numerical problems */
+#  define VOLUME_MAX_STEPS 1024
 /* Number of mantissa bits of floating-point numbers. */
 #  define MANTISSA_BITS 23

@@ -126,36 +126,6 @@ ccl_device_inline bool volume_shader_sample(KernelGlobals kg,
  return true;
 }

-/* Determines the next shading position. */
-struct VolumeStep {
-  /* Shift starting point of all segments by a random amount to avoid banding artifacts due to
-   * biased ray marching with insufficient step size. */
-  float offset;
-
-  /* Step size taken at each marching step. */
-  float size;
-
-  /* Perform shading at this offset within a step, to integrate over the entire step segment. */
-  float shade_offset;
-
-  /* Current step. */
-  int step;
-
-  /* Current active segment. */
-  Interval<float> t;
-};
-
-ccl_device_forceinline void volume_step_init(KernelGlobals kg,
-                                             const ccl_private RNGState *rng_state,
-                                             const float tmin,
-                                             ccl_private VolumeStep *vstep)
-{
-  vstep->step = 0;
-  vstep->t.min = vstep->t.max = tmin;
-  vstep->shade_offset = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_SHADE_OFFSET);
-  vstep->offset = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_OFFSET);
-}
-
 /* -------------------------------------------------------------------- */
 /** \name Hierarchical DDA for ray tracing the volume octree
 *
@@ -417,8 +387,7 @@ ccl_device bool volume_octree_setup(KernelGlobals kg,
                                    const IntegratorGenericState state,
                                    const ccl_private RNGState *rng_state,
                                    const uint32_t path_flag,
-                                    ccl_private OctreeTracing &global,
-                                    ccl_private VolumeStep &vstep)
+                                    ccl_private OctreeTracing &global)
 {
  if (global.no_overlap) {
    /* If the current active octree is already set up. */
@@ -474,26 +443,11 @@ ccl_device bool volume_octree_setup(KernelGlobals kg,
    }
  }

-  if (!global.node) {
-    /* Stack empty. */
-    return false;
-  }
-
-  if (global.t.is_empty()) {
-    return false;
-  }
-
  if (i == 1) {
    global.no_overlap = true;
  }

-  /* Step size should ideally be as small as the active voxel span, but not so small that we
-   * never exit the volume. */
-  const int steps_left = kernel_data.integrator.volume_max_steps - vstep.step;
-  vstep.size = fminf(global.t.length(), 1.0f / global.sigma.range());
-  vstep.size = fmaxf(vstep.size, (ray->tmax - vstep.t.min) / float(steps_left));
-
-  return true;
+  return global.node && !global.t.is_empty();
 }

 /* Advance to the next adjacent leaf node and update the active interval. */
@@ -505,14 +459,14 @@ ccl_device_inline bool volume_octree_advance(KernelGlobals kg,
                                             const ccl_private RNGState *rng_state,
                                             const uint32_t path_flag,
                                             ccl_private OctreeTracing &octree,
-                                             ccl_private VolumeStep &vstep)
+                                             const int step)
 {
  if (octree.t.max >= ray->tmax) {
    /* Reached the last segment. */
    return false;
  }

-  if (vstep.step > kernel_data.integrator.volume_max_steps) {
+  if (step >= VOLUME_MAX_STEPS) {
    /* Exceeds maximal steps. */
    return false;
  }
@@ -548,51 +502,7 @@ ccl_device_inline bool volume_octree_advance(KernelGlobals kg,

  octree.sigma = volume_object_get_extrema<shadow>(
      kg, ray, sd, state, octree, rng_state, path_flag);
-  return volume_octree_setup<shadow>(kg, ray, sd, state, rng_state, path_flag, octree, vstep);
-}
-
-/* Advance to the next interval. If the step size exceeds the current leaf node, find the next leaf
- * node. */
-template<const bool shadow, typename IntegratorGenericState>
-ccl_device bool volume_integrate_advance(KernelGlobals kg,
-                                         const ccl_private Ray *ccl_restrict ray,
-                                         ccl_private ShaderData *ccl_restrict sd,
-                                         const IntegratorGenericState state,
-                                         const ccl_private RNGState *rng_state,
-                                         const uint32_t path_flag,
-                                         ccl_private OctreeTracing &octree,
-                                         ccl_private VolumeStep &vstep)
-{
-  if (vstep.t.max == ray->tmax) {
-    /* Reached the last segment. */
-    return false;
-  }
-
-  /* Advance to new position. */
-  vstep.t.min = vstep.t.max;
-  bool success = true;
-  if (!octree.node) {
-    /* Initialize octree. */
-    success = volume_octree_setup<shadow>(kg, ray, sd, state, rng_state, path_flag, octree, vstep);
-    vstep.t.max = octree.t.min + vstep.offset * vstep.size;
-  }
-  else {
-    float candidate_t_max = vstep.t.max + vstep.size;
-    if (candidate_t_max >= octree.t.max) {
-      /* Advance to next voxel. */
-      volume_octree_advance<shadow>(kg, ray, sd, state, rng_state, path_flag, octree, vstep);
-      candidate_t_max = fminf(candidate_t_max, vstep.t.max + vstep.size);
-    }
-    vstep.t.max = candidate_t_max;
-  }
-
-  /* Clamp to prevent numerical issues. */
-  vstep.t.max = clamp(vstep.t.max, vstep.t.min + OVERLAP_EXP, ray->tmax);
-
-  const float shade_t = mix(vstep.t.min, vstep.t.max, vstep.shade_offset);
-  sd->P = ray->P + ray->D * shade_t;
-
-  return success && vstep.step++ < kernel_data.integrator.volume_max_steps;
+  return volume_octree_setup<shadow>(kg, ray, sd, state, rng_state, path_flag, octree);
 }

 /** \} */
@@ -611,8 +521,7 @@ ccl_device_inline bool volume_octree_advance_shadow(KernelGlobals kg,
                                                    const IntegratorShadowState state,
                                                    ccl_private RNGState *rng_state,
                                                    const uint32_t path_flag,
-                                                    ccl_private OctreeTracing &octree,
-                                                    ccl_private VolumeStep &vstep)
+                                                    ccl_private OctreeTracing &octree)
 {
  /* Advance random number offset. */
  rng_state->rng_offset += PRNG_BOUNCE_NUM;
@@ -621,7 +530,7 @@ ccl_device_inline bool volume_octree_advance_shadow(KernelGlobals kg,
  const float tmin = octree.t.min;

  while (octree.t.is_empty() || sigma.range() * octree.t.length() < 1.0f) {
-    if (!volume_octree_advance<true>(kg, ray, sd, state, rng_state, path_flag, octree, vstep)) {
+    if (!volume_octree_advance<true>(kg, ray, sd, state, rng_state, path_flag, octree, 0)) {
      return !octree.t.is_empty();
    }

@@ -663,7 +572,7 @@ ccl_device Spectrum volume_transmittance(KernelGlobals kg,
  /* Expected number of steps with residual ratio tracking. */
  const float expected_steps = sigma_c * ray_length;
  /* Number of samples for the biased estimator. */
-  const int k = clamp(int(roundf(expected_steps)), 1, kernel_data.integrator.volume_max_steps);
+  const int k = clamp(int(roundf(expected_steps)), 1, VOLUME_MAX_STEPS);

  /* Sample the evaluation order of the debiasing term. */
  /* Use the same random number for all pixels to sync the workload on GPU. */
@@ -748,18 +657,13 @@ ccl_device void volume_shadow_heterogeneous(KernelGlobals kg,

  path_state_rng_scramble(&rng_state, 0x8647ace4);

-  /* Prepare for stepping. */
-  VolumeStep vstep;
-  volume_step_init(kg, &rng_state, ray->tmin, &vstep);
-
-  const uint32_t path_flag = PATH_RAY_SHADOW;
-
  OctreeTracing octree(ray->tmin);
-  if (!volume_octree_setup<true>(kg, ray, sd, state, &rng_state, path_flag, octree, vstep)) {
+  const uint32_t path_flag = PATH_RAY_SHADOW;
+  if (!volume_octree_setup<true>(kg, ray, sd, state, &rng_state, path_flag, octree)) {
    return;
  }

-  while (volume_octree_advance_shadow(kg, ray, sd, state, &rng_state, path_flag, octree, vstep)) {
+  while (volume_octree_advance_shadow(kg, ray, sd, state, &rng_state, path_flag, octree)) {
    const float sigma = octree.sigma.range();
    *throughput *= volume_transmittance<true>(
        kg, state, ray, sd, sigma, octree.t, &rng_state, path_flag);
@@ -859,36 +763,6 @@ ccl_device_inline bool volume_valid_direct_ray_segment(KernelGlobals kg,
  return !t_range->is_empty();
 }

-/* Emission */
-
-ccl_device Spectrum volume_emission_integrate(ccl_private VolumeShaderCoefficients *coeff,
-                                              const int closure_flag,
-                                              Spectrum transmittance,
-                                              const float t)
-{
-  /* integral E * exp(-sigma_t * t) from 0 to t = E * (1 - exp(-sigma_t * t))/sigma_t
-   * this goes to E * t as sigma_t goes to zero
-   *
-   * todo: we should use an epsilon to avoid precision issues near zero sigma_t */
-  Spectrum emission = coeff->emission;
-
-  if (closure_flag & SD_EXTINCTION) {
-    Spectrum sigma_t = coeff->sigma_t;
-
-    FOREACH_SPECTRUM_CHANNEL (i) {
-      GET_SPECTRUM_CHANNEL(emission, i) *= (GET_SPECTRUM_CHANNEL(sigma_t, i) > 0.0f) ?
-                                               (1.0f - GET_SPECTRUM_CHANNEL(transmittance, i)) /
-                                                   GET_SPECTRUM_CHANNEL(sigma_t, i) :
-                                               t;
-    }
-  }
-  else {
-    emission *= t;
-  }
-
-  return emission;
-}
-
 /* Volume Integration */

 struct VolumeIntegrateState {
@@ -899,157 +773,351 @@ struct VolumeIntegrateState {
  /* Multiple importance sampling. */
  VolumeSampleMethod direct_sample_method;
  bool use_mis;
+  /* Probability of sampling the scatter position using null scattering. */
  float distance_pdf;
+  /* Probability of sampling the scatter position using equiangular sampling. */
  float equiangular_pdf;
+
+  /* Ratio tracking estimator of the volume transmittance, with MIS applied. */
+  float transmittance;
+  /* Steps taken while tracking. Should not exceed `VOLUME_MAX_STEPS`. */
+  int step;
+
+  bool stop;
+
+  /* Extra fields for path guiding and denoising. */
+  Spectrum emission;
+#  ifdef __DENOISING_FEATURES__
+  Spectrum albedo;
+#  endif
 };

-ccl_device bool volume_integrate_should_stop(ccl_private VolumeIntegrateResult &result)
+ccl_device bool volume_integrate_should_stop(const ccl_private VolumeIntegrateResult &result,
+                                             const ccl_private VolumeIntegrateState &vstate)
 {
-  /* Stop if nearly all light blocked. */
-  if (!result.indirect_scatter) {
-    if (reduce_max(result.indirect_throughput) < VOLUME_THROUGHPUT_EPSILON) {
-      result.indirect_throughput = zero_spectrum();
-      return true;
-    }
-  }
-  else if (!result.direct_scatter) {
-    if (reduce_max(result.direct_throughput) < VOLUME_THROUGHPUT_EPSILON) {
-      return true;
-    }
+  if (result.indirect_scatter && result.direct_scatter) {
+    return true;
  }

-  /* If we have scattering data for both direct and indirect, we're done. */
-  return (result.direct_scatter && result.indirect_scatter);
+  return vstate.stop;
 }

-/* Returns true if we found the indirect scatter position within the current active ray segment. */
-ccl_device bool volume_sample_indirect_scatter(
-    const Spectrum transmittance,
-    const Spectrum channel_pdf,
-    const int channel,
-    const ccl_private ShaderData *ccl_restrict sd,
-    const ccl_private VolumeShaderCoefficients &ccl_restrict coeff,
-    const ccl_private Interval<float> &t,
+/* -------------------------------------------------------------------- */
+/** \name Null Scattering
+ * \{ */
+
+/* TODO(weizhen): homogeneous volume sample process can be simplified, especially we can take much
+ * less steps. */
+
+/* In a null-scattering framework, we fill the volume with fictitious particles, so that the
+ * density is `sigma_max` everywhere. The null-scattering coefficients `sigma_n` is then defined by
+ * the density of such particles. */
+ccl_device_inline Spectrum volume_null_event_coefficients(
+    KernelGlobals kg, ccl_private VolumeShaderCoefficients &coeff, const float sigma_max)
+{
+  const Spectrum sigma_max_ = make_spectrum(sigma_max);
+
+  /* Clamp if bias is allowed. */
+  if (!kernel_data.integrator.volume_unbiased) {
+    coeff.sigma_t = min(coeff.sigma_t, sigma_max_);
+    coeff.sigma_s = min(coeff.sigma_s, sigma_max_);
+  }
+
+  return sigma_max_ - coeff.sigma_t;
+}
+
+/* The probability of sampling real scattering event at each candidate point of delta tracking. */
+ccl_device_inline float volume_scatter_probability(
+    const ccl_private VolumeShaderCoefficients &coeff,
+    const Spectrum sigma_n,
+    const Spectrum throughput)
+{
+  /* For procedure volumes `sigma_max` might not be the strict upper bound. Use absolute value
+   * to handle negative null scattering coefficient as suggested in the paper. */
+  const Spectrum abs_sigma_n = fabs(sigma_n);
+
+  /* We use `sigma_s` instead of `sigma_t` to skip sampling the absorption event, because it
+   * always returns zero and has high variance. */
+  Spectrum sigma_c = coeff.sigma_s + abs_sigma_n;
+
+  /* Set `albedo` to 1 for the channel where extinction coefficient `sigma_t` is zero, to make
+   * sure that we sample a distance outside the current segment when that channel is picked,
+   * meaning light passes through without attenuation. */
+  const Spectrum albedo = safe_divide_color(coeff.sigma_s, coeff.sigma_t, 1.0f);
+
+  /* Assign weights per channel to pick scattering event based on throughput and single
+   * scattering albedo.  */
+  /* TODO(weizhen): currently the sample distance is the same for each color channel, revisit
+   * the MIS weight when we use Spectral Majorant. */
+  const Spectrum channel_pdf = volume_sample_channel_pdf(albedo, throughput);
+
+  return dot(coeff.sigma_s / sigma_c, channel_pdf);
+}
+
+/**
+ * Sample indirect scatter position along the ray based on weighted delta tracking, from
+ * [Spectral and Decomposition Tracking for Rendering Heterogeneous Volumes]
+ * (https://disneyanimation.com/publications/spectral-and-decomposition-tracking-for-rendering-heterogeneous-volumes)
+ * by Peter Kutz et. al.
+ *
+ * The recursive Monte Carlo estimation of the Radiative Transfer Equation is
+ * <L> = T(x -> y) / p(x -> y) * (L_e + sigma_s * L_s + sigma_n * L),
+ * where T(x -> y) = exp(-sigma_max * dt) is the majorant transmittance between points x and y,
+ * and p(x -> y) = sigma_max * exp(-sigma_max * dt) is the probability of sampling point y from
+ * point x following exponential distribution.
+ * At each recursive step, we randomly pick one of the two events proportional to their weights:
+ * - If ξ < sigma_s / (sigma_s + |sigma_n|), we sample scatter event and evaluate L_s.
+ * - Otherwise, no real collision happens and we continue the recursive process.
+ * The emission L_e is evaluated at each step.
+ *
+ * \param sigma_max: majorant volume density inside the current octree node
+ * \param interval: interval of t along the ray.
+ */
+ccl_device void volume_sample_indirect_scatter(
+    KernelGlobals kg,
+    const IntegratorState state,
+    const ccl_private Ray *ccl_restrict ray,
+    ccl_private ShaderData *ccl_restrict sd,
+    const float sigma_max,
+    const Interval<float> interval,
+    ccl_private RNGState *rng_state,
    ccl_private VolumeIntegrateState &ccl_restrict vstate,
    ccl_private VolumeIntegrateResult &ccl_restrict result)
 {
  if (result.indirect_scatter) {
    /* Already sampled indirect scatter position. */
-    return false;
+    return;
  }

-  /* If sampled distance does not go beyond the current segment, we have found the scatter
-   * position. Otherwise continue searching and accumulate the transmittance along the ray. */
-  const float sample_transmittance = volume_channel_get(transmittance, channel);
-  if (1.0f - vstate.rscatter >= sample_transmittance) {
-    /* Pick `sigma_t` from a random channel. */
-    const float sample_sigma_t = volume_channel_get(coeff.sigma_t, channel);
+  /* Initialization. */
+  float t = interval.min;
+  const float inv_maj = (sigma_max == 0.0f) ? FLT_MAX : 1.0f / sigma_max;
+  const bool segment_has_equiangular = vstate.direct_sample_method == VOLUME_SAMPLE_EQUIANGULAR &&
+                                       interval.contains(result.direct_t) && vstate.use_mis;
+  bool direct_scatter = false;
+  while (vstate.step++ < VOLUME_MAX_STEPS) {
+    if (reduce_max(fabs(result.indirect_throughput)) < VOLUME_THROUGHPUT_EPSILON) {
+      /* TODO(weizhen): terminate using Russian Roulette. */
+      /* TODO(weizhen): deal with negative transmittance. */
+      /* TODO(weizhen): should we stop if direct_scatter not yet found? */
+      vstate.stop = true;
+      result.indirect_throughput = zero_spectrum();
+      return;
+    }

-    /* Generate the next distance using random walk, following exponential distribution
-     * p(dt) = sigma_t * exp(-sigma_t * dt). */
-    const float new_dt = -logf(1.0f - vstate.rscatter) / sample_sigma_t;
-    const float new_t = t.min + new_dt;
+    /* Generate the next distance using random walk. */
+    const float rand = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_SCATTER_DISTANCE);
+    t += sample_exponential_distribution(rand, inv_maj);

-    const Spectrum new_transmittance = volume_color_transmittance(coeff.sigma_t, new_dt);
-    /* PDF for density-based distance sampling is handled implicitly via
-     * transmittance / pdf = exp(-sigma_t * dt) / (sigma_t * exp(-sigma_t * dt)) = 1 / sigma_t. */
-    const float distance_pdf = dot(channel_pdf, coeff.sigma_t * new_transmittance);
+    /* Advance random number offset. */
+    rng_state->rng_offset += PRNG_BOUNCE_NUM;

-    if (vstate.distance_pdf * distance_pdf > VOLUME_SAMPLE_PDF_CUTOFF) {
-      /* Update throughput. */
+    if (segment_has_equiangular && t > result.direct_t && !direct_scatter) {
+      /* Stepped beyond the equiangular scatter position, compute direct throughput. */
+      direct_scatter = true;
+      result.direct_throughput = result.indirect_throughput * vstate.transmittance;
+      vstate.distance_pdf = vstate.transmittance * sigma_max;
+    }
+
+    if (t > interval.max) {
+      break;
+    }
+
+    sd->P = ray->P + ray->D * t;
+    VolumeShaderCoefficients coeff ccl_optional_struct_init;
+    if (!volume_shader_sample(kg, state, sd, &coeff)) {
+      continue;
+    }
+
+    /* Emission. */
+    if (sd->flag & SD_EMISSION) {
+      /* Emission = inv_sigma * (L_e + sigma_n * (inv_sigma * (L_e + sigma_n * ···))). */
+      const Spectrum emission = inv_maj * coeff.emission;
+      vstate.emission += result.indirect_throughput * emission;
+      guiding_record_volume_emission(kg, state, emission);
+    }
+
+    /* Null scattering coefficients. */
+    const Spectrum sigma_n = volume_null_event_coefficients(kg, coeff, sigma_max);
+
+    if (reduce_add(coeff.sigma_s) == 0.0f) {
+      /* Absorption only. Deterministically choose null scattering and estimate the transmittance
+       * of the current ray segment. */
+      result.indirect_throughput *= sigma_n * inv_maj;
+      continue;
+    }
+
+#  ifdef __DENOISING_FEATURES__
+    if (INTEGRATOR_STATE(state, path, flag) & PATH_RAY_DENOISING_FEATURES) {
+      /* Albedo = inv_sigma * (sigma_s + sigma_n * (inv_sigma * (sigma_s + sigma_n * ···))). */
+      vstate.albedo += result.indirect_throughput * coeff.sigma_s * inv_maj;
+    }
+#  endif
+
+    const float prob_s = volume_scatter_probability(coeff, sigma_n, result.indirect_throughput);
+    if (vstate.rchannel < prob_s) {
+      /* Sampled scatter event. */
+      result.indirect_throughput *= coeff.sigma_s * inv_maj / prob_s;
+      result.indirect_t = t;
      result.indirect_scatter = true;
-      result.indirect_t = new_t;
-      result.indirect_throughput *= coeff.sigma_s * new_transmittance / distance_pdf;
-      if (vstate.direct_sample_method == VOLUME_SAMPLE_DISTANCE) {
-        vstate.distance_pdf *= distance_pdf;
-      }
-
      volume_shader_copy_phases(&result.indirect_phases, sd);

-      return true;
-    }
-  }
-  else {
-    /* Update throughput. */
-    const float distance_pdf = dot(channel_pdf, transmittance);
-    result.indirect_throughput *= transmittance / distance_pdf;
-    if (vstate.direct_sample_method == VOLUME_SAMPLE_DISTANCE) {
-      vstate.distance_pdf *= distance_pdf;
+      if (vstate.direct_sample_method == VOLUME_SAMPLE_DISTANCE) {
+        /* If using distance sampling for direct light, just copy parameters of indirect light
+         * since we scatter at the same point. */
+        result.direct_scatter = true;
+        result.direct_t = result.indirect_t;
+        result.direct_throughput = result.indirect_throughput;
+        volume_shader_copy_phases(&result.direct_phases, sd);
+        if (vstate.use_mis) {
+          vstate.distance_pdf = vstate.transmittance * prob_s * sigma_max;
+        }
+      }
+      return;
    }

-    /* Remap rscatter so we can reuse it and keep thing stratified. */
-    vstate.rscatter = 1.0f - (1.0f - vstate.rscatter) / sample_transmittance;
+    /* Null scattering. Accumulate weight and continue. */
+    const float prob_n = 1.0f - prob_s;
+    result.indirect_throughput *= sigma_n * inv_maj / prob_n;
+
+    if (vstate.use_mis) {
+      vstate.transmittance *= prob_n;
+    }
+
+    /* Rescale random number for reusing. */
+    vstate.rchannel = (vstate.rchannel - prob_s) / prob_n;
  }

-  return false;
+  /* No scatter event sampled in the interval. */
 }

-/* Find direct and indirect scatter positions. */
-ccl_device_forceinline void volume_integrate_step_scattering(
-    const ccl_private ShaderData *sd,
-    const ccl_private Ray *ray,
-    const ccl_private EquiangularCoefficients &equiangular_coeffs,
-    const ccl_private VolumeShaderCoefficients &ccl_restrict coeff,
-    const Spectrum transmittance,
+/* Throughput and pdf for equiangular sampling.
+ * If MIS is used with transmittance-based distance sampling, we compute the direct throughput from
+ * the indirect throughput in the function above. Otherwise, we use telescoping for higher quality.
+ */
+ccl_device_inline void volume_equiangular_direct_scatter(
+    KernelGlobals kg,
+    const IntegratorState state,
+    const ccl_private Ray *ccl_restrict ray,
+    const ccl_private Extrema<float> &sigma,
    const ccl_private Interval<float> &t,
+    ccl_private ShaderData *ccl_restrict sd,
+    ccl_private RNGState *rng_state,
+    ccl_private VolumeIntegrateState &vstate,
+    ccl_private VolumeIntegrateResult &ccl_restrict result)
+{
+  if (vstate.direct_sample_method != VOLUME_SAMPLE_EQUIANGULAR) {
+    return;
+  }
+
+  const uint32_t path_flag = INTEGRATOR_STATE(state, path, flag);
+
+  if (t.contains(result.direct_t)) {
+    /* Equiangular scatter position is inside the current segment. */
+    sd->P = ray->P + ray->D * result.direct_t;
+    VolumeShaderCoefficients coeff ccl_optional_struct_init;
+    if (volume_shader_sample(kg, state, sd, &coeff) && (sd->flag & SD_SCATTER)) {
+      volume_shader_copy_phases(&result.direct_phases, sd);
+      result.direct_scatter = true;
+
+      if (vstate.use_mis) {
+        /* Compute distance pdf for multiple importance sampling. */
+        const Spectrum sigma_n = volume_null_event_coefficients(kg, coeff, sigma.max);
+
+        vstate.distance_pdf *= volume_scatter_probability(
+            coeff, sigma_n, result.direct_throughput);
+      }
+      else {
+        /* Compute transmittance until the direct scatter position. */
+        const Interval<float> t_ = {t.min, result.direct_t};
+        result.direct_throughput *= volume_transmittance<false>(
+            kg, state, ray, sd, sigma.range(), t_, rng_state, path_flag);
+      }
+
+      result.direct_throughput *= coeff.sigma_s / vstate.equiangular_pdf;
+    }
+  }
+  else if (result.direct_t > t.max && !vstate.use_mis) {
+    /* Accumulate transmittance. */
+    result.direct_throughput *= volume_transmittance<false>(
+        kg, state, ray, sd, sigma.range(), t, rng_state, path_flag);
+  }
+}
+
+/* Find direct and indirect scatter positions inside the current active octree leaf node. */
+ccl_device void volume_integrate_step_scattering(
+    KernelGlobals kg,
+    const IntegratorState state,
+    const ccl_private Ray *ccl_restrict ray,
+    const ccl_private Extrema<float> &sigma,
+    const ccl_private Interval<float> &interval,
+    ccl_private ShaderData *ccl_restrict sd,
+    ccl_private RNGState *rng_state,
    ccl_private VolumeIntegrateState &ccl_restrict vstate,
    ccl_private VolumeIntegrateResult &ccl_restrict result)
 {
-  /* Pick random color channel for sampling the scatter distance. We use the Veach one-sample model
-   * with balance heuristic for the channels.
-   * Set `albedo` to 1 for the channel where extinction coefficient `sigma_t` is zero, to make sure
-   * that we sample a distance outside the current segment when that channel is picked, meaning
-   * light passes through without attenuation. */
-  const Spectrum albedo = safe_divide_color(coeff.sigma_s, coeff.sigma_t, 1.0f);
-  Spectrum channel_pdf;
-  const int channel = volume_sample_channel(
-      albedo, result.indirect_throughput, &vstate.rchannel, &channel_pdf);
+  /* Distance sampling for indirect and optional direct lighting. */
+  volume_sample_indirect_scatter(
+      kg, state, ray, sd, sigma.max, interval, rng_state, vstate, result);

  /* Equiangular sampling for direct lighting. */
-  if (vstate.direct_sample_method == VOLUME_SAMPLE_EQUIANGULAR && !result.direct_scatter) {
-    if (t.contains(result.direct_t) && vstate.equiangular_pdf > VOLUME_SAMPLE_PDF_CUTOFF) {
-      const float new_dt = result.direct_t - t.min;
-      const Spectrum new_transmittance = volume_color_transmittance(coeff.sigma_t, new_dt);
+  volume_equiangular_direct_scatter(
+      kg, state, ray, sigma, interval, sd, rng_state, vstate, result);
+}

-      result.direct_scatter = true;
-      result.direct_throughput *= coeff.sigma_s * new_transmittance / vstate.equiangular_pdf;
-      volume_shader_copy_phases(&result.direct_phases, sd);
-
-      /* Multiple importance sampling. */
-      if (vstate.use_mis) {
-        const float distance_pdf = vstate.distance_pdf *
-                                   dot(channel_pdf, coeff.sigma_t * new_transmittance);
-        const float mis_weight = 2.0f * power_heuristic(vstate.equiangular_pdf, distance_pdf);
-        result.direct_throughput *= mis_weight;
-      }
-    }
-    else {
-      result.direct_throughput *= transmittance;
-      vstate.distance_pdf *= dot(channel_pdf, transmittance);
-    }
+/* Multiple Importance Sampling between equiangular sampling and distance sampling.
+ *
+ * According to [A null-scattering path integral formulation of light transport]
+ * (https://cs.dartmouth.edu/~wjarosz/publications/miller19null.html), the pdf of sampling a
+ * scattering event at point P using distance sampling is the probability of sampling a series of
+ * null events, and then a scatter event at P, i.e.
+ *
+ *                distance_pdf = (∏p_dist * p_null) * p_dist * p_scatter,
+ *
+ * where `p_dist = sigma_max * exp(-sigma_max * dt)` is the probability of sampling an incremental
+ * distance `dt` following exponential distribution, and `p_null = sigma_n / sigma_c` is the
+ * probability of sampling a null event at a certain point, `p_scatter = sigma_s / sigma_c` the
+ * probability of sampling a scatter event.
+ *
+ * The pdf of sampling a scattering event at point P using equiangular sampling is the probability
+ * of sampling a series of null events deterministically, and then a scatter event at the point of
+ * equiangular sampling, i.e.
+ *
+ *                     equiangular_pdf = (∏p_dist * 1) * T * p_equi,
+ *
+ * where `T = exp(-sigma_max * dt)` is the probability of sampling a distance beyond `dt` following
+ * exponential distribution, `p_equi` is the equiangular pdf. Since the null events are sampled
+ * deterministically, the pdf is 1 instead of `p_null`.
+ *
+ * When performing MIS between distance and equiangular sampling, since we use single-channel
+ * majorant, `p_dist` is shared in both pdfs, therefore we can write
+ *
+ *       distance_pdf / equiangular_pdf = (∏p_null) * sigma_max * p_scatter / p_equi.
+ *
+ * If we want to use multi-channel majorants in the future, the components do not cancel, but we
+ * can divide by the `p_dist` of the hero channel to alleviate numerical issues.
+ */
+ccl_device_inline void volume_direct_scatter_mis(
+    const ccl_private Ray *ccl_restrict ray,
+    const ccl_private VolumeIntegrateState &ccl_restrict vstate,
+    const ccl_private EquiangularCoefficients &equiangular_coeffs,
+    ccl_private VolumeIntegrateResult &ccl_restrict result)
+{
+  if (!vstate.use_mis || vstate.direct_sample_method == VOLUME_SAMPLE_NONE) {
+    return;
  }

-  /* Distance sampling for indirect and optional direct lighting. */
-  if (volume_sample_indirect_scatter(
-          transmittance, channel_pdf, channel, sd, coeff, t, vstate, result))
-  {
-    if (vstate.direct_sample_method == VOLUME_SAMPLE_DISTANCE) {
-      /* If using distance sampling for direct light, just copy parameters of indirect light
-       * since we scatter at the same point. */
-      result.direct_scatter = true;
-      result.direct_t = result.indirect_t;
-      result.direct_throughput = result.indirect_throughput;
-      volume_shader_copy_phases(&result.direct_phases, sd);
-
-      /* Multiple importance sampling. */
-      if (vstate.use_mis) {
-        const float equiangular_pdf = volume_equiangular_pdf(
-            ray, equiangular_coeffs, result.indirect_t);
-        const float mis_weight = power_heuristic(vstate.distance_pdf, equiangular_pdf);
-        result.direct_throughput *= 2.0f * mis_weight;
-      }
-    }
+  float mis_weight;
+  if (vstate.direct_sample_method == VOLUME_SAMPLE_DISTANCE) {
+    const float equiangular_pdf = volume_equiangular_pdf(ray, equiangular_coeffs, result.direct_t);
+    mis_weight = power_heuristic(vstate.distance_pdf, equiangular_pdf);
  }
+  else {
+    kernel_assert(vstate.direct_sample_method == VOLUME_SAMPLE_EQUIANGULAR);
+    mis_weight = power_heuristic(vstate.equiangular_pdf, vstate.distance_pdf);
+  }
+
+  result.direct_throughput *= 2.0f * mis_weight;
 }

 ccl_device_inline void volume_integrate_state_init(KernelGlobals kg,
@@ -1073,8 +1141,42 @@ ccl_device_inline void volume_integrate_state_init(KernelGlobals kg,
      vstate.direct_sample_method = VOLUME_SAMPLE_EQUIANGULAR;
    }
  }
+
+  vstate.distance_pdf = 0.0f;
  vstate.equiangular_pdf = 0.0f;
-  vstate.distance_pdf = 1.0f;
+  vstate.transmittance = 1.0f;
+  vstate.step = 0;
+  vstate.stop = false;
+
+  vstate.emission = zero_spectrum();
+#  ifdef __DENOISING_FEATURES__
+  vstate.albedo = zero_spectrum();
+#  endif
+}
+
+ccl_device_inline void volume_integrate_result_init(
+    const IntegratorState state,
+    const ccl_private Ray *ccl_restrict ray,
+    ccl_private VolumeIntegrateState &vstate,
+    const ccl_private EquiangularCoefficients &equiangular_coeffs,
+    ccl_private VolumeIntegrateResult &result)
+{
+  const Spectrum throughput = INTEGRATOR_STATE(state, path, throughput);
+  result.direct_throughput = (vstate.use_mis ||
+                              (vstate.direct_sample_method == VOLUME_SAMPLE_NONE)) ?
+                                 zero_spectrum() :
+                                 throughput;
+  result.indirect_throughput = throughput;
+
+  /* Equiangular sampling: compute distance and PDF in advance. */
+  if (vstate.direct_sample_method == VOLUME_SAMPLE_EQUIANGULAR) {
+    result.direct_t = volume_equiangular_sample(
+        ray, equiangular_coeffs, vstate.rscatter, &vstate.equiangular_pdf);
+  }
+
+#  if defined(__PATH_GUIDING__)
+  result.direct_sample_method = vstate.direct_sample_method;
+#  endif
 }

 /* Path tracing: sample point on light using equiangular sampling. */
@@ -1154,7 +1256,7 @@ ccl_device_forceinline void volume_integrate_heterogeneous(
    const IntegratorState state,
    const ccl_private Ray *ccl_restrict ray,
    ccl_private ShaderData *ccl_restrict sd,
-    const ccl_private RNGState *ccl_restrict rng_state,
+    ccl_private RNGState *ccl_restrict rng_state,
    ccl_global float *ccl_restrict render_buffer,
    ccl_private LightSample *ls,
    ccl_private VolumeIntegrateResult &result)
@@ -1165,102 +1267,47 @@ ccl_device_forceinline void volume_integrate_heterogeneous(
  const VolumeSampleMethod direct_sample_method = volume_direct_sample_method(
      kg, state, ray, sd, rng_state, &equiangular_coeffs, ls);

-  /* Prepare for stepping. */
-  VolumeStep vstep;
-  volume_step_init(kg, rng_state, ray->tmin, &vstep);
-
-  /* Initialize volume integration state. */
  VolumeIntegrateState vstate ccl_optional_struct_init;
  volume_integrate_state_init(kg, direct_sample_method, rng_state, vstate);

  /* Initialize volume integration result. */
-  const Spectrum throughput = INTEGRATOR_STATE(state, path, throughput);
-  result.direct_throughput = (vstate.direct_sample_method == VOLUME_SAMPLE_NONE) ?
-                                 zero_spectrum() :
-                                 throughput;
-  result.indirect_throughput = throughput;
-
-  /* Equiangular sampling: compute distance and PDF in advance. */
-  if (vstate.direct_sample_method == VOLUME_SAMPLE_EQUIANGULAR) {
-    result.direct_t = volume_equiangular_sample(
-        ray, equiangular_coeffs, vstate.rscatter, &vstate.equiangular_pdf);
-  }
-#  if defined(__PATH_GUIDING__)
-  result.direct_sample_method = vstate.direct_sample_method;
-#  endif
-
-#  ifdef __DENOISING_FEATURES__
-  const bool write_denoising_features = (INTEGRATOR_STATE(state, path, flag) &
-                                         PATH_RAY_DENOISING_FEATURES);
-  Spectrum accum_albedo = zero_spectrum();
-#  endif
-  Spectrum accum_emission = zero_spectrum();
+  volume_integrate_result_init(state, ray, vstate, equiangular_coeffs, result);

  OctreeTracing octree(ray->tmin);
  const uint32_t path_flag = INTEGRATOR_STATE(state, path, flag);
-  for (; volume_integrate_advance<false>(kg, ray, sd, state, rng_state, path_flag, octree, vstep);)
-  {
-    /* compute segment */
-    VolumeShaderCoefficients coeff ccl_optional_struct_init;
-    if (volume_shader_sample(kg, state, sd, &coeff)) {
-      const int closure_flag = sd->flag;
-
-      /* Evaluate transmittance over segment. */
-      const float dt = vstep.t.length();
-      const Spectrum transmittance = (closure_flag & SD_EXTINCTION) ?
-                                         volume_color_transmittance(coeff.sigma_t, dt) :
-                                         one_spectrum();
-
-      /* Emission. */
-      if (closure_flag & SD_EMISSION) {
-        /* Only write emission before indirect light scatter position, since we terminate
-         * stepping at that point if we have already found a direct light scatter position. */
-        if (!result.indirect_scatter) {
-          const Spectrum emission = volume_emission_integrate(
-              &coeff, closure_flag, transmittance, dt);
-          accum_emission += result.indirect_throughput * emission;
-          guiding_record_volume_emission(kg, state, emission);
-        }
-      }
-
-      if (closure_flag & SD_SCATTER) {
-#  ifdef __DENOISING_FEATURES__
-        /* Accumulate albedo for denoising features. */
-        if (write_denoising_features && (closure_flag & SD_SCATTER)) {
-          const Spectrum albedo = safe_divide_color(coeff.sigma_s, coeff.sigma_t);
-          accum_albedo += result.indirect_throughput * albedo * (one_spectrum() - transmittance);
-        }
-#  endif
-
-        /* Scattering and absorption. */
-        volume_integrate_step_scattering(
-            sd, ray, equiangular_coeffs, coeff, transmittance, vstep.t, vstate, result);
-      }
-      else if (closure_flag & SD_EXTINCTION) {
-        /* Absorption only. */
-        result.indirect_throughput *= transmittance;
-        result.direct_throughput *= transmittance;
-      }
-
-      if (volume_integrate_should_stop(result)) {
-        break;
-      }
-    }
+  if (!volume_octree_setup<false>(kg, ray, sd, state, rng_state, path_flag, octree)) {
+    return;
  }

+  /* Scramble for stepping through volume. */
+  path_state_rng_scramble(rng_state, 0xe35fad82);
+
+  do {
+    volume_integrate_step_scattering(
+        kg, state, ray, octree.sigma, octree.t, sd, rng_state, vstate, result);
+
+    if (volume_integrate_should_stop(result, vstate)) {
+      break;
+    }
+
+  } while (
+      volume_octree_advance<false>(kg, ray, sd, state, rng_state, path_flag, octree, vstate.step));
+
+  volume_direct_scatter_mis(ray, vstate, equiangular_coeffs, result);
+
  /* Write accumulated emission. */
-  if (!is_zero(accum_emission)) {
+  if (!is_zero(vstate.emission)) {
    if (light_link_object_match(kg, light_link_receiver_forward(kg, state), sd->object)) {
      film_write_volume_emission(
-          kg, state, accum_emission, render_buffer, object_lightgroup(kg, sd->object));
+          kg, state, vstate.emission, render_buffer, object_lightgroup(kg, sd->object));
    }
  }

 #  ifdef __DENOISING_FEATURES__
  /* Write denoising features. */
-  if (write_denoising_features) {
+  if (INTEGRATOR_STATE(state, path, flag) & PATH_RAY_DENOISING_FEATURES) {
    film_write_denoising_features_volume(
-        kg, state, accum_albedo, result.indirect_scatter, render_buffer);
+        kg, state, vstate.albedo, result.indirect_scatter, render_buffer);
  }
 #  endif /* __DENOISING_FEATURES__ */
 }
@@ -1727,6 +1774,8 @@ ccl_device void integrator_shade_volume(KernelGlobals kg,
  }
 #  endif /* __SHADOW_LINKING__ */

+  kernel_assert(event == VOLUME_PATH_SCATTERED);
+
  /* Queue intersect_closest kernel. */
  integrator_path_next(
      state, DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME, DEVICE_KERNEL_INTEGRATOR_INTERSECT_CLOSEST);
--- a/intern/cycles/kernel/sample/mapping.h
+++ b/intern/cycles/kernel/sample/mapping.h
@@ -221,4 +221,11 @@ ccl_device_inline int sample_geometric_distribution(const float rand,
  return n;
 }

+/* Generate random variable x following exponential distribution p(x) = lambda * exp(-lambda * x),
+ * where lambda > 0 is the rate parameter. */
+ccl_device_inline float sample_exponential_distribution(const float rand, const float inv_lambda)
+{
+  return -logf(1.0f - rand) * inv_lambda;
+}
+
 CCL_NAMESPACE_END
--- a/intern/cycles/kernel/types.h
+++ b/intern/cycles/kernel/types.h
@@ -293,11 +293,10 @@ enum PathTraceDimension {
  PRNG_VOLUME_PHASE = 3,
  PRNG_VOLUME_COLOR_CHANNEL = 4,
  PRNG_VOLUME_SCATTER_DISTANCE = 5,
-  PRNG_VOLUME_OFFSET = 6,
+  PRNG_VOLUME_EXPANSION_ORDER = 6,
  PRNG_VOLUME_SHADE_OFFSET = 7,
  PRNG_VOLUME_PHASE_GUIDING_DISTANCE = 8,
  PRNG_VOLUME_PHASE_GUIDING_EQUIANGULAR = 9,
-  PRNG_VOLUME_EXPANSION_ORDER = 12,

  /* Subsurface random walk bounces */
  PRNG_SUBSURFACE_BSDF = 0,
--- a/intern/cycles/scene/background.cpp
+++ b/intern/cycles/scene/background.cpp
@@ -28,8 +28,6 @@ NODE_DEFINE(Background)
  SOCKET_BOOLEAN(transparent_glass, "Transparent Glass", false);
  SOCKET_FLOAT(transparent_roughness_threshold, "Transparent Roughness Threshold", 0.0f);

-  SOCKET_FLOAT(volume_step_size, "Volume Step Size", 0.1f);
-
  SOCKET_NODE(shader, "Shader", Shader::get_node_type());

  SOCKET_STRING(lightgroup, "Light Group", ustring());
@@ -86,8 +84,6 @@ void Background::device_update(Device *device, DeviceScene *dscene, Scene *scene
    kbackground->volume_shader = SHADER_NONE;
  }

-  kbackground->volume_step_size = volume_step_size * scene->integrator->get_volume_step_rate();
-
  /* No background node, make world shader invisible to all rays, to skip evaluation in kernel. */
  if (bg_shader->graph->nodes.size() <= 1) {
    kbackground->surface_shader |= SHADER_EXCLUDE_ANY;
--- a/intern/cycles/scene/background.h
+++ b/intern/cycles/scene/background.h
@@ -28,8 +28,6 @@ class Background : public Node {
  NODE_SOCKET_API(bool, transparent_glass)
  NODE_SOCKET_API(float, transparent_roughness_threshold)

-  NODE_SOCKET_API(float, volume_step_size)
-
  NODE_SOCKET_API(ustring, lightgroup)

  Background();
--- a/intern/cycles/scene/devicescene.cpp
+++ b/intern/cycles/scene/devicescene.cpp
@@ -30,7 +30,6 @@ DeviceScene::DeviceScene(Device *device)
      object_motion_pass(device, "object_motion_pass", MEM_GLOBAL),
      object_motion(device, "object_motion", MEM_GLOBAL),
      object_flag(device, "object_flag", MEM_GLOBAL),
-      object_volume_step(device, "object_volume_step", MEM_GLOBAL),
      object_prim_offset(device, "object_prim_offset", MEM_GLOBAL),
      camera_motion(device, "camera_motion", MEM_GLOBAL),
      attributes_map(device, "attributes_map", MEM_GLOBAL),
--- a/intern/cycles/scene/devicescene.h
+++ b/intern/cycles/scene/devicescene.h
@@ -42,7 +42,6 @@ class DeviceScene {
  device_vector<Transform> object_motion_pass;
  device_vector<DecomposedTransform> object_motion;
  device_vector<uint> object_flag;
-  device_vector<float> object_volume_step;
  device_vector<uint> object_prim_offset;

  /* cameras */
--- a/intern/cycles/scene/integrator.cpp
+++ b/intern/cycles/scene/integrator.cpp
@@ -57,8 +57,7 @@ NODE_DEFINE(Integrator)
  SOCKET_FLOAT(ao_distance, "AO Distance", FLT_MAX);
  SOCKET_FLOAT(ao_additive_factor, "AO Additive Factor", 0.0f);

-  SOCKET_INT(volume_max_steps, "Volume Max Steps", 1024);
-  SOCKET_FLOAT(volume_step_rate, "Volume Step Rate", 1.0f);
+  SOCKET_BOOLEAN(volume_unbiased, "Unbiased", false);

  static NodeEnum guiding_distribution_enum;
  guiding_distribution_enum.insert("PARALLAX_AWARE_VMM", GUIDING_TYPE_PARALLAX_AWARE_VMM);
@@ -227,8 +226,7 @@ void Integrator::device_update(Device *device, DeviceScene *dscene, Scene *scene
    }
  }

-  kintegrator->volume_max_steps = volume_max_steps;
-  kintegrator->volume_step_rate = volume_step_rate;
+  kintegrator->volume_unbiased = volume_unbiased;

  kintegrator->caustics_reflective = caustics_reflective;
  kintegrator->caustics_refractive = caustics_refractive;
--- a/intern/cycles/scene/integrator.h
+++ b/intern/cycles/scene/integrator.h
@@ -41,8 +41,7 @@ class Integrator : public Node {
  NODE_SOCKET_API(float, ao_distance)
  NODE_SOCKET_API(float, ao_additive_factor)

-  NODE_SOCKET_API(int, volume_max_steps)
-  NODE_SOCKET_API(float, volume_step_rate)
+  NODE_SOCKET_API(bool, volume_unbiased)

  NODE_SOCKET_API(bool, use_guiding);
  NODE_SOCKET_API(bool, deterministic_guiding);
--- a/intern/cycles/scene/object.cpp
+++ b/intern/cycles/scene/object.cpp
@@ -51,7 +51,6 @@ struct UpdateObjectTransformState {
  KernelObject *objects;
  Transform *object_motion_pass;
  DecomposedTransform *object_motion;
-  float *object_volume_step;

  /* Flags which will be synchronized to Integrator. */
  bool have_motion;
@@ -290,86 +289,6 @@ uint Object::visibility_for_tracing() const
  return SHADOW_CATCHER_OBJECT_VISIBILITY(is_shadow_catcher, visibility & PATH_RAY_ALL_VISIBILITY);
 }

-float Object::compute_volume_step_size() const
-{
-  if (geometry->geometry_type != Geometry::MESH && geometry->geometry_type != Geometry::VOLUME) {
-    return FLT_MAX;
-  }
-
-  Mesh *mesh = static_cast<Mesh *>(geometry);
-
-  if (!mesh->has_volume) {
-    return FLT_MAX;
-  }
-
-  /* Compute step rate from shaders. */
-  float step_rate = FLT_MAX;
-
-  for (Node *node : mesh->get_used_shaders()) {
-    Shader *shader = static_cast<Shader *>(node);
-    if (shader->has_volume) {
-      if ((shader->get_heterogeneous_volume() && shader->has_volume_spatial_varying) ||
-          (shader->has_volume_attribute_dependency))
-      {
-        step_rate = fminf(shader->get_volume_step_rate(), step_rate);
-      }
-    }
-  }
-
-  if (step_rate == FLT_MAX) {
-    return FLT_MAX;
-  }
-
-  /* Compute step size from voxel grids. */
-  float step_size = FLT_MAX;
-
-  if (geometry->is_volume()) {
-    Volume *volume = static_cast<Volume *>(geometry);
-
-    for (Attribute &attr : volume->attributes.attributes) {
-      if (attr.element == ATTR_ELEMENT_VOXEL) {
-        ImageHandle &handle = attr.data_voxel();
-        const ImageMetaData &metadata = handle.metadata();
-        if (metadata.byte_size == 0) {
-          continue;
-        }
-
-        /* User specified step size. */
-        float voxel_step_size = volume->get_step_size();
-
-        if (voxel_step_size == 0.0f) {
-          /* Auto detect step size.
-           * Step size is transformed from voxel to world space. */
-          Transform voxel_tfm = tfm;
-          if (metadata.use_transform_3d) {
-            voxel_tfm = tfm * transform_inverse(metadata.transform_3d);
-          }
-          voxel_step_size = reduce_min(fabs(transform_direction(&voxel_tfm, one_float3())));
-        }
-        else if (volume->get_object_space()) {
-          /* User specified step size in object space. */
-          const float3 size = make_float3(voxel_step_size, voxel_step_size, voxel_step_size);
-          voxel_step_size = reduce_min(fabs(transform_direction(&tfm, size)));
-        }
-
-        if (voxel_step_size > 0.0f) {
-          step_size = fminf(voxel_step_size, step_size);
-        }
-      }
-    }
-  }
-
-  if (step_size == FLT_MAX) {
-    /* Fall back to 1/10th of bounds for procedural volumes. */
-    assert(bounds.valid());
-    step_size = 0.1f * average(bounds.size());
-  }
-
-  step_size *= step_rate;
-
-  return step_size;
-}
-
 int Object::get_device_index() const
 {
  return index;
@@ -619,7 +538,6 @@ void ObjectManager::device_update_object_transform(UpdateObjectTransformState *s
    flag |= SD_OBJECT_HOLDOUT_MASK;
  }
  state->object_flag[ob->index] = flag;
-  state->object_volume_step[ob->index] = FLT_MAX;

  /* Have curves. */
  if (geom->is_hair()) {
@@ -688,7 +606,6 @@ void ObjectManager::device_update_transforms(DeviceScene *dscene, Scene *scene,

  state.objects = dscene->objects.alloc(scene->objects.size());
  state.object_flag = dscene->object_flag.alloc(scene->objects.size());
-  state.object_volume_step = dscene->object_volume_step.alloc(scene->objects.size());
  state.object_motion = nullptr;
  state.object_motion_pass = nullptr;

@@ -772,7 +689,6 @@ void ObjectManager::device_update(Device *device,
    dscene->object_motion_pass.tag_realloc();
    dscene->object_motion.tag_realloc();
    dscene->object_flag.tag_realloc();
-    dscene->object_volume_step.tag_realloc();

    /* If objects are added to the scene or deleted, the object indices might change, so we need to
     * update the root indices of the volume octrees. */
@@ -814,7 +730,6 @@ void ObjectManager::device_update(Device *device,
        dscene->object_motion_pass.tag_modified();
        dscene->object_motion.tag_modified();
        dscene->object_flag.tag_modified();
-        dscene->object_volume_step.tag_modified();
      }

      /* Update world object index. */
@@ -879,28 +794,17 @@ void ObjectManager::device_update_flags(Device * /*unused*/,

  /* Object info flag. */
  uint *object_flag = dscene->object_flag.data();
-  float *object_volume_step = dscene->object_volume_step.data();

  /* Object volume intersection. */
  vector<Object *> volume_objects;
  bool has_volume_objects = false;
  for (Object *object : scene->objects) {
    if (object->geometry->has_volume) {
-      /* If the bounds are not valid it is not always possible to calculate the volume step, and
-       * the step size is not needed for the displacement. So, delay calculation of the volume
-       * step size until the final bounds are known. */
      if (bounds_valid) {
        volume_objects.push_back(object);
-        object_volume_step[object->index] = object->compute_volume_step_size();
-      }
-      else {
-        object_volume_step[object->index] = FLT_MAX;
      }
      has_volume_objects = true;
    }
-    else {
-      object_volume_step[object->index] = FLT_MAX;
-    }
  }

  for (Object *object : scene->objects) {
@@ -948,10 +852,8 @@ void ObjectManager::device_update_flags(Device * /*unused*/,

  /* Copy object flag. */
  dscene->object_flag.copy_to_device();
-  dscene->object_volume_step.copy_to_device();

  dscene->object_flag.clear_modified();
-  dscene->object_volume_step.clear_modified();
 }

 void ObjectManager::device_update_geom_offsets(Device * /*unused*/,
@@ -1001,7 +903,6 @@ void ObjectManager::device_free(Device * /*unused*/, DeviceScene *dscene, bool f
  dscene->object_motion_pass.free_if_need_realloc(force_free);
  dscene->object_motion.free_if_need_realloc(force_free);
  dscene->object_flag.free_if_need_realloc(force_free);
-  dscene->object_volume_step.free_if_need_realloc(force_free);
  dscene->object_prim_offset.free_if_need_realloc(force_free);
 }

--- a/intern/cycles/scene/object.h
+++ b/intern/cycles/scene/object.h
@@ -112,9 +112,6 @@ class Object : public Node {
  /* Returns the index that is used in the kernel for this object. */
  int get_device_index() const;

-  /* Compute step size from attributes, shaders, transforms. */
-  float compute_volume_step_size() const;
-
  /* Check whether this object can be used as light-emissive. */
  bool usable_as_light() const;

--- a/intern/cycles/scene/shader.cpp
+++ b/intern/cycles/scene/shader.cpp
@@ -54,7 +54,6 @@ NODE_DEFINE(Shader)

  SOCKET_BOOLEAN(use_transparent_shadow, "Use Transparent Shadow", true);
  SOCKET_BOOLEAN(use_bump_map_correction, "Bump Map Correction", true);
-  SOCKET_BOOLEAN(heterogeneous_volume, "Heterogeneous Volume", true);

  static NodeEnum volume_sampling_method_enum;
  volume_sampling_method_enum.insert("distance", VOLUME_SAMPLING_DISTANCE);
@@ -73,8 +72,6 @@ NODE_DEFINE(Shader)
              volume_interpolation_method_enum,
              VOLUME_INTERPOLATION_LINEAR);

-  SOCKET_FLOAT(volume_step_rate, "Volume Step Rate", 1.0f);
-
  static NodeEnum displacement_method_enum;
  displacement_method_enum.insert("bump", DISPLACE_BUMP);
  displacement_method_enum.insert("true", DISPLACE_TRUE);
@@ -104,7 +101,6 @@ Shader::Shader() : Node(get_node_type())
  has_volume_spatial_varying = false;
  has_volume_attribute_dependency = false;
  has_volume_connected = false;
-  prev_volume_step_rate = 0.0f;
  has_light_path_node = false;

  emission_estimate = zero_float3();
@@ -394,10 +390,9 @@ void Shader::tag_update(Scene *scene)
    scene->procedural_manager->tag_update();
  }

-  if (has_volume != prev_has_volume || volume_step_rate != prev_volume_step_rate) {
+  if (has_volume != prev_has_volume) {
    scene->geometry_manager->need_flags_update = true;
    scene->object_manager->need_flags_update = true;
-    prev_volume_step_rate = volume_step_rate;
  }

  if (has_volume || prev_has_volume) {
@@ -613,10 +608,8 @@ void ShaderManager::device_update_common(Device * /*device*/,
    if (shader->has_volume_connected && !shader->has_surface) {
      flag |= SD_HAS_ONLY_VOLUME;
    }
-    if (shader->has_volume) {
-      if (shader->get_heterogeneous_volume() && shader->has_volume_spatial_varying) {
-        flag |= SD_HETEROGENEOUS_VOLUME;
-      }
+    if (shader->has_volume && shader->has_volume_spatial_varying) {
+      flag |= SD_HETEROGENEOUS_VOLUME;
    }
    if (shader->has_volume_attribute_dependency) {
      flag |= SD_NEED_VOLUME_ATTRIBUTES;
--- a/intern/cycles/scene/shader.h
+++ b/intern/cycles/scene/shader.h
@@ -81,16 +81,12 @@ class Shader : public Node {
  NODE_SOCKET_API(EmissionSampling, emission_sampling_method)
  NODE_SOCKET_API(bool, use_transparent_shadow)
  NODE_SOCKET_API(bool, use_bump_map_correction)
-  NODE_SOCKET_API(bool, heterogeneous_volume)
  NODE_SOCKET_API(VolumeSampling, volume_sampling_method)
  NODE_SOCKET_API(int, volume_interpolation_method)
-  NODE_SOCKET_API(float, volume_step_rate)

  /* displacement */
  NODE_SOCKET_API(DisplacementMethod, displacement_method)

-  float prev_volume_step_rate;
-
  /* synchronization */
  bool need_update_uvs;
  bool need_update_attribute;
--- a/intern/cycles/scene/volume.cpp
+++ b/intern/cycles/scene/volume.cpp
@@ -33,7 +33,6 @@ NODE_DEFINE(Volume)
 {
  NodeType *type = NodeType::add("volume", create, NodeType::NONE, Mesh::get_node_type());

-  SOCKET_FLOAT(step_size, "Step Size", 0.0f);
  SOCKET_BOOLEAN(object_space, "Object Space", false);
  SOCKET_FLOAT(velocity_scale, "Velocity Scale", 1.0f);

@@ -42,7 +41,6 @@ NODE_DEFINE(Volume)

 Volume::Volume() : Mesh(get_node_type(), Geometry::VOLUME)
 {
-  step_size = 0.0f;
  object_space = false;
 }

--- a/intern/cycles/scene/volume.h
+++ b/intern/cycles/scene/volume.h
@@ -23,7 +23,6 @@ class Volume : public Mesh {

  Volume();

-  NODE_SOCKET_API(float, step_size)
  NODE_SOCKET_API(bool, object_space)
  NODE_SOCKET_API(float, velocity_scale)