From 2b0a1cae0610611c0cbb98218a3247712b0e05f5 Mon Sep 17 00:00:00 2001
From: Weizhen Huang <weizhen@blender.org>
Date: Fri, 26 Sep 2025 12:14:45 +0200
Subject: [PATCH] Cycles: Add an option to use ray marching for volume
 rendering

Null Scattering currently has performance and noise issues, and it will
take time to address them. For now add the previous Ray Marching back as
an option.

Co-authored-by: Brecht Van Lommel <brecht@blender.org>
Pull Request: https://projects.blender.org/blender/blender/pulls/146317
---
 intern/cycles/blender/addon/properties.py     |  45 +-
 intern/cycles/blender/addon/ui.py             |  13 +-
 intern/cycles/blender/shader.cpp              |   2 +
 intern/cycles/blender/sync.cpp                |   7 +-
 intern/cycles/blender/volume.cpp              |   1 +
 intern/cycles/device/kernel.cpp               |   3 +
 intern/cycles/device/oneapi/device_impl.cpp   |   1 +
 intern/cycles/device/optix/device_impl.cpp    |   4 +
 intern/cycles/device/optix/device_impl.h      |   1 +
 intern/cycles/device/optix/queue.cpp          |   5 +
 .../cycles/integrator/path_trace_work_gpu.cpp |   3 +-
 intern/cycles/kernel/data_arrays.h            |   1 +
 intern/cycles/kernel/data_template.h          |   4 +-
 intern/cycles/kernel/device/gpu/kernel.h      |  15 +
 intern/cycles/kernel/device/oneapi/kernel.cpp |   9 +
 .../cycles/kernel/device/optix/kernel_osl.cu  |   9 +
 .../kernel/integrator/intersect_closest.h     |   8 +-
 intern/cycles/kernel/integrator/megakernel.h  |   3 +
 .../cycles/kernel/integrator/shade_shadow.h   |   8 +-
 .../cycles/kernel/integrator/shade_volume.h   | 656 ++++++++++++++++--
 .../cycles/kernel/integrator/volume_stack.h   |  22 +
 intern/cycles/kernel/types.h                  |   3 +
 intern/cycles/scene/devicescene.cpp           |   3 +-
 intern/cycles/scene/devicescene.h             |   1 +
 intern/cycles/scene/integrator.cpp            |   7 +-
 intern/cycles/scene/integrator.h              |   4 +-
 intern/cycles/scene/object.cpp                |  97 ++-
 intern/cycles/scene/object.h                  |   3 +
 intern/cycles/scene/shader.cpp                |   6 +-
 intern/cycles/scene/shader.h                  |   3 +
 intern/cycles/scene/volume.cpp                | 147 +++-
 intern/cycles/scene/volume.h                  |  12 +
 .../startup/bl_ui/properties_data_volume.py   |   2 +
 source/blender/makesrna/intern/rna_volume.cc  |  11 +
 .../camera_in_volume_camera_and_volume.png    |   3 +
 .../camera_in_volume_camera_only.png          |   3 +
 .../camera_in_volume_nested.png               |   3 +
 .../camera_in_volume_simple.png               |   3 +
 .../implicit_volume.png                       |   3 +
 .../overlapping_different_anisotropy.png      |   3 +
 .../overlapping_octrees.png                   |   3 +
 .../principled_absorption.png                 |   3 +
 .../principled_bsdf_interior.png              |   3 +
 .../volume_absobtion.png                      |   3 +
 .../volume_deep_stack.png                     |   3 +
 .../volume_edge_fireflies.png                 |   3 +
 .../volume_instance.png                       |   3 +
 .../volume_light_path.png                     |   3 +
 .../volume_output_absorption.png              |   3 +
 .../volume_output_mix_volume.png              |   3 +
 .../volume_output_scatter.png                 |   3 +
 .../volume_output_surface.png                 |   3 +
 .../volume_overlap.png                        |   3 +
 .../volume_scatter.png                        |   3 +
 .../volume_scatter_albedo.png                 |   3 +
 .../volume_scatter_draine.png                 |   3 +
 .../volume_scatter_fournier-forand.png        |   3 +
 .../volume_scatter_mie.png                    |   3 +
 .../volume_scatter_mie_small_particles.png    |   3 +
 .../volume_scatter_rayleigh.png               |   3 +
 .../volume_step_offset.png                    |   3 +
 .../volume_transparent_shadow.png             |   3 +
 .../volume_two_shaders.png                    |   3 +
 .../volume_zero_extinction_channel.png        |   3 +
 .../world_volume.png                          |   3 +
 tests/python/cycles_render_tests.py           |  35 +-
 66 files changed, 1136 insertions(+), 111 deletions(-)
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_and_volume.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_only.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_nested.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_simple.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/implicit_volume.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/overlapping_different_anisotropy.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/overlapping_octrees.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/principled_absorption.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/principled_bsdf_interior.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_absobtion.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_deep_stack.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_edge_fireflies.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_instance.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_light_path.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_output_absorption.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_output_mix_volume.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_output_scatter.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_output_surface.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_overlap.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_scatter.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_albedo.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_draine.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_fournier-forand.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie_small_particles.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_rayleigh.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_step_offset.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_transparent_shadow.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_two_shaders.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/volume_zero_extinction_channel.png
 create mode 100644 tests/files/render/volume/cycles_ray_marching_renders/world_volume.png

diff --git a/intern/cycles/blender/addon/properties.py b/intern/cycles/blender/addon/properties.py
index ca3936c10ad..329ed77424a 100644
--- a/intern/cycles/blender/addon/properties.py
+++ b/intern/cycles/blender/addon/properties.py
@@ -779,10 +779,32 @@ class CyclesRenderSettings(bpy.types.PropertyGroup):
         default=8,
     )
 
-    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",
+    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_biased: BoolProperty(
+        name="Biased",
+        description="Default volume rendering uses null scattering, which is unbiased and has less artifacts, "
+        "but could be noisier. Biased option uses ray marching, with controls for steps size and max steps",
         default=False,
     )
 
@@ -1134,6 +1156,14 @@ 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(
@@ -1228,6 +1258,13 @@ 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 fc181744180..102bfb7cda3 100644
--- a/intern/cycles/blender/addon/ui.py
+++ b/intern/cycles/blender/addon/ui.py
@@ -584,8 +584,13 @@ class CYCLES_RENDER_PT_volumes(CyclesButtonsPanel, Panel):
         scene = context.scene
         cscene = scene.cycles
 
-        col = layout.column()
-        col.prop(cscene, "volume_unbiased", text="Unbiased")
+        col = layout.column(align=True)
+        col.prop(cscene, "volume_biased", text="Biased")
+        if cscene.volume_biased:
+            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")
 
 
 class CYCLES_RENDER_PT_light_paths(CyclesButtonsPanel, Panel):
@@ -1836,6 +1841,8 @@ 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")
+        if context.scene.cycles.volume_biased:
+            col.prop(cworld, "volume_step_size")
 
 
 class CYCLES_WORLD_PT_settings_light_group(CyclesButtonsPanel, Panel):
@@ -2008,6 +2015,8 @@ 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")
+        if context.scene.cycles.volume_biased:
+            col.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 17634d334c3..efb2e7bcaa5 100644
--- a/intern/cycles/blender/shader.cpp
+++ b/intern/cycles/blender/shader.cpp
@@ -1625,6 +1625,7 @@ void BlenderSync::sync_materials(BL::Depsgraph &b_depsgraph, bool update_all)
       shader->set_use_bump_map_correction(get_boolean(cmat, "use_bump_map_correction"));
       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));
