From f3ccf0ac19c9b977ff931443cebe199240ad13b4 Mon Sep 17 00:00:00 2001 From: nicolae marton Date: Wed, 9 Sep 2026 14:22:33 +0200 Subject: [PATCH] New Extremum plugin: Coarse structures of majorant and minorant values over volumetric data. At runtime can be traversed to perform tracking algorithms (e.g. delta tracking or ratio tracking) using the `traverse_extremum` function. Built from sigma_t volumes in media, but intentionally decoupled to allow for different volumetric representations. --- docs/generate_plugin_doc.py | 8 +- docs/src/plugin_reference/section_extrema.rst | 11 + include/mitsuba/core/object.h | 5 +- include/mitsuba/python/docstr.h | 470 +++++++++++++++--- include/mitsuba/render/extremum.h | 137 +++++ include/mitsuba/render/extremum_segment.h | 98 ++++ include/mitsuba/render/fwd.h | 8 + include/mitsuba/render/medium.h | 27 +- include/mitsuba/render/tracking.h | 87 ++++ include/mitsuba/render/volume.h | 40 +- include/mitsuba/render/volume_utils.h | 28 ++ src/CMakeLists.txt | 1 + src/core/parser.cpp | 3 +- src/core/plugin.cpp | 1 + src/extrema/CMakeLists.txt | 6 + src/extrema/extremum_global.cpp | 92 ++++ src/extrema/extremum_grid.cpp | 336 +++++++++++++ src/extrema/tests/test_extremum_global.py | 112 +++++ src/extrema/tests/test_extremum_grid.py | 184 +++++++ src/integrators/volpath.cpp | 351 ++++++++----- src/media/heterogeneous.cpp | 36 +- src/media/homogeneous.cpp | 27 +- src/python/main_v.cpp | 8 + src/render/CMakeLists.txt | 4 + src/render/extremum.cpp | 48 ++ src/render/medium.cpp | 28 ++ src/render/python/CMakeLists.txt | 2 + src/render/python/extremum_v.cpp | 166 +++++++ src/render/python/medium_v.cpp | 9 +- src/render/python/tracking_v.cpp | 21 + src/render/volume.cpp | 6 + src/volumes/grid.cpp | 119 +++++ 32 files changed, 2279 insertions(+), 200 deletions(-) create mode 100644 docs/src/plugin_reference/section_extrema.rst create mode 100644 include/mitsuba/render/extremum.h create mode 100644 include/mitsuba/render/extremum_segment.h create mode 100644 include/mitsuba/render/tracking.h create mode 100644 include/mitsuba/render/volume_utils.h create mode 100644 src/extrema/CMakeLists.txt create mode 100644 src/extrema/extremum_global.cpp create mode 100644 src/extrema/extremum_grid.cpp create mode 100644 src/extrema/tests/test_extremum_global.py create mode 100644 src/extrema/tests/test_extremum_grid.py create mode 100644 src/render/extremum.cpp create mode 100644 src/render/python/extremum_v.cpp create mode 100644 src/render/python/tracking_v.cpp diff --git a/docs/generate_plugin_doc.py b/docs/generate_plugin_doc.py index 76bc6868ec..2e98547410 100644 --- a/docs/generate_plugin_doc.py +++ b/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: diff --git a/docs/src/plugin_reference/section_extrema.rst b/docs/src/plugin_reference/section_extrema.rst new file mode 100644 index 0000000000..4b563b4595 --- /dev/null +++ b/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 `) to perform delta/ratio tracking with +locally-adaptive majorants instead of a single global majorant. + diff --git a/include/mitsuba/core/object.h b/include/mitsuba/core/object.h index c9c68eec17..f59cf8a6a8 100644 --- a/include/mitsuba/core/object.h +++ b/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 }; /** diff --git a/include/mitsuba/python/docstr.h b/include/mitsuba/python/docstr.h index fb35be2037..84025a0aba 100644 --- a/include/mitsuba/python/docstr.h +++ b/include/mitsuba/python/docstr.h @@ -483,6 +483,50 @@ static const char *__doc_OptixProgramGroupSingleModule_entryFunctionName = R"doc static const char *__doc_OptixProgramGroupSingleModule_module = R"doc()doc"; +static const char *__doc_OptixSRTData = R"doc()doc"; + +static const char *__doc_OptixSRTData_a = R"doc()doc"; + +static const char *__doc_OptixSRTData_b = R"doc()doc"; + +static const char *__doc_OptixSRTData_c = R"doc()doc"; + +static const char *__doc_OptixSRTData_pvx = R"doc()doc"; + +static const char *__doc_OptixSRTData_pvy = R"doc()doc"; + +static const char *__doc_OptixSRTData_pvz = R"doc()doc"; + +static const char *__doc_OptixSRTData_qw = R"doc()doc"; + +static const char *__doc_OptixSRTData_qx = R"doc()doc"; + +static const char *__doc_OptixSRTData_qy = R"doc()doc"; + +static const char *__doc_OptixSRTData_qz = R"doc()doc"; + +static const char *__doc_OptixSRTData_sx = R"doc()doc"; + +static const char *__doc_OptixSRTData_sy = R"doc()doc"; + +static const char *__doc_OptixSRTData_sz = R"doc()doc"; + +static const char *__doc_OptixSRTData_tx = R"doc()doc"; + +static const char *__doc_OptixSRTData_ty = R"doc()doc"; + +static const char *__doc_OptixSRTData_tz = R"doc()doc"; + +static const char *__doc_OptixSRTMotionTransform = R"doc()doc"; + +static const char *__doc_OptixSRTMotionTransform_child = R"doc()doc"; + +static const char *__doc_OptixSRTMotionTransform_motionOptions = R"doc()doc"; + +static const char *__doc_OptixSRTMotionTransform_pad = R"doc()doc"; + +static const char *__doc_OptixSRTMotionTransform_srtData = R"doc()doc"; + static const char *__doc_OptixShaderBindingTable = R"doc()doc"; static const char *__doc_OptixShaderBindingTable_callablesRecordBase = R"doc()doc"; @@ -2213,6 +2257,14 @@ static const char *__doc_mitsuba_BoundingSphere_empty = R"doc(Return whether thi static const char *__doc_mitsuba_BoundingSphere_expand = R"doc(Expand the bounding sphere radius to contain another point)doc"; +static const char *__doc_mitsuba_BoundingSphere_fields = R"doc()doc"; + +static const char *__doc_mitsuba_BoundingSphere_fields_2 = R"doc()doc"; + +static const char *__doc_mitsuba_BoundingSphere_labels = R"doc()doc"; + +static const char *__doc_mitsuba_BoundingSphere_name = R"doc()doc"; + static const char *__doc_mitsuba_BoundingSphere_operator_eq = R"doc(Equality test against another bounding sphere)doc"; static const char *__doc_mitsuba_BoundingSphere_operator_ne = R"doc(Inequality test against another bounding sphere)doc"; @@ -2221,8 +2273,6 @@ static const char *__doc_mitsuba_BoundingSphere_radius = R"doc()doc"; static const char *__doc_mitsuba_BoundingSphere_ray_intersect = R"doc(Check if a ray intersects the bounding sphere)doc"; -static const char *__doc_mitsuba_BoundingSphere_traverse_cb = R"doc()doc"; - static const char *__doc_mitsuba_Color = R"doc()doc"; static const char *__doc_mitsuba_Color_Color = R"doc()doc"; @@ -3273,7 +3323,15 @@ static const char *__doc_mitsuba_DiscreteDistribution2D_DiscreteDistribution2D = R"doc(Construct a marginal sample warping scheme for floating point data of resolution ``size``.)doc"; -static const char *__doc_mitsuba_DiscreteDistribution2D_DiscreteDistribution2D_2 = R"doc()doc"; +static const char *__doc_mitsuba_DiscreteDistribution2D_DiscreteDistribution2D_2 = +R"doc(Construct a marginal sample warping scheme from a row-major array of +``prod(size)`` values that lives on the device + +The tables are computed with prefix sums on the device, which avoids a +synchronization with the host. On scalar backends this delegates to +the host constructor.)doc"; + +static const char *__doc_mitsuba_DiscreteDistribution2D_DiscreteDistribution2D_3 = R"doc()doc"; static const char *__doc_mitsuba_DiscreteDistribution2D_eval = R"doc(Evaluate the function value at the given integer position)doc"; @@ -3581,7 +3639,7 @@ static const char *__doc_mitsuba_EmitterFlags_Surface = R"doc(The emitter is att static const char *__doc_mitsuba_Emitter_Emitter = R"doc(This is both a class and the base of various Mitsuba plugins)doc"; -static const char *__doc_mitsuba_Emitter_class_name = R"doc(This is both a class and the base of various Mitsuba plugins)doc"; +static const char *__doc_mitsuba_Emitter_class_name = R"doc()doc"; static const char *__doc_mitsuba_Emitter_dirty = R"doc(Return whether the emitter parameters have changed)doc"; @@ -3616,9 +3674,9 @@ static const char *__doc_mitsuba_Emitter_traverse = R"doc()doc"; static const char *__doc_mitsuba_Emitter_traverse_cb = R"doc()doc"; -static const char *__doc_mitsuba_Emitter_type = R"doc(This is both a class and the base of various Mitsuba plugins)doc"; +static const char *__doc_mitsuba_Emitter_type = R"doc()doc"; -static const char *__doc_mitsuba_Emitter_variant_name = R"doc(This is both a class and the base of various Mitsuba plugins)doc"; +static const char *__doc_mitsuba_Emitter_variant_name = R"doc()doc"; static const char *__doc_mitsuba_Emitter_visibility_mask = R"doc(Return the 8-bit visibility mask (see `RayMask`). Invisible emitters @@ -3930,12 +3988,154 @@ static const char *__doc_mitsuba_Endpoint_traverse_cb_fields = R"doc()doc"; static const char *__doc_mitsuba_Endpoint_world_transform = R"doc(Return the local space to world space transformation)doc"; static const char *__doc_mitsuba_Endpoint_world_transform_scalar = -R"doc(Return the local space to world space transformation as a scalar -transform +R"doc(\brief Return the local space to world space transformation as a +scalar transform This mirrors `world_transform()` without touching the device: reading the JIT representation back would wait for all queued work.)doc"; +static const char *__doc_mitsuba_Extremum = +R"doc(\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.)doc"; + +static const char *__doc_mitsuba_Extremum_2 = R"doc()doc"; + +static const char *__doc_mitsuba_Extremum_3 = R"doc()doc"; + +static const char *__doc_mitsuba_Extremum_4 = R"doc()doc"; + +static const char *__doc_mitsuba_Extremum_5 = R"doc()doc"; + +static const char *__doc_mitsuba_Extremum_6 = R"doc()doc"; + +static const char *__doc_mitsuba_Extremum_7 = R"doc()doc"; + +static const char *__doc_mitsuba_ExtremumSegment = +R"doc(\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.)doc"; + +static const char *__doc_mitsuba_ExtremumSegment_2 = R"doc()doc"; + +static const char *__doc_mitsuba_ExtremumSegment_ExtremumSegment = R"doc(Default constructor — creates an invalid segment via reset())doc"; + +static const char *__doc_mitsuba_ExtremumSegment_ExtremumSegment_2 = R"doc(Construct from entry/exit distances and a combined extremum vector.)doc"; + +static const char *__doc_mitsuba_ExtremumSegment_ExtremumSegment_3 = R"doc(Construct from entry/exit distances and separate minorant/majorant values.)doc"; + +static const char *__doc_mitsuba_ExtremumSegment_fields = R"doc()doc"; + +static const char *__doc_mitsuba_ExtremumSegment_fields_2 = R"doc()doc"; + +static const char *__doc_mitsuba_ExtremumSegment_labels = R"doc()doc"; + +static const char *__doc_mitsuba_ExtremumSegment_majorant = R"doc(Majorant value over the segment. Accessor to the second element of ``value``.)doc"; + +static const char *__doc_mitsuba_ExtremumSegment_maxt = R"doc(Segment exit distance along ray)doc"; + +static const char *__doc_mitsuba_ExtremumSegment_minorant = R"doc(Minorant value over the segment. Accessor to the first element of ``value``.)doc"; + +static const char *__doc_mitsuba_ExtremumSegment_mint = R"doc(Segment entry distance along ray)doc"; + +static const char *__doc_mitsuba_ExtremumSegment_name = R"doc()doc"; + +static const char *__doc_mitsuba_ExtremumSegment_reset = +R"doc(\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.)doc"; + +static const char *__doc_mitsuba_ExtremumSegment_valid = +R"doc(\brief Check whether this is a valid segment + +A segment is considered valid when +\code +segment.mint < segment.maxt +\endcode)doc"; + +static const char *__doc_mitsuba_ExtremumSegment_value = R"doc(Extremum data stored as [minorant, majorant])doc"; + +static const char *__doc_mitsuba_ExtremumSegment_zero = +R"doc(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.)doc"; + +static const char *__doc_mitsuba_Extremum_Extremum = R"doc(@})doc"; + +static const char *__doc_mitsuba_Extremum_Extremum_2 = R"doc()doc"; + +static const char *__doc_mitsuba_Extremum_bbox = R"doc(@{ \name Non-virtual query methods)doc"; + +static const char *__doc_mitsuba_Extremum_build = +R"doc(\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)doc"; + +static const char *__doc_mitsuba_Extremum_class_name = R"doc()doc"; + +static const char *__doc_mitsuba_Extremum_m_bbox = R"doc(The bbox over which the extremum structure must be valid.)doc"; + +static const char *__doc_mitsuba_Extremum_m_scale = R"doc(Scale by which to multiply the extremum values.)doc"; + +static const char *__doc_mitsuba_Extremum_set_bbox = R"doc(Setter for the bbox over which the structure must be valid.)doc"; + +static const char *__doc_mitsuba_Extremum_set_scale = R"doc(Setter for the scale by which to multiply the extremum values.)doc"; + +static const char *__doc_mitsuba_Extremum_traverse_extremum = +R"doc(\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.)doc"; + +static const char *__doc_mitsuba_Extremum_type = R"doc()doc"; + +static const char *__doc_mitsuba_Extremum_update_extremum = +R"doc(\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)doc"; + +static const char *__doc_mitsuba_Extremum_variant_name = R"doc()doc"; + static const char *__doc_mitsuba_FileResolver = R"doc(Simple class for resolving paths on Linux/Windows/Mac OS @@ -5495,6 +5695,8 @@ static const char *__doc_mitsuba_Medium_Medium_2 = R"doc()doc"; static const char *__doc_mitsuba_Medium_class_name = R"doc()doc"; +static const char *__doc_mitsuba_Medium_extremum = R"doc(Returns the extremum structure for local extremum acceleration.)doc"; + static const char *__doc_mitsuba_Medium_get_majorant = R"doc(Returns the medium's majorant used for delta tracking)doc"; static const char *__doc_mitsuba_Medium_get_scattering_coefficients = @@ -5507,6 +5709,8 @@ static const char *__doc_mitsuba_Medium_intersect_aabb = R"doc(Intersects a ray static const char *__doc_mitsuba_Medium_is_homogeneous = R"doc(Returns whether this medium is homogeneous)doc"; +static const char *__doc_mitsuba_Medium_m_extremum = R"doc()doc"; + static const char *__doc_mitsuba_Medium_m_has_spectral_extinction = R"doc()doc"; static const char *__doc_mitsuba_Medium_m_is_homogeneous = R"doc()doc"; @@ -5517,6 +5721,18 @@ static const char *__doc_mitsuba_Medium_m_sample_emitters = R"doc()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_prepare_medium_traversal = +R"doc(\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].)doc"; + static const char *__doc_mitsuba_Medium_sample_interaction = R"doc(Sample a free-flight distance in the medium. @@ -5966,13 +6182,11 @@ static const char *__doc_mitsuba_Mesh_Part_face_offset = R"doc(Index of the firs static const char *__doc_mitsuba_Mesh_Part_id = R"doc(Label of the part, e.g. the id of the source mesh)doc"; static const char *__doc_mitsuba_Mesh_Part_label = -R"doc(Name that the source mesh carries itself +R"doc(Name assigned by the mes source -This names the mesh independently of the scene it appears in: the label -that a :monosp:`serialized` file stores next to the geometry, or the -file name for formats that store no name of their own. The ``id`` field -instead holds the name that the Mitsuba scene knows the shape by, and -the two differ whenever the scene assigns an id.)doc"; +This field records names that a :monosp:`serialized` file stores +along with the geometry (e.g., the mesh name exported from Blender), +or the file name for formats that store no name of their own.)doc"; static const char *__doc_mitsuba_Mesh_add_attribute = R"doc(Add the mesh attribute ``name`` @@ -6245,7 +6459,9 @@ static const char *__doc_mitsuba_Mesh_is_vertex_attribute = R"doc(Does the attri static const char *__doc_mitsuba_Mesh_m_area_pmf = R"doc()doc"; -static const char *__doc_mitsuba_Mesh_m_bbox = R"doc(Bounding box of the mesh positions)doc"; +static const char *__doc_mitsuba_Mesh_m_bbox = R"doc(Bounding box of the mesh positions, computed on demand by `bbox()`)doc"; + +static const char *__doc_mitsuba_Mesh_m_bbox_valid = R"doc(Does `m_bbox` reflect the current positions?)doc"; static const char *__doc_mitsuba_Mesh_m_bsdf_index = R"doc()doc"; @@ -6383,7 +6599,10 @@ packed buffers in a single pass, and ends in ``refresh()``. The ``flip_normals`` flag turns the surface inside out as the records are written, which `from_fields()` uses to bake the property of -the same name. See ``validate_impl()`` for ``updating``.)doc"; +the same name. See ``validate_impl()`` for ``updating``. + +The caller can provide a bounding box (if known), in which case +the implementation does not need to recompute it.)doc"; static const char *__doc_mitsuba_Mesh_packed_face = R"doc(Returns the packed face record of triangle ``index``)doc"; @@ -6450,7 +6669,7 @@ static const char *__doc_mitsuba_Mesh_ray_intersect_triangle_packet_3 = R"doc()d static const char *__doc_mitsuba_Mesh_ray_intersect_triangle_scalar = R"doc()doc"; -static const char *__doc_mitsuba_Mesh_recompute_bbox = R"doc((Re-)compute the bounding box from the packed positions)doc"; +static const char *__doc_mitsuba_Mesh_recompute_bbox = R"doc((Re-)compute the bounding box from the packed positions.)doc"; static const char *__doc_mitsuba_Mesh_recompute_normals = R"doc((Re-) compute smooth interpolated normals from the positions)doc"; @@ -6668,7 +6887,9 @@ static const char *__doc_mitsuba_MicrofacetDistribution_2 = R"doc()doc"; static const char *__doc_mitsuba_MicrofacetDistribution_G = R"doc(Smith's separable shadowing-masking approximation)doc"; -static const char *__doc_mitsuba_MicrofacetDistribution_MicrofacetDistribution = +static const char *__doc_mitsuba_MicrofacetDistribution_MicrofacetDistribution = R"doc(Create an uninitialized distribution, to be assigned later)doc"; + +static const char *__doc_mitsuba_MicrofacetDistribution_MicrofacetDistribution_2 = R"doc(Create an isotropic microfacet distribution of the specified type Args: @@ -6676,7 +6897,7 @@ R"doc(Create an isotropic microfacet distribution of the specified type alpha: The surface roughness)doc"; -static const char *__doc_mitsuba_MicrofacetDistribution_MicrofacetDistribution_2 = +static const char *__doc_mitsuba_MicrofacetDistribution_MicrofacetDistribution_3 = R"doc(Create an anisotropic microfacet distribution of the specified type Args: @@ -6686,7 +6907,7 @@ R"doc(Create an anisotropic microfacet distribution of the specified type alpha_v: The surface roughness in the bitangent direction)doc"; -static const char *__doc_mitsuba_MicrofacetDistribution_MicrofacetDistribution_3 = +static const char *__doc_mitsuba_MicrofacetDistribution_MicrofacetDistribution_4 = R"doc(Create a microfacet distribution from a `Properties` data structure)doc"; @@ -6704,18 +6925,34 @@ R"doc(Evaluate the microfacet distribution function Args: m: The microfacet normal)doc"; +static const char *__doc_mitsuba_MicrofacetDistribution_fields = R"doc()doc"; + +static const char *__doc_mitsuba_MicrofacetDistribution_fields_2 = R"doc()doc"; + static const char *__doc_mitsuba_MicrofacetDistribution_is_anisotropic = R"doc(Is this an anisotropic microfacet distribution?)doc"; static const char *__doc_mitsuba_MicrofacetDistribution_is_isotropic = R"doc(Is this an isotropic microfacet distribution?)doc"; -static const char *__doc_mitsuba_MicrofacetDistribution_m_alpha_u = R"doc()doc"; +static const char *__doc_mitsuba_MicrofacetDistribution_labels = R"doc()doc"; -static const char *__doc_mitsuba_MicrofacetDistribution_m_alpha_v = R"doc()doc"; +static const char *__doc_mitsuba_MicrofacetDistribution_m_alpha_u = R"doc(Roughness values and derived constants, exposed for traversal)doc"; + +static const char *__doc_mitsuba_MicrofacetDistribution_m_alpha_v = R"doc(Roughness values and derived constants, exposed for traversal)doc"; + +static const char *__doc_mitsuba_MicrofacetDistribution_m_inv_alpha_u = R"doc()doc"; + +static const char *__doc_mitsuba_MicrofacetDistribution_m_inv_alpha_v = R"doc()doc"; + +static const char *__doc_mitsuba_MicrofacetDistribution_m_isotropic = R"doc()doc"; + +static const char *__doc_mitsuba_MicrofacetDistribution_m_norm = R"doc()doc"; static const char *__doc_mitsuba_MicrofacetDistribution_m_sample_visible = R"doc()doc"; static const char *__doc_mitsuba_MicrofacetDistribution_m_type = R"doc()doc"; +static const char *__doc_mitsuba_MicrofacetDistribution_name = R"doc()doc"; + static const char *__doc_mitsuba_MicrofacetDistribution_pdf = R"doc(Returns the density function associated with the `sample()` function. @@ -6727,6 +6964,10 @@ the `sample()` function. static const char *__doc_mitsuba_MicrofacetDistribution_project_roughness_2 = R"doc(Compute the squared 1D roughness along direction ``v``)doc"; +static const char *__doc_mitsuba_MicrofacetDistribution_same_alpha = +R"doc(Do both roughness values denote the same quantity? Two separately +differentiable parameters can share a primal, hence the combined index.)doc"; + static const char *__doc_mitsuba_MicrofacetDistribution_sample = R"doc(Draw a sample from the microfacet normal distribution and return the associated probability density @@ -6755,10 +6996,10 @@ R"doc(Smith's shadowing-masking function for a single direction m: The microfacet normal)doc"; -static const char *__doc_mitsuba_MicrofacetDistribution_traverse_cb = R"doc()doc"; - static const char *__doc_mitsuba_MicrofacetDistribution_type = R"doc(Return the distribution type)doc"; +static const char *__doc_mitsuba_MicrofacetDistribution_update_derived = R"doc(Recompute the quantities that depend on the roughness values)doc"; + static const char *__doc_mitsuba_MicrofacetType = R"doc(Supported normal distribution functions)doc"; static const char *__doc_mitsuba_MicrofacetType_Beckmann = R"doc(Beckmann distribution derived from Gaussian random surfaces)doc"; @@ -6885,6 +7126,8 @@ static const char *__doc_mitsuba_ObjectType_BSDF = R"doc(A bidirectional reflect static const char *__doc_mitsuba_ObjectType_Emitter = R"doc(Emits radiance, subclasses `Emitter`)doc"; +static const char *__doc_mitsuba_ObjectType_Extremum = R"doc(A medium acceleration structure.)doc"; + static const char *__doc_mitsuba_ObjectType_Film = R"doc(Storage representation of the sensor)doc"; static const char *__doc_mitsuba_ObjectType_Integrator = R"doc(A rendering algorithm aka. `Integrator`)doc"; @@ -7256,6 +7499,15 @@ turn 3-channel ``*color*`` attributes into sRGB upsampling coefficients while adopting the buffer. Clear ``upsample_srgb`` when the producer already stores coefficients.)doc"; +static const char *__doc_mitsuba_PackedMesh_add_tangents = +R"doc(Generate tangents for a mesh with normals and texture coordinates + +This commit populates the packed mesh data with vertex tangent frames, +matching the behavior of `Mesh::compute_tangents()`. In contrast to +this method, the computation is done on the host machine. +Call this after every face and vertex record has been written and +transformed.)doc"; + static const char *__doc_mitsuba_PackedMesh_attrs = R"doc()doc"; static const char *__doc_mitsuba_PackedMesh_backend = R"doc(Dr.Jit backend of the allocated buffers)doc"; @@ -9099,28 +9351,6 @@ three key abstractions implemented on top of these groups, specifically: static const char *__doc_mitsuba_Scene_2 = R"doc()doc"; -static const char *__doc_mitsuba_Scene_PortalData = -R"doc(Packed light portal records - -Environment emitters gather from this buffer so that their kernels do -not depend on the number of portals. Each record holds 12 entries: the -rectangle's center, its two orthogonal half-edge vectors, and its unit -normal, which points into the region that receives light.)doc"; - -static const char *__doc_mitsuba_Scene_PortalData_count = R"doc()doc"; - -static const char *__doc_mitsuba_Scene_PortalData_fields = R"doc()doc"; - -static const char *__doc_mitsuba_Scene_PortalData_fields_2 = R"doc()doc"; - -static const char *__doc_mitsuba_Scene_PortalData_labels = R"doc()doc"; - -static const char *__doc_mitsuba_Scene_PortalData_name = R"doc()doc"; - -static const char *__doc_mitsuba_Scene_PortalData_records = R"doc()doc"; - -static const char *__doc_mitsuba_Scene_PortalData_weight = R"doc(Probability of sampling the portals rather than the emitter's own strategy)doc"; - static const char *__doc_mitsuba_Scene_3 = R"doc()doc"; static const char *__doc_mitsuba_Scene_4 = R"doc()doc"; @@ -9172,6 +9402,28 @@ static const char *__doc_mitsuba_SceneIR_instances = R"doc(Flattened TLAS/IAS in static const char *__doc_mitsuba_SceneIR_top_blases = R"doc(Indices in ``blases`` that belong to top-level scene geometry.)doc"; +static const char *__doc_mitsuba_Scene_PortalData = +R"doc(Packed light portal records + +Environment emitters gather from this buffer so that their kernels do +not depend on the number of portals. Each record holds 12 entries: the +rectangle's center, its two orthogonal half-edge vectors, and its unit +normal, which points into the region that receives light.)doc"; + +static const char *__doc_mitsuba_Scene_PortalData_count = R"doc()doc"; + +static const char *__doc_mitsuba_Scene_PortalData_fields = R"doc()doc"; + +static const char *__doc_mitsuba_Scene_PortalData_fields_2 = R"doc()doc"; + +static const char *__doc_mitsuba_Scene_PortalData_labels = R"doc()doc"; + +static const char *__doc_mitsuba_Scene_PortalData_name = R"doc()doc"; + +static const char *__doc_mitsuba_Scene_PortalData_records = R"doc()doc"; + +static const char *__doc_mitsuba_Scene_PortalData_weight = R"doc(Probability of sampling the portals rather than the emitter's own strategy)doc"; + static const char *__doc_mitsuba_Scene_Scene = R"doc(Instantiate a scene from a `Properties` object)doc"; static const char *__doc_mitsuba_Scene_bbox = R"doc(Return a bounding box surrounding the scene)doc"; @@ -9300,6 +9552,10 @@ by one, since 0 marks non-instanced intersections.)doc"; static const char *__doc_mitsuba_Scene_m_integrator = R"doc()doc"; +static const char *__doc_mitsuba_Scene_m_portal_data = R"doc()doc"; + +static const char *__doc_mitsuba_Scene_m_portals = R"doc(Light portals, excluded from ``m_emitters`` and never sampled directly)doc"; + static const char *__doc_mitsuba_Scene_m_sensors = R"doc()doc"; static const char *__doc_mitsuba_Scene_m_sensors_dr = R"doc()doc"; @@ -9344,6 +9600,10 @@ point or directional emitters/sensors). Returns: The solid angle density of the sample)doc"; +static const char *__doc_mitsuba_Scene_portal_data = R"doc(Return the light portal records)doc"; + +static const char *__doc_mitsuba_Scene_portals = R"doc(Return the list of light portals)doc"; + static const char *__doc_mitsuba_Scene_ray_intersect = R"doc(Intersect a ray with the shapes comprising the scene and return a detailed data structure describing the intersection, if one is found. @@ -9829,8 +10089,6 @@ sampling weight (see `Shape.silhouette_sampling_weight()`). Returns: Silhouette sample record.)doc"; -static const char *__doc_mitsuba_Scene_portals = R"doc(Return the list of light portals)doc"; - static const char *__doc_mitsuba_Scene_sensors = R"doc(Return the list of sensors)doc"; static const char *__doc_mitsuba_Scene_sensors_2 = R"doc(Return the list of sensors (const version))doc"; @@ -9873,6 +10131,8 @@ static const char *__doc_mitsuba_Scene_update_emitter_sampling_distribution = R" static const char *__doc_mitsuba_Scene_update_instance_transforms = R"doc(Repack the per-instance transform records (see below))doc"; +static const char *__doc_mitsuba_Scene_update_portal_data = R"doc(Build the light portal records from ``m_portals``)doc"; + static const char *__doc_mitsuba_Scene_update_silhouette_sampling_distribution = R"doc(Updates the discrete distribution used to select a shape's silhouette)doc"; static const char *__doc_mitsuba_Scene_variant_name = R"doc()doc"; @@ -9916,7 +10176,7 @@ static const char *__doc_mitsuba_Sensor_7 = R"doc()doc"; static const char *__doc_mitsuba_Sensor_Sensor = R"doc(This is both a class and the base of various Mitsuba plugins)doc"; -static const char *__doc_mitsuba_Sensor_class_name = R"doc(This is both a class and the base of various Mitsuba plugins)doc"; +static const char *__doc_mitsuba_Sensor_class_name = R"doc()doc"; static const char *__doc_mitsuba_Sensor_film = R"doc(Return the `Film` instance associated with this sensor)doc"; @@ -10024,9 +10284,9 @@ static const char *__doc_mitsuba_Sensor_traverse = R"doc()doc"; static const char *__doc_mitsuba_Sensor_traverse_cb = R"doc()doc"; -static const char *__doc_mitsuba_Sensor_type = R"doc(This is both a class and the base of various Mitsuba plugins)doc"; +static const char *__doc_mitsuba_Sensor_type = R"doc()doc"; -static const char *__doc_mitsuba_Sensor_variant_name = R"doc(This is both a class and the base of various Mitsuba plugins)doc"; +static const char *__doc_mitsuba_Sensor_variant_name = R"doc()doc"; static const char *__doc_mitsuba_SerializedFlags = R"doc(Flag word of a ``.serialized`` file. The low bits store the @@ -10970,9 +11230,9 @@ static const char *__doc_mitsuba_Spectrum_operator_assign_2 = R"doc()doc"; static const char *__doc_mitsuba_SphericalRectangle = R"doc(Spherical rectangle -This class represents the projection of a planar rectangle onto the -unit sphere centered at a reference point ``p``. It samples directions -towards the rectangle uniformly in solid angle using the area-preserving +This class represents the projection of a planar rectangle onto the unit +sphere centered at a reference point ``p``. It samples directions towards +the rectangle uniformly in solid angle using the area-preserving parametrization by Ureña, Fajardo and King (EGSR 2013). The constructor builds a frame with ``p`` at its origin so that the @@ -10997,10 +11257,6 @@ static const char *__doc_mitsuba_SphericalRectangle_frame = R"doc()doc"; static const char *__doc_mitsuba_SphericalRectangle_k = R"doc(Parametrization constants used by ``sample()``)doc"; -static const char *__doc_mitsuba_SphericalRectangle_operator_assign = R"doc()doc"; - -static const char *__doc_mitsuba_SphericalRectangle_operator_assign_2 = R"doc()doc"; - static const char *__doc_mitsuba_SphericalRectangle_pdf = R"doc(Solid angle density of ``sample()``, zero for directions that miss the rectangle)doc"; static const char *__doc_mitsuba_SphericalRectangle_rect = R"doc()doc"; @@ -12132,6 +12388,44 @@ static const char *__doc_mitsuba_Timer_start = R"doc()doc"; static const char *__doc_mitsuba_Timer_value = R"doc()doc"; +static const char *__doc_mitsuba_TrackingState = +R"doc(\brief State carried through extremum traversal to accumulate the throughput +and its PDF. + +Can be used to for delta tracking, ratio tracking, and residual ratio tracking. +Note: Since the number of required dimensions is different for all pixel +samples, ``rng`` is used to sample distances and event types.)doc"; + +static const char *__doc_mitsuba_TrackingState_TrackingState = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_TrackingState_2 = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_TrackingState_3 = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_fields = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_fields_2 = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_has_spectral_extinction = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_labels = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_mei = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_name = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_operator_assign = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_operator_assign_2 = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_ray = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_rng = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_target_ot = R"doc()doc"; + +static const char *__doc_mitsuba_TrackingState_throughput = R"doc()doc"; + static const char *__doc_mitsuba_Transform = R"doc(Unified homogeneous coordinate transformation @@ -12318,10 +12612,8 @@ R"doc(Extend the lifetime of a Python object whose contents were reported to `put()` A plugin written in Python can report a value that exists only for the -duration of the call, which the bindings have no way of detecting. They -therefore offer the underlying Python object here, and a callback that -holds on to reported pointers should reference it. The argument is a -borrowed `PyObject *`.)doc"; +duration of the call. It must invoke this funcdtion with the borrowed +``PyObject*`` so that the caller can keep it alive.)doc"; static const char *__doc_mitsuba_TraversalCallback_put = R"doc()doc"; @@ -12464,6 +12756,32 @@ R"doc(Write an encoded form of the volume grid to a stream Args: stream: Target stream that will receive the encoded output)doc"; +static const char *__doc_mitsuba_VolumeParametrization = +R"doc(\brief Frame parameters of a volume's local coordinates in world space. + +.. note:: + Currently only holds a ``to_world`` matrix but becomes useful in the + presence of volume that use different coordinate systems.)doc"; + +static const char *__doc_mitsuba_VolumeParametrization_VolumeParametrization = R"doc()doc"; + +static const char *__doc_mitsuba_VolumeParametrization_VolumeParametrization_2 = R"doc()doc"; + +static const char *__doc_mitsuba_VolumeParametrization_to_world = R"doc()doc"; + +static const char *__doc_mitsuba_Volume_PinGuard = +R"doc(\brief A scoped guard that pins the reference count of the volume. + +Use for bulk operations in scalar mode.)doc"; + +static const char *__doc_mitsuba_Volume_PinGuard_PinGuard = R"doc()doc"; + +static const char *__doc_mitsuba_Volume_PinGuard_PinGuard_2 = R"doc()doc"; + +static const char *__doc_mitsuba_Volume_PinGuard_operator_assign = R"doc()doc"; + +static const char *__doc_mitsuba_Volume_PinGuard_volume = R"doc()doc"; + static const char *__doc_mitsuba_Volume_Volume = R"doc()doc"; static const char *__doc_mitsuba_Volume_bbox = R"doc(Returns the bounding box of the volume)doc"; @@ -12496,6 +12814,17 @@ R"doc(Evaluate this volume as a n-channel float quantity This interface is specifically intended to encode a variable number of parameters. Pointer allocation/deallocation must be performed by the caller.)doc"; +static const char *__doc_mitsuba_Volume_extremum = +R"doc(\brief Compute local extrema over a spatial region + +Returns the minorant (minimum) and majorant (maximum) value over the +specified bounding box region. Only fully implemented for grid-based +volumes; the default falls back to a placeholder minorant of 0 and the +volume's global majorant. + +\param bbox Bounding box defining the query region in local space +\return (minorant, majorant) pair)doc"; + static const char *__doc_mitsuba_Volume_m_bbox = R"doc(Bounding box)doc"; static const char *__doc_mitsuba_Volume_m_channel_count = R"doc(Number of channels stored in the volume)doc"; @@ -12510,6 +12839,12 @@ the maximum value for each channel. Pointer allocation/deallocation must be performed by the caller.)doc"; +static const char *__doc_mitsuba_Volume_parametrization = R"doc(Returns the world-space parametrization of the volume's local coordinates.)doc"; + +static const char *__doc_mitsuba_Volume_pin = R"doc()doc"; + +static const char *__doc_mitsuba_Volume_pin_ref_count = R"doc(Pin the reference count of the data that constitutes the volume, e.g. a Texture.)doc"; + static const char *__doc_mitsuba_Volume_resolution = R"doc(Returns the resolution of the volume, assuming that it is based on a discrete representation. @@ -12522,6 +12857,8 @@ static const char *__doc_mitsuba_Volume_traverse_cb = R"doc()doc"; static const char *__doc_mitsuba_Volume_type = R"doc()doc"; +static const char *__doc_mitsuba_Volume_unpin_ref_count = R"doc(Unpin the reference count.)doc"; + static const char *__doc_mitsuba_Volume_update_bbox = R"doc()doc"; static const char *__doc_mitsuba_Volume_variant_name = R"doc()doc"; @@ -13100,6 +13437,9 @@ at a planar interface between two dielectrics means that the surface normal is pointing into the region of lower density. + inv_eta: Reciprocal of ``eta``. Callers that store this value can pass + it to avoid a division per evaluation. + Returns: A tuple ``(F, cos_theta_t, eta_it, eta_ti)`` where @@ -13116,6 +13456,8 @@ at a planar interface between two dielectrics factor that must be applied to the X and Y component of the refracted direction.)doc"; +static const char *__doc_mitsuba_fresnel_2 = R"doc()doc"; + static const char *__doc_mitsuba_fresnel_conductor = R"doc(Calculates the unpolarized Fresnel reflection coefficient at a planar interface of a conductor, i.e. a surface with a complex-valued relative index @@ -13297,6 +13639,8 @@ static const char *__doc_mitsuba_hasher = R"doc()doc"; static const char *__doc_mitsuba_hasher_operator_call = R"doc()doc"; +static const char *__doc_mitsuba_index_spectrum = R"doc(Helper function to index the channel of an ``UnpolarizedSpectrum``.)doc"; + static const char *__doc_mitsuba_ior_from_file = R"doc()doc"; static const char *__doc_mitsuba_key_iterator = R"doc(Forward declaration of iterator)doc"; @@ -13944,6 +14288,8 @@ static const char *__doc_mitsuba_operator_lshift_31 = R"doc()doc"; static const char *__doc_mitsuba_operator_lshift_32 = R"doc()doc"; +static const char *__doc_mitsuba_operator_lshift_33 = R"doc()doc"; + static const char *__doc_mitsuba_operator_sub = R"doc(Subtracting two points should always yield a vector)doc"; static const char *__doc_mitsuba_operator_sub_2 = R"doc(Subtracting a vector from a point should always yield a point)doc"; diff --git a/include/mitsuba/render/extremum.h b/include/mitsuba/render/extremum.h new file mode 100644 index 0000000000..10e8ccdc47 --- /dev/null +++ b/include/mitsuba/render/extremum.h @@ -0,0 +1,137 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#include + +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 +class MI_EXPORT_LIB Extremum : public JitObject> { +public: + MI_IMPORT_TYPES(Medium, Sampler, Volume) + + using TrackingStateType = TrackingState; + using TrackingFunctionType = TrackingFunction; + + /// 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 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() + +//! @} +// ----------------------------------------------------------------------- diff --git a/include/mitsuba/render/extremum_segment.h b/include/mitsuba/render/extremum_segment.h new file mode 100644 index 0000000000..49d94bcd7e --- /dev/null +++ b/include/mitsuba/render/extremum_segment.h @@ -0,0 +1,98 @@ +#pragma once + +#include + +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 +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(dr::Infinity, size); + maxt = dr::full(-dr::Infinity, size); + value = dr::zeros(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; + maxt = -dr::Infinity; + } + + /// 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) diff --git a/include/mitsuba/render/fwd.h b/include/mitsuba/render/fwd.h index 55029a6532..4299bc4bdb 100644 --- a/include/mitsuba/render/fwd.h +++ b/include/mitsuba/render/fwd.h @@ -13,6 +13,7 @@ template class DirectedEdge; template class OptixDenoiser; template class Emitter; template class Endpoint; +template class Extremum; template class Film; template class ImageBlock; template class Integrator; @@ -41,6 +42,7 @@ template struct PositionSample; template struct BSDFSample3; template struct SilhouetteSample; template struct PhaseFunctionContext; +template struct ExtremumSegment; template struct Interaction; template struct MediumInteraction; template struct SurfaceInteraction; @@ -132,6 +134,8 @@ template struct RenderAliases { using Emitter = mitsuba::Emitter; using Endpoint = mitsuba::Endpoint; using Medium = mitsuba::Medium; + using Extremum = mitsuba::Extremum; + using ExtremumSegment = mitsuba::ExtremumSegment; using PhaseFunction = mitsuba::PhaseFunction; using Film = mitsuba::Film; using ImageBlock = mitsuba::ImageBlock; @@ -145,6 +149,7 @@ template struct RenderAliases { using ObjectPtr = dr::replace_scalar_t; using BSDFPtr = dr::replace_scalar_t; using MediumPtr = dr::replace_scalar_t; + using ExtremumPtr = dr::replace_scalar_t; using PhaseFunctionPtr = dr::replace_scalar_t; using ShapePtr = dr::replace_scalar_t; using MeshPtr = dr::replace_scalar_t; @@ -196,6 +201,7 @@ template 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 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 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; \ diff --git a/include/mitsuba/render/medium.h b/include/mitsuba/render/medium.h index 833e5e0fc3..f33e1d78d3 100644 --- a/include/mitsuba/render/medium.h +++ b/include/mitsuba/render/medium.h @@ -11,7 +11,7 @@ NAMESPACE_BEGIN(mitsuba) template class MI_EXPORT_LIB Medium : public JitObject> { 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> { 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 + 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> { bool m_sample_emitters; bool m_is_homogeneous; bool m_has_spectral_extinction; + ref 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() // ----------------------------------------------------------------------- diff --git a/include/mitsuba/render/tracking.h b/include/mitsuba/render/tracking.h new file mode 100644 index 0000000000..baf012e627 --- /dev/null +++ b/include/mitsuba/render/tracking.h @@ -0,0 +1,87 @@ +#pragma once + +#include +#include +#include +#include +#include + +NAMESPACE_BEGIN(mitsuba) + +/** + * \brief State carried through extremum traversal to accumulate the throughput + * and its PDF. + * + * Can be used to for delta tracking, ratio tracking, and residual ratio tracking. + * Note: Since the number of required dimensions is different for all pixel + * samples, ``rng`` is used to sample distances and event types. + * + */ +template< typename Float, typename Spectrum > +struct TrackingState { + MI_IMPORT_TYPES() + + Ray3f ray; + dr::PCG32 rng; + MediumInteraction3f mei; + Float target_ot; + Mask has_spectral_extinction; + + // Note that ``throughput`` is shared between algorithm and should be accumulated + // accordingly. If used for volpathmis, new members and data types will need to + // be introduced. + UnpolarizedSpectrum throughput; + + DRJIT_STRUCT(TrackingState, ray, rng, mei, target_ot, \ + has_spectral_extinction, throughput) +}; + +/** + * \brief Signature of the tracking function callback accepted by + * ``Extremum::traverse_extremum``. + * + * \param segment + * An extremum segment along a ray. + * \param state + * Pointer to the tracking state that holds interaction information and + * accumulates throughput and pdfs. + * \param channel + * The channel to use for sampling. + * \param active + * Represents the active lanes. + * + * + * \return A pair (advance, active): + * advance: If true, tracking has exited the segment and requires a + * new one. If false, repeat the loop with the same segment. + * active: Represent active lanes. Lanes that have sampled a real + * interaction or terminated for other reasons will return + * ``false``, prompting the termination of the traversal. + */ +template< typename Float, typename Spectrum > +using TrackingFunction = std::function< + std::pair, dr::mask_t>( + const ExtremumSegment& /*segment*/, + TrackingState* /*state*/, + const dr::uint32_array_t& /*channel*/, + dr::mask_t /*active*/ +)>; + +/// Helper function to index the channel of an ``UnpolarizedSpectrum``. +template< typename Float, typename Spectrum > +MI_INLINE +Float index_spectrum( + const unpolarized_spectrum_t &spec, + const dr::uint32_array_t &idx +) { + Float m = spec[0]; + if constexpr (is_rgb_v) { // Handle RGB rendering + dr::masked(m, idx == 1u) = spec[1]; + dr::masked(m, idx == 2u) = spec[2]; + } else { + DRJIT_MARK_USED(idx); + } + return m; +} + +NAMESPACE_END(mitsuba) diff --git a/include/mitsuba/render/volume.h b/include/mitsuba/render/volume.h index bf16fccaf2..2020a31469 100644 --- a/include/mitsuba/render/volume.h +++ b/include/mitsuba/render/volume.h @@ -6,6 +6,7 @@ #include #include #include +#include NAMESPACE_BEGIN(mitsuba) @@ -13,7 +14,7 @@ NAMESPACE_BEGIN(mitsuba) template class MI_EXPORT_LIB Volume : public JitObject> { public: - MI_IMPORT_TYPES(Texture) + MI_IMPORT_TYPES(Texture, Extremum) // ====================================================================== // Volume interface @@ -60,6 +61,26 @@ class MI_EXPORT_LIB Volume : public JitObject> { */ virtual void max_per_channel(ScalarFloat *out) const; + /** + * \brief Compute local extrema over a spatial region + * + * Returns the minorant (minimum) and majorant (maximum) value over the + * specified bounding box region. + * + * \param bbox Bounding box defining the query region, expected to be + * in the volume's local coordinate frame + * \return (minorant, majorant) pair + */ + virtual std::pair + extremum(BoundingBox3f bbox) const; + + /// Returns the world-space parametrization of the volume's local coordinates. + virtual VolumeParametrization parametrization() const { + VolumeParametrization param; + param.to_world = m_to_local.inverse(); + return param; + } + /// Returns the bounding box of the volume ScalarBoundingBox3f bbox() const { return m_bbox; } @@ -92,6 +113,18 @@ class MI_EXPORT_LIB Volume : public JitObject> { MI_DECLARE_PLUGIN_BASE_CLASS(Volume) + /// A scoped guard that pins the reference count for bulk operations in scalar mode. + struct PinGuard { + const Volume *volume; + explicit PinGuard(const Volume *v) : volume(v) { volume->pin_ref_count(); } + ~PinGuard() { volume->unpin_ref_count(); } + + PinGuard(const PinGuard &) = delete; + PinGuard &operator=(const PinGuard &) = delete; + }; + + virtual PinGuard pin() const { return PinGuard(this); } + protected: Volume(const Properties &props); @@ -108,6 +141,11 @@ class MI_EXPORT_LIB Volume : public JitObject> { m_bbox.expand(to_world * ScalarPoint3f(1.f, 1.f, 1.f)); } + /// Pin the reference count of the data that constitutes the volume, e.g. a Texture. + virtual void pin_ref_count() const {} + /// Unpin the reference count. + virtual void unpin_ref_count() const {} + protected: /// Used to bring points in world coordinates to local coordinates. ScalarAffineTransform4f m_to_local; diff --git a/include/mitsuba/render/volume_utils.h b/include/mitsuba/render/volume_utils.h new file mode 100644 index 0000000000..03b7667fb2 --- /dev/null +++ b/include/mitsuba/render/volume_utils.h @@ -0,0 +1,28 @@ +#pragma once + +#include +#include +#include + +NAMESPACE_BEGIN(mitsuba) + +/** + * \brief Frame parameters of a volume's local coordinates in world space. + * + * .. note:: + * Currently only holds a ``to_world`` matrix but becomes useful in the + * presence of volume that use different coordinate systems. + */ +template +struct VolumeParametrization { + + using AffineTransform4f = Transform, true>; + + AffineTransform4f to_world; + + VolumeParametrization(): to_world(AffineTransform4f()) {} + + VolumeParametrization(AffineTransform4f to_world): to_world(to_world) {} +}; + +NAMESPACE_END(mitsuba) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 248e548f50..d497c6e367 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -32,6 +32,7 @@ add_subdirectory(mitsuba) add_subdirectory(bsdfs) add_subdirectory(emitters) +add_subdirectory(extrema) add_subdirectory(films) add_subdirectory(integrators) add_subdirectory(media) diff --git a/src/core/parser.cpp b/src/core/parser.cpp index 8910da41d9..455a351050 100644 --- a/src/core/parser.cpp +++ b/src/core/parser.cpp @@ -1349,7 +1349,8 @@ void transform_merge_equivalent(const ParserConfig &/*config*/, ParserState &sta // Skip merging for emitters and shapes if (state[repr].type == ObjectType::Emitter || - state[repr].type == ObjectType::Shape) + state[repr].type == ObjectType::Shape || + state[repr].type == ObjectType::Extremum) continue; // Try to find an equivalent node diff --git a/src/core/plugin.cpp b/src/core/plugin.cpp index 8b1f6f8f56..11195b23fe 100644 --- a/src/core/plugin.cpp +++ b/src/core/plugin.cpp @@ -229,6 +229,7 @@ std::string_view plugin_type_name(ObjectType ot) { case ObjectType::Integrator: return "integrator"; case ObjectType::PhaseFunction: return "phase"; case ObjectType::ReconstructionFilter: return "rfilter"; + case ObjectType::Extremum: return "extremum"; } return "invalid"; // (to avoid a compiler warning; this should never happen) } diff --git a/src/extrema/CMakeLists.txt b/src/extrema/CMakeLists.txt new file mode 100644 index 0000000000..67b5fb2651 --- /dev/null +++ b/src/extrema/CMakeLists.txt @@ -0,0 +1,6 @@ +set(MI_PLUGIN_PREFIX "extremum") + +add_plugin(extremum_global extremum_global.cpp) +add_plugin(extremum_grid extremum_grid.cpp) + +set(MI_PLUGIN_TARGETS "${MI_PLUGIN_TARGETS}" PARENT_SCOPE) diff --git a/src/extrema/extremum_global.cpp b/src/extrema/extremum_global.cpp new file mode 100644 index 0000000000..3820be01de --- /dev/null +++ b/src/extrema/extremum_global.cpp @@ -0,0 +1,92 @@ +#include +#include +#include +#include +#include + +NAMESPACE_BEGIN(mitsuba) + +/**! +.. _extremum-extremum_global: + +Extremum global structure (:monosp:`extremum_global`) +----------------------------------------------------- + +This plugin holds the global minorant and majorant values of a volume. +At runtime, traversal is performed via a single segment determined by the +passed ``mint`` and ``maxt`` values. +*/ + +template +class ExtremumGlobal final : public Extremum { +public: + MI_IMPORT_BASE(Extremum, m_bbox, m_scale) + MI_IMPORT_TYPES(Volume) + + using TrackingStateType = TrackingState; + using TrackingFunctionType = TrackingFunction; + + ExtremumGlobal(const Properties &props) : Base(props) {} + + void build(const Volume *volume) override { + // placeholder minorant value + m_minorant = 0.f; + m_majorant = volume->max(); + } + + TrackingStateType traverse_extremum( + const Ray3f &/*ray*/, + Float mint, + Float maxt, + UInt32 channel, + TrackingStateType state, + const TrackingFunctionType &func, + Mask active + ) const override { + ExtremumSegment segment(mint, maxt, m_scale*m_minorant, m_scale*m_majorant); + + struct LoopState { + ExtremumSegment segment; + TrackingStateType state; + Mask advance; + Mask active; + + DRJIT_STRUCT(LoopState, segment, state, advance, active) + } ls { + segment, + state, + /*advance =*/true, + active + }; + + dr::tie(ls) = dr::while_loop( + dr::make_tuple(ls), + [](const LoopState &ls){ return ls.active; }, + [func, channel](LoopState &ls){ + std::tie(ls.advance, ls.active) = + func(ls.segment, &ls.state, channel, ls.active); + ls.active &= !ls.advance; + }); + + return ls.state; + } + + std::string to_string() const override { + std::ostringstream oss; + oss << "ExtremumGlobal[" << std::endl + << " minorant = " << m_minorant << "," << std::endl + << " majorant = " << m_majorant << "," << std::endl + << " scale = " << m_scale << "," << std::endl + << "]"; + return oss.str(); + } + + MI_DECLARE_CLASS(ExtremumGlobal) + +private: + ScalarFloat m_minorant; + ScalarFloat m_majorant; +}; + +MI_EXPORT_PLUGIN(ExtremumGlobal) +NAMESPACE_END(mitsuba) diff --git a/src/extrema/extremum_grid.cpp b/src/extrema/extremum_grid.cpp new file mode 100644 index 0000000000..565e5e2ec8 --- /dev/null +++ b/src/extrema/extremum_grid.cpp @@ -0,0 +1,336 @@ +#include +#include +#include +#include +#include +#include + +NAMESPACE_BEGIN(mitsuba) + +/**! +.. _extremum-extremum_grid: + +Extremum grid structure (:monosp:`extremum_grid`) +------------------------------------------------- + +.. pluginparameters:: + + * - resolution + - |vector| + - Grid resolution along the XYZ axis. Does not have to be a multiple of + the underlying volume. Default: [1,1,1] + +This plugin creates a regular grid structure storing local extrema value. +The grid is constructed by querying a dense grid of bounding boxes corresponding +to the volumes cells. +``traverse_extremum`` uses a DDA (Digital Differential Analyzer) to efficiently +track through the grid. +*/ + +template +class ExtremumGrid final : public Extremum { +public: + MI_IMPORT_BASE(Extremum, m_bbox, m_scale) + MI_IMPORT_TYPES(Volume) + + using TrackingStateType = TrackingState; + using TrackingFunctionType = TrackingFunction; + using FloatStorage = DynamicBuffer; + + ExtremumGrid(const Properties &props) : Base(props) { + + // Resolution Parameters + m_resolution = props.get("resolution", ScalarVector3i(1,1,1)); + } + + void build(const Volume *volume) override { + + VolumeParametrization volume_param = volume->parametrization(); + + m_to_local = volume_param.to_world.inverse(); + + ScalarVector3i clamped = dr::clip(m_resolution, 1, volume->resolution()); + if (dr::any(clamped != m_resolution)) { + Log(Info, + "ExtremumGrid: requested resolution %s is finer than the " + "underlying volume's resolution %s; clamping to %s.", + m_resolution, volume->resolution(), clamped); + m_resolution = clamped; + } + + build_grid(volume, m_resolution); + } + + TrackingStateType traverse_extremum( + const Ray3f &ray, + Float mint, + Float maxt, + UInt32 channel, + TrackingStateType state, + const TrackingFunctionType &func, + Mask active + ) const override { + return traverse_dda( + func, + state, + ray, + mint, + maxt, + channel, + active + ); + } + + void traverse(TraversalCallback *cb) override { + cb->put("resolution", m_resolution, ParamFlags::NonDifferentiable); + Base::traverse(cb); + } + + std::string to_string() const override { + std::ostringstream oss; + oss << "ExtremumGrid[" << std::endl + << " resolution = " << m_resolution << "," << std::endl + << " bbox = " << m_bbox << "," << std::endl + << " scale = " << m_scale << "," << std::endl + << "]"; + return oss.str(); + } + + MI_DECLARE_CLASS(ExtremumGrid) + +private: + + /** + * \brief Build the extremum grid from a volume + * + * This method constructs a lower-resolution grid where each cell stores + * the majorant (maximum) extinction value over the corresponding region + * of the high-resolution volume. + */ + void build_grid(const Volume *volume, ScalarVector3i resolution) { + + // local space supergrid cell size + const ScalarVector3f cell_size = dr::rcp(ScalarVector3f(resolution)); + + ScalarVector2f safety_factor(1.f - dr::Epsilon, + 1.f + dr::Epsilon); + + // Allocate extremum grid data + size_t n = dr::prod(resolution); + + size_t n_threads = pool_size() + 1; + size_t grain_size = std::max(n / (4 * n_threads), (size_t) 1); + + m_extremum_grid = dr::empty(n * 2); + + // Early return if using the global majorant. + if (n == 1) { + ScalarFloat max = volume->max(); + // Global minorant tracking is not yet implemented on the volume + // side; use a placeholder minorant until it lands. + ScalarFloat min = 0.f; + dr::scatter(m_extremum_grid, + Vector2f(min, max) * safety_factor, + UInt32(0)); + return; + } + + if constexpr (!dr::is_jit_v) { + auto guard = volume->pin(); + + dr::parallel_for( + dr::blocked_range(0, n, grain_size), + [&](const dr::blocked_range &range) { + // Recover x, y, z from block start (one-time div/mod per + // block) + + for (auto idx = range.begin(); idx != range.end(); ++idx) { + // Store in linear array (Z-slowest, X-fastest) + int32_t x = idx % resolution.x(); + int32_t y = (idx / resolution.x()) % resolution.y(); + int32_t z = idx / (resolution.x() * resolution.y()); + + ScalarPoint3f cell_min = + ScalarVector3f(x, y, z) * cell_size; + ScalarPoint3f cell_max = cell_min + cell_size; + ScalarBoundingBox3f cell_bounds( + cell_min + math::RayEpsilon, + cell_max - math::RayEpsilon); + + // Query volume for local extremum, currently assume + // local bounds. + auto [min, maj] = volume->extremum(cell_bounds); + + dr::scatter(m_extremum_grid, + Vector2f(min, maj) * safety_factor, + UInt32(idx)); + } + }); + } else { + + UInt32 idx = dr::arange((uint32_t) n); + + UInt32 x = idx % resolution.x() ; + UInt32 y = (idx / resolution.x()) % resolution.y(); + UInt32 z = idx / (resolution.x() * resolution.y()); + + Point3f cell_min = Vector3f(x, y, z) * cell_size; + Point3f cell_max = cell_min + cell_size; + BoundingBox3f cell_bounds( + cell_min + math::RayEpsilon, + cell_max - math::RayEpsilon + ); + + auto [min, maj]= volume->extremum(cell_bounds); + + dr::scatter(m_extremum_grid, min * safety_factor.x(), idx*2); + dr::scatter(m_extremum_grid, maj * safety_factor.y(), idx*2+1); + dr::sync_thread(); + } + + Log(Info, "Extremum grid constructed successfully"); + } + + /** \brief General regular grid DDA traversal algorithm. + * + * This method traverses the regular grid along the provided ray using the + * DDA algorithm. + * + * \param func Function to be called at each step of the traversal. + * Must be of type \ref TrackingFunction (tracking.h). + * \param state The payload passed to ``func``. + * \param ray The ray along which the structure is traversed. + * \param mint The minimum distance along the ray. + * \param maxt The maximum distance along the ray. + * \param active + * + * \return + * Returns the final state at the end of the traversal. + */ + template + std::decay_t traverse_dda( + FuncT&& func, + StateT&& state, + const Ray3f& ray, + Float mint, + Float maxt, + UInt32 channel, + Mask active + ) const { + using StateD = std::decay_t; + + ExtremumSegment segment = dr::zeros(); + + // Currently assuming that the majorant aligns perfectly with the + // volume and that values outside the bbox cannot be evaluated. + // Transform ray to local grid coordinates [0,res]³ + Vector3f res = Vector3f(m_resolution); + Ray3f local_ray((m_to_local * ray.o) * res, // Normalize origin + (m_to_local * ray.d) * res, // Normalize direction + ray.time, ray.wavelengths); + Vector3f rcp_d = dr::rcp(local_ray.d); + auto inf_t = local_ray.d == 0.f; + auto d_pos = local_ray.d >= 0.f; + + Float t_min = mint; + Float t_max = maxt; + + active &= t_max > t_min && dr::isfinite(t_max); + + // Advance the ray to the start of the interval + local_ray.o = dr::fmadd(local_ray.d, t_min, local_ray.o); + t_max = t_max - t_min; + t_min = 0.f; + + // Compute the integer step direction + Vector3i step = dr::select(d_pos, 1, -1); + Vector3f abs_rcp_d = abs(rcp_d); + + // Integer grid coordinates + Vector3i pi = dr::floor2int(local_ray.o); + + // Fractional entry position + Vector3f p0 = local_ray.o - Vector3f(pi); + // Step size to next interaction + Vector3f dt_v = + dr::select(d_pos, dr::fmadd(-p0, rcp_d, rcp_d), -p0 * rcp_d); + dr::masked(dt_v, inf_t) = dr::Infinity; + + struct LoopState { + ExtremumSegment segment; + StateD state; + Mask advance; + Mask active; + Vector3f dt_v; + Vector3i pi; + Float t_rem; + + DRJIT_STRUCT(LoopState, segment, state, advance, active, dt_v, \ + pi, t_rem) + } ls = { + segment, + state, + /*advance=*/active, + active, + dt_v, + pi, + t_max + }; + + dr::tie(ls) = dr::while_loop( + dr::make_tuple(ls), + [](const LoopState& ls) { return ls.active; }, + [this, func, step, abs_rcp_d, t_max, mint, channel](LoopState& ls) { + + ExtremumSegment& segment = ls.segment; + StateD& state = ls.state; + Mask& advance = ls.advance; + Mask& active = ls.active; + Vector3f& dt_v = ls.dt_v; + Vector3i& pi = ls.pi; + Float& t_rem = ls.t_rem; + + // Select the smallest step. + Float dt = dr::minimum(dr::min(dt_v), t_rem); + auto mask = dt_v == dt; + + // Note: not multiplying the index by 2 because we gather using + // Vector2f. + Vector3i piw = dr::clip(pi, 0, m_resolution - 1); + UInt32 idx = dr::fmadd(dr::fmadd(piw.z(), m_resolution.y(), piw.y()), + m_resolution.x(), piw.x()); + + const Vector2f extremum = + m_scale * dr::gather(m_extremum_grid, idx); + + // Store segment for lanes that reached target + Float t_curr = mint + t_max - t_rem; + dr::masked(segment, active) = ExtremumSegment( + t_curr, + t_curr + dt, + extremum + ); + + std::tie(advance, active) = func(segment, &state, channel, active); + + // Advance + dr::masked(dt_v, advance) = dr::select(mask, abs_rcp_d, dt_v - dt); + dr::masked(pi, mask && advance) += step; + dr::masked(t_rem, advance) -= dt; + + active &= t_rem > 0.f; + }, + "DDA Traversal"); + + return ls.state; + }; + +private: + /// Grid storing pre-computed local majorants + FloatStorage m_extremum_grid; + ScalarVector3i m_resolution; + ScalarAffineTransform4f m_to_local; +}; + +MI_EXPORT_PLUGIN(ExtremumGrid) +NAMESPACE_END(mitsuba) diff --git a/src/extrema/tests/test_extremum_global.py b/src/extrema/tests/test_extremum_global.py new file mode 100644 index 0000000000..25b92006da --- /dev/null +++ b/src/extrema/tests/test_extremum_global.py @@ -0,0 +1,112 @@ +import mitsuba as mi +import numpy as np +import pytest + + +def _make_grid_volume(values, n): + data = np.array(values, dtype=float).reshape(n, 1, 1) + return mi.load_dict( + { + "type": "gridvolume", + "grid": mi.VolumeGrid(data), + "filter_type": "nearest", + "accel": False, + } + ) + + +def test_build(variant_scalar_rgb): + volume = _make_grid_volume([1.0, 2.0, 3.0, 4.0], 4) + extremum = mi.load_dict({"type": "extremum_global"}) + extremum.update_extremum(volume.bbox(), volume) + + ray = mi.Ray3f(o=[0, 0, 0], d=[0, 0, 1]) + distance, leftover_ot = extremum.sample_test(ray, 0.0, 1.0, target_ot=1.0) + + assert np.allclose(distance, 0.25) + assert np.allclose(leftover_ot, 1.0) + + +def test_set_scale(variant_scalar_rgb): + volume = _make_grid_volume([1.0, 2.0, 3.0, 4.0], 4) + extremum = mi.load_dict({"type": "extremum_global"}) + extremum.update_extremum(volume.bbox(), volume, 2.0) + + ray = mi.Ray3f(o=[0, 0, 0], d=[0, 0, 1]) + distance, leftover_ot = extremum.sample_test(ray, 0.0, 1.0, target_ot=1.0) + + assert np.allclose(distance, 0.125) + assert np.allclose(leftover_ot, 1.0) + + +@pytest.mark.parametrize("medium_type", ["heterogeneous", "homogeneous"]) +def test_update_on_sigma_t_change(variant_scalar_rgb, medium_type): + n = 4 + before = [1.0, 2.0, 3.0, 4.0] + after = [5.0, 6.0, 7.0, 8.0] + + volume = _make_grid_volume(before, n) + medium = mi.load_dict({"type": medium_type, "sigma_t": volume, "albedo": 0.5}) + + # Ground truth: an extremum structure built directly from the "after" + # data, independently of the update mechanism under test. + ref_volume = _make_grid_volume(after, n) + ref_extremum = mi.load_dict({"type": "extremum_global"}) + ref_extremum.update_extremum(ref_volume.bbox(), ref_volume) + + params = mi.traverse(medium) + params["sigma_t.data"] = mi.TensorXf(np.array(after).reshape(n, 1, 1, 1)) + params.update() + + ray = mi.Ray3f(o=[0, 0, 0], d=[0, 0, 1]) + expected = ref_extremum.sample_test(ray, 0.0, 1.0, target_ot=1.0) + got = medium.extremum().sample_test(ray, 0.0, 1.0, target_ot=1.0) + assert np.allclose(got, expected) + + +@pytest.mark.parametrize("medium_type", ["heterogeneous", "homogeneous"]) +def test_update_on_scale_change(variant_scalar_rgb, medium_type): + # A scale-only update must go through `set_scale()` without building + # the extremum structure from `sigma_t`. + volume = _make_grid_volume([1.0, 2.0, 3.0, 4.0], 4) + medium = mi.load_dict( + {"type": medium_type, "sigma_t": volume, "albedo": 0.5, "scale": 1.0} + ) + + params = mi.traverse(medium) + params["scale"] = mi.Float(2.0) + params.update() + + ray = mi.Ray3f(o=[0, 0, 0], d=[0, 0, 1]) + distance, leftover_ot = medium.extremum().sample_test(ray, 0.0, 1.0, target_ot=1.0) + + assert np.allclose(distance, 0.125) + assert np.allclose(leftover_ot, 1.0) + + +def test_sample_test_sampled(variant_scalar_rgb): + # Homogeneous majorant of 3.0 over [1, 4]: segment_ot = 3 * 3 = 9. + volume = _make_grid_volume([3.0, 3.0], 2) + extremum = mi.load_dict({"type": "extremum_global"}) + extremum.update_extremum(volume.bbox(), volume) + + ray = mi.Ray3f(o=[0.5, 0.5, 0], d=[0, 0, 1]) + distance, leftover_ot = extremum.sample_test(ray, 1.0, 4.0, target_ot=6.0) + + # target_ot < segment_ot: interaction sampled inside the segment. + assert np.allclose(distance, 3.0) + assert np.allclose(leftover_ot, 6.0) + + +def test_sample_test_escapes(variant_scalar_rgb): + # Homogeneous majorant of 3.0 over [1, 4]: segment_ot = 3 * 3 = 9. + volume = _make_grid_volume([3.0, 3.0], 2) + extremum = mi.load_dict({"type": "extremum_global"}) + extremum.update_extremum(volume.bbox(), volume) + + ray = mi.Ray3f(o=[0.5, 0.5, 0], d=[0, 0, 1]) + distance, leftover_ot = extremum.sample_test(ray, 1.0, 4.0, target_ot=12.0) + + # target_ot > segment_ot: ray exits the medium before sampling. + assert np.isinf(distance) + assert np.allclose(leftover_ot, 3.0) diff --git a/src/extrema/tests/test_extremum_grid.py b/src/extrema/tests/test_extremum_grid.py new file mode 100644 index 0000000000..cf8770b93b --- /dev/null +++ b/src/extrema/tests/test_extremum_grid.py @@ -0,0 +1,184 @@ +import mitsuba as mi +import numpy as np +import pytest + + +def test_build_resolution_clamped_to_volume(variant_scalar_rgb): + # A fixed resolution requested finer than the volume's own resolution + # along some axis is not a valid extremum grid (it wouldn't be coarser + # than what it summarizes), so that axis must be clamped down to match + # the volume's resolution. + n_x, n_y, n_z = 8, 8, 4 + n_prod = n_x * n_y * n_z + data = np.linspace(1, n_prod, n_prod).reshape(n_x, n_y, n_z) + volume_grid = mi.VolumeGrid(data.transpose(2, 1, 0)) + + volume = mi.load_dict( + { + "type": "gridvolume", + "grid": volume_grid, + "filter_type": "nearest", + "accel": False, + } + ) + + extremum_struct = mi.load_dict( + { + "type": "extremum_grid", + "resolution": mi.ScalarVector3i(n_x, n_y * 2, n_z), + } + ) + extremum_struct.update_extremum(volume.bbox(), volume) + + resolution = np.array(mi.traverse(extremum_struct)["resolution"]) + assert np.array_equal(resolution, [n_x, n_y, n_z]) + + +def _make_grid_volume(values, n): + data = np.array(values, dtype=float).reshape(n, 1, 1) + return mi.load_dict( + { + "type": "gridvolume", + "grid": mi.VolumeGrid(data), + "filter_type": "nearest", + "accel": False, + } + ) + + +def _make_medium(medium_type, volume, extremum): + return mi.load_dict( + { + "type": medium_type, + "sigma_t": volume, + "albedo": 0.5, + "extremum": extremum, + } + ) + + +def test_update_on_sigma_t_change(variant_scalar_rgb): + n = 4 + resolution = mi.ScalarVector3i(1, 1, n) + before = [1.0, 2.0, 3.0, 4.0] + after = [5.0, 6.0, 7.0, 8.0] + + volume = _make_grid_volume(before, n) + extremum = mi.load_dict({"type": "extremum_grid", "resolution": resolution}) + medium = _make_medium("heterogeneous", volume, extremum) + + # Build the ground truth directly from the "after" data. + ref_volume = _make_grid_volume(after, n) + ref_extremum = mi.load_dict({"type": "extremum_grid", "resolution": resolution}) + ref_extremum.update_extremum(ref_volume.bbox(), ref_volume) + + params = mi.traverse(medium) + params["sigma_t.data"] = mi.TensorXf(np.array(after).reshape(n, 1, 1, 1)) + params.update() + + ray = mi.Ray3f(o=[0.5, 0.5, 0], d=[0, 0, 1]) + expected = ref_extremum.sample_test(ray, 0.0, 1.0, target_ot=6.0) + got = medium.extremum().sample_test(ray, 0.0, 1.0, target_ot=6.0) + assert np.allclose(got, expected) + + +def _make_x_volume(values, n, to_world=None): + # Radial/X resolution must be the grid's fastest (last) axis. + data = np.array(values, dtype=float).reshape(1, 1, n) + d = { + "type": "gridvolume", + "grid": mi.VolumeGrid(data), + "filter_type": "nearest", + "accel": False, + } + if to_world is not None: + d["to_world"] = to_world + return mi.load_dict(d) + + +def _make_x_extremum(bbox, volume, n): + extremum = mi.load_dict( + { + "type": "extremum_grid", + "resolution": mi.ScalarVector3i(n, 1, 1), + } + ) + extremum.update_extremum(bbox, volume) + return extremum + + +def test_sample_tight_sampled(variant_scalar_rgb): + # 4 cells of width 0.25 along x, values [1, 2, 3, 4], domain == volume + # bbox. Sampled inside the 4th cell: ot(cell1..3) = 0.25 + 0.5 + 0.75, + # leaving 0.1 of the 1.6 target. + volume = _make_x_volume([1.0, 2.0, 3.0, 4.0], 4) + extremum = _make_x_extremum(volume.bbox(), volume, 4) + + ray = mi.Ray3f(o=[0, 0.5, 0.5], d=[1, 0, 0]) + distance, leftover_ot = extremum.sample_test(ray, 0.0, 1.0, target_ot=1.6) + + assert np.allclose(distance, 0.775) + assert np.allclose(leftover_ot, 0.1) + + +def test_sample_tight_escapes(variant_scalar_rgb): + # Same setup as `test_sample_tight_sampled`, + # total ot = 0.25 * (1 + 2 + 3 + 4) = 2.5, never reached. + volume = _make_x_volume([1.0, 2.0, 3.0, 4.0], 4) + extremum = _make_x_extremum(volume.bbox(), volume, 4) + + ray = mi.Ray3f(o=[0, 0.5, 0.5], d=[1, 0, 0]) + distance, leftover_ot = extremum.sample_test(ray, 0.0, 1.0, target_ot=3.0) + + assert np.isinf(distance) + assert np.allclose(leftover_ot, 0.5) + + +def test_sample_non_tight_sampled(variant_scalar_rgb): + # Domain twice as wide as the volume's own [0, 1] bbox along x: the + # second half is only reachable through the extremum's edge-clamped + # indexing, which is the only supported behavior past the volume's data. + volume = _make_x_volume([1.0, 2.0, 3.0, 4.0], 4) + domain = mi.BoundingBox3f([0, 0, 0], [2, 2, 2]) + extremum = _make_x_extremum(domain, volume, 4) + + ray = mi.Ray3f(o=[0, 0.5, 0.5], d=[1, 0, 0]) + distance, leftover_ot = extremum.sample_test(ray, 0.0, 2.0, target_ot=4.1) + + assert np.allclose(distance, 1.4) + assert np.allclose(leftover_ot, 0.6) + + +def test_sample_rotated_axis_aligned(variant_scalar_rgb): + # Volume rotated 90 degrees about z: its local x-variation now runs + # along world y. Domain is volume.bbox(), so the ray never leaves it. + to_world = mi.ScalarAffineTransform4f.rotate([0, 0, 1], 90) + volume = _make_x_volume([1.0, 2.0, 3.0, 4.0], 4, to_world=to_world) + extremum = _make_x_extremum(volume.bbox(), volume, 4) + + ray = mi.Ray3f(o=[-0.5, 0, 0.5], d=[0, 1, 0]) + distance, leftover_ot = extremum.sample_test(ray, 0.0, 1.0, target_ot=1.6) + + assert np.allclose(distance, 0.775) + assert np.allclose(leftover_ot, 0.1) + + +def test_sample_rotated(variant_scalar_rgb): + # Rotated 45 degrees about z, so the local x-variation runs diagonally + # in world space. The domain is the tightest axis-aligned box around the + # volume's actual footprint, leaving a triangular gap between the box + # and the volume. Start on the domain edge, at local x=-0.5, 2 cells + # before the real data starts at x=0 heading inward. + to_world = mi.ScalarAffineTransform4f.rotate([0, 0, 1], 45) + volume = _make_x_volume([1.0, 2.0, 3.0, 4.0], 4, to_world=to_world) + extremum = _make_x_extremum(volume.bbox(), volume, 4) + + ray = mi.Ray3f( + o=to_world @ mi.ScalarPoint3f(-0.5, 0.5, 0.5), + d=to_world @ mi.ScalarVector3f(1, 0, 0), + ) + distance, leftover_ot = extremum.sample_test(ray, 0.0, 1.5, target_ot=1.6) + + # The 2 gap cells both clip to cell 0 into one 0.75-long segment. + assert np.allclose(distance, 1.0 + 0.35 / 3.0) + assert np.allclose(leftover_ot, 0.35) diff --git a/src/integrators/volpath.cpp b/src/integrators/volpath.cpp index 8840b16856..eec823763b 100644 --- a/src/integrators/volpath.cpp +++ b/src/integrators/volpath.cpp @@ -1,11 +1,15 @@ +#include #include #include #include #include +#include #include #include #include #include +#include +#include NAMESPACE_BEGIN(mitsuba) @@ -69,21 +73,22 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator { public: MI_IMPORT_BASE(MonteCarloIntegrator, m_max_depth, m_rr_depth) MI_IMPORT_TYPES(Scene, Sampler, Emitter, EmitterPtr, BSDF, BSDFPtr, - Medium, MediumPtr, PhaseFunctionContext) + Medium, MediumPtr, PhaseFunctionContext, Extremum) + + using TrackingStateType = TrackingState; VolumetricPathIntegrator(const Properties &props) : Base(props) { } + /// Create seed and offsets that can be used to generate a new PCG32 rng. MI_INLINE - Float index_spectrum(const UnpolarizedSpectrum &spec, const UInt32 &idx) const { - Float m = spec[0]; - if constexpr (is_rgb_v) { // Handle RGB rendering - dr::masked(m, idx == 1u) = spec[1]; - dr::masked(m, idx == 2u) = spec[2]; - } else { - DRJIT_MARK_USED(idx); - } - return m; + std::pair new_seed_offset(Float sample1, Float sample2) const { + UInt32 s0 = UInt32(sample1 * 4294967296.f); // [0,1) -> [0, 2^32) + UInt32 s1 = UInt32(sample2 * 4294967296.f); + UInt64 seed, offset; + seed = sample_tea_64(s0, s1); + offset = sample_tea_64(s1, s0); + return {seed, offset}; } std::pair sample(const Scene *scene, @@ -115,7 +120,6 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator { } SurfaceInteraction3f si = dr::zeros(); - Mask needs_intersection = true; Interaction3f last_scatter_event = dr::zeros(); Float last_scatter_direction_pdf = 1.f; @@ -134,14 +138,13 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator { Float eta; Interaction3f last_scatter_event; Float last_scatter_direction_pdf; - Mask needs_intersection; Mask specular_chain; Mask valid_ray; Sampler* sampler; DRJIT_STRUCT(LoopState, active, depth, ray, throughput, result, \ si, mei, medium, eta, last_scatter_event, \ - last_scatter_direction_pdf, needs_intersection, \ + last_scatter_direction_pdf, \ specular_chain, valid_ray, sampler) } ls = { active, @@ -155,7 +158,6 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator { eta, last_scatter_event, last_scatter_direction_pdf, - needs_intersection, specular_chain, valid_ray, sampler @@ -176,7 +178,6 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator { Float& eta = ls.eta; Interaction3f& last_scatter_event = ls.last_scatter_event; Float& last_scatter_direction_pdf = ls.last_scatter_direction_pdf; - Mask& needs_intersection = ls.needs_intersection; Mask& specular_chain = ls.specular_chain; Mask& valid_ray = ls.valid_ray; Sampler* sampler = ls.sampler; @@ -206,45 +207,132 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator { Mask act_null_scatter = false, act_medium_scatter = false, escaped_medium = false; - // If the medium does not have a spectrally varying extinction, - // we can perform a few optimizations to speed up rendering - Mask is_spectral = active_medium; - Mask not_spectral = false; + dr::masked(si, active) = scene->ray_intersect( + ray, +RayFlags::Default, false, active, ray_mask); + dr::masked(ray.maxt, active) = si.t; + + Float mint, maxt; if (dr::any_or(active_medium)) { - is_spectral &= medium->has_spectral_extinction(); - not_spectral = !is_spectral && active_medium; + std::tie(mei, mint, maxt) = medium->prepare_medium_traversal(ray, active_medium); + escaped_medium = active_medium && !dr::isfinite(maxt); + active_medium &= !escaped_medium; } if (dr::any_or(active_medium)) { - mei = medium->sample_interaction(ray, sampler->next_1d(active_medium), channel, active_medium); - dr::masked(ray.maxt, active_medium && medium->is_homogeneous() && mei.is_valid()) = mei.t; - Mask intersect = needs_intersection && active_medium; - if (dr::any_or(intersect)) - dr::masked(si, intersect) = scene->ray_intersect( - ray, +RayFlags::Default, false, intersect, ray_mask); - needs_intersection &= !active_medium; - - dr::masked(mei.t, active_medium && (si.t < mei.t)) = dr::Infinity; - if (dr::any_or(is_spectral)) { - auto [tr, free_flight_pdf] = medium->transmittance_eval_pdf(mei, si, is_spectral); - Float tr_pdf = index_spectrum(free_flight_pdf, channel); - dr::masked(throughput, is_spectral) *= dr::select(tr_pdf > 0.f, tr / tr_pdf, 0.f); - } - - escaped_medium = active_medium && !mei.is_valid(); + auto extremum = medium->extremum(); + + Float sample1 = sampler->next_1d(active_medium); + Float sample2 = sampler->next_1d(active_medium); + auto [seed, offset] = new_seed_offset(sample1, sample2); + dr::PCG32 rng; + rng.seed(seed, offset); + + Float target_ot = -dr::log(1.f - sampler->next_1d(active_medium)); + + TrackingStateType state { + ray, rng, mei, target_ot, + medium->has_spectral_extinction(), + /*throughput=*/UnpolarizedSpectrum(1.f), + }; + + // Traverse extremum segments and perform delta tracking + state = extremum->traverse_extremum( + ray, mint, maxt, channel, state, + [](const ExtremumSegment &segment, TrackingStateType *state, + const UInt32 &channel, Mask active) { + UnpolarizedSpectrum &throughput = state->throughput; + dr::PCG32 &rng = state->rng; + MediumInteraction3f &mei = state->mei; + MediumPtr medium = mei.medium; + Mask act_spectral = state->has_spectral_extinction && active; + Mask act_not_spectral = !state->has_spectral_extinction && active; + + // Check if the last iteration had a valid interaction + // within the segment + Float mint = dr::select( + mei.is_valid(), dr::maximum(segment.mint, mei.t), + segment.mint); + + Float segment_ot = (segment.maxt - mint) * segment.majorant(); + Mask sampled = (state->target_ot < segment_ot) && active; + Float maxt = segment.maxt; + + if (dr::any_or(sampled)) { + dr::masked(maxt, sampled) = + mint + state->target_ot / + dr::maximum(segment.majorant(), dr::Epsilon); + } + + Float dt = maxt - mint; + + if (dr::any_or(act_spectral)) { + // Accumulate transmittance in the throughput and + // pdf (spectral only). + UnpolarizedSpectrum tr = dr::exp(-dt * segment.majorant()); + Float pdf = index_spectrum( + dr::select(sampled, tr * segment.majorant(), tr), channel); + dr::masked(throughput, act_spectral) *= tr / pdf; + } + + if (dr::any_or(sampled)) { + mei.t = maxt; + mei.p = state->ray(maxt); + + // Retrieve scattering coefficients at position. + UnpolarizedSpectrum sigma_s, sigma_n, sigma_t; + std::tie(sigma_s, std::ignore, sigma_t) = + medium->get_scattering_coefficients(mei, sampled); + sigma_n = segment.majorant() - sigma_t; + + // Sample event type + Float null_scatter_prob = dr::mean(sigma_n / segment.majorant()); + Mask null_scatter = + (rng.template next_float(sampled) < null_scatter_prob) + && sampled; + Mask real_scatter = !null_scatter && sampled; + + // Accumulate throughput and pdf given the event + // type and is_spectral. + if (dr::any_or(null_scatter && act_spectral)) { + dr::masked(throughput, null_scatter && act_spectral) *= + sigma_n / null_scatter_prob; + } + + if (dr::any_or(real_scatter)) { + if (dr::any_or(act_spectral)) { + dr::masked(throughput, real_scatter && act_spectral) *= + sigma_s / (1.0f - null_scatter_prob); + } + + if (dr::any_or(act_not_spectral)) { + dr::masked(throughput, real_scatter && act_not_spectral) *= + sigma_s / sigma_t; + } + + // disable the loop once we encounter a real + // scattering interaction + active &= !real_scatter; + } + + dr::masked(state->target_ot, sampled) = + -dr::log(1.f - state->rng.template next_float(sampled)); + } + + dr::masked(mei.t, !sampled) = dr::Infinity; + dr::masked(state->target_ot, !sampled && active) -= segment_ot; + + Mask step = !sampled; + return std::pair(step, active); + }, active_medium); + + // Update throughput by the transmittance and pdf weight + dr::masked(throughput, active_medium) *= state.throughput; + dr::masked(mei, active_medium) = state.mei; + + escaped_medium |= active_medium && !mei.is_valid(); active_medium &= mei.is_valid(); - // Handle null and real scatter events - Float null_scatter_prob = dr::mean(mei.sigma_n / mei.combined_extinction); - Mask null_scatter = sampler->next_1d(active_medium) < null_scatter_prob; - - act_null_scatter |= null_scatter && active_medium; - act_medium_scatter |= !act_null_scatter && active_medium; - - if (dr::any_or(is_spectral && act_null_scatter)) - dr::masked(throughput, is_spectral && act_null_scatter) *= - mei.sigma_n / null_scatter_prob; - + act_medium_scatter = !escaped_medium && active_medium; dr::masked(depth, act_medium_scatter) += 1; dr::masked(last_scatter_event, act_medium_scatter) = mei; } @@ -253,18 +341,7 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator { active &= depth < (uint32_t) m_max_depth; act_medium_scatter &= active; - if (dr::any_or(act_null_scatter)) { - dr::masked(ray.o, act_null_scatter) = mei.p; - dr::masked(si.t, act_null_scatter) = si.t - mei.t; - } - if (dr::any_or(act_medium_scatter)) { - if (dr::any_or(is_spectral)) - dr::masked(throughput, is_spectral && act_medium_scatter) *= - mei.sigma_s / dr::mean(mei.sigma_t / mei.combined_extinction); - if (dr::any_or(not_spectral)) - dr::masked(throughput, not_spectral && act_medium_scatter) *= mei.sigma_s / mei.sigma_t; - PhaseFunctionContext phase_ctx(sampler); auto phase = mei.medium->phase_function(); @@ -291,17 +368,12 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator { act_medium_scatter &= phase_pdf > 0.f; Ray3f new_ray = mei.spawn_ray(wo); dr::masked(ray, act_medium_scatter) = new_ray; - needs_intersection |= act_medium_scatter; dr::masked(last_scatter_direction_pdf, act_medium_scatter) = phase_pdf; dr::masked(throughput, act_medium_scatter) *= phase_weight; } // --------------------- Surface Interactions --------------------- active_surface |= escaped_medium; - Mask intersect = active_surface && needs_intersection; - if (dr::any_or(intersect)) - dr::masked(si, intersect) = scene->ray_intersect( - ray, +RayFlags::Default, false, intersect, ray_mask); if (dr::any_or(active_surface)) { // ---------------- Intersection with emitters ---------------- @@ -359,7 +431,6 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator { Ray3f bsdf_ray = si.spawn_ray(si.to_world(bs.wo)); dr::masked(ray, active_surface) = bsdf_ray; - needs_intersection |= active_surface; Mask non_null_bsdf = active_surface && !has_flag(bs.sampled_type, BSDFFlags::Null); dr::masked(depth, non_null_bsdf) += 1; @@ -409,32 +480,26 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator { Float total_dist = 0.f; SurfaceInteraction3f si = dr::zeros(); - Mask needs_intersection = true; DirectionSample3f dir_sample = ds; struct LoopState { Mask active; Ray3f ray; Float total_dist; - Mask needs_intersection; MediumPtr medium; SurfaceInteraction3f si; Spectrum transmittance; - DirectionSample3f dir_sample; Sampler* sampler; DRJIT_STRUCT(LoopState, active, ray, total_dist, \ - needs_intersection, medium, si, transmittance, \ - dir_sample, sampler) + medium, si, transmittance, sampler) } ls = { active, ray, total_dist, - needs_intersection, medium, si, transmittance, - dir_sample, sampler }; @@ -445,11 +510,9 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator { Mask& active = ls.active; Ray3f& ray = ls.ray; Float& total_dist = ls.total_dist; - Mask& needs_intersection = ls.needs_intersection; MediumPtr& medium = ls.medium; SurfaceInteraction3f& si = ls.si; Spectrum& transmittance = ls.transmittance; - DirectionSample3f& dir_sample = ls.dir_sample; Sampler* sampler = ls.sampler; Float remaining_dist = max_dist - total_dist; @@ -458,61 +521,102 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator { if (dr::none_or(active)) return; - Mask escaped_medium = false; Mask active_medium = active && (medium != nullptr); - Mask active_surface = active && !active_medium; - if (dr::any_or(active_medium)) { - auto mei = medium->sample_interaction(ray, sampler->next_1d(active_medium), channel, active_medium); - dr::masked(ray.maxt, active_medium && medium->is_homogeneous() && mei.is_valid()) = dr::minimum(mei.t, remaining_dist); - Mask intersect = needs_intersection && active_medium; - if (dr::any_or(intersect)) - dr::masked(si, intersect) = scene->ray_intersect(ray, intersect); - - dr::masked(mei.t, active_medium && (si.t < mei.t)) = dr::Infinity; - needs_intersection &= !active_medium; - - Mask is_spectral = medium->has_spectral_extinction() && active_medium; - Mask not_spectral = !is_spectral && active_medium; - if (dr::any_or(is_spectral)) { - Float t = dr::minimum(remaining_dist, dr::minimum(mei.t, si.t)) - mei.mint; - UnpolarizedSpectrum tr = dr::exp(-t * mei.combined_extinction); - UnpolarizedSpectrum free_flight_pdf = dr::select(si.t < mei.t || mei.t > remaining_dist, tr, tr * mei.combined_extinction); - Float tr_pdf = index_spectrum(free_flight_pdf, channel); - dr::masked(transmittance, is_spectral) *= dr::select(tr_pdf > 0.f, tr / tr_pdf, 0.f); - } - - // Handle exceeding the maximum distance by medium sampling - dr::masked(total_dist, active_medium && (mei.t > remaining_dist) && mei.is_valid()) = dir_sample.dist; - dr::masked(mei.t, active_medium && (mei.t > remaining_dist)) = dr::Infinity; - - escaped_medium = active_medium && !mei.is_valid(); - active_medium &= mei.is_valid(); - is_spectral &= active_medium; - not_spectral &= active_medium; - - dr::masked(total_dist, active_medium) += mei.t; + dr::masked(si, active) = scene->ray_intersect(ray, active); + dr::masked(ray.maxt, active) = dr::minimum(si.t, remaining_dist); + dr::masked(total_dist, active) += ray.maxt; - if (dr::any_or(active_medium)) { - dr::masked(ray.o, active_medium) = mei.p; - dr::masked(si.t, active_medium) = si.t - mei.t; + MediumInteraction3f mei; + Float mint, maxt; + if (dr::any_or(active_medium)) { + std::tie(mei, mint, maxt) = medium->prepare_medium_traversal(ray, active_medium); + active_medium &= dr::isfinite(maxt); + } - if (dr::any_or(is_spectral)) - dr::masked(transmittance, is_spectral) *= mei.sigma_n; - if (dr::any_or(not_spectral)) - dr::masked(transmittance, not_spectral) *= mei.sigma_n / mei.combined_extinction; - } + if (dr::any_or(active_medium)) { + auto extremum = medium->extremum(); + + Float sample1 = sampler->next_1d(active_medium); + Float sample2 = sampler->next_1d(active_medium); + auto [seed, offset] = new_seed_offset(sample1, sample2); + dr::PCG32 rng; + rng.seed(seed, offset); + + Float target_ot = -dr::log(1.f - sampler->next_1d(active_medium)); + + TrackingStateType state { + ray, rng, mei, target_ot, + medium->has_spectral_extinction(), + /*throughput=*/UnpolarizedSpectrum(1.f), + }; + + // Ratio tracking through the medium + state = extremum->traverse_extremum( + ray, mint, maxt, channel, state, + [](const ExtremumSegment &segment, TrackingStateType *state, + const UInt32 &channel, Mask active) { + UnpolarizedSpectrum &throughput = state->throughput; + dr::PCG32 &rng = state->rng; + MediumInteraction3f &mei = state->mei; + MediumPtr medium = mei.medium; + Mask act_spectral = state->has_spectral_extinction && active; + Mask act_not_spectral = !state->has_spectral_extinction && active; + + Float majorant = segment.majorant(); + + Float mint = dr::select( + mei.is_valid(), dr::maximum(segment.mint, mei.t), + segment.mint); + Float segment_ot = (segment.maxt - mint) * majorant; + Mask sampled = (state->target_ot < segment_ot) && active; + Float maxt = segment.maxt; + + if (dr::any_or(sampled)) + dr::masked(maxt, sampled) = + mint + state->target_ot / + dr::maximum(majorant, dr::Epsilon); + + Float dt = maxt - mint; + + if (dr::any_or(act_spectral)) { + UnpolarizedSpectrum tr = dr::exp(-dt * majorant); + Float pdf = index_spectrum( + dr::select(sampled, tr * majorant, tr), channel); + dr::masked(throughput, act_spectral) *= tr / pdf; + } + + if (dr::any_or(sampled)) { + mei.t = maxt; + mei.p = state->ray(maxt); + + UnpolarizedSpectrum sigma_t; + std::tie(std::ignore, std::ignore, sigma_t) = + medium->get_scattering_coefficients(mei, sampled); + UnpolarizedSpectrum sigma_n = majorant - sigma_t; + + if (dr::any_or(act_spectral)) + dr::masked(throughput, sampled && act_spectral) *= sigma_n; + + if (dr::any_or(act_not_spectral)) + dr::masked(throughput, sampled && act_not_spectral) *= + dr::maximum(1.f - sigma_t / majorant, 0.f); + + dr::masked(state->target_ot, sampled) = + -dr::log(1.f - rng.template next_float(active)); + } + dr::masked(mei.t, !sampled) = dr::Infinity; + dr::masked(state->target_ot, !sampled && active) -= segment_ot; + + Mask step = !sampled; + return std::pair(step, active); + }, active_medium); + + dr::masked(transmittance, active_medium) *= state.throughput; } // Handle interactions with surfaces - Mask intersect = active_surface && needs_intersection; - if (dr::any_or(intersect)) - dr::masked(si, intersect) = scene->ray_intersect(ray, intersect); - needs_intersection &= !intersect; - active_surface |= escaped_medium; - dr::masked(total_dist, active_surface) += si.t; - - active_surface &= si.is_valid() && active && !active_medium; + Mask active_surface = si.is_valid() && active; if (dr::any_or(active_surface)) { auto bsdf = si.bsdf(ray); Spectrum bsdf_val = bsdf->eval_null_transmission(si, active_surface); @@ -523,10 +627,9 @@ class VolumetricPathIntegrator : public MonteCarloIntegrator { // Update the ray with new origin & t parameter dr::masked(ray, active_surface) = si.spawn_ray(ray.d); ray.maxt = remaining_dist; - needs_intersection |= active_surface; // Continue tracing through scene if non-zero weights exist - active &= (active_medium || active_surface) && + active &= active_surface && dr::any(unpolarized_spectrum(transmittance) != 0.f); // If a medium transition is taking place: Update the medium pointer diff --git a/src/media/heterogeneous.cpp b/src/media/heterogeneous.cpp index f7d6540a19..382665dad9 100644 --- a/src/media/heterogeneous.cpp +++ b/src/media/heterogeneous.cpp @@ -1,9 +1,12 @@ #include +#include #include #include +#include #include #include #include +#include #include #include #include @@ -149,8 +152,8 @@ template class HeterogeneousMedium final : public Medium { public: MI_IMPORT_BASE(Medium, m_is_homogeneous, m_has_spectral_extinction, - m_phase_function) - MI_IMPORT_TYPES(Scene, Sampler, Texture, Volume) + m_phase_function, m_extremum) + MI_IMPORT_TYPES(Scene, Sampler, Texture, Volume, Extremum) HeterogeneousMedium(const Properties &props) : Base(props) { m_is_homogeneous = false; @@ -161,6 +164,23 @@ class HeterogeneousMedium final : public Medium { m_has_spectral_extinction = props.get("has_spectral_extinction", true); m_max_density = dr::opaque(m_scale * m_sigmat->max()); + + for (auto &prop : props.objects()) { + if (auto *extremum = prop.try_get()) { + if (m_extremum) + Throw("Only a single extremum structure can be specified per medium"); + m_extremum = extremum; + } + } + + if (!m_extremum) { + // Create a default global extremum structure. + m_extremum = + PluginManager::instance()->create_object(Properties("extremum_global")); + } + + m_extremum->update_extremum( + m_sigmat->bbox(), m_sigmat.get(), m_scale); } void traverse(TraversalCallback *cb) override { @@ -170,8 +190,15 @@ class HeterogeneousMedium final : public Medium { Base::traverse(cb); } - void parameters_changed(const std::vector &/*keys*/ = {}) override { + void parameters_changed(const std::vector &keys = {}) override { m_max_density = dr::opaque(m_scale * m_sigmat->max()); + + if (string::contains(keys, "sigma_t")) + m_extremum->update_extremum( + m_sigmat->bbox(), m_sigmat.get(), std::nullopt); + + if (string::contains(keys, "scale")) + m_extremum->set_scale(m_scale); } UnpolarizedSpectrum @@ -205,7 +232,8 @@ class HeterogeneousMedium final : public Medium { oss << "HeterogeneousMedium[" << std::endl << " albedo = " << string::indent(m_albedo) << std::endl << " sigma_t = " << string::indent(m_sigmat) << std::endl - << " scale = " << string::indent(m_scale) << std::endl + << " scale = " << string::indent(m_scale) << "," << std::endl + << " extremum = " << string::indent(m_extremum) << std::endl << "]"; return oss.str(); } diff --git a/src/media/homogeneous.cpp b/src/media/homogeneous.cpp index cbb9a8d61c..05c3d171cc 100644 --- a/src/media/homogeneous.cpp +++ b/src/media/homogeneous.cpp @@ -1,9 +1,12 @@ #include +#include #include #include +#include #include #include #include +#include #include #include #include @@ -131,8 +134,8 @@ However, it supports the use of a spatially varying albedo. template class HomogeneousMedium final : public Medium { public: - MI_IMPORT_BASE(Medium, m_is_homogeneous, m_has_spectral_extinction, m_phase_function) - MI_IMPORT_TYPES(Scene, Sampler, Texture, Volume) + MI_IMPORT_BASE(Medium, m_is_homogeneous, m_has_spectral_extinction, m_phase_function, m_extremum) + MI_IMPORT_TYPES(Scene, Sampler, Texture, Volume, Extremum) HomogeneousMedium(const Properties &props) : Base(props) { m_is_homogeneous = true; @@ -141,6 +144,15 @@ class HomogeneousMedium final : public Medium { m_scale = props.get("scale", 1.0f); m_has_spectral_extinction = props.get("has_spectral_extinction", true); + + // Create a default global extremum structure + m_extremum = + PluginManager::instance()->create_object(Properties("extremum_global")); + + m_extremum->update_extremum( + ScalarBoundingBox3f(-dr::Infinity, + dr::Infinity), + m_sigmat.get(), m_scale); } void traverse(TraversalCallback *cb) override { @@ -150,6 +162,17 @@ class HomogeneousMedium final : public Medium { Base::traverse(cb); } + void parameters_changed(const std::vector &keys = {}) override { + if (string::contains(keys, "sigma_t")) + m_extremum->update_extremum( + ScalarBoundingBox3f(-dr::Infinity, + dr::Infinity), + m_sigmat.get(), std::nullopt); + + if (string::contains(keys, "scale")) + m_extremum->set_scale(m_scale); + } + MI_INLINE auto eval_sigmat(const MediumInteraction3f &mi, Mask active) const { auto sigmat = m_sigmat->eval(mi) * m_scale; if (has_flag(m_phase_function->flags(), PhaseFunctionFlags::Microflake)) diff --git a/src/python/main_v.cpp b/src/python/main_v.cpp index bdfbaa3788..af03216034 100644 --- a/src/python/main_v.cpp +++ b/src/python/main_v.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -49,6 +50,7 @@ static nb::object caster(Object *o) { PY_TRY_CAST(PhaseFunction); PY_TRY_CAST(Medium); + PY_TRY_CAST(Extremum); return nb::object(); } @@ -96,6 +98,9 @@ MI_PY_DECLARE(MediumInteraction); MI_PY_DECLARE(PreliminaryIntersection); MI_PY_DECLARE(DirectedEdge); MI_PY_DECLARE(Medium); +MI_PY_DECLARE(ExtremumSegment); +MI_PY_DECLARE(TrackingState); +MI_PY_DECLARE(Extremum); MI_PY_DECLARE(mueller); MI_PY_DECLARE(MicrofacetDistribution); MI_PY_DECLARE(MicroflakeDistribution); @@ -212,11 +217,14 @@ NB_MODULE(MI_VARIANT_NAME, m) { MI_PY_IMPORT(Shape); MI_PY_IMPORT(DirectedEdge); MI_PY_IMPORT(Medium); + MI_PY_IMPORT(ExtremumSegment); + MI_PY_IMPORT(Extremum); MI_PY_IMPORT(Endpoint); MI_PY_IMPORT(Emitter); MI_PY_IMPORT(Interaction); MI_PY_IMPORT(SurfaceInteraction); MI_PY_IMPORT(MediumInteraction); + MI_PY_IMPORT(TrackingState); MI_PY_IMPORT(PreliminaryIntersection); MI_PY_IMPORT(PositionSample); MI_PY_IMPORT(SilhouetteSample); diff --git a/src/render/CMakeLists.txt b/src/render/CMakeLists.txt index 1ca07dcbc6..7cfe4049a1 100644 --- a/src/render/CMakeLists.txt +++ b/src/render/CMakeLists.txt @@ -76,6 +76,10 @@ add_library(mitsuba-render OBJECT dedge.cpp ${INC_DIR}/dedge.h emitter.cpp ${INC_DIR}/emitter.h endpoint.cpp ${INC_DIR}/endpoint.h + extremum.cpp ${INC_DIR}/extremum.h + ${INC_DIR}/extremum_segment.h + ${INC_DIR}/tracking.h + ${INC_DIR}/volume_utils.h film.cpp ${INC_DIR}/film.h ${INC_DIR}/fresnel.h imageblock.cpp ${INC_DIR}/imageblock.h diff --git a/src/render/extremum.cpp b/src/render/extremum.cpp new file mode 100644 index 0000000000..00ac021fae --- /dev/null +++ b/src/render/extremum.cpp @@ -0,0 +1,48 @@ +#include +#include + +NAMESPACE_BEGIN(mitsuba) + +MI_VARIANT Extremum::Extremum() + : JitObject(""), m_scale(1.f) { +} + +MI_VARIANT Extremum::Extremum(const Properties &props) + : JitObject(props.id()), m_scale(1.f) { +} + +MI_VARIANT Extremum::~Extremum() { +} + +MI_VARIANT void Extremum::update_extremum( + const ScalarBoundingBox3f &bbox, const Volume *volume, + std::optional scale) { + // set scale if provided + if (scale) + set_scale(scale.value()); + // set validity bbox + set_bbox(bbox); + + if (!m_bbox.valid()) + Throw("Extremum::update_extremum() called with an invalid bbox."); + + // rebuild the extremum structure + build(volume); +} + +MI_VARIANT +TrackingState +Extremum::traverse_extremum( + const Ray3f &/*ray*/, + Float /*mint*/, + Float /*maxt*/, + UInt32 /*channel*/, + TrackingStateType /*state*/, + const TrackingFunctionType & /*func*/, + Mask /*active*/ +) const { + NotImplementedError("traverse_extremum"); +} + +MI_INSTANTIATE_CLASS(Extremum) +NAMESPACE_END(mitsuba) diff --git a/src/render/medium.cpp b/src/render/medium.cpp index 9987ba978b..800e01d72a 100644 --- a/src/render/medium.cpp +++ b/src/render/medium.cpp @@ -1,6 +1,8 @@ #include #include #include +#include +#include #include #include #include @@ -35,6 +37,7 @@ MI_VARIANT Medium::~Medium() { } MI_VARIANT void Medium::traverse(TraversalCallback *cb) { cb->put("phase_function", m_phase_function, ParamFlags::Differentiable); + cb->put("extremum", m_extremum, ParamFlags::NonDifferentiable); } MI_VARIANT @@ -95,6 +98,31 @@ Medium::transmittance_eval_pdf(const MediumInteraction3f &mi, return { tr, pdf }; } +MI_VARIANT +std::tuple::MediumInteraction3f, Float, Float> +Medium::prepare_medium_traversal(const Ray3f &ray, Mask active) const { + // Initialize basic medium interaction fields + MediumInteraction3f mei = dr::zeros(); + mei.wi = -ray.d; + mei.sh_frame = Frame3f(mei.wi); + mei.time = ray.time; + mei.wavelengths = ray.wavelengths; + mei.medium = this; + + // Intersect AABB + auto [aabb_its, mint, maxt] = intersect_aabb(ray); + aabb_its &= (dr::isfinite(mint) || dr::isfinite(maxt)); + active &= aabb_its; + dr::masked(mint, !active) = 0.f; + dr::masked(maxt, !active) = dr::Infinity; + + dr::masked(mint, active) = dr::maximum(0.f, mint); + dr::masked(maxt, active) = dr::minimum(ray.maxt, maxt); + mei.mint = mint; + + return {mei, mint, maxt}; +} + MI_IMPLEMENT_TRAVERSE_CB(Medium, Object) MI_INSTANTIATE_CLASS(Medium) NAMESPACE_END(mitsuba) diff --git a/src/render/python/CMakeLists.txt b/src/render/python/CMakeLists.txt index 5ff34eba52..dd9c2dce24 100644 --- a/src/render/python/CMakeLists.txt +++ b/src/render/python/CMakeLists.txt @@ -3,6 +3,8 @@ set(RENDER_PY_V_SRC ${CMAKE_CURRENT_SOURCE_DIR}/dedge_v.cpp ${CMAKE_CURRENT_SOURCE_DIR}/emitter_v.cpp ${CMAKE_CURRENT_SOURCE_DIR}/endpoint_v.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/extremum_v.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/tracking_v.cpp ${CMAKE_CURRENT_SOURCE_DIR}/film_v.cpp ${CMAKE_CURRENT_SOURCE_DIR}/fresnel_v.cpp ${CMAKE_CURRENT_SOURCE_DIR}/imageblock_v.cpp diff --git a/src/render/python/extremum_v.cpp b/src/render/python/extremum_v.cpp new file mode 100644 index 0000000000..3d4915f505 --- /dev/null +++ b/src/render/python/extremum_v.cpp @@ -0,0 +1,166 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +MI_PY_EXPORT(ExtremumSegment) { + MI_PY_IMPORT_TYPES() + + auto es = nb::class_(m, "ExtremumSegment", D(ExtremumSegment)) + .def(nb::init<>()) + .def(nb::init(), "other"_a, "Copy constructor") + .def(nb::init(), + D(ExtremumSegment, ExtremumSegment, 2), + "mint"_a, "maxt"_a, "minorant"_a, "majorant"_a) + .def(nb::init(), + D(ExtremumSegment, ExtremumSegment, 2), + "mint"_a, "maxt"_a, "value"_a) + .def("valid", &ExtremumSegment::valid, D(ExtremumSegment, valid)) + .def("reset", &ExtremumSegment::reset, D(ExtremumSegment, reset)) + .def("zero_", &ExtremumSegment::zero_, "size"_a = 1) + .def("zero_", &ExtremumSegment::zero_, D(ExtremumSegment, zero)) + .def("minorant", &ExtremumSegment::minorant, D(ExtremumSegment, minorant)) + .def("majorant", &ExtremumSegment::majorant, D(ExtremumSegment, majorant)) + .def_field(ExtremumSegment, mint, D(ExtremumSegment, mint)) + .def_field(ExtremumSegment, maxt, D(ExtremumSegment, maxt)) + .def_field(ExtremumSegment, value, D(ExtremumSegment, value)) + .def_repr(ExtremumSegment); + + MI_PY_DRJIT_STRUCT(es, ExtremumSegment, mint, maxt, value); +} + +/// Trampoline for derived types implemented in Python +MI_VARIANT class PyExtremum : public Extremum { +public: + MI_IMPORT_TYPES(Extremum, Volume) + NB_TRAMPOLINE(Extremum); + + PyExtremum(const Properties &props) : Extremum(props) {} + + void build(const Volume * volume) override { + NB_OVERRIDE_PURE(build, volume); + } + + std::string to_string() const override { + NB_OVERRIDE(to_string); + } + + void traverse(TraversalCallback *cb) override { + NB_OVERRIDE(traverse, cb); + } + + void parameters_changed(const std::vector &keys) override { + NB_OVERRIDE(parameters_changed, keys); + } +}; + +template void bind_extremum_generic(Cls &cls) { + MI_PY_IMPORT_TYPES(Extremum, Medium) + using TrackingStateType = TrackingState; + using TrackingFunctionType = TrackingFunction; + + cls.def("traverse_extremum", + [](Ptr ptr, const Ray3f &ray, Float mint, Float maxt, + UInt32 channel, TrackingStateType state, + const TrackingFunctionType &func, Mask active) { + return ptr->traverse_extremum(ray, mint, maxt, channel, + state, func, active); + }, + "ray"_a, "mint"_a, "maxt"_a, "channel"_a, "state"_a, "func"_a, + "active"_a = true, + D(Extremum, traverse_extremum)); + + + // Test utility: deterministic delta tracking driven by a fixed target + // optical thickness. + cls.def("sample_test", + [](Ptr ptr, const Ray3f &ray, Float mint, Float maxt, + Float target_ot, UInt32 channel, Mask active) { + using TrackingStateType = TrackingState; + + TrackingStateType state = dr::zeros(); + state.ray = ray; + state.target_ot = target_ot; + state.mei = dr::zeros(); + + state = ptr->traverse_extremum( + ray, mint, maxt, channel, state, + [](const ExtremumSegment &segment, TrackingStateType *state, + const UInt32 &, Mask active) { + Float mint = dr::select( + state->mei.is_valid(), + dr::maximum(segment.mint, state->mei.t), + segment.mint); + Float segment_ot = + (segment.maxt - mint) * segment.majorant(); + Mask sampled = (state->target_ot < segment_ot) && active; + + Float maxt = dr::select( + sampled, + mint + state->target_ot / + dr::maximum(segment.majorant(), + dr::Epsilon), + segment.maxt); + + dr::masked(state->mei.t, sampled) = maxt; + dr::masked(state->mei.t, !sampled) = dr::Infinity; + dr::masked(state->target_ot, !sampled && active) -= + segment_ot; + + return std::pair(/*advance=*/!sampled, + active && !sampled); + }, + active); + + return std::make_tuple(state.mei.t, state.target_ot); + }, + "ray"_a, "mint"_a, "maxt"_a, "target_ot"_a, "channel"_a = 0u, + "active"_a = true, + "Deterministic delta-tracking test utility. Traverses the extremum " + "structure's segments, until an interaction is sampled based on " + "`target_ot`. Returns (distance, leftover_ot); `distance` is " + "infinite if `target_ot` is not reached before `maxt`."); +} + + +MI_PY_EXPORT(Extremum) { + MI_PY_IMPORT_TYPES(Extremum, ExtremumPtr) + using PyExtremum = PyExtremum; + using Properties = mitsuba::Properties; + + auto extremum = MI_PY_TRAMPOLINE_CLASS(PyExtremum, Extremum, Object) + .def(nb::init(), "props"_a) + .def("__repr__", &Extremum::to_string) + .def("set_bbox", &Extremum::set_bbox, + "bbox"_a, D(Extremum, set_bbox)) + .def("set_scale", &Extremum::set_scale, + "scale"_a, D(Extremum, set_scale)) + .def("update_extremum", &Extremum::update_extremum, + "bbox"_a, "volume"_a, "scale"_a = nb::none(), + D(Extremum, update_extremum)) + .def("build", &Extremum::build, + "volume"_a, D(Extremum, build)) + .def("bbox", &Extremum::bbox, D(Extremum, bbox)); + + drjit::bind_traverse(extremum); + + bind_extremum_generic(extremum); + + if constexpr (dr::is_array_v) { + dr::ArrayBinding b; + auto extremum_ptr = dr::bind_array_t(b, m, "ExtremumPtr"); + bind_extremum_generic(extremum_ptr); + extremum_ptr.freeze(); + } +} diff --git a/src/render/python/medium_v.cpp b/src/render/python/medium_v.cpp index 957c6022ba..0ac1ace9a1 100644 --- a/src/render/python/medium_v.cpp +++ b/src/render/python/medium_v.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -50,9 +51,10 @@ MI_VARIANT class PyMedium : public Medium { }; template void bind_medium_generic(Cls &cls) { - MI_PY_IMPORT_TYPES(PhaseFunctionContext) + MI_PY_IMPORT_TYPES(PhaseFunctionContext, Extremum, ExtremumPtr) using RetPhaseFunction = std::conditional_t, PhaseFunctionPtr, drjit::scalar_t>; + using RetExtremum = std::conditional_t, ExtremumPtr, drjit::scalar_t>; cls.def("phase_function", [](Ptr ptr) -> RetPhaseFunction { return ptr->phase_function(); }, @@ -91,7 +93,10 @@ template void bind_medium_generic(Cls &cls) { [](Ptr ptr, const MediumInteraction3f &mi, Mask active = true) { return ptr->get_scattering_coefficients(mi, active); }, "mi"_a, "active"_a=true, - D(Medium, get_scattering_coefficients)); + D(Medium, get_scattering_coefficients)) + .def("extremum", + [](Ptr ptr) -> RetExtremum { return ptr->extremum(); }, + D(Medium, extremum)); } MI_PY_EXPORT(Medium) { diff --git a/src/render/python/tracking_v.cpp b/src/render/python/tracking_v.cpp new file mode 100644 index 0000000000..198828730c --- /dev/null +++ b/src/render/python/tracking_v.cpp @@ -0,0 +1,21 @@ +#include +#include +#include + +MI_PY_EXPORT(TrackingState) { + MI_PY_IMPORT_TYPES() + using TrackingState = mitsuba::TrackingState; + + auto ts = nb::class_(m, "TrackingState", D(TrackingState)) + .def(nb::init<>()) + .def(nb::init(), "other"_a, "Copy constructor") + .def_field(TrackingState, ray, D(TrackingState, ray)) + .def_field(TrackingState, rng, D(TrackingState, rng)) + .def_field(TrackingState, mei, D(TrackingState, mei)) + .def_field(TrackingState, target_ot, D(TrackingState, target_ot)) + .def_field(TrackingState, has_spectral_extinction, D(TrackingState, has_spectral_extinction)) + .def_field(TrackingState, throughput, D(TrackingState, throughput)); + + MI_PY_DRJIT_STRUCT(ts, TrackingState, ray, rng, mei, target_ot, + has_spectral_extinction, throughput); +} diff --git a/src/render/volume.cpp b/src/render/volume.cpp index d6973887a8..a887c86731 100644 --- a/src/render/volume.cpp +++ b/src/render/volume.cpp @@ -57,6 +57,12 @@ Volume::max_per_channel(ScalarFloat * /*out*/) const { NotImplementedError("max_per_channel"); } +MI_VARIANT std::pair +Volume::extremum(BoundingBox3f /*bbox*/) const { + // placeholder minorant + return { 0.f, max() }; +} + MI_VARIANT typename Volume::ScalarVector3i Volume::resolution() const { return ScalarVector3i(1, 1, 1); diff --git a/src/volumes/grid.cpp b/src/volumes/grid.cpp index 98b1869d08..4480566667 100644 --- a/src/volumes/grid.cpp +++ b/src/volumes/grid.cpp @@ -427,6 +427,116 @@ class GridVolume final : public Volume { out[i] = m_max_per_channel[i]; } + std::pair + extremum(BoundingBox3f bbox) const override { + + if (m_texture.shape()[3] != 1) + NotImplementedError("extremum() only supported for single-channel volumes"); + + bbox.clip(BoundingBox3f(Point3f(0.f), Point3f(1.f))); + + // early exit in scalar mode + if (dr::any_or(!bbox.valid())) + return { 0.f, 0.f }; + + Mask active = bbox.valid(); + const Vector3i res = resolution(); + + // Convert to voxel indices with proper padding for interpolation + int32_t padding = + (m_texture.filter_mode() == dr::FilterMode::Linear) ? 1 : 0; + + Vector3i voxel_min = + dr::maximum(dr::floor(bbox.min * Vector3f(res)) - Vector3i(padding), + Vector3i(0)); + Vector3i voxel_max = dr::minimum( + dr::floor(bbox.max * Vector3f(res)) + Vector3i(padding), res - 1); + + UInt32 n = dr::prod((voxel_max - voxel_min) + 1); + Vector3i range = (voxel_max - voxel_min) + 1; + + // Scan voxels in bounds and find min/max + Float max_val = -dr::Infinity; + Float min_val = dr::Infinity; + + if constexpr ( !dr::is_jit_v){ + // If possible use pinned data to avoid ref count issues. + const ScalarFloat *data = m_pinned_data + ? m_pinned_data + : m_texture.tensor().data(); + + for (int32_t z = voxel_min.z(); z <= voxel_max.z(); ++z) { + for (int32_t y = voxel_min.y(); y <= voxel_max.y(); ++y) { + for (int32_t x = voxel_min.x(); x <= voxel_max.x(); ++x) { + size_t idx = ( x + + y * res.x() + + z * res.x() * res.y() ); + ScalarFloat val = data[idx]; + max_val = dr::maximum(max_val, val); + min_val = dr::minimum(min_val, val); + } + } + } + + } else { + struct LoopState { + UInt32 x; + UInt32 y; + UInt32 z; + Float min_val; + Float max_val; + Mask active; + + DRJIT_STRUCT(LoopState, x, y, z, min_val, max_val, active) + } ls = { + UInt32(0), UInt32(0), UInt32(0), + min_val, + max_val, + active + }; + + auto array = m_texture.tensor().array(); + dr::tie(ls) = dr::while_loop( + dr::make_tuple(ls), + [](const LoopState &ls) { return ls.active; }, + [array, res, n, voxel_min, range](LoopState &ls) { + Float &min_val = ls.min_val; + Float &max_val = ls.max_val; + Mask &active = ls.active; + + // volume wide indices + UInt32 x = voxel_min.x() + ls.x; + UInt32 y = voxel_min.y() + ls.y; + UInt32 z = voxel_min.z() + ls.z; + + // serial index + UInt32 tex_idx = x + y * res.x() + z * res.x() * res.y(); + + Float val = dr::gather(array, tex_idx, active); + dr::masked(max_val, active) = dr::maximum(max_val, val); + dr::masked(min_val, active) = dr::minimum(min_val, val); + + // This approach avoids modulo and division which are + // detrimental to performance. + ls.x += 1; + + Mask carry_x = ls.x >= range.x(); + ls.x = dr::select(carry_x, UInt32(0), ls.x); + ls.y += dr::select(carry_x, UInt32(1), UInt32(0)); + + Mask carry_y = carry_x && (ls.y >= range.y()); + ls.y = dr::select(carry_y, UInt32(0), ls.y); + ls.z += dr::select(carry_y, UInt32(1), UInt32(0)); + + ls.active &= ls.z < range.z(); + }); + max_val = ls.max_val; + min_val = ls.min_val; + } + + return { min_val, max_val }; + } + ScalarVector3i resolution() const override { const size_t *shape = m_texture.shape(); return { (int) shape[2], (int) shape[1], (int) shape[0] }; @@ -589,6 +699,14 @@ class GridVolume final : public Volume { out[i] = result.entry(i); } + void pin_ref_count() const override { + m_pinned_data = m_texture.tensor().data(); + } + + void unpin_ref_count() const override { + m_pinned_data = nullptr; + } + protected: Texture3f m_texture; bool m_accel; @@ -596,6 +714,7 @@ class GridVolume final : public Volume { bool m_fixed_max = false; ScalarFloat m_max; std::vector m_max_per_channel; + mutable const ScalarFloat* m_pinned_data = nullptr; MI_TRAVERSE_CB(Base, m_texture) };