Get everything working with the latest Mitsuba master branch. Fully correct implementation of integrators, volpath, volpathmis and prbvolpath. Comparison with PBRT-V4 in a homogeneous medium bounded by the sphere shape with emission, absorption and scattering gives identical results.

TODO:
  * Test mesh sampling and directional PDFs are correct
  * Validate `heterogeneous` media results
  * Update reference images in `resources` submodule

Caveats:
  * PRBVolpath gives incorrect gradients (correct sign, incorrect magnitude) for the medium density `sigma_t` and the phase function `phase`, but this exists in the original implementation as well.

Author: Microno95

include/mitsuba/core/warp.h

@@ -282,19 +282,19 @@ MI_INLINE Vector<Value, 3> cube_to_uniform_sphere(const Point<Value, 3> &sample)
     Value radius = dr::safe_cbrt(sample.x()),
           phi = 2.f * dr::Pi<Value> * sample.z(),
           u = dr::fmsub(2.f, sample.y(), 1.f);
-    radius = dr::select(dr::neq(sample.x(), 0.f), radius, sample.x());
+    radius = dr::select(sample.x() == 0.f, sample.x(), radius);
     auto [s, c] = dr::sincos(phi);
-    return { radius * c * dr::safe_sqrt(1 - dr::sqr(u)), radius * s * dr::safe_sqrt(1 - dr::sqr(u)), radius * u };
+    return { radius * c * dr::safe_sqrt(1 - dr::square(u)), radius * s * dr::safe_sqrt(1 - dr::square(u)), radius * u };
 }
 
 /// Inverse of the mapping \ref square_to_uniform_sphere
 template <typename Value>
 MI_INLINE Point<Value, 3> uniform_sphere_to_cube(const Vector<Value, 3> &p) {
     Value phi = dr::atan2(p.y(), p.x()) * dr::InvTwoPi<Value>,
-          radius = dr::safe_sqrt(dr::sqr(p.x()) + dr::sqr(p.y()) + dr::sqr(p.z())),
+          radius = dr::safe_sqrt(dr::square(p.x()) + dr::square(p.y()) + dr::square(p.z())),
           u = (p.z() / radius + 1.f) / 2.f;
     return {
-        dr::sqr(radius) * radius,
+        dr::square(radius) * radius,
         u,
         dr::select(phi < 0.f, phi + 1.f, phi),
     };

include/mitsuba/python/docstr.h

@@ -2549,6 +2549,8 @@ static const char *__doc_mitsuba_Endpoint_m_needs_sample_2 = R"doc()doc";
 
 static const char *__doc_mitsuba_Endpoint_m_needs_sample_3 = R"doc()doc";
 
+static const char *__doc_mitsuba_Endpoint_m_needs_sample_2_3d = R"doc()doc";
+
 static const char *__doc_mitsuba_Endpoint_m_shape = R"doc()doc";
 
 static const char *__doc_mitsuba_Endpoint_m_to_world = R"doc()doc";
@@ -4259,6 +4261,8 @@ static const char *__doc_mitsuba_MediumInteraction_sigma_s = R"doc()doc";
 
 static const char *__doc_mitsuba_MediumInteraction_sigma_t = R"doc()doc";
 
+static const char *__doc_mitsuba_MediumInteraction_radiance = R"doc()doc";
+
 static const char *__doc_mitsuba_MediumInteraction_to_local =
 R"doc(Convert a world-space vector into local shading coordinates (defined
 by ``wi``))doc";
@@ -4291,6 +4295,18 @@ static const char *__doc_mitsuba_Medium_get_interaction_probabilities =
     R"doc(Returns the real scatter event probability, and the
 weights of a real and null scattering event at a given MediumInteraction mi)doc";
 
+static const char *__doc_mitsuba_Medium_medium_probabilities_analog =
+    R"doc(Computes the probabilities of interacting with a medium based on the point-wise
+density of each event)doc";
+
+static const char *__doc_mitsuba_Medium_medium_probabilities_max =
+    R"doc(Computes the probabilities of interacting with a medium based on the maximum
+density of each event multiplied by the throughput)doc";
+
+static const char *__doc_mitsuba_Medium_medium_probabilities_mean =
+    R"doc(Computes the probabilities of interacting with a medium based on the mean
+density of each event multiplied by the throughput)doc";
+
 static const char *__doc_mitsuba_Medium_get_radiance = R"doc(Returns the medium's radiance used for emissive media)doc";
 
 static const char *__doc_mitsuba_Medium_has_spectral_extinction = R"doc(Returns whether this medium has a spectrally varying extinction)doc";
