Shadow terminator improvements

docs/references.bib

@@ -675,3 +675,16 @@ @article{Schuessler2017Microfacet
 	articleno    = {205},
 	numpages     = {12},
 }
+
+@inbook{Hanika2021,
+	title        = {Hacking the Shadow Terminator},
+	author       = {Hanika, Johannes},
+	year         = {2021},
+	booktitle    = {Ray Tracing Gems II: Next Generation Real-Time Rendering with DXR, Vulkan, and OptiX},
+	publisher    = {Apress},
+	address      = {Berkeley, CA},
+	pages        = {65--76},
+	doi          = {10.1007/978-1-4842-7185-8_4},
+	isbn         = {978-1-4842-7185-8},
+	editor       = {Marrs, Adam and Shirley, Peter and Wald, Ingo}
+}

include/mitsuba/python/docstr.h

@@ -4476,6 +4476,8 @@ static const char *__doc_mitsuba_Interaction_operator_assign_2 = R"doc(//! @})do
 
 static const char *__doc_mitsuba_Interaction_p = R"doc(Position of the interaction in world coordinates)doc";
 
+static const char *__doc_mitsuba_Interaction_ray_offset = R"doc((Optional) offset to be used to spawn rays from this interaction)doc";
+
 static const char *__doc_mitsuba_Interaction_spawn_ray = R"doc(Spawn a semi-infinite ray towards the given direction)doc";
 
 static const char *__doc_mitsuba_Interaction_spawn_ray_to = R"doc(Spawn a finite ray towards the given position)doc";

include/mitsuba/render/interaction.h

