Filled out methods in volumelight.cpp and added appropriate tests for both shapes and volumelights.

Author: Microno95

include/mitsuba/python/docstr.h

@@ -4082,6 +4082,9 @@ 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_emitter =
+    R"doc(Return the emitter associated with the intersection (if any))doc";
+
 static const char *__doc_mitsuba_MediumInteraction_to_local = R"doc(Convert a world-space vector into local shading coordinates)doc";
 
 static const char *__doc_mitsuba_MediumInteraction_to_world = R"doc(Convert a local shading-space vector into world space)doc";

include/mitsuba/render/interaction.h

@@ -552,6 +552,10 @@ struct MediumInteraction : Interaction<Float_, Spectrum_> {
     //! @{ \name Methods
     // =============================================================
 
+    /// Get emitter attached to the medium associated with this interaction
+    /// \note Defined in scene.h
+    EmitterPtr emitter(Mask active = true) const;
+
 
     /// Convert a local shading-space vector into world space
     Vector3f to_world(const Vector3f &v) const {

include/mitsuba/render/medium.h

@@ -18,7 +18,7 @@ MI_DECLARE_ENUM_OPERATORS(MediumEventSamplingMode)
 template <typename Float, typename Spectrum>
 class MI_EXPORT_LIB Medium : public Object {
 public:
-    MI_IMPORT_TYPES(PhaseFunction, Sampler, Scene, Texture, Emitter, Volume);
+    MI_IMPORT_TYPES(PhaseFunction, Sampler, Scene, Texture, Emitter, Volume, EmitterPtr);
 
     /// Intersects a ray with the medium's bounding box
     virtual std::tuple<Mask, Float, Float>
@@ -123,7 +123,12 @@ class MI_EXPORT_LIB Medium : public Object {
     }
 
     /// Return the emitter of this medium
-    MI_INLINE const Emitter *emitter() const {
+    const EmitterPtr emitter(Mask /*unused*/ = true) const {
+        return m_emitter.get();
+    }
+
+    /// Return the emitter of this medium
+    EmitterPtr emitter(Mask /*unused*/ = true) {
         return m_emitter.get();
     }
 
@@ -185,7 +190,6 @@ DRJIT_VCALL_TEMPLATE_BEGIN(mitsuba::Medium)
     DRJIT_VCALL_GETTER(emitter, const typename Class::Emitter*)
     DRJIT_VCALL_GETTER(use_emitter_sampling, bool)
     DRJIT_VCALL_GETTER(is_homogeneous, bool)
-    DRJIT_VCALL_GETTER(is_emitter, bool)
     DRJIT_VCALL_GETTER(has_spectral_extinction, bool)
     DRJIT_VCALL_METHOD(get_majorant)
     DRJIT_VCALL_METHOD(get_radiance)
@@ -197,6 +201,7 @@ DRJIT_VCALL_TEMPLATE_BEGIN(mitsuba::Medium)
     DRJIT_VCALL_METHOD(medium_probabilities_analog)
     DRJIT_VCALL_METHOD(medium_probabilities_max)
     DRJIT_VCALL_METHOD(medium_probabilities_mean)
+    auto is_emitter() const { return neq(emitter(), nullptr); }
 DRJIT_VCALL_TEMPLATE_END(mitsuba::Medium)
 
 //! @}

include/mitsuba/render/mesh.h

@@ -334,6 +334,15 @@ class MI_EXPORT_LIB Mesh : public Shape<Float, Spectrum> {
      */
     void build_parameterization();
 
+    /**
+     * \brief Initialize the \c m_volume_parameterization field for running
+     * ray/object intersections.
+     *
+     * Internally, the function creates a nested scene to leverage optimized
+     * ray tracing functionality in \ref pdf_position_3d()
+     */
+    void build_volume_parameterization();
+
     // Ensures that the sampling table are ready.
     DRJIT_INLINE void ensure_pmf_built() const {
         if (unlikely(m_area_pmf.empty()))
@@ -487,7 +496,7 @@ class MI_EXPORT_LIB Mesh : public Shape<Float, Spectrum> {
     Float m_inv_volume;
 
     /// Optional: used in eval_parameterization()
-    ref<Scene<Float, Spectrum>> m_parameterization;
+    ref<Scene<Float, Spectrum>> m_parameterization, m_volume_parameterization;
 
     /// Pointer to the scene that owns this mesh
     Scene<Float, Spectrum>* m_scene = nullptr;

include/mitsuba/render/records.h

@@ -26,6 +26,7 @@ struct PositionSample {
     using Spectrum = Spectrum_;
     MI_IMPORT_RENDER_BASIC_TYPES()
     using SurfaceInteraction3f = typename RenderAliases::SurfaceInteraction3f;
+    using MediumInteraction3f  = typename RenderAliases::MediumInteraction3f;
 
     //! @}
     // =============================================================
@@ -77,6 +78,17 @@ struct PositionSample {
         : p(si.p), n(si.sh_frame.n), uv(si.uv), time(si.time), pdf(0.f),
           delta(false) { }
 
+    /**
+     * \brief Create a position sampling record from a medium interaction
+     *
+     * This is useful to determine the hypothetical sampling density in a
+     * medium after hitting it using standard ray tracing. This happens for
+     * instance in volumetric path tracing with multiple importance sampling.
+     */
+    PositionSample(const MediumInteraction3f &mei)
+        : p(mei.p), n(mei.sh_frame.n), uv(0.f), time(mei.time), pdf(0.f),
+          delta(false) { }
+
     /// Basic field constructor
     PositionSample(const Point3f &p, const Normal3f &n, const Point2f &uv,
                    Float time, Float pdf, Mask delta)
@@ -119,7 +131,9 @@ struct DirectionSample : public PositionSample<Float_, Spectrum_> {
 
     using Interaction3f        = typename RenderAliases::Interaction3f;
     using SurfaceInteraction3f = typename RenderAliases::SurfaceInteraction3f;
+    using MediumInteraction3f  = typename RenderAliases::MediumInteraction3f;
     using EmitterPtr           = typename RenderAliases::EmitterPtr;
+    using MediumPtr            = typename RenderAliases::MediumPtr;
 
     //! @}
     // =============================================================
@@ -179,6 +193,29 @@ struct DirectionSample : public PositionSample<Float_, Spectrum_> {
         emitter = si.emitter(scene);
     }
 
+    /**
+     * \brief Create a direct sampling record, which can be used to \a query
+     * the density of a medium position with respect to a given reference
+     * position.
+     *
+     * Direction `s` is set so that it points from the reference interaction to
+     * the interacted medium, as required when using e.g. the \ref Endpoint
+     * interface to compute PDF values.
+     *
+     * \param it
+     *     Medium interaction
+     *
+     * \param ref
+     *     Reference position
+     */
+    DirectionSample(const MediumInteraction3f &mei,
+                    const Interaction3f &ref) : Base(mei) {
+        Vector3f rel = mei.p - ref.p;
+        dist = dr::norm(rel);
+        d = select(mei.is_valid(), rel / dist, -mei.wi);
+        emitter = mei.emitter();
+    }
+
     /// Element-by-element constructor
     DirectionSample(const Point3f &p, const Normal3f &n, const Point2f &uv,
                     const Float &time, const Float &pdf, const Mask &delta,

include/mitsuba/render/scene.h

@@ -599,5 +599,18 @@ SurfaceInteraction<Float, Spectrum>::emitter(const Scene *scene, Mask active) co
     }
 }
 
+// See interaction.h
+template <typename Float, typename Spectrum>
+typename MediumInteraction<Float, Spectrum>::EmitterPtr
+MediumInteraction<Float, Spectrum>::emitter(Mask active) const {
+    if constexpr (!dr::is_jit_v<Float>) {
+        DRJIT_MARK_USED(active);
+        return is_valid() ? medium->emitter() : nullptr;
+    } else {
+        EmitterPtr emitter = medium->emitter(active & is_valid());
+        return emitter;
+    }
+}
+
 MI_EXTERN_CLASS(Scene)
 NAMESPACE_END(mitsuba)

include/mitsuba/render/shape.h

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

src/emitters/tests/test_volumelight.py

@@ -0,0 +1,132 @@
+import pytest
+import drjit as dr
+import mitsuba as mi
+
+from mitsuba.scalar_rgb.test.util import fresolver_append_path
+
+
+spectrum_dicts = {
+    'd65': {
+        "type": "d65",
+    },
+    'regular': {
+        "type": "regular",
+        "wavelength_min": 500,
+        "wavelength_max": 600,
+        "values": "1, 2"
+    }
+}
+
+
+@fresolver_append_path
+def create_emitter_and_spectrum(s_key='d65'):
+    emitter = mi.load_dict({
+        "type": "obj",
+        "filename": "resources/data/tests/obj/cbox_smallbox.obj",
+        "emitter" : { "type": "volumelight", "radiance" : spectrum_dicts[s_key] }
+    })
+    spectrum = mi.load_dict(spectrum_dicts[s_key])
+    expanded = spectrum.expand()
+    if len(expanded) == 1:
+        spectrum = expanded[0]
+
+    return emitter, spectrum
+
+
+def test01_constructor(variant_scalar_rgb):
+    # Check that the shape is properly bound to the emitter
+    shape, spectrum = create_emitter_and_spectrum()
+    assert shape.emitter().bbox() == shape.bbox()
+
+    # Check that we are not allowed to specify a to_world transform directly in the emitter.
+    with pytest.raises(RuntimeError):
+        e = mi.load_dict({
+            "type" : "volumelight",
+            "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()
+
+    it = dr.zeros(mi.SurfaceInteraction3f, 3)
+    assert dr.allclose(emitter.eval(it), spectrum.eval(it))
+
+    # Check that eval returns 0.0 when the sample point is outside the shape
+    it.p = mi.ScalarPoint3f([0.0, 0.0, 0.0])
+    assert dr.allclose(emitter.eval(it), 0.0)
+
+
+@pytest.mark.parametrize("spectrum_key", spectrum_dicts.keys())
+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()
+
+    time = 0.5
+    wavelength_sample = [0.5, 0.33, 0.1]
+    pos_sample = [[0.2, 0.1, 0.2], [0.6, 0.9, 0.2], [0.5]*3]
+    dir_sample = [[0.4, 0.5, 0.3], [0.1, 0.4, 0.9]]
+
+    # Sample a ray (position, direction, wavelengths) on the emitter
+    ray, res = emitter.sample_ray(time, wavelength_sample, pos_sample, dir_sample)
+
+    # Sample wavelengths on the spectrum
+    it = dr.zeros(mi.SurfaceInteraction3f, 3)
+
+    # Sample a position in the shape
+    ps = shape.sample_position_3d(time, pos_sample)
+    pdf = mi.warp.square_to_uniform_sphere_pdf(mi.warp.square_to_uniform_sphere(dir_sample))
+    it.p = ps.p
+    it.n = ps.n
+    it.time = time
+    it.wavelengths = mi.Float(wavelength_sample) * (mi.MI_CIE_MAX - mi.MI_CIE_MIN) + mi.MI_CIE_MIN
+
+    spec = dr.select(ps.pdf > 0.0, emitter.eval(it) * (mi.MI_CIE_MAX - mi.MI_CIE_MIN) / (ps.pdf * pdf), 0.0)
+    assert dr.allclose(res, spec)
+    assert dr.allclose(ray.time, time)
+    assert dr.allclose(ray.o, ps.p, atol=2e-2)
+    assert dr.allclose(ray.d, mi.Frame3f(ps.n).to_world(mi.warp.square_to_uniform_sphere(dir_sample)))
+
+
+@pytest.mark.parametrize("spectrum_key", spectrum_dicts.keys())
+def test04_sample_direction(variants_vec_spectral, spectrum_key):
+    # Check the correctness of the sample_direction(), pdf_direction(), and eval_direction() methods
+
+    shape, spectrum = create_emitter_and_spectrum(spectrum_key)
+    emitter = shape.emitter()
+
+    # Direction sampling is conditioned on a sampled position
+    it = dr.zeros(mi.SurfaceInteraction3f, 3)
+    it.p = [[0.2, 0.1, 0.2], [0.6, -0.9, 0.2],
+            [0.4, 0.9, -0.2]]  # Some positions
+    it.time = 1.0
+
+    # Sample direction on the emitter
+    samples = [[0.4, 0.5, 0.3], [0.1, 0.4, 0.9], [0.5]*3]
+    ds, res = emitter.sample_direction(it, samples)
+
+    # Sample direction on the shape
+    sphere_ds = dr.zeros(mi.DirectionSample3f)
+    sphere_ds.p = it.p
+    ps = shape.sample_position_3d(it.time, samples)
+    sphere_ds.d = ps.p - it.p
+    dist2 = dr.squared_norm(sphere_ds.d)
+    sphere_ds.d = sphere_ds.d / dr.sqrt(dist2)
+    sphere_ds.pdf = ps.pdf * dist2
+
+    assert dr.allclose(ds.pdf, sphere_ds.pdf)
+    assert dr.allclose(ds.pdf, emitter.pdf_direction(it, ds))
+    assert dr.allclose(ds.d, sphere_ds.d, atol=1e-3)
+    assert dr.allclose(ds.time, it.time)
+
+    # Evaluate the spectrum (divide by the pdf)
+    spec = dr.select(ds.pdf > 0.0, emitter.eval(it) / ds.pdf, 0.0)
+    assert dr.allclose(res, spec)
+
+    assert dr.allclose(emitter.eval_direction(it, ds), spec)