@@ -4313,6 +4329,16 @@ static const char *__doc_mitsuba_Medium_m_phase_function = R"doc()doc";
 
 static const char *__doc_mitsuba_Medium_m_sample_emitters = R"doc()doc";
 
+static const char *__doc_mitsuba_Medium_m_medium_sampling_mode =
+R"doc(Determine how free-flight distances are sampled in the medium
+
+Takes on one of three values --- Analogue, Maximum and Mean --- which
+consider the probability of interacting with the medium based on the
+point-wise densities of each event, the maximum density of each event
+convolved with the throughput of the path, and the mean density of
+each event convolved with the throughput of the path.
+)doc";
+
 static const char *__doc_mitsuba_Medium_phase_function = R"doc(Return the phase function of this medium)doc";
 
 static const char *__doc_mitsuba_Medium_emitter = R"doc(Return the emitter associated with this medium)doc";

include/mitsuba/render/interaction.h

@@ -564,7 +564,7 @@ struct MediumInteraction : Interaction<Float_, Spectrum_> {
     /// Incident direction in world frame
     Vector3f wi;
 
-    UnpolarizedSpectrum sigma_s, sigma_n, sigma_t, combined_extinction;
+    UnpolarizedSpectrum sigma_s, sigma_n, sigma_t, radiance, combined_extinction;
 
     /// mint used when sampling the given distance ``t``
     Float mint;
@@ -587,6 +587,7 @@ struct MediumInteraction : Interaction<Float_, Spectrum_> {
         sigma_s             = dr::zeros<UnpolarizedSpectrum>(size);
         sigma_n             = dr::zeros<UnpolarizedSpectrum>(size);
         sigma_t             = dr::zeros<UnpolarizedSpectrum>(size);
+        radiance            = dr::zeros<UnpolarizedSpectrum>(size);
         combined_extinction = dr::zeros<UnpolarizedSpectrum>(size);
         mint                = dr::zeros<Float>(size);
 
@@ -615,7 +616,7 @@ struct MediumInteraction : Interaction<Float_, Spectrum_> {
     // =============================================================
 
     DRJIT_STRUCT(MediumInteraction, t, time, wavelengths, p, n, medium,
-                 sh_frame, wi, sigma_s, sigma_n, sigma_t,
+                 sh_frame, wi, sigma_s, sigma_n, sigma_t, radiance,
                  combined_extinction, mint)
 };
 

include/mitsuba/render/medium.h

@@ -41,15 +41,16 @@ class MI_EXPORT_LIB Medium : public Object {
 
     /// Returns the radiance, the probability of a scatter event, and
     /// the weights associated with real and null scattering events
-    std::pair<std::pair<UnpolarizedSpectrum, UnpolarizedSpectrum>,
+    virtual std::pair<std::pair<UnpolarizedSpectrum, UnpolarizedSpectrum>,
                std::pair<UnpolarizedSpectrum, UnpolarizedSpectrum>>
     get_interaction_probabilities(const Spectrum &radiance,
                                   const MediumInteraction3f &mi,
                                   const Spectrum &throughput) const;
 
     MI_INLINE std::pair<UnpolarizedSpectrum, UnpolarizedSpectrum>
     medium_probabilities_analog(const UnpolarizedSpectrum &radiance,
-                                const MediumInteraction3f &mi) const {
+                                const MediumInteraction3f &mi,
+                                const UnpolarizedSpectrum &/*throughput*/) const {
         auto prob_s = mi.sigma_t;
         auto prob_n = mi.sigma_n + dr::maximum(radiance, dr::mean(dr::abs(radiance)));
         return std::make_pair( prob_s, prob_n );
@@ -192,6 +193,7 @@ MI_CALL_TEMPLATE_BEGIN(Medium)
     DRJIT_CALL_GETTER(emitter)
     DRJIT_CALL_GETTER(use_emitter_sampling)
     DRJIT_CALL_GETTER(is_homogeneous)
+    DRJIT_CALL_GETTER(is_emitter)
     DRJIT_CALL_GETTER(has_spectral_extinction)
     DRJIT_CALL_METHOD(get_majorant)
     DRJIT_CALL_METHOD(get_radiance)
@@ -203,7 +205,6 @@ MI_CALL_TEMPLATE_BEGIN(Medium)
     DRJIT_CALL_METHOD(medium_probabilities_analog)
     DRJIT_CALL_METHOD(medium_probabilities_max)
     DRJIT_CALL_METHOD(medium_probabilities_mean)
-    auto is_emitter() const { return emitter() != nullptr; }
 MI_CALL_TEMPLATE_END(Medium)
 
 //! @}

include/mitsuba/render/scene.h

@@ -647,9 +647,11 @@ typename SurfaceInteraction<Float, Spectrum>::EmitterPtr
 SurfaceInteraction<Float, Spectrum>::emitter(const Scene *scene, Mask active) const {
     if constexpr (!dr::is_jit_v<Float>) {
         DRJIT_MARK_USED(active);
-        return is_valid() ? shape->emitter() : scene->environment();
+        EmitterPtr emitter = is_valid() ? shape->emitter() : scene->environment();
+        return (emitter != nullptr) && has_flag(emitter->flags(), EmitterFlags::Medium) ? nullptr : emitter;
     } else {
         EmitterPtr emitter = shape->emitter(active);
+        emitter = dr::select(active & !has_flag(emitter->flags(), EmitterFlags::Medium), emitter, nullptr);
         if (scene && scene->environment())
             emitter = dr::select(is_valid(), emitter, scene->environment() & active);
         return emitter;
@@ -662,9 +664,9 @@ typename MediumInteraction<Float, Spectrum>::EmitterPtr
 MediumInteraction<Float, Spectrum>::emitter(Mask active) const {
     if constexpr (!dr::is_jit_v<Float>) {
         DRJIT_MARK_USED(active);
-        return dr::neq(medium, nullptr) ? medium->emitter() : nullptr;
+        return medium != nullptr ? medium->emitter() : nullptr;
     } else {
-        EmitterPtr emitter = medium->emitter(active && dr::neq(medium, nullptr));
+        EmitterPtr emitter = medium->emitter(active && medium != nullptr);
         return emitter;
     }
 }

include/mitsuba/render/shape.h

@@ -221,7 +221,7 @@ struct SilhouetteSample : public PositionSample<Float_, Spectrum_> {
 template <typename Float, typename Spectrum>
 class MI_EXPORT_LIB Shape : public Object {
 public:
-    MI_IMPORT_TYPES(BSDF, Medium, Volume, Emitter, Sensor, MeshAttribute, Texture)
+    MI_IMPORT_TYPES(BSDF, Medium, Volume, Emitter, Sensor, MeshAttribute, Texture, EmitterPtr)
 
     // Use 32 bit indices to keep track of indices to conserve memory
     using ScalarIndex = uint32_t;

src/conftest.py

@@ -9,6 +9,7 @@
 import drjit as dr
 import mitsuba as mi
 
+
 re1 = re.compile(r'<built-in method (\w*) of PyCapsule object at 0x[0-9a-f]*>')
 re2 = re.compile(r'<bound method PyCapsule.(\w*)[^>]*>')
 

src/core/tests/test_logger.py

@@ -29,11 +29,12 @@ def append(self, level, text):
         logger.set_formatter(MyFormatter())
         logger.add_appender(MyAppender())
 
-        mi.Log(mi.LogLevel.Warn, "This is a test message")
-        assert len(messages) == 1
-        assert messages[0].startswith(
-                '300: class=None, thread=main, text=test01_custom(): This is a'
-                ' test message, filename=')
+        with mi.scoped_log_level(mi.LogLevel.Warn):
+            mi.Log(mi.LogLevel.Warn, "This is a test message")
+            assert len(messages) == 1
+            assert messages[0].startswith(
+                    '300: class=None, thread=main, text=test01_custom(): This is a'
+                    ' test message, filename=')
     finally:
         logger.clear_appenders()
         for app in appenders:

src/core/tests/test_warp.py

@@ -124,6 +124,8 @@ def test_cube_to_uniform_sphere_vec(variant_scalar_rgb):
     check_inverse(mi.warp.cube_to_uniform_sphere, mi.warp.uniform_sphere_to_cube, is_3d=True)
     check_warp_vectorization("cube_to_uniform_sphere", is_3d=True)
 
+    dr.set_flag(dr.JitFlag.Debug, False)
+
 
 def test_square_to_uniform_hemisphere(variant_scalar_rgb):
     assert(dr.allclose(mi.warp.square_to_uniform_hemisphere([0.5, 0.5]), [0, 0, 1]))

src/emitters/tests/test_volumelight.py

@@ -1,3 +1,5 @@
+import gc
+
 import pytest
 import drjit as dr
 import mitsuba as mi
@@ -23,6 +25,12 @@ def create_emitter_and_spectrum(s_key='d65'):
     emitter = mi.load_dict({
         "type": "obj",
         "filename": "resources/data/tests/obj/cbox_smallbox.obj",
+        "interior": {
+            "type": "homogeneous", "sigma_t": {
+                "type": "rgb",
+                "value": [1.0, 1.0, 1.0]
+            }
+        },
         "emitter" : { "type": "volumelight", "radiance" : spectrum_dicts[s_key] }
     })
     spectrum = mi.load_dict(spectrum_dicts[s_key])
@@ -43,7 +51,7 @@ def create_emitter_rgb():
         "emitter" : { "type": "volumelight", "radiance" : {
             "type": "rgb",
             "value": [1.0, 2.0, 1.0]
-        } }
+        }}
     })
 
     return r, c, emitter
@@ -58,14 +66,13 @@ def test01_constructor(variant_scalar_rgb):
     with pytest.raises(RuntimeError):
         e = mi.load_dict({
             "type" : "volumelight",
-            "to_world" : mi.ScalarTransform4f.translate([5, 0, 0])
+            "to_world" : mi.ScalarTransform4f().translate([5, 0, 0])
         })
 
 
 @pytest.mark.parametrize("spectrum_key", spectrum_dicts.keys())
 def test02_eval(variants_vec_spectral, spectrum_key):
     # Check that eval() return the same values as the 'radiance' spectrum
-
     shape, spectrum = create_emitter_and_spectrum(spectrum_key)
     emitter = shape.emitter()
 
@@ -81,7 +88,6 @@ def test02_eval(variants_vec_spectral, spectrum_key):
 def test03_sample_ray(variants_vec_spectral, spectrum_key):
     # Check the correctness of the sample_ray() method
 
-
     shape, spectrum = create_emitter_and_spectrum(spectrum_key)
     emitter = shape.emitter()
 
@@ -188,7 +194,7 @@ def test06_sample_direction_rgb(variants_vec_rgb):
 
     # Direction sampling is conditioned on a sampled position
     it = dr.zeros(mi.SurfaceInteraction3f, 3)
-    it.p = [[0.0, 0.48, 0.19, 10, -2], [0.0, 0.9018, 2.19, 2, 1.02], [0.0, 0.5, 3.291, 3.1, 9.1]]  # Some positions
+    it.p = mi.Point3f([[0.0, 0.48, 0.19, 10, -2], [0.0, 0.9018, 2.19, 2, 1.02], [0.0, 0.5, 3.291, 3.1, 9.1]])  # Some positions
     it.time = 1.0
 
     # Sample direction on the emitter
@@ -211,12 +217,12 @@ def test06_sample_direction_rgb(variants_vec_rgb):
     A = 1
     dist_from_c = dr.norm(it.p - center)
     B = 2*dr.dot(shape_ds.d, it.p - center + dist_from_c*shape_ds.d)
-    C = dr.squared_norm(it.p - center + dist_from_c*shape_ds.d) - dr.sqr(radius)
+    C = dr.squared_norm(it.p - center + dist_from_c*shape_ds.d) - dr.square(radius)
 
-    r0 = 0.5 * (-B - dr.sqrt(dr.sqr(B) - 4*A*C))+dist_from_c
-    r1 = 0.5 * (-B + dr.sqrt(dr.sqr(B) - 4*A*C))+dist_from_c
-    r0 = dr.clamp(r0, 0.0, dr.inf)
-    r1 = dr.clamp(r1, 0.0, dr.inf)
+    r0 = 0.5 * (-B - dr.sqrt(dr.square(B) - 4*A*C))+dist_from_c
+    r1 = 0.5 * (-B + dr.sqrt(dr.square(B) - 4*A*C))+dist_from_c
+    r0 = dr.clip(r0, 0.0, dr.inf)
+    r1 = dr.clip(r1, 0.0, dr.inf)
 
     analytic_pdf = dr.select(
         dr.isfinite(r0) & dr.isfinite(r1),
@@ -231,3 +237,17 @@ def test06_sample_direction_rgb(variants_vec_rgb):
     assert dr.allclose(res, spec)
 
     assert dr.allclose(emitter.eval_direction(it, ds) / ds.pdf, spec)
+
+
+def test07_shape_accessors(variants_vec_rgb):
+    shape, _ = create_emitter_and_spectrum()
+    shape_ptr = mi.ShapePtr(shape)
+
+    assert type(shape.emitter()) == mi.Emitter
+    assert type(shape_ptr.emitter()) == mi.EmitterPtr
+
+    emitter = shape.emitter()
+    emitter_ptr = mi.EmitterPtr(emitter)
+
+    assert type(emitter.get_shape()) == mi.Mesh
+    assert type(emitter_ptr.get_shape()) == mi.ShapePtr