@@ -100,6 +100,9 @@ struct Interaction {
     /// Position of the interaction in world coordinates
     Point3f p;
 
+    /// (Optional) offset to be used to spawn rays from this interaction
+    Vector3f ray_offset;
+
     /// Geometric normal (only valid for \c SurfaceInteraction)
     Normal3f n;
 
@@ -113,7 +116,7 @@ struct Interaction {
     /// Constructor
     Interaction(Float t, Float time, const Wavelength &wavelengths,
                 const Point3f &p, const Normal3f &n = 0.f)
-        : t(t), time(time), wavelengths(wavelengths), p(p), n(n) { }
+        : t(t), time(time), wavelengths(wavelengths), p(p), ray_offset(0.f), n(n) { }
 
     /// Virtual destructor
     virtual ~Interaction() = default;
@@ -124,11 +127,12 @@ struct Interaction {
      * field should be set to an infinite value to mark invalid intersection records.
      */
     virtual void zero_(size_t size = 1) {
-        t           = dr::full<Float>(dr::Infinity<Float>, size);
-        time        = dr::zeros<Float>(size);
-        wavelengths = dr::zeros<Wavelength>(size);
-        p           = dr::zeros<Point3f>(size);
-        n           = dr::zeros<Normal3f>(size);
+        t             = dr::full<Float>(dr::Infinity<Float>, size);
+        time          = dr::zeros<Float>(size);
+        wavelengths   = dr::zeros<Wavelength>(size);
+        p             = dr::zeros<Point3f>(size);
+        ray_offset    = dr::zeros<Vector3f>(size);
+        n             = dr::zeros<Normal3f>(size);
     }
 
     /// Is the current interaction valid?
@@ -154,7 +158,7 @@ struct Interaction {
     //! @}
     // =============================================================
 
-    DRJIT_STRUCT(Interaction, t, time, wavelengths, p, n);
+    DRJIT_STRUCT(Interaction, t, time, wavelengths, p, ray_offset, n);
 
 private:
     /**
@@ -165,7 +169,8 @@ struct Interaction {
     Point3f offset_p(const Vector3f &d) const {
         Float mag = (1.f + dr::max(dr::abs(p))) * math::RayEpsilon<Float>;
         mag = dr::detach(dr::mulsign(mag, dr::dot(n, d)));
-        return dr::fmadd(mag, dr::detach(n), p);
+        Vector3f shading_offset = dr::select(mag > 0.f, dr::detach(ray_offset), 0.f);
+        return dr::fmadd(mag, dr::detach(n), p + shading_offset);
     }
 };
 
@@ -186,7 +191,7 @@ struct SurfaceInteraction : Interaction<Float_, Spectrum_> {
     using Spectrum = Spectrum_;
 
     // Make parent fields/functions visible
-    MI_IMPORT_BASE(Interaction, t, time, wavelengths, p, n, is_valid)
+    MI_IMPORT_BASE(Interaction, t, time, wavelengths, p, ray_offset, n, is_valid)
 
     MI_IMPORT_RENDER_BASIC_TYPES()
     MI_IMPORT_OBJECT_TYPES()
@@ -524,8 +529,8 @@ struct SurfaceInteraction : Interaction<Float_, Spectrum_> {
     //! @}
     // =============================================================
 
-    DRJIT_STRUCT(SurfaceInteraction, t, time, wavelengths, p, n, shape, uv,
-                 sh_frame, dp_du, dp_dv, dn_du, dn_dv, duv_dx,
+    DRJIT_STRUCT(SurfaceInteraction, t, time, wavelengths, p, ray_offset, n,
+                 shape, uv, sh_frame, dp_du, dp_dv, dn_du, dn_dv, duv_dx,
                  duv_dy, wi, prim_index, instance)
 };
 
@@ -545,7 +550,7 @@ struct MediumInteraction : Interaction<Float_, Spectrum_> {
     using Index = typename CoreAliases::UInt32;
 
     // Make parent fields/functions visible
-    MI_IMPORT_BASE(Interaction, t, time, wavelengths, p, n, is_valid)
+    MI_IMPORT_BASE(Interaction, t, time, wavelengths, p, ray_offset, n, is_valid)
     //! @}
     // =============================================================
 
@@ -609,7 +614,7 @@ struct MediumInteraction : Interaction<Float_, Spectrum_> {
     //! @}
     // =============================================================
 
-    DRJIT_STRUCT(MediumInteraction, t, time, wavelengths, p, n, medium,
+    DRJIT_STRUCT(MediumInteraction, t, time, wavelengths, p, ray_offset, n, medium,
                  sh_frame, wi, sigma_s, sigma_n, sigma_t,
                  combined_extinction, mint)
 };
@@ -755,7 +760,9 @@ std::ostream &operator<<(std::ostream &os, const Interaction<Float, Spectrum> &i
            << "  t = " << it.t << "," << std::endl
            << "  time = " << it.time << "," << std::endl
            << "  wavelengths = " << it.wavelengths << "," << std::endl
-           << "  p = " << string::indent(it.p, 6) << std::endl
+           << "  p = " << string::indent(it.p, 6) << "," << std::endl
+           << "  ray_offset = " << string::indent(it.ray_offset, 15) << "," << std::endl
+           << "  n = " << string::indent(it.n, 6) << std::endl
            << "]";
     }
     return os;

include/mitsuba/render/mesh.h

@@ -604,6 +604,7 @@ class MI_EXPORT_LIB Mesh : public Shape<Float, Spectrum> {
     /// Flag that can be set by the user to disable loading/computation of vertex normals
     bool m_face_normals = false;
     bool m_flip_normals = false;
+    ScalarFloat m_ray_offset_scale = 1.0f;
 
     /* Surface area distribution -- generated on demand when \ref
        prepare_area_pmf() is first called. */

src/render/mesh.cpp

@@ -28,6 +28,7 @@ MI_VARIANT Mesh<Float, Spectrum>::Mesh(const Properties &props) : Base(props) {
        appearance. Default: ``false`` */
     m_face_normals = props.get<bool>("face_normals", false);
     m_flip_normals = props.get<bool>("flip_normals", false);
+    m_ray_offset_scale = props.get<ScalarFloat>("ray_offset_scale", 1.0f);
 
     m_discontinuity_types = (uint32_t) DiscontinuityFlags::PerimeterType;
 
@@ -1553,6 +1554,22 @@ Mesh<Float, Spectrum>::compute_surface_interaction(const Ray3f &ray,
         si.sh_frame.n = si.n;
     }
 
+    if (has_vertex_normals() &&
+        likely(has_flag(ray_flags, RayFlags::ShadingFrame) ||
+               has_flag(ray_flags, RayFlags::dNSdUV)) &&
+        (m_ray_offset_scale > 0.f)) {
+        // Implements "Hacking the shadow terminator" by Johannes Hanika (2021).
+        // The code matches the original implementation, but was slightly
+        // simplified since b0 * p0 + b1 * p1 + b2 * p2 = si.p, and b0 + b1 + b2 = 1.
+        Vector3f tmp0 = p0 - si.p, tmp1 = p1 - si.p, tmp2 = p2 - si.p;
+        Float dot0 = dr::maximum(dr::dot(tmp0, n0), 0.f),
+              dot1 = dr::maximum(dr::dot(tmp1, n1), 0.f),
+              dot2 = dr::maximum(dr::dot(tmp2, n2), 0.f);
+        si.ray_offset =
+            m_ray_offset_scale *
+            dr::detach(dr::fmadd(b0, dot0 * n0, dr::fmadd(b1, dot1 * n1, b2 * dot2 * n2)));
+    }
+
     if (m_flip_normals) {
         si.n = -si.n;
         si.sh_frame.n = -si.sh_frame.n;

src/render/python/interaction_v.cpp

@@ -16,6 +16,7 @@ MI_PY_EXPORT(Interaction) {
         .def_field(Interaction3f, time,        D(Interaction, time))
         .def_field(Interaction3f, wavelengths, D(Interaction, wavelengths))
         .def_field(Interaction3f, p,           D(Interaction, p))
+        .def_field(Interaction3f, ray_offset,  D(Interaction, ray_offset))
         .def_field(Interaction3f, n,           D(Interaction, n))
         // Methods
         .def(nb::init<>(), D(Interaction, Interaction))
@@ -30,7 +31,7 @@ MI_PY_EXPORT(Interaction) {
         .def("zero_",        &Interaction3f::zero_, D(Interaction, zero))
         .def_repr(Interaction3f);
 
-    MI_PY_DRJIT_STRUCT(it, Interaction3f, t, time, wavelengths, p, n)
+    MI_PY_DRJIT_STRUCT(it, Interaction3f, t, time, wavelengths, p, ray_offset, n)
 }
 
 MI_PY_EXPORT(SurfaceInteraction) {
@@ -91,7 +92,7 @@ MI_PY_EXPORT(SurfaceInteraction) {
             D(SurfaceInteraction, has_n_partials))
         .def_repr(SurfaceInteraction3f);
 
-    MI_PY_DRJIT_STRUCT(si, SurfaceInteraction3f, t, time, wavelengths, p, n,
+    MI_PY_DRJIT_STRUCT(si, SurfaceInteraction3f, t, time, wavelengths, p, ray_offset, n,
                        shape, uv, sh_frame, dp_du, dp_dv, dn_du, dn_dv, duv_dx,
                        duv_dy, wi, prim_index, instance)
 }
@@ -118,7 +119,7 @@ MI_PY_EXPORT(MediumInteraction) {
         .def("to_local", &MediumInteraction3f::to_local, "v"_a, D(MediumInteraction, to_local))
         .def_repr(MediumInteraction3f);
 
-    MI_PY_DRJIT_STRUCT(mi, MediumInteraction3f, t, time, wavelengths, p, n,
+    MI_PY_DRJIT_STRUCT(mi, MediumInteraction3f, t, time, wavelengths, p, ray_offset, n,
                        medium, sh_frame, wi, sigma_s, sigma_n, sigma_t,
                        combined_extinction, mint)
 }

src/render/tests/test_interaction.py

@@ -42,6 +42,7 @@ def test02_intersection_construction(variant_scalar_rgb):
   time=2,
   wavelengths=[],
   p=[1, 2, 3],
+  ray_offset=[0, 0, 0],
   n=[4, 5, 6],
   shape=0x0,
   uv=[7, 8],

src/render/tests/test_mesh.py

@@ -1410,3 +1410,29 @@ def test36_mesh_vcalls_with_directed_edges(variants_vec_rgb):
 
     result = mesh_ptr.opposite_dedge(mi.UInt32([2, 3, 2]))
     assert dr.all(result == mi.UInt32([3, 2, 3]))
+
+
+@pytest.mark.parametrize("ray_offset_scale", [0.0, 1.0])
+def test37_mesh_ray_offset(variants_vec_rgb, ray_offset_scale):
+    props = mi.Properties()
+    props["ray_offset_scale"] = ray_offset_scale
+    mesh = mi.Mesh(props)
+    params = mi.traverse(mesh)
+    params['vertex_positions'] = [-1, -1, 0, 1, -1, 0, 0,  1, 0]
+    params['faces'] = [0, 1, 2]
+    params.update()
+    params['vertex_normals'] = dr.ravel(dr.normalize(
+        mi.Vector3f([-1, 1, 0], [-1, -1, 1], [1, 1, 1])))
+    params.update()
+    scene = mi.load_dict({
+        "type": "scene",
+        "mesh": mesh,
+    })
+    si = scene.ray_intersect(mi.Ray3f([0, 0, 1], [0, 0, -1]))
+    ray = si.spawn_ray(dr.normalize(mi.Vector3f(1, -2, 3)))
+    if ray_offset_scale == 0.0:
+        assert dr.allclose(si.ray_offset, 0.0)
+        assert dr.allclose(ray.o.z, 0.0, atol=1e-4)
+    else:
+        assert dr.allclose(si.ray_offset.x, 0.0) # Due to symmetry, x offset is 0.
+        assert ray.o.z > 0.5 # Offset along z axis should be significant.

src/shapes/obj.cpp

@@ -20,7 +20,7 @@ Wavefront OBJ mesh loader (:monosp:`obj`)
 -----------------------------------------
 
 .. pluginparameters::
- :extra-rows: 5
+ :extra-rows: 6
 
  * - filename
    - |string|
@@ -41,6 +41,15 @@ Wavefront OBJ mesh loader (:monosp:`obj`)
    - Is the mesh inverted, i.e. should the normal vectors be flipped? (Default:|false|, i.e.
      the normals point outside)
 
+ * - ray_offset_scale
+   - |float|
+   - Scale factor used to offset ray origins when tracing secondary rays
+     (e.g., shadow rays) to avoid self-intersections. Rays are offset
+     according to the mismatch between geometric and shading normals, following
+     :cite:`Hanika2021`. This offset can be reduced by using a factor < 1,
+     and disabled entirely by setting it to 0 (e.g., to prevent light leaks).
+     (Default: 1.0)
+
  * - to_world
    - |transform|
    - Specifies an optional linear object-to-world transformation.

src/shapes/ply.cpp

@@ -23,7 +23,7 @@ PLY (Stanford Triangle Format) mesh loader (:monosp:`ply`)
 ----------------------------------------------------------
 
 .. pluginparameters::
- :extra-rows: 4
+ :extra-rows: 5
 
  * - filename
    - |string|
@@ -44,6 +44,15 @@ PLY (Stanford Triangle Format) mesh loader (:monosp:`ply`)
    - Is the mesh inverted, i.e. should the normal vectors be flipped? (Default:|false|, i.e.
      the normals point outside)
 
+ * - ray_offset_scale
+    - |float|
+    - Scale factor used to offset ray origins when tracing secondary rays
+      (e.g., shadow rays) to avoid self-intersections. Rays are offset
+      according to the mismatch between geometric and shading normals, following
+      :cite:`Hanika2021`. This offset can be reduced by using a factor < 1,
+      and disabled entirely by setting it to 0 (e.g., to prevent light leaks).
+      (Default: 1.0)
+
  * - to_world
    - |transform|
    - Specifies an optional linear object-to-world transformation.