From d1f2950921ff1e5534ee5785b9dbbac2de236c56 Mon Sep 17 00:00:00 2001 From: Vincent Leroy Date: Tue, 19 May 2020 17:21:42 +0200 Subject: [PATCH] Added distant sensor plugin Co-authored-by: Sebastian Schunke --- src/librender/scene.cpp | 4 + src/sensors/CMakeLists.txt | 1 + src/sensors/distant.cpp | 191 ++++++++++++++++++++++++++ src/sensors/tests/test_distant.py | 219 ++++++++++++++++++++++++++++++ 4 files changed, 415 insertions(+) create mode 100644 src/sensors/distant.cpp create mode 100644 src/sensors/tests/test_distant.py diff --git a/src/librender/scene.cpp b/src/librender/scene.cpp index 5275f9112..339d0e795 100644 --- a/src/librender/scene.cpp +++ b/src/librender/scene.cpp @@ -96,6 +96,10 @@ MTS_VARIANT Scene::Scene(const Properties &props) { // Create emitters' shapes (environment luminaires) for (Emitter *emitter: m_emitters) emitter->set_scene(this); + + // Create sensors' shapes (environment sensors) + for (Sensor *sensor: m_sensors) + sensor->set_scene(this); } MTS_VARIANT Scene::~Scene() { diff --git a/src/sensors/CMakeLists.txt b/src/sensors/CMakeLists.txt index 6c643a55c..ba9579b58 100644 --- a/src/sensors/CMakeLists.txt +++ b/src/sensors/CMakeLists.txt @@ -4,6 +4,7 @@ add_plugin(perspective perspective.cpp) add_plugin(radiancemeter radiancemeter.cpp) add_plugin(thinlens thinlens.cpp) add_plugin(irradiancemeter irradiancemeter.cpp) +add_plugin(distant distant.cpp) # Register the test directory add_tests(${CMAKE_CURRENT_SOURCE_DIR}/tests) diff --git a/src/sensors/distant.cpp b/src/sensors/distant.cpp new file mode 100644 index 000000000..0a71ae40b --- /dev/null +++ b/src/sensors/distant.cpp @@ -0,0 +1,191 @@ +#include +#include +#include +#include +#include +#include +#include +#include + +NAMESPACE_BEGIN(mitsuba) + +/**! + +.. _sensor-distant: + +Distant directional sensor (:monosp:`distant`) +---------------------------------------------- + +.. pluginparameters:: + + * - to_world + - |transform| + - Sensor-to-world transformation matrix. + * - direction + - |vector| + - Alternative (and exclusive) to `to_world`. Direction from which the + sensor will be recording in world coordinates. + * - target + - |point| + - *Optional.* Point (in world coordinates) to which sampled rays will be + cast. Useful for one-dimensional scenes. If unset, rays will be cast + uniformly over the entire scene. + +This sensor plugin implements a distant directional sensor which records +radiation leaving the scene in a given direction. If the ``target`` parameter +is not set, rays cast by the sensor will be distributed uniformly on the cross +section of the scene's bounding sphere. + +*/ + +MTS_VARIANT class DistantSensor final : public Sensor { +public: + MTS_IMPORT_BASE(Sensor, m_world_transform, m_needs_sample_3, m_film) + MTS_IMPORT_TYPES(Scene) + + DistantSensor(const Properties &props) : Base(props) { + /* Until `set_scene` is called, we have no information + about the scene and default to the unit bounding sphere. */ + m_bsphere = ScalarBoundingSphere3f(ScalarPoint3f(0.f), 1.f); + + if (props.has_property("direction")) { + if (props.has_property("to_world")) + Throw("Only one of the parameters 'direction' and 'to_world' " + "can be specified at the same time!'"); + + ScalarVector3f direction(normalize(props.vector3f("direction"))); + auto [up, unused] = coordinate_system(direction); + + m_world_transform = + new AnimatedTransform(ScalarTransform4f::look_at( + ScalarPoint3f(0.0f), ScalarPoint3f(direction), up)); + } + + if (props.has_property("target")) { + m_target = props.point3f("target"); + m_has_target = true; + Log(Debug, "Targeting point %s", m_target); + } else { + m_has_target = false; + } + + if (m_film->size() != ScalarPoint2i(1, 1)) + Throw("This sensor only supports films of size 1x1 Pixels!"); + + if (m_film->reconstruction_filter()->radius() > + 0.5f + math::RayEpsilon) + Log(Warn, "This sensor should be used with a reconstruction filter " + "with a radius of 0.5 or lower (e.g. default box)"); + + m_needs_sample_3 = false; + } + + void set_scene(const Scene *scene) override { + m_bsphere = scene->bbox().bounding_sphere(); + m_bsphere.radius = + max(math::RayEpsilon, + m_bsphere.radius * (1.f + math::RayEpsilon) ); + } + + std::pair sample_ray(Float time, Float wavelength_sample, + const Point2f &spatial_sample, + const Point2f & /*direction_sample*/, + Mask active) const override { + MTS_MASKED_FUNCTION(ProfilerPhase::EndpointSampleRay, active); + Ray3f ray; + ray.time = time; + + // 1. Sample spectrum + auto [wavelengths, wav_weight] = + sample_wavelength(wavelength_sample); + ray.wavelengths = wavelengths; + + // 2. Set ray direction + auto trafo = m_world_transform->eval(time, active); + ray.d = trafo.transform_affine(Vector3f{ 0.f, 0.f, 1.f }); + + // 3. Sample ray origin + if (!m_has_target) { + // If no target point is defined, sample a target point on the + // bounding sphere's cross section + Point2f offset = + warp::square_to_uniform_disk_concentric(spatial_sample); + Vector3f perp_offset = + trafo.transform_affine(Vector3f{ offset.x(), offset.y(), 0.f }); + ray.o = m_bsphere.center + (perp_offset - ray.d) * m_bsphere.radius; + } else { + ray.o = m_target - 2.f * ray.d * m_bsphere.radius; + } + + ray.update(); + return std::make_pair( + ray, m_has_target + ? wav_weight + : wav_weight * (math::Pi * sqr(m_bsphere.radius))); + } + + std::pair sample_ray_differential( + Float time, Float wavelength_sample, const Point2f &spatial_sample, + const Point2f & /*direction_sample*/, Mask active) const override { + MTS_MASKED_FUNCTION(ProfilerPhase::EndpointSampleRay, active); + RayDifferential3f ray; + ray.time = time; + + // 1. Sample spectrum + auto [wavelengths, wav_weight] = + sample_wavelength(wavelength_sample); + ray.wavelengths = wavelengths; + + // 2. Set ray direction + auto trafo = m_world_transform->eval(time, active); + ray.d = trafo.transform_affine(Vector3f{ 0.f, 0.f, 1.f }); + + // 3. Sample ray origin + if (!m_has_target) { + // If no target point is defined, sample a target point on the + // bounding sphere's cross section + Point2f offset = + warp::square_to_uniform_disk_concentric(spatial_sample); + Vector3f perp_offset = + trafo.transform_affine(Vector3f{ offset.x(), offset.y(), 0.f }); + ray.o = m_bsphere.center + (perp_offset - ray.d) * m_bsphere.radius; + } else { + ray.o = m_target - 2.f * ray.d * m_bsphere.radius; + } + + // 4. Set differentials; since the film size is always 1x1, we don't + // have differentials + ray.has_differentials = false; + + ray.update(); + return std::make_pair( + ray, m_has_target + ? wav_weight + : wav_weight * (math::Pi * sqr(m_bsphere.radius))); + } + + /// 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 << "DistantSensor[" << std::endl + << " world_transform = " << m_world_transform << "," << std::endl + << " bsphere = " << m_bsphere << "," << std::endl + << " film = " << m_film << "," << std::endl + << "]"; + return oss.str(); + } + + MTS_DECLARE_CLASS() + +protected: + ScalarBoundingSphere3f m_bsphere; + ScalarPoint3f m_target; + bool m_has_target; +}; + +MTS_IMPLEMENT_CLASS_VARIANT(DistantSensor, Sensor) +MTS_EXPORT_PLUGIN(DistantSensor, "DistantSensor"); +NAMESPACE_END(mitsuba) diff --git a/src/sensors/tests/test_distant.py b/src/sensors/tests/test_distant.py new file mode 100644 index 000000000..4329dea10 --- /dev/null +++ b/src/sensors/tests/test_distant.py @@ -0,0 +1,219 @@ +import numpy as np +import pytest + +import enoki as ek +import mitsuba + + +def xml_sensor(direction=None, target=None, xpixel=1): + if direction is None: + xml_direction = "" + else: + if type(direction) is not str: + direction = ",".join([str(x) for x in direction]) + xml_direction = f"""""" + + if target is None: + xml_target = "" + else: + if type(target) is not str: + target = ",".join([str(x) for x in target]) + xml_target = f"""""" + + xml_film = f""" + + + + """ + + return f""" + {xml_direction} + {xml_target} + {xml_film} + """ + + +def make_sensor(direction=None, target=None, xpixel=1): + from mitsuba.core.xml import load_string + return load_string(xml_sensor(direction, target, xpixel)) + + +def test_construct(variant_scalar_rgb): + from mitsuba.core.xml import load_string + + # Construct without parameters (wrong film size) + with pytest.raises(RuntimeError): + sensor = load_string("""""") + + # Construct with wrong film size + with pytest.raises(RuntimeError): + sensor = make_sensor(xpixel=2) + + # Construct with minimal parameters + sensor = make_sensor() + assert sensor is not None + assert not sensor.bbox().valid() # Degenerate bounding box + + # Construct with direction, check transform setup correctness + world_reference = [[0, 1, 0, 0], + [1, 0, 0, 0], + [0, 0, -1, 0], + [0, 0, 0, 1]] + sensor = make_sensor(direction=[0, 0, -1]) + assert ek.allclose( + sensor.world_transform().eval(0.).matrix, + world_reference + ) + + sensor = make_sensor(direction=[0, 0, -2]) + assert ek.allclose( + sensor.world_transform().eval(0.).matrix, + world_reference + ) + + +@pytest.mark.parametrize("target", [ + None, + [0.0, 0.0, 0.0], + [4.0, 1.0, 0.0], + [1.0, 1.0, 0.0], + [0.5, -0.3, 0.0], + +]) +@pytest.mark.parametrize("direction", [ + [0.0, 0.0, 1.0], + [-1.0, -1.0, 0.0], + [2.0, 0.0, 0.0] +]) +@pytest.mark.parametrize("ray_kind", ["regular", "differential"]) +def test_sample_ray(variant_scalar_rgb, direction, target, ray_kind): + sensor = make_sensor(direction, target) + + for (sample1, sample2) in [[[0.32, 0.87], [0.16, 0.44]], + [[0.17, 0.44], [0.22, 0.81]], + [[0.12, 0.82], [0.99, 0.42]], + [[0.72, 0.40], [0.01, 0.61]]]: + + if ray_kind == "regular": + ray, _ = sensor.sample_ray(1., 1., sample1, sample2, True) + elif ray_kind == "differential": + ray, _ = sensor.sample_ray_differential( + 1., 1., sample1, sample2, True) + assert not ray.has_differentials + + # Check that ray direction is what is expected + assert ek.allclose(ray.d, ek.normalize(direction)) + + # Check that ray origin is outside of bounding sphere + # Bounding sphere is centered at world origin and has radius 1 without scene + assert ek.norm(ray.o) >= 1. + + +def make_scene(direction=[0, 0, -1], target=None): + from mitsuba.core.xml import load_string + + scene_xml = f""" + + {xml_sensor(direction, target)} + + + """ + + return load_string(scene_xml) + + +@pytest.mark.parametrize("target", [[0, 0, 0], [0.5, 0, 1]]) +def test_target(variant_scalar_rgb, target): + # Check if the sensor correctly targets the point it is supposed to + scene = make_scene(direction=[0, 0, -1], target=target) + sensor = scene.sensors()[0] + sampler = sensor.sampler() + + ray, _ = sensor.sample_ray( + sampler.next_1d(), + sampler.next_1d(), + sampler.next_2d(), + sampler.next_2d() + ) + si = scene.ray_intersect(ray) + assert si.is_valid() + assert ek.allclose(si.p, [target[0], target[1], 0.], atol=1e-6) + + +@pytest.mark.parametrize("direction", [[0, 0, -1], [0.5, 0.5, -1]]) +def test_intersection(variant_scalar_rgb, direction): + # Check if the sensor correctly casts rays spread uniformly in the scene + direction = ek.normalize(direction) + scene = make_scene(direction=direction) + sensor = scene.sensors()[0] + sampler = sensor.sampler() + + n_rays = 10000 + isect = np.empty((n_rays, 3)) + + for i in range(n_rays): + ray, _ = sensor.sample_ray( + sampler.next_1d(), + sampler.next_1d(), + sampler.next_2d(), + sampler.next_2d() + ) + si = scene.ray_intersect(ray) + + if not si.is_valid(): + isect[i, :] = np.nan + else: + isect[i, :] = si.p[:] + + # Average intersection locations should be (in average) centered + # around (0, 0, 0) + isect_valid = isect[~np.isnan(isect).all(axis=1)] + assert np.allclose(isect_valid[:, :2].mean(axis=0), 0., atol=1e-2) + assert np.allclose(isect_valid[:, 2], 0., atol=1e-5) + + # Check number of invalid intersections + # We expect a ratio of invalid interactions equal to the square's area + # divided by the bounding sphere's cross section, weighted by the surface's + # slanting factor (cos theta) w.r.t the sensor's direction + n_invalid = np.count_nonzero(np.isnan(isect).all(axis=1)) + assert np.allclose(n_invalid / n_rays, 1. - 2. / np.pi * + ek.dot(direction, [0, 0, -1]), atol=1e-2) + + +@pytest.mark.parametrize("radiance", [10**x for x in range(-3, 4)]) +def test_render(variant_scalar_rgb, radiance): + # Test render results with a simple scene + from mitsuba.core.xml import load_string + import numpy as np + + scene_xml = """ + + + + + + + + + + + + + + + + + + + + + + + + """ + + scene = load_string(scene_xml, spp=1, radiance=radiance) + sensor = scene.sensors()[0] + scene.integrator().render(scene, sensor) + img = sensor.film().bitmap() + assert np.allclose(np.array(img), radiance)