diff --git a/src/integrators/volpath.cpp b/src/integrators/volpath.cpp
index 99e3232ea..37cc76c45 100644
--- a/src/integrators/volpath.cpp
+++ b/src/integrators/volpath.cpp
@@ -341,7 +341,8 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator<Float, Spectrum> {
             return { emitter_val, ds };
         }
 
-        Ray3f ray = ref_interaction.spawn_ray(ds.d);
+        Ray3f ray = ref_interaction.spawn_ray_to(ds.p);
+        Float max_dist = ray.maxt;
 
         // Potentially escaping the medium if this is the current medium's boundary
         if constexpr (std::is_convertible_v<Interaction, SurfaceInteraction3f>)
@@ -358,7 +359,7 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator<Float, Spectrum> {
         sampler->loop_put(loop);
         loop.init();
         while (loop(dr::detach(active))) {
-            Float remaining_dist = ds.dist * (1.f - math::ShadowEpsilon<Float>) - total_dist;
+            Float remaining_dist = max_dist - total_dist;
             ray.maxt = remaining_dist;
             active &= remaining_dist > 0.f;
             if (dr::none_or<false>(active))
diff --git a/src/integrators/volpathmis.cpp b/src/integrators/volpathmis.cpp
index 976d3ef06..4f6239f38 100644
--- a/src/integrators/volpathmis.cpp
+++ b/src/integrators/volpathmis.cpp
@@ -392,7 +392,8 @@ class VolpathMisIntegratorImpl final : public MonteCarloIntegrator<Float, Spectr
             return { p_over_f_nee, p_over_f_uni, emitter_val, ds};
         }
 
-        Ray3f ray = ref_interaction.spawn_ray(ds.d);
+        Ray3f ray = ref_interaction.spawn_ray_to(ds.p);
+        Float max_dist = ray.maxt;
 
         // Potentially escaping the medium if this is the current medium's boundary
         if constexpr (std::is_convertible_v<Interaction, SurfaceInteraction3f>)
@@ -409,8 +410,7 @@ class VolpathMisIntegratorImpl final : public MonteCarloIntegrator<Float, Spectr
         sampler->loop_put(loop);
         loop.init();
         while (loop(dr::detach(active))) {
-
-            Float remaining_dist = ds.dist * (1.f - math::ShadowEpsilon<Float>) - total_dist;
+            Float remaining_dist = max_dist - total_dist;
             ray.maxt = remaining_dist;
             active &= remaining_dist > 0.f;
             if (dr::none_or<false>(active))
diff --git a/src/python/python/ad/integrators/prbvolpath.py b/src/python/python/ad/integrators/prbvolpath.py
index ab391dfc7..73bf3550b 100644
--- a/src/python/python/ad/integrators/prbvolpath.py
+++ b/src/python/python/ad/integrators/prbvolpath.py
@@ -332,7 +332,8 @@ def sample_emitter(self, mei, si, active_medium, active_surface, scene, sampler,
         medium = dr.select(active, medium, dr.zeros(mi.MediumPtr))
         medium[(active_surface & si.is_medium_transition())] = si.target_medium(ds.d)
 
-        ray = ref_interaction.spawn_ray(ds.d)
+        ray = ref_interaction.spawn_ray_to(ds.p)
+        max_dist = mi.Float(ray.maxt)
         total_dist = mi.Float(0.0)
         si = dr.zeros(mi.SurfaceInteraction3f)
         needs_intersection = mi.Bool(True)
@@ -341,7 +342,7 @@ def sample_emitter(self, mei, si, active_medium, active_surface, scene, sampler,
                        state=lambda: (sampler, active, medium, ray, total_dist,
                                       needs_intersection, si, transmittance))
         while loop(active):
-            remaining_dist = ds.dist * (1.0 - mi.math.ShadowEpsilon) - total_dist
+            remaining_dist = max_dist - total_dist
             ray.maxt = dr.detach(remaining_dist)
             active &= remaining_dist > 0.0