diff --git a/resources/data b/resources/data index cba1bf767e..c7537606a7 160000 --- a/resources/data +++ b/resources/data @@ -1 +1 @@ -Subproject commit cba1bf767efb6df6d7fcfc9e19270884d0386dd4 +Subproject commit c7537606a7384806b3ae0380b2e956006c0f4361 diff --git a/src/sensors/CMakeLists.txt b/src/sensors/CMakeLists.txt index d247196aef..7215934e1b 100644 --- a/src/sensors/CMakeLists.txt +++ b/src/sensors/CMakeLists.txt @@ -5,3 +5,4 @@ add_plugin(radiancemeter radiancemeter.cpp) add_plugin(thinlens thinlens.cpp) add_plugin(irradiancemeter irradiancemeter.cpp) add_plugin(distant distant.cpp) +add_plugin(hdistant hdistant.cpp) diff --git a/src/sensors/distant.cpp b/src/sensors/distant.cpp index e42cda1408..cd41cc3b4b 100644 --- a/src/sensors/distant.cpp +++ b/src/sensors/distant.cpp @@ -30,8 +30,7 @@ Distant radiancemeter sensor (:monosp:`distant`) - Sensor-to-world transformation matrix. * - direction - |vector| - - Alternative (and exclusive) to ``to_world``. Direction orienting the - sensor's reference hemisphere. + - Alternative (and exclusive) to ``to_world``. Direction orienting the sensor. * - target - |point| or nested :paramtype:`shape` plugin - *Optional.* Define the ray target sampling strategy. diff --git a/src/sensors/hdistant.cpp b/src/sensors/hdistant.cpp new file mode 100644 index 0000000000..f78ef4c62f --- /dev/null +++ b/src/sensors/hdistant.cpp @@ -0,0 +1,346 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +NAMESPACE_BEGIN(mitsuba) + +enum class RayTargetType { Shape, Point, None }; + +// Forward declaration of specialized HemisphericalDistantSensor +template +class HemisphericalDistantSensorImpl; + +/**! + +.. _sensor-hdistant: + +Hemispherical distant radiancemeter sensor (:monosp:`hdistant`) +--------------------------------------------------------------- + +.. pluginparameters:: + + * - to_world + - |transform| + - Sensor-to-world transformation matrix. + * - target + - |point| or nested :paramtype:`shape` plugin + - *Optional.* Define the ray target sampling strategy. + If this parameter is unset, ray target points are sampled uniformly on + the cross section of the scene's bounding sphere. + If a |point| is passed, rays will target it. + If a shape plugin is passed, ray target points will be sampled from its + surface. + +This sensor plugin implements a distant directional sensor which records +radiation leaving the scene. It records the spectral radiance leaving the scene +in directions covering a hemisphere defined by its ``to_world`` parameter and +mapped to film coordinates. To some extent, it can be seen as the adjoint to +the ``envmap`` emitter. + +The ``to_world`` transform is best set using a +:py:meth:`~mitsuba.core.Transform4f.look_at`. The default orientation covers a +hemisphere defined by the [0, 0, 1] direction, and the ``up`` film direction is +set to [0, 1, 0]. + +The following XML snippet creates a scene with a ``roughconductor`` +surface illuminated by three ``directional`` emitter, each emitting in +a single RGB channel. A ``hdistant`` plugin with default orientation is +defined. + +.. subfigstart:: +.. subfigure:: ../../resources/data/docs/images/sensor/sensor_hdistant_illumination_optimized.svg + :caption: Example scene illumination setup +.. subfigend:: + :label: fig-hdistant-illumination + +.. code:: xml + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +The following figures show the recorded exitant radiance with the default film +orientation (left, ``up = [0,1,0]``) and with a rotated film (right, +``up = [1,1,0]``). Colored dots on the plots materialize emitter directions. +The orange arrow represents the ``up`` direction on the film. +Note that on the plots, the origin of pixel coordinates is taken at the bottom +left. + +.. subfigstart:: +.. subfigure:: ../../resources/data/docs/images/sensor/sensor_hdistant_film_default_optimized.svg + :caption: Default film orientation +.. subfigure:: ../../resources/data/docs/images/sensor/sensor_hdistant_film_rotated_optimized.svg + :caption: Rotated film +.. subfigure:: ../../resources/data/docs/images/sensor/sensor_hdistant_default.svg + :caption: Exitant radiance +.. subfigure:: ../../resources/data/docs/images/sensor/sensor_hdistant_rotated.svg + :caption: Exitant radiance +.. subfigend:: + :label: fig-hdistant-film + +By default, ray target points are sampled from the cross section of the scene's +bounding sphere. The ``target`` parameter can be set to restrict ray target +sampling to a specific subregion of the scene. The recorded radiance is averaged +over the targeted geometry. + +Ray origins are positioned outside of the scene's geometry, such that it is +as if the sensor would be located at an infinite distance from the scene. + +By default, ray target points are sampled from the cross section of the scene's +bounding sphere. The ``target`` parameter should be set to restrict ray target +sampling to a specific subregion of the scene using a flat surface. The recorded +radiance is averaged over the targeted geometry. + +.. warning:: + + * While setting ``target`` using any shape plugin is possible, only specific + configurations will produce meaningful results. This is due to ray sampling + method: when ``target`` is a shape, a point is sampled at its surface, + then shifted along the ``-direction`` vector by the diameter of the scene's + bounding sphere, effectively positioning the ray origin outside of the + geometry. The ray's weight is set to :math:`\frac{1}{A \, p}`, where + :math:`A` is the shape's surface area and :math:`p` is the shape's position + sampling PDF value. This weight definition is irrelevant when the sampled + origin may corresponds to multiple points on the shape, *i.e.* when the + sampled ray can intersect the target shape multiple times. From this + follows that only flat surfaces should be used to set the ``target`` + parameter. Typically, one will rather use a ``rectangle`` or ``disk`` + shape. + * If this sensor is used with a targeting strategy leading to rays not + hitting the scene's geometry (*e.g.* default targeting strategy), it will + pick up ambient emitter radiance samples (or zero values if no ambient + emitter is defined). Therefore, it is almost always preferrable to use a + nondefault targeting strategy. +*/ + +template +class HemisphericalDistantSensor final : public Sensor { +public: + MTS_IMPORT_BASE(Sensor, m_to_world, m_film) + MTS_IMPORT_TYPES(Scene, Shape) + + HemisphericalDistantSensor(const Properties &props) : Base(props) { + // Check reconstruction filter radius + if (m_film->reconstruction_filter()->radius() > + 0.5f + math::RayEpsilon) { + Log(Warn, "This sensor is best used with a reconstruction filter " + "with a radius of 0.5 or lower (e.g. default box)"); + } + + // Store film sample location spacing for performance + m_d.x() = 1.0f / m_film->size().x(); + m_d.y() = 1.0f / m_film->size().y(); + + // Set ray target if relevant + // Get target + if (props.has_property("target")) { + if (props.type("target") == Properties::Type::Array3f) { + m_target_type = RayTargetType::Point; + m_target_point = props.get("target"); + } else if (props.type("target") == Properties::Type::Object) { + // We assume it's a shape + m_target_type = RayTargetType::Shape; + auto obj = props.object("target"); + m_target_shape = dynamic_cast(obj.get()); + + if (!m_target_shape) + Throw("Invalid parameter target, must be a Point3f or a " + "Shape."); + } else { + Throw("Unsupported 'target' parameter type"); + } + } else { + m_target_type = RayTargetType::None; + Log(Debug, "No target specified."); + } + } + + void set_scene(const Scene *scene) override { + m_bsphere = scene->bbox().bounding_sphere(); + m_bsphere.radius = + ek::max(math::RayEpsilon, + m_bsphere.radius * (1.f + math::RayEpsilon) ); + } + + std::pair sample_ray(Float time, Float wavelength_sample, + const Point2f &film_sample, + const Point2f &aperture_sample, + Mask active) const override { + MTS_MASKED_FUNCTION(ProfilerPhase::EndpointSampleRay, active); + + Ray3f ray; + ray.time = time; + + // Sample spectrum + auto [wavelengths, wav_weight] = + sample_wavelength(wavelength_sample); + ray.wavelengths = wavelengths; + + // Sample ray origin + Spectrum ray_weight = 0.f; + + // Sample ray direction + ray.d = -m_to_world.value().transform_affine( + warp::square_to_uniform_hemisphere(film_sample)); + + // Sample target point and position ray origin + if (m_target_type == RayTargetType::Point) { + ray.o = m_target_point - 2.f * ray.d * m_bsphere.radius; + ray_weight = wav_weight; + } else if (m_target_type == RayTargetType::Shape) { + // Use area-based sampling of shape + PositionSample3f ps = + m_target_shape->sample_position(time, aperture_sample, active); + ray.o = ps.p - 2.f * ray.d * m_bsphere.radius; + ray_weight = wav_weight / (ps.pdf * m_target_shape->surface_area()); + } else { // if (m_target_type == RayTargetType::None) { + // Sample target uniformly on bounding sphere cross section + Point2f offset = + warp::square_to_uniform_disk_concentric(aperture_sample); + Vector3f perp_offset = m_to_world.value().transform_affine( + Vector3f(offset.x(), offset.y(), 0.f)); + ray.o = m_bsphere.center + perp_offset * m_bsphere.radius - + ray.d * m_bsphere.radius; + ray_weight = wav_weight; + } + + return { ray, ray_weight & active }; + } + + std::pair sample_ray_differential( + Float time, Float wavelength_sample, const Point2f &film_sample, + const Point2f &aperture_sample, Mask active) const override { + MTS_MASKED_FUNCTION(ProfilerPhase::EndpointSampleRay, active); + + RayDifferential3f ray; + ray.has_differentials = true; + ray.time = time; + + // Sample spectrum + auto [wavelengths, wav_weight] = + sample_wavelength(wavelength_sample); + ray.wavelengths = wavelengths; + + // Sample ray origin + Spectrum ray_weight = 0.f; + + // Sample ray direction + ray.d = -m_to_world.value().transform_affine( + warp::square_to_uniform_hemisphere(film_sample)); + ray.d_x = -m_to_world.value().transform_affine( + warp::square_to_uniform_hemisphere( + Point2f{ film_sample.x() + m_d.x(), film_sample.y() })); + ray.d_y = -m_to_world.value().transform_affine( + warp::square_to_uniform_hemisphere( + Point2f{ film_sample.x(), film_sample.y() + m_d.y() })); + + // Sample target point and position ray origin + if (m_target_type == RayTargetType::Point) { + ray.o = m_target_point - 2.f * ray.d * m_bsphere.radius; + ray.o_x = m_target_point - 2.f * ray.d_x * m_bsphere.radius; + ray.o_y = m_target_point - 2.f * ray.d_y * m_bsphere.radius; + ray_weight = wav_weight; + } else if (m_target_type == RayTargetType::Shape) { + // Use area-based sampling of shape + PositionSample3f ps = + m_target_shape->sample_position(time, aperture_sample, active); + ray.o = ps.p - 2.f * ray.d * m_bsphere.radius; + ray.o_x = ps.p - 2.f * ray.d_x * m_bsphere.radius; + ray.o_y = ps.p - 2.f * ray.d_y * m_bsphere.radius; + ray_weight = wav_weight / (ps.pdf * m_target_shape->surface_area()); + } else { // if (m_target_type == RayTargetType::None) { + // Sample target uniformly on bounding sphere cross section + Point2f offset = + warp::square_to_uniform_disk_concentric(aperture_sample); + Vector3f perp_offset = m_to_world.value().transform_affine( + Vector3f(offset.x(), offset.y(), 0.f)); + ray.o = m_bsphere.center + perp_offset * m_bsphere.radius - + ray.d * m_bsphere.radius; + ray.o_x = m_bsphere.center + perp_offset * m_bsphere.radius - + ray.d_x * m_bsphere.radius; + ray.o_y = m_bsphere.center + perp_offset * m_bsphere.radius - + ray.d_y * m_bsphere.radius; + ray_weight = wav_weight; + } + + return { ray, ray_weight & active }; + } + + // This sensor does not occupy any particular region of space, return an + // invalid bounding box + ScalarBoundingBox3f bbox() const override { return ScalarBoundingBox3f(); } + + std::string to_string() const override { + std::ostringstream oss; + oss << "HemisphericalDistantSensor[" << std::endl + << " to_world = " << string::indent(m_to_world) << "," << std::endl + << " film = " << string::indent(m_film) << "," << std::endl; + + if (m_target_type == RayTargetType::Point) + oss << " target = " << m_target_point << std::endl; + else if (m_target_type == RayTargetType::Shape) + oss << " target = " << string::indent(m_target_shape) << std::endl; + else // if (m_target_type == RayTargetType::None) + oss << " target = None" << std::endl; + + oss << "]"; + + return oss.str(); + } + + MTS_DECLARE_CLASS() + +protected: + // Scene bounding sphere + ScalarBoundingSphere3f m_bsphere; + // Ray target type + RayTargetType m_target_type; + // Target shape if any + ref m_target_shape; + // Target point if any + Point3f m_target_point; + // Spacing between two pixels in film coordinates + ScalarPoint2f m_d; +}; + +MTS_IMPLEMENT_CLASS_VARIANT(HemisphericalDistantSensor, Sensor) +MTS_EXPORT_PLUGIN(HemisphericalDistantSensor, "HemisphericalDistantSensor") +NAMESPACE_END(mitsuba) diff --git a/src/sensors/tests/test_hdistant.py b/src/sensors/tests/test_hdistant.py new file mode 100644 index 0000000000..0d9bd82ed0 --- /dev/null +++ b/src/sensors/tests/test_hdistant.py @@ -0,0 +1,359 @@ +import enoki as ek +import numpy as np +import pytest + + +def sensor_dict(target=None, to_world=None): + result = {"type": "hdistant"} + + if to_world is not None: + result["to_world"] = to_world + + if target == "point": + result.update({"target": [0, 0, 0]}) + + elif target == "shape": + result.update({"target": {"type": "rectangle"}}) + + elif isinstance(target, dict): + result.update({"target": target}) + + return result + + +def make_sensor(d): + from mitsuba.core.xml import load_dict + + return load_dict(d) + + +def test_construct(variant_scalar_rgb): + from mitsuba.core import ScalarTransform4f + + # Construct without parameters + sensor = make_sensor({"type": "hdistant"}) + assert sensor is not None + assert not sensor.bbox().valid() # Degenerate bounding box + + # Construct with transform + sensor = make_sensor( + sensor_dict( + to_world=ScalarTransform4f.look_at( + origin=[0, 0, 0], target=[0, 0, 1], up=[1, 0, 0] + ) + ) + ) + + # Test different target values + # -- No target, + sensor = make_sensor(sensor_dict()) + assert sensor is not None + + # -- Point target + sensor = make_sensor(sensor_dict(target="point")) + assert sensor is not None + + # -- Shape target + sensor = make_sensor(sensor_dict(target="shape")) + assert sensor is not None + + # -- Random object target (we expect to raise) + with pytest.raises(RuntimeError): + make_sensor(sensor_dict(target={"type": "constant"})) + + +def test_sample_ray_direction(variant_scalar_rgb): + sensor = make_sensor(sensor_dict()) + + # Check that directions are appropriately set + for (sample1, sample2, expected) in [ + [[0.5, 0.5], [0.16, 0.44], [0, 0, -1]], + [[0.0, 0.0], [0.23, 0.40], [0.707107, 0.707107, 0]], + [[1.0, 0.0], [0.22, 0.81], [-0.707107, 0.707107, 0]], + [[0.0, 1.0], [0.99, 0.42], [0.707107, -0.707107, 0]], + [[1.0, 1.0], [0.52, 0.31], [-0.707107, -0.707107, 0]], + ]: + ray, _ = sensor.sample_ray(1.0, 1.0, sample1, sample2, True) + + # Check that ray direction is what is expected + assert ek.allclose(ray.d, expected, atol=1e-7) + + +@pytest.mark.parametrize( + "sensor_setup", + [ + "default", + "target_square", + "target_square_small", + "target_square_large", + "target_disk", + "target_point", + ], +) +@pytest.mark.parametrize("w_e", [[0, 0, -1], [0, 1, -1]]) +def test_sample_target(variant_scalar_rgb, sensor_setup, w_e): + # Check if targeting works as intended by rendering a basic scene + from mitsuba.core import Bitmap, ScalarTransform4f, Struct + from mitsuba.core.xml import load_dict + + # Basic illumination and sensing parameters + l_e = 1.0 # Emitted radiance + w_e = list(w_e / np.linalg.norm(w_e)) # Emitter direction + cos_theta_e = abs(np.dot(w_e, [0, 0, 1])) + + # Reflecting surface specification + surface_scale = 1.0 + rho = 1.0 # Surface reflectance + + # Sensor definitions + sensors = { + "default": { # No target, origin projected to bounding sphere + "type": "hdistant", + "sampler": { + "type": "independent", + "sample_count": 100000, + }, + "film": { + "type": "hdrfilm", + "height": 4, + "width": 4, + "rfilter": {"type": "box"}, + }, + }, + "target_square": { # Targeting square, origin projected to bounding sphere + "type": "hdistant", + "target": { + "type": "rectangle", + "to_world": ScalarTransform4f.scale(surface_scale), + }, + "sampler": { + "type": "independent", + "sample_count": 1000, + }, + "film": { + "type": "hdrfilm", + "height": 16, + "width": 16, + "rfilter": {"type": "box"}, + }, + }, + "target_square_small": { # Targeting small square, origin projected to bounding sphere + "type": "hdistant", + "target": { + "type": "rectangle", + "to_world": ScalarTransform4f.scale(0.5 * surface_scale), + }, + "sampler": { + "type": "independent", + "sample_count": 1000, + }, + "film": { + "type": "hdrfilm", + "height": 16, + "width": 16, + "rfilter": {"type": "box"}, + }, + }, + "target_square_large": { # Targeting large square, origin projected to bounding sphere + "type": "hdistant", + "target": { + "type": "rectangle", + "to_world": ScalarTransform4f.scale(2.0 * surface_scale), + }, + "sampler": { + "type": "independent", + "sample_count": 1000000, + }, + "film": { + "type": "hdrfilm", + "height": 4, + "width": 4, + "rfilter": {"type": "box"}, + }, + }, + "target_point": { # Targeting point, origin projected to bounding sphere + "type": "hdistant", + "target": [0, 0, 0], + "sampler": { + "type": "independent", + "sample_count": 1000, + }, + "film": { + "type": "hdrfilm", + "height": 16, + "width": 16, + "rfilter": {"type": "box"}, + }, + }, + "target_disk": { # Targeting disk, origin projected to bounding sphere + "type": "hdistant", + "target": { + "type": "disk", + "to_world": ScalarTransform4f.scale(surface_scale), + }, + "sampler": { + "type": "independent", + "sample_count": 1000, + }, + "film": { + "type": "hdrfilm", + "height": 16, + "width": 16, + "rfilter": {"type": "box"}, + }, + }, + } + + # Scene setup + scene_dict = { + "type": "scene", + "shape": { + "type": "rectangle", + "to_world": ScalarTransform4f.scale(surface_scale), + "bsdf": { + "type": "diffuse", + "reflectance": rho, + }, + }, + "emitter": {"type": "directional", "direction": w_e, "irradiance": l_e}, + "integrator": {"type": "path"}, + } + + scene = load_dict({**scene_dict, "sensor": sensors[sensor_setup]}) + + # Run simulation + scene.render() + result = np.array( + scene.sensors()[0].film().bitmap().convert( + Bitmap.PixelFormat.RGB, Struct.Type.Float32, False + ) + ).squeeze() + + l_o = l_e * cos_theta_e * rho / np.pi # Outgoing radiance + expected = { # Special expected values for some cases + "default": l_o * 2.0 / ek.Pi, + "target_square_large": l_o * 0.25, + } + expected_value = expected.get(sensor_setup, l_o) + + rtol = { # Special tolerance values for some cases + "default": 1e-2, + "target_square_large": 1e-2, + } + rtol_value = rtol.get(sensor_setup, 5e-3) + + assert np.allclose(expected_value, result, rtol=rtol_value) + + +@pytest.mark.parametrize("target", ("point", "shape")) +def test_sample_ray_differential(variant_scalar_rgb, target): + from mitsuba.core.warp import uniform_hemisphere_to_square + + # We set odd and even values on purpose + n_x = 5 + n_y = 6 + + d = sensor_dict(target=target) + d.update( + { + "film": { + "type": "hdrfilm", + "width": n_x, + "height": n_y, + "rfilter": {"type": "box"}, + } + } + ) + sensor = make_sensor(d) + + sample1 = [0.5, 0.5] + sample2 = [0.5, 0.5] + ray, _ = sensor.sample_ray_differential(1.0, 1.0, sample1, sample2, True) + + # Sampled ray differential directions are expected to map to film + # coordinates shifted by one pixel + sample_dx = [0.5 + 1.0 / n_x, 0.5] + expected_ray_dx, _ = sensor.sample_ray(1.0, 1.0, sample_dx, sample2, True) + assert ek.allclose(sample_dx, uniform_hemisphere_to_square(-ray.d_x)) + assert ek.allclose(ray.d_x, expected_ray_dx.d) + assert ek.allclose(ray.o_x, expected_ray_dx.o) + + sample_dy = [0.5, 0.5 + 1.0 / n_y] + expected_ray_dy, _ = sensor.sample_ray(1.0, 1.0, sample_dy, sample2, True) + assert ek.allclose(sample_dy, uniform_hemisphere_to_square(-ray.d_y)) + assert ek.allclose(ray.d_y, expected_ray_dy.d) + assert ek.allclose(ray.o_y, expected_ray_dy.o) + + +def test_checkerboard(variants_all_rgb): + """ + Very basic render test with checkerboard texture and square target. + """ + from mitsuba.core import Bitmap, ScalarTransform4f, Struct + from mitsuba.core.xml import load_dict + + l_e = 1.0 # Emitted radiance + rho0 = 0.5 + rho1 = 1.0 + + # Scene setup + scene_dict = { + "type": "scene", + "shape": { + "type": "rectangle", + "bsdf": { + "type": "diffuse", + "reflectance": { + "type": "checkerboard", + "color0": rho0, + "color1": rho1, + "to_uv": ScalarTransform4f.scale(2), + }, + }, + }, + "emitter": {"type": "directional", "direction": [0, 0, -1], "irradiance": 1.0}, + "sensor0": { + "type": "hdistant", + "target": {"type": "rectangle"}, + "sampler": { + "type": "independent", + "sample_count": 50000, + }, + "film": { + "type": "hdrfilm", + "height": 4, + "width": 4, + "pixel_format": "luminance", + "component_format": "float32", + "rfilter": {"type": "box"}, + }, + }, + # "sensor1": { # In case one would like to check what the scene looks like + # "type": "perspective", + # "to_world": ScalarTransform4f.look_at(origin=[0, 0, 5], target=[0, 0, 0], up=[0, 1, 0]), + # "sampler": { + # "type": "independent", + # "sample_count": 10000, + # }, + # "film": { + # "type": "hdrfilm", + # "height": 256, + # "width": 256, + # "pixel_format": "luminance", + # "component_format": "float32", + # "rfilter": {"type": "box"}, + # }, + # }, + "integrator": {"type": "path"}, + } + + scene = load_dict(scene_dict) + scene.render() + result = np.array( + scene.sensors()[0].film().bitmap().convert( + Bitmap.PixelFormat.RGB, Struct.Type.Float32, False + ) + ).squeeze() + + expected = l_e * 0.5 * (rho0 + rho1) / ek.Pi + assert np.allclose(expected, result, atol=1e-3)