@@ -1691,6 +1692,7 @@ void BlenderSync::sync_world(BL::Depsgraph &b_depsgraph, BL::SpaceView3D &b_v3d,
       PointerRNA cworld = RNA_pointer_get(&b_world.ptr, "cycles");
       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 0e24821409e..239b234dea0 100644
--- a/intern/cycles/blender/sync.cpp
+++ b/intern/cycles/blender/sync.cpp
@@ -348,10 +348,15 @@ 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_ray_marching(get_boolean(cscene, "volume_biased"));
+  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 1e19baa7c92..7b3e52fbd76 100644
--- a/intern/cycles/blender/volume.cpp
+++ b/intern/cycles/blender/volume.cpp
@@ -284,6 +284,7 @@ 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/device/kernel.cpp b/intern/cycles/device/kernel.cpp
index bc8481789e6..ec59bbfd7ca 100644
--- a/intern/cycles/device/kernel.cpp
+++ b/intern/cycles/device/kernel.cpp
@@ -18,6 +18,7 @@ bool device_kernel_has_shading(DeviceKernel kernel)
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_MNEE ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME ||
+          kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_SHADOW ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT ||
           kernel == DEVICE_KERNEL_SHADER_EVAL_DISPLACE ||
@@ -69,6 +70,8 @@ const char *device_kernel_as_string(DeviceKernel kernel)
       return "integrator_shade_surface_mnee";
     case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME:
       return "integrator_shade_volume";
+    case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING:
+      return "integrator_shade_volume_ray_marching";
     case DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT:
       return "integrator_shade_dedicated_light";
     case DEVICE_KERNEL_INTEGRATOR_MEGAKERNEL:
diff --git a/intern/cycles/device/oneapi/device_impl.cpp b/intern/cycles/device/oneapi/device_impl.cpp
index 77536dae6dc..8363cb8bb2f 100644
--- a/intern/cycles/device/oneapi/device_impl.cpp
+++ b/intern/cycles/device/oneapi/device_impl.cpp
@@ -1284,6 +1284,7 @@ void OneapiDevice::get_adjusted_global_and_local_sizes(SyclQueue *queue,
     case DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_MNEE:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME:
+    case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_SHADOW:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT: {
       const bool device_is_simd8 =
diff --git a/intern/cycles/device/optix/device_impl.cpp b/intern/cycles/device/optix/device_impl.cpp
index 228761ea2fe..d28afbdd2c7 100644
--- a/intern/cycles/device/optix/device_impl.cpp
+++ b/intern/cycles/device/optix/device_impl.cpp
@@ -574,6 +574,10 @@ bool OptiXDevice::load_kernels(const uint kernel_features)
     group_descs[PG_RGEN_SHADE_VOLUME].raygen.module = optix_module;
     group_descs[PG_RGEN_SHADE_VOLUME].raygen.entryFunctionName =
         "__raygen__kernel_optix_integrator_shade_volume";
+    group_descs[PG_RGEN_SHADE_VOLUME_RAY_MARCHING].kind = OPTIX_PROGRAM_GROUP_KIND_RAYGEN;
+    group_descs[PG_RGEN_SHADE_VOLUME_RAY_MARCHING].raygen.module = optix_module;
+    group_descs[PG_RGEN_SHADE_VOLUME_RAY_MARCHING].raygen.entryFunctionName =
+        "__raygen__kernel_optix_integrator_shade_volume_ray_marching";
     group_descs[PG_RGEN_SHADE_SHADOW].kind = OPTIX_PROGRAM_GROUP_KIND_RAYGEN;
     group_descs[PG_RGEN_SHADE_SHADOW].raygen.module = optix_module;
     group_descs[PG_RGEN_SHADE_SHADOW].raygen.entryFunctionName =
diff --git a/intern/cycles/device/optix/device_impl.h b/intern/cycles/device/optix/device_impl.h
index f1db63feb9b..35626c1f940 100644
--- a/intern/cycles/device/optix/device_impl.h
+++ b/intern/cycles/device/optix/device_impl.h
@@ -32,6 +32,7 @@ enum {
   PG_RGEN_SHADE_SURFACE_RAYTRACE,
   PG_RGEN_SHADE_SURFACE_MNEE,
   PG_RGEN_SHADE_VOLUME,
+  PG_RGEN_SHADE_VOLUME_RAY_MARCHING,
   PG_RGEN_SHADE_SHADOW,
   PG_RGEN_SHADE_DEDICATED_LIGHT,
   PG_RGEN_EVAL_DISPLACE,
diff --git a/intern/cycles/device/optix/queue.cpp b/intern/cycles/device/optix/queue.cpp
index d86667e5230..be7cb6f51f8 100644
--- a/intern/cycles/device/optix/queue.cpp
+++ b/intern/cycles/device/optix/queue.cpp
@@ -124,6 +124,11 @@ bool OptiXDeviceQueue::enqueue(DeviceKernel kernel,
       pipeline = optix_device->pipelines[PIP_SHADE];
       sbt_params.raygenRecord = sbt_data_ptr + PG_RGEN_SHADE_VOLUME * sizeof(SbtRecord);
       break;
+    case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING:
+      pipeline = optix_device->pipelines[PIP_SHADE];
+      sbt_params.raygenRecord = sbt_data_ptr +
+                                PG_RGEN_SHADE_VOLUME_RAY_MARCHING * sizeof(SbtRecord);
+      break;
     case DEVICE_KERNEL_INTEGRATOR_SHADE_SHADOW:
       pipeline = optix_device->pipelines[PIP_SHADE];
       sbt_params.raygenRecord = sbt_data_ptr + PG_RGEN_SHADE_SHADOW * sizeof(SbtRecord);
diff --git a/intern/cycles/integrator/path_trace_work_gpu.cpp b/intern/cycles/integrator/path_trace_work_gpu.cpp
index af029e8ac42..9f5b147b52a 100644
--- a/intern/cycles/integrator/path_trace_work_gpu.cpp
+++ b/intern/cycles/integrator/path_trace_work_gpu.cpp
@@ -562,6 +562,7 @@ void PathTraceWorkGPU::enqueue_path_iteration(DeviceKernel kernel, const int num
     case DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_MNEE:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME:
+    case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING:
     case DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT: {
       /* Shading kernels with integrator state and render buffer. */
       const DeviceKernelArguments args(
@@ -570,7 +571,6 @@ void PathTraceWorkGPU::enqueue_path_iteration(DeviceKernel kernel, const int num
       queue_->enqueue(kernel, work_size, args);
       break;
     }
-
     default:
       LOG_FATAL << "Unhandled kernel " << device_kernel_as_string(kernel)
                 << " used for path iteration, should never happen.";
@@ -1272,6 +1272,7 @@ bool PathTraceWorkGPU::kernel_creates_shadow_paths(DeviceKernel kernel)
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_MNEE ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME ||
+          kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING ||
           kernel == DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT);
 }
 
diff --git a/intern/cycles/kernel/data_arrays.h b/intern/cycles/kernel/data_arrays.h
index 43c3d57e707..f4ad2abd8de 100644
--- a/intern/cycles/kernel/data_arrays.h
+++ b/intern/cycles/kernel/data_arrays.h
@@ -88,5 +88,6 @@ KERNEL_DATA_ARRAY(float, ies)
 KERNEL_DATA_ARRAY(KernelOctreeNode, volume_tree_nodes)
 KERNEL_DATA_ARRAY(KernelOctreeRoot, volume_tree_roots)
 KERNEL_DATA_ARRAY(int, volume_tree_root_ids)
+KERNEL_DATA_ARRAY(float, volume_step_size)
 
 #undef KERNEL_DATA_ARRAY
diff --git a/intern/cycles/kernel/data_template.h b/intern/cycles/kernel/data_template.h
index 74ad79f5db8..13eb1a6153d 100644
--- a/intern/cycles/kernel/data_template.h
+++ b/intern/cycles/kernel/data_template.h
@@ -208,7 +208,8 @@ 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_unbiased)
+KERNEL_STRUCT_MEMBER(integrator, int, volume_ray_marching)
+KERNEL_STRUCT_MEMBER(integrator, int, volume_max_steps)
 /* Shadow catcher. */
 KERNEL_STRUCT_MEMBER(integrator, int, has_shadow_catcher)
 /* Closure filter. */
@@ -231,7 +232,6 @@ KERNEL_STRUCT_MEMBER(integrator, int, use_guiding_mis_weights)
 /* Padding. */
 KERNEL_STRUCT_MEMBER(integrator, int, pad1)
 KERNEL_STRUCT_MEMBER(integrator, int, pad2)
-KERNEL_STRUCT_MEMBER(integrator, int, pad3)
 KERNEL_STRUCT_END(KernelIntegrator)
 
 /* SVM. For shader specialization. */
diff --git a/intern/cycles/kernel/device/gpu/kernel.h b/intern/cycles/kernel/device/gpu/kernel.h
index f237258925d..efcb06fb7f9 100644
--- a/intern/cycles/kernel/device/gpu/kernel.h
+++ b/intern/cycles/kernel/device/gpu/kernel.h
@@ -354,6 +354,21 @@ ccl_gpu_kernel(GPU_KERNEL_BLOCK_NUM_THREADS, GPU_KERNEL_MAX_REGISTERS)
 }
 ccl_gpu_kernel_postfix
 
+ccl_gpu_kernel(GPU_KERNEL_BLOCK_NUM_THREADS, GPU_KERNEL_MAX_REGISTERS)
+    ccl_gpu_kernel_signature(integrator_shade_volume_ray_marching,
+                             const ccl_global int *path_index_array,
+                             ccl_global float *render_buffer,
+                             const int work_size)
+{
+  const int global_index = ccl_gpu_global_id_x();
+
+  if (ccl_gpu_kernel_within_bounds(global_index, work_size)) {
+    const int state = (path_index_array) ? path_index_array[global_index] : global_index;
+    ccl_gpu_kernel_call(integrator_shade_volume_ray_marching(nullptr, state, render_buffer));
+  }
+}
+ccl_gpu_kernel_postfix
+
 ccl_gpu_kernel(GPU_KERNEL_BLOCK_NUM_THREADS, GPU_KERNEL_MAX_REGISTERS)
     ccl_gpu_kernel_signature(integrator_shade_dedicated_light,
                              const ccl_global int *path_index_array,
diff --git a/intern/cycles/kernel/device/oneapi/kernel.cpp b/intern/cycles/kernel/device/oneapi/kernel.cpp
index fbf2d0940c3..df477246a3a 100644
--- a/intern/cycles/kernel/device/oneapi/kernel.cpp
+++ b/intern/cycles/kernel/device/oneapi/kernel.cpp
@@ -460,6 +460,15 @@ bool oneapi_enqueue_kernel(KernelContext *kernel_context,
               kg, cgh, global_size, local_size, args, oneapi_kernel_integrator_shade_volume);
           break;
         }
+        case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING: {
+          oneapi_call(kg,
+                      cgh,
+                      global_size,
+                      local_size,
+                      args,
+                      oneapi_kernel_integrator_shade_volume_ray_marching);
+          break;
+        }
         case DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT: {
           oneapi_call(kg,
                       cgh,
diff --git a/intern/cycles/kernel/device/optix/kernel_osl.cu b/intern/cycles/kernel/device/optix/kernel_osl.cu
index d32845ae3bc..faea22dda83 100644
--- a/intern/cycles/kernel/device/optix/kernel_osl.cu
+++ b/intern/cycles/kernel/device/optix/kernel_osl.cu
@@ -53,6 +53,15 @@ extern "C" __global__ void __raygen__kernel_optix_integrator_shade_volume()
   integrator_shade_volume(nullptr, path_index, kernel_params.render_buffer);
 }
 
+extern "C" __global__ void __raygen__kernel_optix_integrator_shade_volume_ray_marching()
+{
+  const int global_index = optixGetLaunchIndex().x;
+  const int path_index = (kernel_params.path_index_array) ?
+                             kernel_params.path_index_array[global_index] :
+                             global_index;
+  integrator_shade_volume_ray_marching(nullptr, path_index, kernel_params.render_buffer);
+}
+
 extern "C" __global__ void __raygen__kernel_optix_integrator_shade_shadow()
 {
   const int global_index = optixGetLaunchIndex().x;
diff --git a/intern/cycles/kernel/integrator/intersect_closest.h b/intern/cycles/kernel/integrator/intersect_closest.h
index 4977766d4bb..633abb4334c 100644
--- a/intern/cycles/kernel/integrator/intersect_closest.h
+++ b/intern/cycles/kernel/integrator/intersect_closest.h
@@ -233,7 +233,13 @@ ccl_device_forceinline void integrator_intersect_next_kernel(
     const int flags = (hit_surface) ? kernel_data_fetch(shaders, shader).flags : 0;
 
     if (!integrator_intersect_terminate(kg, state, flags)) {
-      integrator_path_next(state, current_kernel, DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME);
+      if (kernel_data.integrator.volume_ray_marching) {
+        integrator_path_next(
+            state, current_kernel, DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING);
+      }
+      else {
+        integrator_path_next(state, current_kernel, DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME);
+      }
     }
     else {
       integrator_path_terminate(kg, state, render_buffer, current_kernel);
diff --git a/intern/cycles/kernel/integrator/megakernel.h b/intern/cycles/kernel/integrator/megakernel.h
index 3e843cb0172..96b9791e134 100644
--- a/intern/cycles/kernel/integrator/megakernel.h
+++ b/intern/cycles/kernel/integrator/megakernel.h
@@ -76,6 +76,9 @@ ccl_device void integrator_megakernel(KernelGlobals kg,
         case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME:
           integrator_shade_volume(kg, state, render_buffer);
           break;
+        case DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING:
+          integrator_shade_volume_ray_marching(kg, state, render_buffer);
+          break;
         case DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE:
           integrator_shade_surface_raytrace(kg, state, render_buffer);
           break;
diff --git a/intern/cycles/kernel/integrator/shade_shadow.h b/intern/cycles/kernel/integrator/shade_shadow.h
index af5c0972f2d..daf2dbf78ed 100644
--- a/intern/cycles/kernel/integrator/shade_shadow.h
+++ b/intern/cycles/kernel/integrator/shade_shadow.h
@@ -94,7 +94,13 @@ ccl_device_inline void integrate_transparent_volume_shadow(KernelGlobals kg,
   /* `object` is only needed for light tree with light linking, it is irrelevant for shadow. */
   shader_setup_from_volume(shadow_sd, &ray, OBJECT_NONE);
 
-  volume_shadow_heterogeneous(kg, state, &ray, shadow_sd, throughput);
+  if (kernel_data.integrator.volume_ray_marching) {
+    const float step_size = volume_stack_step_size<true>(kg, state);
+    volume_shadow_ray_marching(kg, state, &ray, shadow_sd, throughput, step_size);
+  }
+  else {
+    volume_shadow_null_scattering(kg, state, &ray, shadow_sd, throughput);
+  }
 }
 #  endif
 
diff --git a/intern/cycles/kernel/integrator/shade_volume.h b/intern/cycles/kernel/integrator/shade_volume.h
index d4dc26ea72e..b12be528cf7 100644
--- a/intern/cycles/kernel/integrator/shade_volume.h
+++ b/intern/cycles/kernel/integrator/shade_volume.h
@@ -636,11 +636,11 @@ ccl_device Spectrum volume_transmittance(KernelGlobals kg,
 
 /* Compute the volumetric transmittance of the segment [ray->tmin, ray->tmax],
  * used for the shadow ray throughput. */
-ccl_device void volume_shadow_heterogeneous(KernelGlobals kg,
-                                            IntegratorShadowState state,
-                                            ccl_private Ray *ccl_restrict ray,
-                                            ccl_private ShaderData *ccl_restrict sd,
-                                            ccl_private Spectrum *ccl_restrict throughput)
+ccl_device void volume_shadow_null_scattering(KernelGlobals kg,
+                                              IntegratorShadowState state,
+                                              ccl_private Ray *ccl_restrict ray,
+                                              ccl_private ShaderData *ccl_restrict sd,
+                                              ccl_private Spectrum *ccl_restrict throughput)
 {
   /* Load random number state. */
   RNGState rng_state;
@@ -1891,14 +1891,14 @@ volume_direct_sample_method(KernelGlobals kg,
 }
 
 /* Shared function of integrating homogeneous and heterogeneous volume. */
-ccl_device void volume_integrate_shared(KernelGlobals kg,
-                                        const IntegratorState state,
-                                        const ccl_private Ray *ccl_restrict ray,
-                                        ccl_private ShaderData *ccl_restrict sd,
-                                        const ccl_private RNGState *rng_state,
-                                        ccl_global float *ccl_restrict render_buffer,
-                                        ccl_private LightSample *ls,
-                                        ccl_private VolumeIntegrateResult &result)
+ccl_device void volume_integrate_null_scattering(KernelGlobals kg,
+                                                 const IntegratorState state,
+                                                 const ccl_private Ray *ccl_restrict ray,
+                                                 ccl_private ShaderData *ccl_restrict sd,
+                                                 const ccl_private RNGState *rng_state,
+                                                 ccl_global float *ccl_restrict render_buffer,
+                                                 ccl_private LightSample *ls,
+                                                 ccl_private VolumeIntegrateResult &result)
 {
   PROFILING_INIT(kg, PROFILING_SHADE_VOLUME_INTEGRATE);
 
@@ -1944,6 +1944,426 @@ ccl_device void volume_integrate_shared(KernelGlobals kg,
   }
 }
 
+/* -------------------------------------------------------------------- */
+/** \name Ray Marching
+ * \{ */
+
+/* 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;
+
+  /* Maximal steps allowed between `ray->tmin` and `ray->tmax`. */
+  int max_steps;
+
+  /* Current active segment. */
+  Interval<float> t;
+};
+
+template<const bool shadow>
+ccl_device_forceinline void volume_step_init(KernelGlobals kg,
+                                             const ccl_private RNGState *rng_state,
+                                             const float object_step_size,
+                                             const float tmin,
+                                             const float tmax,
+                                             ccl_private VolumeStep *vstep)
+{
+  vstep->t.min = vstep->t.max = tmin;
+
+  if (object_step_size == FLT_MAX) {
+    /* Homogeneous volume. */
+    vstep->size = tmax - tmin;
+    vstep->shade_offset = 0.0f;
+    vstep->offset = 1.0f;
+    vstep->max_steps = 1;
+  }
+  else {
+    /* Heterogeneous volume. */
+    vstep->max_steps = kernel_data.integrator.volume_max_steps;
+    const float t = tmax - tmin;
+    float step_size = min(object_step_size, t);
+
+    if (t > vstep->max_steps * step_size) {
+      /* Increase step size to cover the whole ray segment. */
+      step_size = t / (float)vstep->max_steps;
+    }
+
+    vstep->size = step_size;
+    vstep->shade_offset = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_SHADE_OFFSET);
+
+    if (shadow) {
+      /* For shadows we do not offset all segments, since the starting point is already a random
+       * distance inside the volume. It also appears to create banding artifacts for unknown
+       * reasons. */
+      vstep->offset = 1.0f;
+    }
+    else {
+      vstep->offset = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_OFFSET);
+    }
+  }
+}
+
+ccl_device_inline bool volume_ray_marching_advance(const int step,
+                                                   const ccl_private Ray *ccl_restrict ray,
+                                                   ccl_private float3 *shade_P,
+                                                   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;
+  vstep.t.max = min(ray->tmax, ray->tmin + (step + vstep.offset) * vstep.size);
+  const float shade_t = mix(vstep.t.min, vstep.t.max, vstep.shade_offset);
+  *shade_P = ray->P + ray->D * shade_t;
+
+  return step < vstep.max_steps;
+}
+
+ccl_device void volume_shadow_ray_marching(KernelGlobals kg,
+                                           IntegratorShadowState state,
+                                           ccl_private Ray *ccl_restrict ray,
+                                           ccl_private ShaderData *ccl_restrict sd,
+                                           ccl_private Spectrum *ccl_restrict throughput,
+                                           const float object_step_size)
+{
+  /* Load random number state. */
+  RNGState rng_state;
+  shadow_path_state_rng_load(state, &rng_state);
+
+  /* For stochastic texture sampling. */
+  sd->lcg_state = lcg_state_init(
+      rng_state.rng_pixel, rng_state.rng_offset, rng_state.sample, 0xd9111870);
+
+  Spectrum tp = *throughput;
+
+  /* Prepare for stepping. */
+  VolumeStep vstep;
+  volume_step_init<true>(kg, &rng_state, object_step_size, ray->tmin, ray->tmax, &vstep);
+
+  /* compute extinction at the start */
+  Spectrum sum = zero_spectrum();
+  for (int step = 0; volume_ray_marching_advance(step, ray, &sd->P, vstep); step++) {
+    /* compute attenuation over segment */
+    const Spectrum sigma_t = volume_shader_eval_extinction<true>(kg, state, sd, PATH_RAY_SHADOW);
+    /* Compute `expf()` only for every Nth step, to save some calculations
+     * because `exp(a)*exp(b) = exp(a+b)`, also do a quick #VOLUME_THROUGHPUT_EPSILON
+     * check then. */
+    sum += (-sigma_t * vstep.t.length());
+    if ((step & 0x07) == 0) { /* TODO: Other interval? */
+      tp = *throughput * exp(sum);
+
+      /* stop if nearly all light is blocked */
+      if (reduce_max(tp) < VOLUME_THROUGHPUT_EPSILON) {
+        break;
+      }
+    }
+  }
+
+  if (vstep.t.max == ray->tmax) {
+    /* Update throughput in case we haven't done it above. */
+    tp = *throughput * exp(sum);
+  }
+
+  *throughput = tp;
+}
+
+struct VolumeRayMarchingState {
+  /* Random numbers for scattering. */
+  float rscatter;
+  float rchannel;
+
+  /* Multiple importance sampling. */
+  VolumeSampleMethod direct_sample_method;
+  bool use_mis;
+  float distance_pdf;
+  float equiangular_pdf;
+};
+
+ccl_device_inline void volume_ray_marching_state_init(
+    KernelGlobals kg,
+    const ccl_private RNGState *rng_state,
+    const VolumeSampleMethod direct_sample_method,
+    ccl_private VolumeRayMarchingState &vstate)
+{
+  vstate.rscatter = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_SCATTER_DISTANCE);
+  vstate.rchannel = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_COLOR_CHANNEL);
+
+  /* Multiple importance sampling: pick between equiangular and distance sampling strategy. */
+  vstate.direct_sample_method = direct_sample_method;
+  vstate.use_mis = (direct_sample_method == VOLUME_SAMPLE_MIS);
+  if (vstate.use_mis) {
+    if (vstate.rscatter < 0.5f) {
+      vstate.rscatter *= 2.0f;
+      vstate.direct_sample_method = VOLUME_SAMPLE_DISTANCE;
+    }
+    else {
+      vstate.rscatter = (vstate.rscatter - 0.5f) * 2.0f;
+      vstate.direct_sample_method = VOLUME_SAMPLE_EQUIANGULAR;
+    }
+  }
+  vstate.equiangular_pdf = 0.0f;
+  vstate.distance_pdf = 1.0f;
+}
+
+/* Returns true if we found the indirect scatter position within the current active ray segment. */
+ccl_device bool volume_sample_indirect_scatter_ray_marching(
+    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,
+    ccl_private VolumeRayMarchingState &ccl_restrict vstate,
+    ccl_private VolumeIntegrateResult &ccl_restrict result)
+{
+  if (result.indirect_scatter) {
+    /* Already sampled indirect scatter position. */
+    return false;
+  }
+
+  /* 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);
+
+    /* 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;
+
+    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);
+
+    if (vstate.distance_pdf * distance_pdf > VOLUME_SAMPLE_PDF_CUTOFF) {
+      /* Update throughput. */
+      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;
+    }
+
+    /* Remap rscatter so we can reuse it and keep thing stratified. */
+    vstate.rscatter = 1.0f - (1.0f - vstate.rscatter) / sample_transmittance;
+  }
+
+  return false;
+}
+
+/* Find direct and indirect scatter positions. */
+ccl_device_forceinline void volume_ray_marching_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,
+    const ccl_private Interval<float> &t,
+    ccl_private VolumeRayMarchingState &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);
+
+  /* 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);
+
+      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);
+    }
+  }
+
+  /* Distance sampling for indirect and optional direct lighting. */
+  if (volume_sample_indirect_scatter_ray_marching(
+          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;
+      }
+    }
+  }
+}
+
+/* heterogeneous volume distance sampling: integrate stepping through the
+ * volume until we reach the end, get absorbed entirely, or run out of
+ * iterations. this does probabilistically scatter or get transmitted through
+ * for path tracing where we don't want to branch. */
+ccl_device_forceinline void volume_integrate_ray_marching(
+    KernelGlobals kg,
+    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_global float *ccl_restrict render_buffer,
+    const float object_step_size,
+    ccl_private LightSample *ls,
+    ccl_private VolumeIntegrateResult &result)
+{
+  PROFILING_INIT(kg, PROFILING_SHADE_VOLUME_INTEGRATE);
+
+  EquiangularCoefficients equiangular_coeffs = {zero_float3(), {ray->tmin, ray->tmax}};
+  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<false>(kg, rng_state, object_step_size, ray->tmin, ray->tmax, &vstep);
+
+  /* Initialize volume integration state. */
+  VolumeRayMarchingState vstate ccl_optional_struct_init;
+  volume_ray_marching_state_init(kg, rng_state, direct_sample_method, 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();
+
+  for (int step = 0; volume_ray_marching_advance(step, ray, &sd->P, vstep); step++) {
+    /* 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, 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_ray_marching_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;
+      }
+    }
+  }
+
+  /* Write accumulated emission. */
+  if (!is_zero(accum_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));
+    }
+  }
+
+#  ifdef __DENOISING_FEATURES__
+  /* Write denoising features. */
+  if (write_denoising_features) {
+    film_write_denoising_features_volume(
+        kg, state, accum_albedo, result.indirect_scatter, render_buffer);
+  }
+#  endif /* __DENOISING_FEATURES__ */
+}
+
+/** \} */
+
 /* Path tracing: sample point on light and evaluate light shader, then
  * queue shadow ray to be traced. */
 ccl_device_forceinline void integrate_volume_direct_light(
@@ -2199,38 +2619,15 @@ ccl_device_forceinline bool integrate_volume_phase_scatter(
   return true;
 }
 
-/* get the volume attenuation and emission over line segment defined by
- * ray, with the assumption that there are no surfaces blocking light
- * between the endpoints. distance sampling is used to decide if we will
- * scatter or not. */
-ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
-                                                 IntegratorState state,
-                                                 ccl_private Ray *ccl_restrict ray,
-                                                 ccl_global float *ccl_restrict render_buffer)
+ccl_device_inline VolumeIntegrateEvent
+volume_integrate_event(KernelGlobals kg,
+                       IntegratorState state,
+                       const ccl_private Ray *ccl_restrict ray,
+                       ccl_private ShaderData *sd,
+                       const ccl_private RNGState *rng_state,
+                       ccl_private LightSample &ls,
+                       ccl_private VolumeIntegrateResult &result)
 {
-  if (integrator_state_volume_stack_is_empty(kg, state)) {
-    return VOLUME_PATH_ATTENUATED;
-  }
-
-  ShaderData sd;
-  /* FIXME: `object` is used for light linking. We read the bottom of the stack for simplicity, but
-   * this does not work for overlapping volumes. */
-  shader_setup_from_volume(&sd, ray, INTEGRATOR_STATE_ARRAY(state, volume_stack, 0, object));
-
-  /* Load random number state. */
-  RNGState rng_state;
-  path_state_rng_load(state, &rng_state);
-
-  /* For stochastic texture sampling. */
-  sd.lcg_state = lcg_state_init(
-      rng_state.rng_pixel, rng_state.rng_offset, rng_state.sample, 0x15b4f88d);
-
-  LightSample ls ccl_optional_struct_init;
-
-  /* TODO: expensive to zero closures? */
-  VolumeIntegrateResult result = {};
-  volume_integrate_shared(kg, state, ray, &sd, &rng_state, render_buffer, &ls, result);
-
 #  if defined(__PATH_GUIDING__) && PATH_GUIDING_LEVEL >= 1
   /* The current path throughput which is used later to calculate per-segment throughput. */
   const float3 initial_throughput = INTEGRATOR_STATE(state, path, throughput);
@@ -2267,10 +2664,10 @@ ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
         const float3 transmittance_weight = spectrum_to_rgb(
             safe_divide_color(result.indirect_throughput, initial_throughput));
         guiding_record_volume_transmission(kg, state, transmittance_weight);
-        guiding_record_volume_segment(kg, state, direct_P, sd.wi);
+        guiding_record_volume_segment(kg, state, direct_P, sd->wi);
         guiding_generated_new_segment = true;
         unlit_throughput = result.indirect_throughput / continuation_probability;
-        rand_phase_guiding = path_state_rng_1D(kg, &rng_state, PRNG_VOLUME_PHASE_GUIDING_DISTANCE);
+        rand_phase_guiding = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_PHASE_GUIDING_DISTANCE);
       }
       else if (result.direct_sample_method == VOLUME_SAMPLE_EQUIANGULAR) {
         /* If the direct scatter event is generated using VOLUME_SAMPLE_EQUIANGULAR the direct
@@ -2287,7 +2684,7 @@ ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
         unlit_throughput /= scatterEval;
         unlit_throughput *= continuation_probability;
         rand_phase_guiding = path_state_rng_1D(
-            kg, &rng_state, PRNG_VOLUME_PHASE_GUIDING_EQUIANGULAR);
+            kg, rng_state, PRNG_VOLUME_PHASE_GUIDING_EQUIANGULAR);
       }
 #    endif
 #    if PATH_GUIDING_LEVEL >= 4
@@ -2302,8 +2699,8 @@ ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
     result.direct_throughput /= continuation_probability;
     integrate_volume_direct_light(kg,
                                   state,
-                                  &sd,
-                                  &rng_state,
+                                  sd,
+                                  rng_state,
                                   direct_P,
                                   &result.direct_phases,
 #  if defined(__PATH_GUIDING__)
@@ -2330,28 +2727,28 @@ ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
   INTEGRATOR_STATE_WRITE(state, path, throughput) = result.indirect_throughput;
 
   if (result.indirect_scatter) {
-    sd.P = ray->P + result.indirect_t * ray->D;
+    sd->P = ray->P + result.indirect_t * ray->D;
 
 #  if defined(__PATH_GUIDING__)
     if ((kernel_data.kernel_features & KERNEL_FEATURE_PATH_GUIDING)) {
 #    if PATH_GUIDING_LEVEL >= 1
       if (!guiding_generated_new_segment) {
-        guiding_record_volume_segment(kg, state, sd.P, sd.wi);
+        guiding_record_volume_segment(kg, state, sd->P, sd->wi);
       }
 #    endif
 #    if PATH_GUIDING_LEVEL >= 4
       /* If the direct scatter event was generated using VOLUME_SAMPLE_EQUIANGULAR we need to
        * initialize the guiding distribution at the indirect scatter position. */
       if (result.direct_sample_method == VOLUME_SAMPLE_EQUIANGULAR) {
-        rand_phase_guiding = path_state_rng_1D(kg, &rng_state, PRNG_VOLUME_PHASE_GUIDING_DISTANCE);
+        rand_phase_guiding = path_state_rng_1D(kg, rng_state, PRNG_VOLUME_PHASE_GUIDING_DISTANCE);
         volume_shader_prepare_guiding(
-            kg, state, rand_phase_guiding, sd.P, ray->D, &result.indirect_phases);
+            kg, state, rand_phase_guiding, sd->P, ray->D, &result.indirect_phases);
       }
 #    endif
     }
 #  endif
 
