Extremum: new extremum plugin family for coarse majorant and minorant superstructures.

docs/generate_plugin_doc.py

@@ -138,6 +138,11 @@
     'gridvolume'
 ]
 
+EXTREMUM_ORDERING = [
+    'extremum_global',
+    'extremum_grid'
+]
+
 
 def find_order_id(filename, ordering):
     f = os.path.split(filename)[-1].split('.')[0]
@@ -246,7 +251,8 @@ def generate(build_dir):
         ('samplers',    SAMPLER_ORDERING),
         ('films',       FILM_ORDERING),
         ('rfilters',    RFILTER_ORDERING),
-        ('volumes',     VOLUME_ORDERING)
+        ('volumes',     VOLUME_ORDERING),
+        ('extrema',     EXTREMUM_ORDERING)
     ]
 
     for section, ordering in sections:

docs/src/plugin_reference/section_extrema.rst

@@ -0,0 +1,11 @@
+.. _sec-extremum:
+
+Extremum Structures
+===================
+
+This section covers the different types of extremum structures included with
+Mitsuba. These plugins store local majorant/minorant bounds of a medium's
+extinction coefficient and are used by tracking-based integrators (e.g.
+:ref:`volpath <integrator-volpath>`) to perform delta/ratio tracking with
+locally-adaptive majorants instead of a single global majorant.
+

include/mitsuba/core/object.h

@@ -69,7 +69,10 @@ enum class ObjectType : uint32_t {
     PhaseFunction,
 
     /// A rendering algorithm aka. `Integrator`
-    Integrator
+    Integrator,
+
+    /// A medium acceleration structure.
+    Extremum
 };
 
 /**

include/mitsuba/render/extremum.h

@@ -0,0 +1,137 @@
+#pragma once
+
+#include <mitsuba/core/object.h>
+#include <mitsuba/render/interaction.h>
+#include <mitsuba/render/volume.h>
+#include <mitsuba/render/extremum_segment.h>
+#include <mitsuba/render/tracking.h>
+#include <drjit/call.h>
+
+#include <optional>
+
+NAMESPACE_BEGIN(mitsuba)
+
+/**
+ * \brief Abstract base class for extremum structures
+ *
+ * This class provides an interface for spatial data structures that store
+ * coarse volumetric local extrema (majorant/minorant). This enables efficient
+ * use of tracking algorithms with locally-adaptive majorants and minorants.
+ *
+ * The extremum structure needs to be built using the ``update_extremum``
+ * function, it is **not** called automatically in the constructor. It is the
+ * caller's responsability to pass the ``Volume`` plugin the extremum is
+ * derived from.
+ */
+template <typename Float, typename Spectrum>
+class MI_EXPORT_LIB Extremum : public JitObject<Extremum<Float, Spectrum>> {
+public:
+    MI_IMPORT_TYPES(Medium, Sampler, Volume)
+
+    using TrackingStateType    = TrackingState<Float, Spectrum>;
+    using TrackingFunctionType = TrackingFunction<Float, Spectrum>;
+
+    /// Destructor
+    ~Extremum();
+
+    /// Setter for the bbox over which the structure must be valid.
+    MI_INLINE void set_bbox(ScalarBoundingBox3f bbox) { m_bbox = bbox; };
+
+    /// Setter for the scale by which to multiply the extremum values.
+    MI_INLINE void set_scale(ScalarFloat scale) { m_scale = scale; }
+
+    /**
+     * \brief Update the bbox and scale, and rebuild the structure.
+     *
+     * The \c bbox parameters indicates the domain over which the extremum
+     * can be queried. It can be larger or smaller than the underlying
+     * volume bbox. It is the extremum's responsibility to be valid over this
+     * area. The building implementation is handled in ``build``.
+     *
+     * \param bbox      The validity bbox of the extremum structure
+     * \param volume    The volume from which to derive the extremum structure
+     * \param scale     The scale by which to multiply the extremum values
+     */
+    void update_extremum(const ScalarBoundingBox3f &bbox,
+                         const Volume *volume,
+                         std::optional<ScalarFloat> scale);
+
+    /**
+     * \brief Build the extremum structure of \c volume.
+     *
+     * Implements the logic that constructs the extremum structure from a
+     * \c volume. Called by ``update_extremum`` which is itself called by
+     * the owning ``Medium``
+     *
+     * \param volume  Volume to compute extremum values from
+     */
+    virtual void build(const Volume *volume) = 0;
+
+
+    /**
+     * \brief Traverse the extremum along a ray and applies a callback at each
+     * encountered segment.
+     *
+     * This method traverses the extremum structure segment by segment. At each
+     * segment, the callback ``func`` is called to advance the ``state``. This
+     * is useful for example to implement Delta Tracking, Ratio Tracking, and
+     * Residual Ratio Tracking. The callback is typically defined in the
+     * integrator.
+     *
+     * \param ray           Ray along which to sample
+     * \param mint          Minimum distance to consider
+     * \param maxt          Maximum distance to consider
+     * \param channel       Channel from which to sample
+     * \param state         Mutable tracking state carried through the traversal loop
+     * \param func          Callback function called at every segment.
+     * \param active        Mask for active lanes
+     *
+     * \return
+     *      The final tracking state, that includes the medium interaction if
+     *      a real scattering event was sampled, and the throughput and pdfs
+     *      accumulated throughout the traversal.
+     */
+    virtual TrackingStateType traverse_extremum(
+        const Ray3f &ray,
+        Float mint,
+        Float maxt,
+        UInt32 channel,
+        TrackingStateType state,
+        const TrackingFunctionType &func,
+        Mask active = true
+    ) const;
+
+    // =============================================================
+    //! @{ \name Non-virtual query methods
+    // =============================================================
+
+    ScalarBoundingBox3f bbox() const { return m_bbox; }
+    //! @}
+    // =============================================================
+
+    MI_DECLARE_PLUGIN_BASE_CLASS(Extremum)
+
+protected:
+    Extremum();
+    Extremum(const Properties &props);
+
+protected:
+    /// The bbox over which the extremum structure must be valid.
+    ScalarBoundingBox3f m_bbox;
+    /// Scale by which to multiply the extremum values.
+    ScalarFloat m_scale;
+};
+
+MI_EXTERN_CLASS(Extremum)
+NAMESPACE_END(mitsuba)
+
+// -----------------------------------------------------------------------
+//! @{ \name Enables vectorized method calls on Dr.Jit medium arrays
+// -----------------------------------------------------------------------
+
+DRJIT_CALL_TEMPLATE_BEGIN(mitsuba::Extremum)
+    DRJIT_CALL_METHOD(traverse_extremum)
+DRJIT_CALL_END()
+
+//! @}
+// -----------------------------------------------------------------------

