From ebe56bc6c3a5b6b1d20a1c2fcf13d0189574d95c Mon Sep 17 00:00:00 2001
From: gcem <34423323+gcem@users.noreply.github.com>
Date: Wed, 14 Apr 2021 00:32:59 +0300
Subject: [PATCH] add faster intersection tests for shadow rays

---
 src/Objects/Surface/Mesh.cpp     | 11 +++++++++++
 src/Objects/Surface/Mesh.hpp     | 16 ++++++++++++++++
 src/Objects/Surface/Sphere.cpp   | 21 ++++++++++++++++++++
 src/Objects/Surface/Sphere.hpp   | 16 ++++++++++++++++
 src/Objects/Surface/Surface.hpp  | 16 ++++++++++++++++
 src/Objects/Surface/Triangle.cpp | 33 ++++++++++++++++++++++++++++++++
 src/Objects/Surface/Triangle.hpp | 18 +++++++++++++++++
 src/PathTracer/PathTracer.cpp    |  4 +---
 8 files changed, 132 insertions(+), 3 deletions(-)

diff --git a/src/Objects/Surface/Mesh.cpp b/src/Objects/Surface/Mesh.cpp
index 8a76c03..e69dd43 100644
--- a/src/Objects/Surface/Mesh.cpp
+++ b/src/Objects/Surface/Mesh.cpp
@@ -30,4 +30,15 @@ Mesh::intersect(const Ray& ray,
     }
     return minT < std::numeric_limits<FloatT>::infinity() ? minT : -1;
 }
+
+bool
+Mesh::intersectsBefore(const Ray& ray, FloatT maxT, FloatT epsilon) const
+{
+    for (auto& triangle : triangles) {
+        if (triangle.intersectsBefore(ray, maxT, epsilon))
+            return true;
+    }
+
+    return false;
+}
 }
\ No newline at end of file
diff --git a/src/Objects/Surface/Mesh.hpp b/src/Objects/Surface/Mesh.hpp
index 6057713..72d6982 100644
--- a/src/Objects/Surface/Mesh.hpp
+++ b/src/Objects/Surface/Mesh.hpp
@@ -63,6 +63,22 @@ class Mesh : public Surface
                      LinearAlgebra::Vec3& normalOut,
                      FloatT epsilon = 0) const override;
 