-    if (integrate_volume_phase_scatter(kg, state, &sd, ray, &rng_state, &result.indirect_phases)) {
+    if (integrate_volume_phase_scatter(kg, state, sd, ray, rng_state, &result.indirect_phases)) {
       return VOLUME_PATH_SCATTERED;
     }
     return VOLUME_PATH_MISSED;
@@ -2365,27 +2762,99 @@ ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
   return VOLUME_PATH_ATTENUATED;
 }
 
+/* get the volume attenuation and emission over line segment defined by
+ * ray, with the assumption that there are no surfaces blocking light
+ * between the endpoints. distance sampling is used to decide if we will
+ * scatter or not. */
+ccl_device VolumeIntegrateEvent volume_integrate(KernelGlobals kg,
+                                                 IntegratorState state,
+                                                 ccl_private Ray *ccl_restrict ray,
+                                                 ccl_global float *ccl_restrict render_buffer)
+{
+  kernel_assert(!kernel_data.integrator.volume_ray_marching);
+
+  if (integrator_state_volume_stack_is_empty(kg, state)) {
+    return VOLUME_PATH_ATTENUATED;
+  }
+
+  ShaderData sd;
+  /* FIXME: `object` is used for light linking. We read the bottom of the stack for simplicity, but
+   * this does not work for overlapping volumes. */
+  shader_setup_from_volume(&sd, ray, INTEGRATOR_STATE_ARRAY(state, volume_stack, 0, object));
+
+  /* Load random number state. */
+  RNGState rng_state;
+  path_state_rng_load(state, &rng_state);
+
+  /* For stochastic texture sampling. */
+  sd.lcg_state = lcg_state_init(
+      rng_state.rng_pixel, rng_state.rng_offset, rng_state.sample, 0x15b4f88d);
+
+  LightSample ls ccl_optional_struct_init;
+
+  /* TODO: expensive to zero closures? */
+  VolumeIntegrateResult result = {};
+  volume_integrate_null_scattering(kg, state, ray, &sd, &rng_state, render_buffer, &ls, result);
+
+  return volume_integrate_event(kg, state, ray, &sd, &rng_state, ls, result);
+}
+
+ccl_device VolumeIntegrateEvent
+volume_integrate_ray_marching(KernelGlobals kg,
+                              IntegratorState state,
+                              ccl_private Ray *ccl_restrict ray,
+                              ccl_global float *ccl_restrict render_buffer)
+{
+  kernel_assert(kernel_data.integrator.volume_ray_marching);
+
+  if (integrator_state_volume_stack_is_empty(kg, state)) {
+    return VOLUME_PATH_ATTENUATED;
+  }
+
+  ShaderData sd;
+  /* FIXME: `object` is used for light linking. We read the bottom of the stack for simplicity, but
+   * this does not work for overlapping volumes. */
+  shader_setup_from_volume(&sd, ray, INTEGRATOR_STATE_ARRAY(state, volume_stack, 0, object));
+
+  /* Load random number state. */
+  RNGState rng_state;
+  path_state_rng_load(state, &rng_state);
+
+  /* For stochastic texture sampling. */
+  sd.lcg_state = lcg_state_init(
+      rng_state.rng_pixel, rng_state.rng_offset, rng_state.sample, 0x15b4f88d);
+
+  LightSample ls ccl_optional_struct_init;
+
+  /* TODO: expensive to zero closures? */
+  VolumeIntegrateResult result = {};
+
+  const float step_size = volume_stack_step_size<false>(kg, state);
+  volume_integrate_ray_marching(
+      kg, state, ray, &sd, &rng_state, render_buffer, step_size, &ls, result);
+
+  return volume_integrate_event(kg, state, ray, &sd, &rng_state, ls, result);
+}
+
 #endif
 