include/mitsuba/render/extremum_segment.h

@@ -0,0 +1,98 @@
+#pragma once
+
+#include <mitsuba/core/fwd.h>
+
+NAMESPACE_BEGIN(mitsuba)
+
+/**
+ * \brief Stores the extremum (minorant/majorant) data for a ray segment.
+ *
+ * Used as the output type of Extremum traversal. Tracks the
+ * segment's entry/exit distances and the local extinction coefficient
+ * bounds within that interval.
+ */
+template<typename Float, typename Spectrum>
+struct ExtremumSegment {
+    MI_IMPORT_CORE_TYPES()                                                                         \
+
+    /// Segment entry distance along ray
+    Float mint;
+    /// Segment exit distance along ray
+    Float maxt;
+    /// Extremum data stored as [minorant, majorant]
+    Vector2f value;
+
+    /// Default constructor — creates an invalid segment via reset()
+    ExtremumSegment(){ reset(); };
+
+
+    /// Construct from entry/exit distances and a combined extremum vector.
+
+    ExtremumSegment(
+        Float mint,
+        Float maxt,
+        Vector2f value
+    ) : mint(mint),
+        maxt(maxt),
+        value(value) {}
+
+    /// Construct from entry/exit distances and separate minorant/majorant values.
+    ExtremumSegment(
+        const Float& mint,
+        const Float& maxt,
+        const Float& minorant,
+        const Float& majorant
+    ) : mint(mint),
+        maxt(maxt),
+        value(Vector2f(minorant, majorant)) {}
+
+    /**
+     * This callback method is invoked by dr::zeros<>, and takes care of fields
+     * that deviate from the standard zero-initialization convention. In
+     * ExtremumSegment, the ``mint`` and ``maxt`` fields are set to  + and -
+     * infinity respectively to to mark invalid intersection records.
+     */
+    void zero_(size_t size = 1) {
+        mint        = dr::full<Float>(dr::Infinity<Float>, size);
+        maxt        = dr::full<Float>(-dr::Infinity<Float>, size);
+        value       = dr::zeros<Vector2f>(size);
+    }
+
+    /**
+     * \brief Check whether this is a valid segment
+     *
+     * A segment is considered valid when
+     * \code
+     * segment.mint < segment.maxt
+     * \endcode
+     */
+    Mask valid() const {
+        return mint < maxt;
+    }
+
+    /**
+     * \brief Mark the extremum segment as invalid.
+     *
+     * This operation sets segment's minimum
+     * and maximum distances to \f$\infty\f$ and \f$-\infty\f$,
+     * respectively.
+     */
+    void reset() {
+        mint = dr::Infinity<Float>;
+        maxt = -dr::Infinity<Float>;
+    }
+
+    /// Minorant value over the segment. Accessor to the first element of ``value``.
+    MI_INLINE Float minorant() const {
+        return value.x();
+    }
+
+    /// Majorant value over the segment. Accessor to the second element of ``value``.
+    MI_INLINE Float majorant() const {
+        return value.y();
+    }
+
+    DRJIT_TRAVERSE(ExtremumSegment, mint, maxt, value)
+};
+
+NAMESPACE_END(mitsuba)