+    /**
+     * @brief Finds if ray intersects this mesh with a t value lower than
+     * maxT
+     *
+     * Looks for an intersection with t in range (0, maxT)
+     *
+     * @param ray
+     * @param maxT Upper limit of t value for intersection test to succeed
+     * (exclusive)
+     * @return true There is a triangle closer than t
+     * @return false There is no such triangle
+     */
+    bool intersectsBefore(const Ray& ray,
+                          FloatT maxT,
+                          FloatT epsilon = 0) const override;
+
 protected:
     /**
      * @brief Triangles in this mesh
diff --git a/src/Objects/Surface/Sphere.cpp b/src/Objects/Surface/Sphere.cpp
index b99280f..8171d4c 100644
--- a/src/Objects/Surface/Sphere.cpp
+++ b/src/Objects/Surface/Sphere.cpp
@@ -47,4 +47,25 @@ Sphere::intersect(const Ray& ray,
     normalOut = (ray.origin + ray.direction * t - center).normalize();
     return t;
 }
+
+bool
+Sphere::intersectsBefore(const Ray& ray, FloatT maxT, FloatT epsilon) const
+{
+    // see intersect() for an explanation of this calculation
+
+    auto a = ray.direction.squaredNorm();
+    auto b = 2 * ray.direction.dot(ray.origin - center);
+    auto c = (ray.origin - center).squaredNorm() - radius * radius;
+
+    auto delta = b * b - 4 * a * c;
+
+    if (delta < 0)
+        return false;
+
+    auto t = (-b - sqrt(delta)) / (2 * a);
+    if (t <= 0)
+        t = (-b + sqrt(delta)) / (2 * a);
+
+    return t > 0 && t < maxT;
+}
 }
\ No newline at end of file
diff --git a/src/Objects/Surface/Sphere.hpp b/src/Objects/Surface/Sphere.hpp
index 83b9e7b..da139a3 100644
--- a/src/Objects/Surface/Sphere.hpp
+++ b/src/Objects/Surface/Sphere.hpp
@@ -54,6 +54,22 @@ class Sphere : public Surface
                      LinearAlgebra::Vec3& normalOut,
                      FloatT epsilon = 0) const override;
 
+    /**
+     * @brief Finds if ray intersects this sphere with a t value lower than
+     * maxT
+     *
+     * Looks for an intersection with t in range (0, maxT)
+     *
+     * @param ray
+     * @param maxT Upper limit of t value for intersection test to succeed
+     * (exclusive)
+     * @return true There is an intersection closer than t
+     * @return false There is no such intersection
+     */
+    bool intersectsBefore(const Ray& ray,
+                          FloatT maxT,
+                          FloatT epsilon = 0) const override;
+
 protected:
     /**
      * @brief Center
diff --git a/src/Objects/Surface/Surface.hpp b/src/Objects/Surface/Surface.hpp
index 485700b..4971c41 100644
--- a/src/Objects/Surface/Surface.hpp
+++ b/src/Objects/Surface/Surface.hpp
@@ -45,6 +45,22 @@ class Surface
                              LinearAlgebra::Vec3& normalOut,
                              FloatT epsilon = 0) const = 0;
 
+    /**
+     * @brief Finds if ray intersects this surface with a t value lower than
+     * maxT
+     *
+     * Looks for an intersection with t in range (0, maxT)
+     *
+     * @param ray
+     * @param maxT Upper limit of t value for intersection test to succeed
+     * (exclusive)
+     * @return true There is an intersection closer than t
+     * @return false There is no such intersection
+     */
+    virtual bool intersectsBefore(const Ray& ray,
+                                  FloatT maxT,
+                                  FloatT epsilon = 0) const = 0;
+
     /**
      * @brief Material of this Surface
      *
diff --git a/src/Objects/Surface/Triangle.cpp b/src/Objects/Surface/Triangle.cpp
index 5660da5..79aeaac 100644
--- a/src/Objects/Surface/Triangle.cpp
+++ b/src/Objects/Surface/Triangle.cpp
@@ -54,6 +54,39 @@ Triangle::intersect(const Ray& ray, FloatT epsilon) const
     return t;
 }
 
+bool
+Triangle::intersectsBefore(const Ray& ray, FloatT maxT, FloatT epsilon) const
+{
+    // see intersect() for an explanation of this calculation
+
+    auto edge21 = v1 - v2;
+    auto edge31 = v1 - v3;
+    auto rhs = v1 - ray.origin;
+
+    auto detA =
+      LinearAlgebra::Mat3(ray.direction, edge21, edge31).determinant();
+
+    auto detT = LinearAlgebra::Mat3(rhs, edge21, edge31).determinant();
+    auto t = detT / detA;
+
+    if (t <= 0 || t >= maxT)
+        return false;
+
+    auto detB = LinearAlgebra::Mat3(ray.direction, rhs, edge31).determinant();
+    auto beta = detB / detA;
+
+    if (beta <= -epsilon || beta >= 1 + epsilon)
+        return false;
+
+    auto detC = LinearAlgebra::Mat3(ray.direction, edge21, rhs).determinant();
+    auto gamma = detC / detA;
+
+    if (gamma <= -epsilon || beta + gamma >= 1 + epsilon)
+        return false;
+
+    return true;
+}
+
 LinearAlgebra::Vec3
 Triangle::getNormal() const
 {
diff --git a/src/Objects/Surface/Triangle.hpp b/src/Objects/Surface/Triangle.hpp
index 4cffebc..6b2b4f3 100644
--- a/src/Objects/Surface/Triangle.hpp
+++ b/src/Objects/Surface/Triangle.hpp
@@ -52,6 +52,24 @@ class Triangle
      */
     FloatT intersect(const Ray& ray, FloatT epsilon = 0) const;
 
+    /**
+     * @brief Finds if ray intersects this triangle with a t value lower than
+     * maxT
+     *
+     * Looks for an intersection with t in range (0, maxT)
+     *
+     * @param ray
+     * @param maxT Upper limit of t value for intersection test to succeed
+     * (exclusive)
+     * @param epsilon Baricentric coordinates are allowed to be in the range
+     * (-epsilon, 1 + epsilon) instead of (0, 1)
+     * @return true This triangle is closer than t
+     * @return false There is no intersection in this range
+     */
+    bool intersectsBefore(const Ray& ray,
+                          FloatT maxT,
+                          FloatT epsilon = 0) const;
+
     /**
      * @brief Normal vector
      *
diff --git a/src/PathTracer/PathTracer.cpp b/src/PathTracer/PathTracer.cpp
index 90a5595..d5a878f 100644
--- a/src/PathTracer/PathTracer.cpp
+++ b/src/PathTracer/PathTracer.cpp
@@ -125,9 +125,7 @@ PathTracer::lightVisible(const LinearAlgebra::Vec3& point,
     // than the light, t > 1 means the surface is behind the light
 
     for (auto surface : scene.surfaces) {
-        LinearAlgebra::Vec3 normal;
-        auto t = surface->intersect(ray, normal, scene.intersectionTestEpsilon);
-        if (t != -1 && t < 1)
+        if (surface->intersectsBefore(ray, 1))
             return false;
     }
     return true;