From a7283fc1d5334149fb142e67805a59383abb273e Mon Sep 17 00:00:00 2001
From: Weizhen Huang <weizhen@blender.org>
Date: Tue, 11 Feb 2025 19:29:51 +0100
Subject: [PATCH] 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.
---
 intern/cycles/blender/addon/properties.py     |  54 +-
 intern/cycles/blender/addon/ui.py             |  15 +-
 intern/cycles/blender/shader.cpp              |   4 -
 intern/cycles/blender/sync.cpp                |   7 +-
 intern/cycles/blender/volume.cpp              |   1 -
 intern/cycles/kernel/closure/volume.h         |   6 +-
 intern/cycles/kernel/data_arrays.h            |   1 -
 intern/cycles/kernel/data_template.h          |   5 +-
 intern/cycles/kernel/geom/object.h            |  11 +-
 .../cycles/kernel/integrator/shade_volume.h   | 719 ++++++++++--------
 intern/cycles/kernel/sample/mapping.h         |   7 +
 intern/cycles/kernel/types.h                  |   3 +-
 intern/cycles/scene/background.cpp            |   4 -
 intern/cycles/scene/background.h              |   2 -
 intern/cycles/scene/devicescene.cpp           |   1 -
 intern/cycles/scene/devicescene.h             |   1 -
 intern/cycles/scene/integrator.cpp            |   6 +-
 intern/cycles/scene/integrator.h              |   3 +-
 intern/cycles/scene/object.cpp                |  99 ---
 intern/cycles/scene/object.h                  |   3 -
 intern/cycles/scene/shader.cpp                |  13 +-
 intern/cycles/scene/shader.h                  |   4 -
 intern/cycles/scene/volume.cpp                |   2 -
 intern/cycles/scene/volume.h                  |   1 -
 .../startup/bl_ui/properties_data_volume.py   |   2 -
 source/blender/makesrna/intern/rna_volume.cc  |  11 -
 26 files changed, 413 insertions(+), 572 deletions(-)

diff --git a/intern/cycles/blender/addon/properties.py b/intern/cycles/blender/addon/properties.py
index 8426c35df19..d43cae2cf39 100644
--- 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):
diff --git a/intern/cycles/blender/addon/ui.py b/intern/cycles/blender/addon/ui.py
index c22a2c4c008..da47762bdae 100644
--- 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)
diff --git a/intern/cycles/blender/shader.cpp b/intern/cycles/blender/shader.cpp
index f77558cc6a2..3a0b21a0af2 100644
--- 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>();
diff --git a/intern/cycles/blender/sync.cpp b/intern/cycles/blender/sync.cpp
index e6207c71310..f843696181b 100644
--- 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"));
diff --git a/intern/cycles/blender/volume.cpp b/intern/cycles/blender/volume.cpp
index 7b3e52fbd76..1e19baa7c92 100644
--- 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();
diff --git a/intern/cycles/kernel/closure/volume.h b/intern/cycles/kernel/closure/volume.h
index c969f56cf49..56f68da40bc 100644
--- 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)
diff --git a/intern/cycles/kernel/data_arrays.h b/intern/cycles/kernel/data_arrays.h
index bea3f2320aa..43c3d57e707 100644
--- 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 */
diff --git a/intern/cycles/kernel/data_template.h b/intern/cycles/kernel/data_template.h
index 67d87e4aeaa..772ecfcf768 100644
--- 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. */
diff --git a/intern/cycles/kernel/geom/object.h b/intern/cycles/kernel/geom/object.h
index 2c656d90c3c..78aa7f2bdab 100644
--- 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)
diff --git a/intern/cycles/kernel/integrator/shade_volume.h b/intern/cycles/kernel/integrator/shade_volume.h
index f5ffe864a3d..4c2c5954d14 100644
--- 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);
diff --git a/intern/cycles/kernel/sample/mapping.h b/intern/cycles/kernel/sample/mapping.h
index 603b0d0d4eb..f25e5cbd22f 100644
--- 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
diff --git a/intern/cycles/kernel/types.h b/intern/cycles/kernel/types.h
index c81fe4daca3..9efaebec28e 100644
--- 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,
diff --git a/intern/cycles/scene/background.cpp b/intern/cycles/scene/background.cpp
index 0983a02e36b..56d32e047d9 100644
--- 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;
diff --git a/intern/cycles/scene/background.h b/intern/cycles/scene/background.h
index af3d297a4ee..ba2f4713ea2 100644
--- 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();
diff --git a/intern/cycles/scene/devicescene.cpp b/intern/cycles/scene/devicescene.cpp
index 1ad1cee9276..27aa20ba1d1 100644
--- 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),
diff --git a/intern/cycles/scene/devicescene.h b/intern/cycles/scene/devicescene.h
index fbb45c970da..9aaa7bbc22e 100644
--- 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 */
diff --git a/intern/cycles/scene/integrator.cpp b/intern/cycles/scene/integrator.cpp
index d915af809cb..2610b4dafe2 100644
--- 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;
diff --git a/intern/cycles/scene/integrator.h b/intern/cycles/scene/integrator.h
index 8860f52e240..bb2ed42c792 100644
--- 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);
diff --git a/intern/cycles/scene/object.cpp b/intern/cycles/scene/object.cpp
index 3d75a8ba76d..f96920a6220 100644
--- 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);
 }
 