-ccl_device void integrator_shade_volume(KernelGlobals kg,
-                                        IntegratorState state,
-                                        ccl_global float *ccl_restrict render_buffer)
+#ifdef __VOLUME__
+ccl_device_inline void integrator_shade_volume_setup(KernelGlobals kg,
+                                                     IntegratorState state,
+                                                     ccl_private Ray *ccl_restrict ray,
+                                                     ccl_private Intersection *ccl_restrict isect)
 {
   PROFILING_INIT(kg, PROFILING_SHADE_VOLUME_SETUP);
 
-#ifdef __VOLUME__
   /* Setup shader data. */
-  Ray ray ccl_optional_struct_init;
-  integrator_state_read_ray(state, &ray);
-
-  Intersection isect ccl_optional_struct_init;
-  integrator_state_read_isect(state, &isect);
+  integrator_state_read_ray(state, ray);
+  integrator_state_read_isect(state, isect);
 
   /* Set ray length to current segment. */
-  ray.tmax = (isect.prim != PRIM_NONE) ? isect.t : FLT_MAX;
+  ray->tmax = (isect->prim != PRIM_NONE) ? isect->t : FLT_MAX;
 
   /* Clean volume stack for background rays. */
-  if (isect.prim == PRIM_NONE) {
+  if (isect->prim == PRIM_NONE) {
     volume_stack_clean(kg, state);
   }
 
@@ -2393,35 +2862,70 @@ ccl_device void integrator_shade_volume(KernelGlobals kg,
   if (INTEGRATOR_STATE(state, path, bounce) == 0) {
     INTEGRATOR_STATE_WRITE(state, path, flag) |= PATH_RAY_VOLUME_PRIMARY_TRANSMIT;
   }
+}
+#endif
 
-  const VolumeIntegrateEvent event = volume_integrate(kg, state, &ray, render_buffer);
+template<DeviceKernel volume_kernel>
+ccl_device_inline void integrator_next_kernel_after_shade_volume(
+    KernelGlobals kg,
+    const IntegratorState state,
+    ccl_global float *ccl_restrict render_buffer,
+    const ccl_private Intersection *ccl_restrict isect,
+    const VolumeIntegrateEvent event)
+{
   if (event == VOLUME_PATH_MISSED) {
     /* End path. */
-    integrator_path_terminate(kg, state, render_buffer, DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME);
+    integrator_path_terminate(kg, state, render_buffer, volume_kernel);
     return;
   }
 
   if (event == VOLUME_PATH_ATTENUATED) {
     /* Continue to background, light or surface. */
-    integrator_intersect_next_kernel_after_volume<DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME>(
-        kg, state, &isect, render_buffer);
+    integrator_intersect_next_kernel_after_volume<volume_kernel>(kg, state, isect, render_buffer);
     return;
   }
 
-#  ifdef __SHADOW_LINKING__
-  if (shadow_linking_schedule_intersection_kernel<DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME>(kg,
-                                                                                         state))
-  {
+#ifdef __SHADOW_LINKING__
+  if (shadow_linking_schedule_intersection_kernel<volume_kernel>(kg, state)) {
     return;
   }
-#  endif /* __SHADOW_LINKING__ */
+#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);
+  integrator_path_next(state, volume_kernel, DEVICE_KERNEL_INTEGRATOR_INTERSECT_CLOSEST);
+}
+
+ccl_device void integrator_shade_volume(KernelGlobals kg,
+                                        IntegratorState state,
+                                        ccl_global float *ccl_restrict render_buffer)
+{
+#ifdef __VOLUME__
+  Ray ray ccl_optional_struct_init;
+  Intersection isect ccl_optional_struct_init;
+  integrator_shade_volume_setup(kg, state, &ray, &isect);
+
+  const VolumeIntegrateEvent event = volume_integrate(kg, state, &ray, render_buffer);
+  integrator_next_kernel_after_shade_volume<DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME>(
+      kg, state, render_buffer, &isect, event);
+
 #endif /* __VOLUME__ */
 }
 
+ccl_device void integrator_shade_volume_ray_marching(KernelGlobals kg,
+                                                     IntegratorState state,
+                                                     ccl_global float *ccl_restrict render_buffer)
+{
+#ifdef __VOLUME__
+  Ray ray ccl_optional_struct_init;
+  Intersection isect ccl_optional_struct_init;
+  integrator_shade_volume_setup(kg, state, &ray, &isect);
+
+  const VolumeIntegrateEvent event = volume_integrate_ray_marching(kg, state, &ray, render_buffer);
+  integrator_next_kernel_after_shade_volume<DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING>(
+      kg, state, render_buffer, &isect, event);
+
+#endif /* __VOLUME__ */
+}
 CCL_NAMESPACE_END
diff --git a/intern/cycles/kernel/integrator/volume_stack.h b/intern/cycles/kernel/integrator/volume_stack.h
index f45b6dc37da..35b0cc9f462 100644
--- a/intern/cycles/kernel/integrator/volume_stack.h
+++ b/intern/cycles/kernel/integrator/volume_stack.h
@@ -185,6 +185,28 @@ ccl_device_inline bool volume_is_homogeneous(KernelGlobals kg, const IntegratorG
   return false;
 }
 
+template<const bool shadow, typename IntegratorGenericState>
+ccl_device float volume_stack_step_size(KernelGlobals kg, const IntegratorGenericState state)
+{
+  kernel_assert(kernel_data.integrator.volume_ray_marching);
+
+  float step_size = FLT_MAX;
+
+  for (int i = 0;; i++) {
+    const VolumeStack entry = volume_stack_read<shadow>(state, i);
+    if (entry.shader == SHADER_NONE) {
+      break;
+    }
+
+    if (!volume_is_homogeneous(kg, entry)) {
+      const float object_step_size = kernel_data_fetch(volume_step_size, entry.object);
+      step_size = fminf(object_step_size, step_size);
+    }
+  }
+
+  return step_size;
+}
+
 enum VolumeSampleMethod {
   VOLUME_SAMPLE_NONE = 0,
   VOLUME_SAMPLE_DISTANCE = (1 << 0),
diff --git a/intern/cycles/kernel/types.h b/intern/cycles/kernel/types.h
index 55de43cedb5..629eeaa97b0 100644
--- a/intern/cycles/kernel/types.h
+++ b/intern/cycles/kernel/types.h
@@ -297,6 +297,8 @@ enum PathTraceDimension {
   PRNG_VOLUME_SHADE_OFFSET = 7,
   PRNG_VOLUME_PHASE_GUIDING_DISTANCE = 8,
   PRNG_VOLUME_PHASE_GUIDING_EQUIANGULAR = 9,
+  PRNG_VOLUME_COLOR_CHANNEL = 4,
+  PRNG_VOLUME_OFFSET = 6,
 
   /* Subsurface random walk bounces */
   PRNG_SUBSURFACE_BSDF = 0,
@@ -1853,6 +1855,7 @@ enum DeviceKernel : int {
   DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_RAYTRACE,
   DEVICE_KERNEL_INTEGRATOR_SHADE_SURFACE_MNEE,
   DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME,
+  DEVICE_KERNEL_INTEGRATOR_SHADE_VOLUME_RAY_MARCHING,
   DEVICE_KERNEL_INTEGRATOR_SHADE_SHADOW,
   DEVICE_KERNEL_INTEGRATOR_SHADE_DEDICATED_LIGHT,
   DEVICE_KERNEL_INTEGRATOR_MEGAKERNEL,
diff --git a/intern/cycles/scene/devicescene.cpp b/intern/cycles/scene/devicescene.cpp
index 27aa20ba1d1..d3630014d15 100644
--- a/intern/cycles/scene/devicescene.cpp
+++ b/intern/cycles/scene/devicescene.cpp
@@ -56,7 +56,8 @@ DeviceScene::DeviceScene(Device *device)
       ies_lights(device, "ies", MEM_GLOBAL),
       volume_tree_nodes(device, "volume_tree_nodes", MEM_GLOBAL),
       volume_tree_roots(device, "volume_tree_roots", MEM_GLOBAL),
-      volume_tree_root_ids(device, "volume_tree_root_ids", MEM_GLOBAL)
+      volume_tree_root_ids(device, "volume_tree_root_ids", MEM_GLOBAL),
+      volume_step_size(device, "volume_step_size", MEM_GLOBAL)
 {
   memset((void *)&data, 0, sizeof(data));
 }
diff --git a/intern/cycles/scene/devicescene.h b/intern/cycles/scene/devicescene.h
index 9aaa7bbc22e..9dbd2ecba62 100644
--- a/intern/cycles/scene/devicescene.h
+++ b/intern/cycles/scene/devicescene.h
@@ -89,6 +89,7 @@ class DeviceScene {
   device_vector<KernelOctreeNode> volume_tree_nodes;
   device_vector<KernelOctreeRoot> volume_tree_roots;
   device_vector<int> volume_tree_root_ids;
+  device_vector<float> volume_step_size;
 
   KernelData data;
 
diff --git a/intern/cycles/scene/integrator.cpp b/intern/cycles/scene/integrator.cpp
index 2610b4dafe2..eef15d0757a 100644
--- a/intern/cycles/scene/integrator.cpp
+++ b/intern/cycles/scene/integrator.cpp
@@ -57,7 +57,9 @@ NODE_DEFINE(Integrator)
   SOCKET_FLOAT(ao_distance, "AO Distance", FLT_MAX);
   SOCKET_FLOAT(ao_additive_factor, "AO Additive Factor", 0.0f);
 
-  SOCKET_BOOLEAN(volume_unbiased, "Unbiased", false);
+  SOCKET_BOOLEAN(volume_ray_marching, "Biased", false);
+  SOCKET_INT(volume_max_steps, "Volume Max Steps", 1024);
+  SOCKET_FLOAT(volume_step_rate, "Volume Step Rate", 1.0f);
 
   static NodeEnum guiding_distribution_enum;
   guiding_distribution_enum.insert("PARALLAX_AWARE_VMM", GUIDING_TYPE_PARALLAX_AWARE_VMM);
@@ -226,7 +228,8 @@ void Integrator::device_update(Device *device, DeviceScene *dscene, Scene *scene
     }
   }
 
-  kintegrator->volume_unbiased = volume_unbiased;
+  kintegrator->volume_ray_marching = volume_ray_marching;
+  kintegrator->volume_max_steps = volume_max_steps;
 
   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 bb2ed42c792..0672250a53b 100644
--- a/intern/cycles/scene/integrator.h
+++ b/intern/cycles/scene/integrator.h
@@ -41,7 +41,9 @@ class Integrator : public Node {
   NODE_SOCKET_API(float, ao_distance)
   NODE_SOCKET_API(float, ao_additive_factor)
 
-  NODE_SOCKET_API(bool, volume_unbiased)
+  NODE_SOCKET_API(bool, volume_ray_marching)
+  NODE_SOCKET_API(int, volume_max_steps)
+  NODE_SOCKET_API(float, volume_step_rate)
 
   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 f96920a6220..5112270fb07 100644
--- a/intern/cycles/scene/object.cpp
+++ b/intern/cycles/scene/object.cpp
@@ -289,6 +289,97 @@ 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->is_light()) {
+    /* World volume. */
+    assert(static_cast<const Light *>(geometry)->get_light_type() == LIGHT_BACKGROUND);
+    for (const Node *node : geometry->get_used_shaders()) {
+      const Shader *shader = static_cast<const Shader *>(node);
+      if (shader->has_volume) {
+        return shader->get_volume_step_rate();
+      }
+    }
+    assert(false);
+    return FLT_MAX;
+  }
+
+  if (!geometry->is_mesh() && !geometry->is_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->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;
@@ -689,6 +780,7 @@ void ObjectManager::device_update(Device *device,
     dscene->object_motion_pass.tag_realloc();
     dscene->object_motion.tag_realloc();
     dscene->object_flag.tag_realloc();
+    dscene->volume_step_size.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. */
@@ -730,6 +822,7 @@ void ObjectManager::device_update(Device *device,
         dscene->object_motion_pass.tag_modified();
         dscene->object_motion.tag_modified();
         dscene->object_flag.tag_modified();
+        dscene->volume_step_size.tag_modified();
       }
 
       /* Update world object index. */
@@ -800,6 +893,9 @@ void ObjectManager::device_update_flags(Device * /*unused*/,
   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);
       }
@@ -852,7 +948,6 @@ void ObjectManager::device_update_flags(Device * /*unused*/,
 
   /* Copy object flag. */
   dscene->object_flag.copy_to_device();
-
   dscene->object_flag.clear_modified();
 }
 
diff --git a/intern/cycles/scene/object.h b/intern/cycles/scene/object.h
index 170fcef1e6b..026bc8dbcc2 100644
--- a/intern/cycles/scene/object.h
+++ b/intern/cycles/scene/object.h
@@ -112,6 +112,9 @@ 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 69c6b1f3dd3..5fc1b0d6f3f 100644
--- a/intern/cycles/scene/shader.cpp
+++ b/intern/cycles/scene/shader.cpp
@@ -72,6 +72,8 @@ 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);
@@ -101,6 +103,7 @@ 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();
@@ -390,9 +393,10 @@ void Shader::tag_update(Scene *scene)
     scene->procedural_manager->tag_update();
   }
 
-  if (has_volume != prev_has_volume) {
+  if (has_volume != prev_has_volume || volume_step_rate != prev_volume_step_rate) {
     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) {
diff --git a/intern/cycles/scene/shader.h b/intern/cycles/scene/shader.h
index 7b323c30119..f14ddaa4c00 100644
--- a/intern/cycles/scene/shader.h
+++ b/intern/cycles/scene/shader.h
@@ -83,10 +83,13 @@ class Shader : public Node {
   NODE_SOCKET_API(bool, use_bump_map_correction)
   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 8aa38bf1710..f7b602c235b 100644
--- a/intern/cycles/scene/volume.cpp
+++ b/intern/cycles/scene/volume.cpp
@@ -6,6 +6,7 @@
 #include "scene/attribute.h"
 #include "scene/background.h"
 #include "scene/image_vdb.h"
+#include "scene/integrator.h"
 #include "scene/light.h"
 #include "scene/object.h"
 #include "scene/scene.h"
@@ -34,6 +35,7 @@ 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,6 +44,7 @@ NODE_DEFINE(Volume)
 
 Volume::Volume() : Mesh(get_node_type(), Geometry::VOLUME)
 {
+  step_size = 0.0f;
   object_space = false;
 }
 
@@ -156,9 +159,12 @@ class VolumeMeshBuilder {
 
   void create_mesh(vector<float3> &vertices,
                    vector<int> &indices,
-                   const float face_overlap_avoidance);
+                   const float face_overlap_avoidance,
+                   const bool ray_marching);
 
-  void generate_vertices_and_quads(vector<int3> &vertices_is, vector<QuadData> &quads);
+  void generate_vertices_and_quads(vector<int3> &vertices_is,
+                                   vector<QuadData> &quads,
+                                   const bool ray_marching);
 
   void convert_object_space(const vector<int3> &vertices,
                             vector<float3> &out_vertices,
@@ -207,23 +213,96 @@ void VolumeMeshBuilder::add_padding(const int pad_size)
 
 void VolumeMeshBuilder::create_mesh(vector<float3> &vertices,
                                     vector<int> &indices,
-                                    const float face_overlap_avoidance)
+                                    const float face_overlap_avoidance,
+                                    const bool ray_marching)
 {
   /* We create vertices in index space (is), and only convert them to object
    * space when done. */
   vector<int3> vertices_is;
   vector<QuadData> quads;
 
-  generate_vertices_and_quads(vertices_is, quads);
+  generate_vertices_and_quads(vertices_is, quads, ray_marching);
 
   convert_object_space(vertices_is, vertices, face_overlap_avoidance);
 
   convert_quads_to_tris(quads, indices);
 }
 
-void VolumeMeshBuilder::generate_vertices_and_quads(vector<ccl::int3> &vertices_is,
-                                                    vector<QuadData> &quads)
+static bool is_non_empty_leaf(const openvdb::MaskGrid::TreeType &tree, const openvdb::Coord coord)
 {
+  const auto *leaf_node = tree.probeLeaf(coord);
+  return (leaf_node && !leaf_node->isEmpty());
+}
+
+void VolumeMeshBuilder::generate_vertices_and_quads(vector<ccl::int3> &vertices_is,
+                                                    vector<QuadData> &quads,
+                                                    const bool ray_marching)
+{
+  if (ray_marching) {
+    /* Make sure we only have leaf nodes in the tree, as tiles are not handled by this algorithm */
+    topology_grid->tree().voxelizeActiveTiles();
+
+    const openvdb::MaskGrid::TreeType &tree = topology_grid->tree();
+    tree.evalLeafBoundingBox(bbox);
+
+    const int3 resolution = make_int3(bbox.dim().x(), bbox.dim().y(), bbox.dim().z());
+
+    unordered_map<size_t, int> used_verts;
+    for (auto iter = tree.cbeginLeaf(); iter; ++iter) {
+      if (iter->isEmpty()) {
+        continue;
+      }
+      openvdb::CoordBBox leaf_bbox = iter->getNodeBoundingBox();
+      /* +1 to convert from exclusive to include bounds. */
+      leaf_bbox.max() = leaf_bbox.max().offsetBy(1);
+      int3 min = make_int3(leaf_bbox.min().x(), leaf_bbox.min().y(), leaf_bbox.min().z());
+      int3 max = make_int3(leaf_bbox.max().x(), leaf_bbox.max().y(), leaf_bbox.max().z());
+      int3 corners[8] = {
+          make_int3(min[0], min[1], min[2]),
+          make_int3(max[0], min[1], min[2]),
+          make_int3(max[0], max[1], min[2]),
+          make_int3(min[0], max[1], min[2]),
+          make_int3(min[0], min[1], max[2]),
+          make_int3(max[0], min[1], max[2]),
+          make_int3(max[0], max[1], max[2]),
+          make_int3(min[0], max[1], max[2]),
+      };
+      /* Only create a quad if on the border between an active and an inactive leaf.
+       *
+       * We verify that a leaf exists by probing a coordinate that is at its center,
+       * to do so we compute the center of the current leaf and offset this coordinate
+       * by the size of a leaf in each direction.
+       */
+      static const int LEAF_DIM = openvdb::MaskGrid::TreeType::LeafNodeType::DIM;
+      auto center = leaf_bbox.min() + openvdb::Coord(LEAF_DIM / 2);
+      if (!is_non_empty_leaf(tree, openvdb::Coord(center.x() - LEAF_DIM, center.y(), center.z())))
+      {
+        create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_X_MIN);
+      }
+      if (!is_non_empty_leaf(tree, openvdb::Coord(center.x() + LEAF_DIM, center.y(), center.z())))
+      {
+        create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_X_MAX);
+      }
+      if (!is_non_empty_leaf(tree, openvdb::Coord(center.x(), center.y() - LEAF_DIM, center.z())))
+      {
+        create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_Y_MIN);
+      }
+      if (!is_non_empty_leaf(tree, openvdb::Coord(center.x(), center.y() + LEAF_DIM, center.z())))
+      {
+        create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_Y_MAX);
+      }
+      if (!is_non_empty_leaf(tree, openvdb::Coord(center.x(), center.y(), center.z() - LEAF_DIM)))
+      {
+        create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_Z_MIN);
+      }
+      if (!is_non_empty_leaf(tree, openvdb::Coord(center.x(), center.y(), center.z() + LEAF_DIM)))
+      {
+        create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_Z_MAX);
+      }
+    }
+    return;
+  }
+
   bbox = topology_grid->evalActiveVoxelBoundingBox();
 
   const int3 resolution = make_int3(bbox.dim().x(), bbox.dim().y(), bbox.dim().z());
@@ -524,7 +603,8 @@ void GeometryManager::create_volume_mesh(const Scene *scene, Volume *volume, Pro
   /* Create mesh. */
   vector<float3> vertices;
   vector<int> indices;
-  builder.create_mesh(vertices, indices, face_overlap_avoidance);
+  const bool ray_marching = scene->integrator->get_volume_ray_marching();
+  builder.create_mesh(vertices, indices, face_overlap_avoidance, ray_marching);
 
   volume->reserve_mesh(vertices.size(), indices.size() / 3);
   volume->used_shaders.clear();
@@ -570,6 +650,11 @@ void VolumeManager::tag_update()
 /* Remove changed object from the list of octrees and tag for rebuild. */
 void VolumeManager::tag_update(const Object *object, uint32_t flag)
 {
+  if (object_octrees_.empty()) {
+    /* Volume object is not in the octree, can happen when using ray marching. */
+    return;
+  }
+
   if (flag & ObjectManager::VISIBILITY_MODIFIED) {
     tag_update();
   }
@@ -1013,12 +1098,52 @@ std::string VolumeManager::visualize_octree(const char *filename) const
   return filename_full;
 }
 
+void VolumeManager::update_step_size(const Scene *scene, DeviceScene *dscene) const
+{
+  assert(scene->integrator->get_volume_ray_marching());
+
+  if (!dscene->volume_step_size.is_modified() && last_algorithm == RAY_MARCHING) {
+    return;
+  }
+
+  if (dscene->volume_step_size.need_realloc()) {
+    dscene->volume_step_size.alloc(scene->objects.size());
+  }
+
+  float *volume_step_size = dscene->volume_step_size.data();
+
+  for (const Object *object : scene->objects) {
+    const Geometry *geom = object->get_geometry();
+    if (!geom->has_volume) {
+      continue;
+    }
+
+    volume_step_size[object->index] = scene->integrator->get_volume_step_rate() *
+                                      object->compute_volume_step_size();
+  }
+
+  dscene->volume_step_size.copy_to_device();
+  dscene->volume_step_size.clear_modified();
+}
+
 void VolumeManager::device_update(Device *device,
                                   DeviceScene *dscene,
                                   const Scene *scene,
                                   Progress &progress)
 {
-  if (need_rebuild_) {
+  if (scene->integrator->get_volume_ray_marching()) {
+    /* No need to update octree for ray marching. */
+    if (last_algorithm == NULL_SCATTERING) {
+      dscene->volume_tree_nodes.free();
+      dscene->volume_tree_roots.free();
+      dscene->volume_tree_root_ids.free();
+    }
+    update_step_size(scene, dscene);
+    last_algorithm = RAY_MARCHING;
+    return;
+  }
+
+  if (need_rebuild_ || last_algorithm == RAY_MARCHING) {
     /* Data needed for volume shader evaluation. */
     device->const_copy_to("data", &dscene->data, sizeof(dscene->data));
 
@@ -1039,6 +1164,11 @@ void VolumeManager::device_update(Device *device,
     LOG_DEBUG << "Octree visualization has been written to " << visualize_octree("octree.py");
     update_visualization_ = false;
   }
+
+  if (last_algorithm == RAY_MARCHING) {
+    dscene->volume_step_size.free();
+  }
+  last_algorithm = NULL_SCATTERING;
 }
 
 void VolumeManager::device_free(DeviceScene *dscene)
@@ -1046,6 +1176,7 @@ void VolumeManager::device_free(DeviceScene *dscene)
   dscene->volume_tree_nodes.free();
   dscene->volume_tree_roots.free();
   dscene->volume_tree_root_ids.free();
+  dscene->volume_step_size.free();
 }
 
 VolumeManager::~VolumeManager()
diff --git a/intern/cycles/scene/volume.h b/intern/cycles/scene/volume.h
index 46fdfbfe7fb..984ffa7d416 100644
--- a/intern/cycles/scene/volume.h
+++ b/intern/cycles/scene/volume.h
@@ -17,12 +17,19 @@ CCL_NAMESPACE_BEGIN
 class Object;
 class Octree;
 
+enum VolumeRenderingAlgorithm {
+  NULL_SCATTERING,
+  RAY_MARCHING,
+  NONE,
+};
+
 class Volume : public Mesh {
  public:
   NODE_DECLARE
 
   Volume();
 
+  NODE_SOCKET_API(float, step_size)
   NODE_SOCKET_API(bool, object_space)
   NODE_SOCKET_API(float, velocity_scale)
 
@@ -72,6 +79,9 @@ class VolumeManager {
    * `filename`, which is a Python script that can be run inside Blender. */
   std::string visualize_octree(const char *filename) const;
 
+  /* Step size for ray marching. */
+  void update_step_size(const Scene *, DeviceScene *) const;
+
   /* One octree per object per shader. */
   std::map<std::pair<const Object *, const Shader *>, std::shared_ptr<Octree>> object_octrees_;
 
@@ -81,6 +91,8 @@ class VolumeManager {
   int num_octree_nodes_;
   int num_octree_roots_;
 
+  VolumeRenderingAlgorithm last_algorithm = NONE;
+
 #ifdef WITH_OPENVDB
   /* Create SDF grid for mesh volumes, to determine whether a certain point is in the
    * interior of the mesh. This reduces evaluation time needed for heterogeneous volume. */
diff --git a/scripts/startup/bl_ui/properties_data_volume.py b/scripts/startup/bl_ui/properties_data_volume.py
index d6ddf6958e1..00f4af7792a 100644
--- a/scripts/startup/bl_ui/properties_data_volume.py
+++ b/scripts/startup/bl_ui/properties_data_volume.py
@@ -128,6 +128,8 @@ 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 33ea7c28e00..f9708564f59 100644
--- a/source/blender/makesrna/intern/rna_volume.cc
+++ b/source/blender/makesrna/intern/rna_volume.cc
@@ -541,6 +541,17 @@ 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);
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_and_volume.png b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_and_volume.png
new file mode 100644
index 00000000000..ea7d1369900
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_and_volume.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7572b66f9ca7f20a5ac1a4697d58eb73b7d26c1df0ed8dc7232de64c17f8642e
+size 31066
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_only.png b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_only.png
new file mode 100644
index 00000000000..ce95ba9dd26
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_camera_only.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:adbfb1fa2d3aa39ab56977cefa648c96c2ee1161c5975c7786296fb26c70634d
+size 32198
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_nested.png b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_nested.png
new file mode 100644
index 00000000000..6f63af7cca3
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_nested.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b19317985ecbeab0101313e38d897db1a0a7f5163b21ff97a564277cc12672f3
+size 34080
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_simple.png b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_simple.png
new file mode 100644
index 00000000000..3e2fd788c08
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/camera_in_volume_simple.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:dcc9ae25fb61bccc21a2cd18eda8011bff44f59bf886d01f0b16ff44f38d8225
+size 32407
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/implicit_volume.png b/tests/files/render/volume/cycles_ray_marching_renders/implicit_volume.png
new file mode 100644
index 00000000000..dd615aefbdf
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/implicit_volume.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:97eae519b8285147f3f9a799886e4973d9d432a50b45497502d023c714b8f7de
+size 45924
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/overlapping_different_anisotropy.png b/tests/files/render/volume/cycles_ray_marching_renders/overlapping_different_anisotropy.png
new file mode 100644
index 00000000000..d7db939d43b
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/overlapping_different_anisotropy.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:02441326ba412e6fa6b103cd18dd7429203fbdc70758ccde585844070c2ba4d1
+size 37744
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/overlapping_octrees.png b/tests/files/render/volume/cycles_ray_marching_renders/overlapping_octrees.png
new file mode 100644
index 00000000000..5a43614a408
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/overlapping_octrees.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c3308fcdf3b240a2fad27d08a7164fffc3c6e1a2e9616fb401af8ce406d7d15d
+size 25479
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/principled_absorption.png b/tests/files/render/volume/cycles_ray_marching_renders/principled_absorption.png
new file mode 100644
index 00000000000..6f3bf47367c
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/principled_absorption.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b0a51cbd9bd51945acd917b137c67f5a04eb4c49ebc34b2a9321ff673833563d
+size 29659
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/principled_bsdf_interior.png b/tests/files/render/volume/cycles_ray_marching_renders/principled_bsdf_interior.png
new file mode 100644
index 00000000000..3cce76d91a5
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/principled_bsdf_interior.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a786c1532d647a03fb809fd18649bca09ee6e3537a1792ef76672723d215247d
+size 29335
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_absobtion.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_absobtion.png
new file mode 100644
index 00000000000..189a252ede2
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_absobtion.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:1731b6e31e4e0210c38921dabb96e800563d24d209dcaebab73396fb59ec22be
+size 23589
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_deep_stack.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_deep_stack.png
new file mode 100644
index 00000000000..605b33f9836
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_deep_stack.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a0a0d4666a11ff6af51a1204220d2b0fdf7c13c0c4bb4c438ed57cc3af695070
+size 4141
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_edge_fireflies.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_edge_fireflies.png
new file mode 100644
index 00000000000..29f6858a82b
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_edge_fireflies.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:dcfb133dbc7e7240c2e1fd787186e5935d4cd7e93c00d7e4a027519ea5f5f713
+size 3298
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_instance.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_instance.png
new file mode 100644
index 00000000000..d03d6ddd2c9
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_instance.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7d993a35b298d5a21616387c84edb0c194c354e7cd58f0a71be44af8329c4ebe
+size 28354
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_light_path.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_light_path.png
new file mode 100644
index 00000000000..52587087999
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_light_path.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ab7c1a2289edb13d746b740a6c001ca482da7d38fb3a7d4b65c9abb269ecaeaf
+size 23088
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_output_absorption.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_absorption.png
new file mode 100644
index 00000000000..0e5c87f4869
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_absorption.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:cc15cddd5027cf590328a6d21a9f9dd4cfb52ead19fc2a7052e6928a90024fd7
+size 23903
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_output_mix_volume.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_mix_volume.png
new file mode 100644
index 00000000000..a38ab36b5b5
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_mix_volume.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a528037e96290f5306330d1fe656b11859339818e999ebc51b979feef2fbb9c9
+size 25603
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_output_scatter.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_scatter.png
new file mode 100644
index 00000000000..5ef5959c128
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_scatter.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:de2cf1e9d414b4eb18c6895c774356d411b33c857935b3a91e32769d9395bda6
+size 26116
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_output_surface.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_surface.png
new file mode 100644
index 00000000000..6d076157aa2
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_output_surface.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f3d2801e019ca3a2136124b5c62b6486227514e7e78a34e4d58a20a4f7b3de31
+size 25025
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_overlap.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_overlap.png
new file mode 100644
index 00000000000..b212ee7aefa
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_overlap.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:65ffd7efcabedd956c5dfb51eb3abb148804c3731181274148920e46e9d043aa
+size 25998
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter.png
new file mode 100644
index 00000000000..6c4adf9b9a7
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0ff643c9c84624544a2d05b2ed070adbab245ce9f00dc83adf718389c4cfce51
+size 26109
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_albedo.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_albedo.png
new file mode 100644
index 00000000000..d4c976f5694
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_albedo.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:11d937b2f5724e88dcac3555c63dfc8e28c51fc5545dc862d971a7c6c8814708
+size 32745
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_draine.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_draine.png
new file mode 100644
index 00000000000..dbfc1a4431e
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_draine.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:592e2d7c80b0ebc4db8985fcd22d18ab4123b841a799ba7d95679e4277e60fb8
+size 26681
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_fournier-forand.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_fournier-forand.png
new file mode 100644
index 00000000000..d57b83e832b
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_fournier-forand.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2dc3e38e17e7edb7719de4d2a5c5808a3eaceaebe98cb5a64cd96688960d9549
+size 36306
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie.png
new file mode 100644
index 00000000000..93d21eca68c
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e59f9a1176a5d278406ee1cab45ec04f791163743de859a85091a0c50e75ac7f
+size 33926
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie_small_particles.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie_small_particles.png
new file mode 100644
index 00000000000..500ef4b5576
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_mie_small_particles.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:78810cc672185f94c1cd83586d605a8ebb05bd52f96cb0043470cac469799a1d
+size 34232
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_rayleigh.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_rayleigh.png
new file mode 100644
index 00000000000..f332ced7af7
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_scatter_rayleigh.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fe4b1b20c4c812af686614898e90a8a40e9270152a8e86e14590e71e5f8bcefc
+size 43458
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_step_offset.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_step_offset.png
new file mode 100644
index 00000000000..28b2059fc5e
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_step_offset.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c9046550e6fb8ebcaa36d4ee4e5e8648a2a4ce36edda564e321652a086e0bde3
+size 9410
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_transparent_shadow.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_transparent_shadow.png
new file mode 100644
index 00000000000..36085dfe5c0
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_transparent_shadow.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3b49e7828249be2f02ba7033999919ac22d01534cdab1b6adcabf9f9832c0009
+size 32327
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_two_shaders.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_two_shaders.png
new file mode 100644
index 00000000000..29a43e25b08
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_two_shaders.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:86c1954c00369c5177a1b033bbfb71c1f601eef4b8e5bbd2ccf9f50a3bc76209
+size 27522
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/volume_zero_extinction_channel.png b/tests/files/render/volume/cycles_ray_marching_renders/volume_zero_extinction_channel.png
new file mode 100644
index 00000000000..f16cce2f009
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/volume_zero_extinction_channel.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:88247095729b93c07dac899809633d8333e5b8bf52523321e75047e90cbd8ad4
+size 42210
diff --git a/tests/files/render/volume/cycles_ray_marching_renders/world_volume.png b/tests/files/render/volume/cycles_ray_marching_renders/world_volume.png
new file mode 100644
index 00000000000..040c9fb6b8c
--- /dev/null
+++ b/tests/files/render/volume/cycles_ray_marching_renders/world_volume.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9ebb133a43bf16fc85199ee45d7dfc702a85e5b7d7f5255536340779ff9a8a6f
+size 31393
diff --git a/tests/python/cycles_render_tests.py b/tests/python/cycles_render_tests.py
index 1f8b69007c0..d5b8ef72d80 100644
--- a/tests/python/cycles_render_tests.py
+++ b/tests/python/cycles_render_tests.py
@@ -128,29 +128,39 @@ BLOCKLIST_GPU = [
 
 
 class CyclesReport(render_report.Report):
-    def __init__(self, title, output_dir, oiiotool, device=None, blocklist=[], osl=False):
+    def __init__(self, title, output_dir, oiiotool, device=None, blocklist=[], osl=False, ray_marching=False):
         # Split device name in format "<device_type>[-<RT>]" into individual
         # tokens, setting the RT suffix to an empty string if its not specified.
         self.device, suffix = (device.split("-") + [""])[:2]
         self.use_hwrt = (suffix == "RT")
         self.osl = osl
+        self.ray_marching = ray_marching
 
         variation = self.device
         if suffix:
             variation += ' ' + suffix
         if self.osl:
             variation += ' OSL'
+        if ray_marching:
+            variation += ' Ray Marching'
 
         super().__init__(title, output_dir, oiiotool, variation, blocklist)
 
+        self.set_pixelated(True)
+        self.set_reference_dir("cycles_renders")
+        if device == 'CPU':
+            self.set_compare_engine('eevee')
+        else:
+            self.set_compare_engine('cycles', 'CPU')
+
     def _get_render_arguments(self, arguments_cb, filepath, base_output_filepath):
-        return arguments_cb(filepath, base_output_filepath, self.use_hwrt, self.osl)
+        return arguments_cb(filepath, base_output_filepath, self.use_hwrt, self.osl, self.ray_marching)
 
     def _get_arguments_suffix(self):
         return ['--', '--cycles-device', self.device] if self.device else []
 
 
-def get_arguments(filepath, output_filepath, use_hwrt=False, osl=False):
+def get_arguments(filepath, output_filepath, use_hwrt=False, osl=False, ray_marching=False):
     dirname = os.path.dirname(filepath)
     basedir = os.path.dirname(dirname)
     subject = os.path.basename(dirname)
@@ -191,6 +201,9 @@ def get_arguments(filepath, output_filepath, use_hwrt=False, osl=False):
     if osl:
         args.extend(["--python-expr", "import bpy; bpy.context.scene.cycles.shading_system = True"])
 
+    if ray_marching:
+        args.extend(["--python-expr", "import bpy; bpy.context.scene.cycles.volume_biased = True"])
+
     if subject == 'bake':
         args.extend(['--python', os.path.join(basedir, "util", "render_bake.py")])
     elif subject == 'denoise_animation':
@@ -215,6 +228,13 @@ def create_argparse():
     return parser
 
 
+def test_volume_ray_marching(args, device, blocklist):
+    # Default volume rendering algorithm is null scattering, but we also want to test ray marching
+    report = CyclesReport('Cycles', args.outdir, args.oiiotool, device, blocklist, args.osl == 'all', ray_marching=True)
+    report.set_reference_dir("cycles_ray_marching_renders")
+    return report.run(args.testdir, args.blender, get_arguments, batch=args.batch)
+
+
 def main():
     parser = create_argparse()
     args = parser.parse_args()
@@ -243,12 +263,6 @@ def main():
         blocklist += BLOCKLIST_METAL
 
     report = CyclesReport('Cycles', args.outdir, args.oiiotool, device, blocklist, args.osl == 'all')
-    report.set_pixelated(True)
-    report.set_reference_dir("cycles_renders")
-    if device == 'CPU':
-        report.set_compare_engine('eevee')
-    else:
-        report.set_compare_engine('cycles', 'CPU')
 
     # Increase threshold for motion blur, see #78777.
     #
@@ -286,6 +300,9 @@ def main():
 
     ok = report.run(args.testdir, args.blender, get_arguments, batch=args.batch)
 
+    if (test_dir_name == 'volume'):
+        ok = ok and test_volume_ray_marching(args, device, blocklist)
+
     sys.exit(not ok)