include/mitsuba/render/fwd.h

@@ -13,6 +13,7 @@ template <typename Float, typename Spectrum> class DirectedEdge;
 template <typename Float, typename Spectrum> class OptixDenoiser;
 template <typename Float, typename Spectrum> class Emitter;
 template <typename Float, typename Spectrum> class Endpoint;
+template <typename Float, typename Spectrum> class Extremum;
 template <typename Float, typename Spectrum> class Film;
 template <typename Float, typename Spectrum> class ImageBlock;
 template <typename Float, typename Spectrum> class Integrator;
@@ -41,6 +42,7 @@ template <typename Float, typename Spectrum> struct PositionSample;
 template <typename Float, typename Spectrum> struct BSDFSample3;
 template <typename Float, typename Spectrum> struct SilhouetteSample;
 template <typename Float, typename Spectrum> struct PhaseFunctionContext;
+template <typename Float, typename Spectrum> struct ExtremumSegment;
 template <typename Float, typename Spectrum> struct Interaction;
 template <typename Float, typename Spectrum> struct MediumInteraction;
 template <typename Float, typename Spectrum> struct SurfaceInteraction;
@@ -132,6 +134,8 @@ template <typename Float_, typename Spectrum_> struct RenderAliases {
     using Emitter                = mitsuba::Emitter<Float, Spectrum>;
     using Endpoint               = mitsuba::Endpoint<Float, Spectrum>;
     using Medium                 = mitsuba::Medium<Float, Spectrum>;
+    using Extremum               = mitsuba::Extremum<Float, Spectrum>;
+    using ExtremumSegment        = mitsuba::ExtremumSegment<Float, Spectrum>;
     using PhaseFunction          = mitsuba::PhaseFunction<Float, Spectrum>;
     using Film                   = mitsuba::Film<Float, Spectrum>;
     using ImageBlock             = mitsuba::ImageBlock<Float, Spectrum>;
@@ -145,6 +149,7 @@ template <typename Float_, typename Spectrum_> struct RenderAliases {
     using ObjectPtr              = dr::replace_scalar_t<Float, const Object *>;
     using BSDFPtr                = dr::replace_scalar_t<Float, const BSDF *>;
     using MediumPtr              = dr::replace_scalar_t<Float, const Medium *>;
+    using ExtremumPtr            = dr::replace_scalar_t<Float, const Extremum *>;
     using PhaseFunctionPtr       = dr::replace_scalar_t<Float, const PhaseFunction *>;
     using ShapePtr               = dr::replace_scalar_t<Float, const Shape *>;
     using MeshPtr                = dr::replace_scalar_t<Float, const Mesh *>;
@@ -196,6 +201,7 @@ template <typename Float_, typename Spectrum_> struct RenderAliases {
     using Interaction3f             = typename RenderAliases::Interaction3f;                       \
     using SurfaceInteraction3f      = typename RenderAliases::SurfaceInteraction3f;                \
     using MediumInteraction3f       = typename RenderAliases::MediumInteraction3f;                 \
+    using ExtremumSegment           = typename RenderAliases::ExtremumSegment;                     \
     using PreliminaryIntersection3f = typename RenderAliases::PreliminaryIntersection3f;           \
     using BSDFSample3f              = typename RenderAliases::BSDFSample3f;                        \
     using SilhouetteSample3f        = typename RenderAliases::SilhouetteSample3f;                  \
@@ -219,6 +225,7 @@ template <typename Float_, typename Spectrum_> struct RenderAliases {
     using Emitter                = typename RenderAliases::Emitter;                                \
     using Endpoint               = typename RenderAliases::Endpoint;                               \
     using Medium                 = typename RenderAliases::Medium;                                 \
+    using Extremum               = typename RenderAliases::Extremum;                               \
     using PhaseFunction          = typename RenderAliases::PhaseFunction;                          \
     using Film                   = typename RenderAliases::Film;                                   \
     using ImageBlock             = typename RenderAliases::ImageBlock;                             \
@@ -228,6 +235,7 @@ template <typename Float_, typename Spectrum_> struct RenderAliases {
     using ObjectPtr              = typename RenderAliases::ObjectPtr;                              \
     using BSDFPtr                = typename RenderAliases::BSDFPtr;                                \
     using MediumPtr              = typename RenderAliases::MediumPtr;                              \
+    using ExtremumPtr            = typename RenderAliases::ExtremumPtr;                            \
     using PhaseFunctionPtr       = typename RenderAliases::PhaseFunctionPtr;                       \
     using ShapePtr               = typename RenderAliases::ShapePtr;                               \
     using MeshPtr                = typename RenderAliases::MeshPtr;                                \

include/mitsuba/render/medium.h

@@ -11,7 +11,7 @@ NAMESPACE_BEGIN(mitsuba)
 template <typename Float, typename Spectrum>
 class MI_EXPORT_LIB Medium : public JitObject<Medium<Float, Spectrum>> {
 public:
-    MI_IMPORT_TYPES(PhaseFunction, Sampler, Scene, Texture);
+    MI_IMPORT_TYPES(PhaseFunction, Sampler, Scene, Texture, Extremum);
 
     /// Destructor
     ~Medium();
@@ -93,6 +93,26 @@ class MI_EXPORT_LIB Medium : public JitObject<Medium<Float, Spectrum>> {
         return m_has_spectral_extinction;
     }
 
+    /**
+     * \brief Intersects ray with the medium bbox and creates a medium interaction.
+     *
+     * \param ray   The ray that is used to test the medium bbox.
+     *
+     * \return
+     *      A tuple (mei, mint, maxt): ``mei`` is a  ``MediumInteraction3f``
+     *      object initialized with the current ray and medium data. ``mint``
+     *      and ``maxt`` represent the minimum and maximum intersection
+     *      distances of the ray with the medium's bbox. In case there are no
+     *      valid intersection, the range defaults to [0, +Inf].
+     */
+    std::tuple<MediumInteraction3f, Float, Float>
+    prepare_medium_traversal(const Ray3f &ray, Mask active) const;
+
+    /// Returns the extremum structure for local extremum acceleration.
+    MI_INLINE const Extremum *extremum() const {
+        return m_extremum.get();
+    }
+
     void traverse(TraversalCallback *callback) override;
 
     /// Return a human-readable representation of the Medium
@@ -109,8 +129,9 @@ class MI_EXPORT_LIB Medium : public JitObject<Medium<Float, Spectrum>> {
     bool m_sample_emitters;
     bool m_is_homogeneous;
     bool m_has_spectral_extinction;
+    ref<Extremum> m_extremum;
 
-    MI_DECLARE_TRAVERSE_CB(m_phase_function)
+    MI_DECLARE_TRAVERSE_CB(m_phase_function, m_extremum)
 };
 
 MI_EXTERN_CLASS(Medium)
@@ -130,6 +151,8 @@ DRJIT_CALL_TEMPLATE_BEGIN(mitsuba::Medium)
     DRJIT_CALL_METHOD(sample_interaction)
     DRJIT_CALL_METHOD(transmittance_eval_pdf)
     DRJIT_CALL_METHOD(get_scattering_coefficients)
+    DRJIT_CALL_GETTER(extremum)
+    DRJIT_CALL_METHOD(prepare_medium_traversal)
 DRJIT_CALL_END()
 
 // -----------------------------------------------------------------------