diff --git a/intern/cycles/scene/object.h b/intern/cycles/scene/object.h
index 026bc8dbcc2..170fcef1e6b 100644
--- 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;
 
diff --git a/intern/cycles/scene/shader.cpp b/intern/cycles/scene/shader.cpp
index 4bdd5908dc4..9678623e78c 100644
--- 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;
diff --git a/intern/cycles/scene/shader.h b/intern/cycles/scene/shader.h
index d74237f0ca1..ecaef5cc827 100644
--- 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;
diff --git a/intern/cycles/scene/volume.cpp b/intern/cycles/scene/volume.cpp
index 7465acb54b8..25ffe8e7884 100644
--- 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;
 }
 
diff --git a/intern/cycles/scene/volume.h b/intern/cycles/scene/volume.h
index 396e1dceb0c..46fdfbfe7fb 100644
--- 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)
 
diff --git a/scripts/startup/bl_ui/properties_data_volume.py b/scripts/startup/bl_ui/properties_data_volume.py
index 00f4af7792a..d6ddf6958e1 100644
--- a/scripts/startup/bl_ui/properties_data_volume.py
+++ b/scripts/startup/bl_ui/properties_data_volume.py
@@ -128,8 +128,6 @@ class DATA_PT_volume_render(DataButtonsPanel, Panel):
         col.prop(render, "space")
 
         if scene.render.engine == 'CYCLES':
-            col.prop(render, "step_size")
-
             col = layout.column(align=True)
             col.prop(render, "clipping")
 
diff --git a/source/blender/makesrna/intern/rna_volume.cc b/source/blender/makesrna/intern/rna_volume.cc
index 1cd735c5df8..00b371e2a82 100644
--- a/source/blender/makesrna/intern/rna_volume.cc
+++ b/source/blender/makesrna/intern/rna_volume.cc
@@ -539,17 +539,6 @@ static void rna_def_volume_render(BlenderRNA *brna)
       prop, "Space", "Specify volume density and step size in object or world space");
   RNA_def_property_update(prop, 0, "rna_Volume_update_display");
 
-  prop = RNA_def_property(srna, "step_size", PROP_FLOAT, PROP_DISTANCE);
-  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
-  RNA_def_property_range(prop, 0.0, FLT_MAX);
-  RNA_def_property_ui_range(prop, 0.0, 100.0, 1, 3);
-  RNA_def_property_ui_text(prop,
-                           "Step Size",
-                           "Distance between volume samples. Lower values render more detail at "
-                           "the cost of performance. If set to zero, the step size is "
-                           "automatically determined based on voxel size.");
-  RNA_def_property_update(prop, 0, "rna_Volume_update_display");
-
   prop = RNA_def_property(srna, "clipping", PROP_FLOAT, PROP_NONE);
   RNA_def_property_float_sdna(prop, nullptr, "clipping");
   RNA_def_property_range(prop, 0.0, 1.0);