Intersection tests

src/bvh.rs

@@ -5,9 +5,26 @@
 
 use crate::geometry::Vec3;
 use crate::geometry::Axis;
+use crate::geometry::Ray;
 use crate::mesh::Face;
 use crate::mesh::Mesh;
 
+pub struct Intersection<'a> {
+    pub t: f32,
+    pub fr: FaceReference<'a>
+}
+
+impl PartialOrd for Intersection<'_> {
+    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
+        return self.t.partial_cmp(&other.t);
+    }
+}
+
+impl PartialEq for Intersection<'_> {
+    fn eq(&self, other: &Self) -> bool {
+        self.t == other.t
+    }
+}
 
 /**
  * Axis-Aligned Bounding Box.
@@ -50,6 +67,33 @@ impl AABB {
 
         return aabb;
     }
+    
+    /**
+     * Test for intersection of ABB with Ray r.
+     * t is the maximum distance of the ray from its origin
+     */
+    pub fn intersect(&self, r: &Ray, t: f32) -> bool {
+        // See [Tavian Barnes' site](https://tavianator.com/2011/ray_box.html) for details.
+        let tx1 = (self.min.x - r.origin.x) * r.n_inv.x;
+        let tx2 = (self.max.x - r.origin.x) * r.n_inv.x;
+
+        let mut tmin = f32::min(tx1, tx2);
+        let mut tmax = f32::max(tx1, tx2);
+
+        let ty1 = (self.min.y - r.origin.y) * r.n_inv.y;
+        let ty2 = (self.max.y - r.origin.y) * r.n_inv.y;
+
+        tmin = f32::max(tmin, f32::min(ty1, ty2));
+        tmax = f32::min(tmax, f32::max(ty1, ty2));
+
+        let tz1 = (self.min.z - r.origin.z) * r.n_inv.z;
+        let tz2 = (self.max.z - r.origin.z) * r.n_inv.z;
+
+        tmin = f32::max(tmin, f32::min(tz1, tz2));
+        tmax = f32::min(tmax, f32::max(tz1, tz2));
+
+        return tmax >= f32::max(0.0, tmin) && tmin < t;
+    }
 }
 
 
@@ -63,7 +107,7 @@ pub enum BVH<'a> {
 }
 
 
-#[derive(Clone)]
+#[derive(Copy, Clone)]
 pub struct FaceReference<'a> {
     m: &'a Mesh,
     f: &'a Face
@@ -119,6 +163,81 @@ impl BVH<'_> {
 
         return bvh;
     }
+
+    fn intersect_ray_triangle<'a>(ray: &Ray, fr: &'a FaceReference) -> Option<Intersection<'a>> {
+        // See Akenine-Möller et al. Real-Time Rendering (2018), 22.8 Ray/Triangle Intersection
+        let p = fr.m.face_coords(fr.f);
+        let e1 = p[1] - p[0];
+        let e2 = p[2] - p[0];
+        let q = ray.direction.cross(e2);
+
+        let a = e1 * q;
+        if f32::abs(a) < 1e-9 {
+            return None;
+        }
+
+        let f = 1.0 / a;
+        let s = &ray.origin - p[0];
+
+        let u = f * (s * q);
+        if u < 0.0 {
+            return None
+        };
+
+        let r = s.cross(e1);
+        let v = f * (ray.direction * r);
+        if v < 0.0 || u + v > 1.0 {
+            return None;
+        }
+
+        let t = f * (e2 * r);
+
+        return Some(Intersection {
+            t: t,
+            fr: *fr
+        });
+
+    }
+
+    pub fn intersect<'a>(&'a self, r : &Ray) -> Option<Intersection<'a>> {
+        match self {
+            BVH::Leaf(aabb, frs) => {
+                // Early-out if the ray does not intersect the AABB
+                if aabb.intersect(r, f32::INFINITY) {
+                    return None;
+                }
+
+                // Make intersection tests with individual triangles
+                let closest_intersection : Option<Intersection> = frs.iter()
+                    .map(|fr| Self::intersect_ray_triangle(r, fr))
+                    .filter(|i| i.is_some())
+                    .map(|i| i.unwrap())
+                    .filter(|i| i.t >= 0.0)
+                    .reduce(|a, b| if a < b { a } else { b});
+
+                return closest_intersection;
+            },
+            BVH::Node(aabb, left, right) => {
+                // Early-out if the ray does not intersect the AABB
+                if aabb.intersect(r, f32::INFINITY) {
+                    return None;
+                }
+
+                // Recurse down the tree
+                let left_intersection = left.intersect(r);
+                let right_intersection = right.intersect(r);
+
+                if left_intersection.is_some() && right_intersection.is_some() {
+                    let (li, ri) = (left_intersection.unwrap(), right_intersection.unwrap());
+                    return if li < ri { Some(li) } else {Some(ri)};
+                } else {
+                    // At least one intersection is None.
+                    return left_intersection.or(right_intersection);
+                }
+            }
+        }
+    }
+
 }
 
 impl std::fmt::Display for AABB {

src/geometry.rs

@@ -103,6 +103,14 @@ impl std::ops::Sub for Vec3 {
 
 }
 
+impl std::ops::Sub for &Vec3 {
+    type Output = Vec3;
+
+    fn sub(self, _rhs: &Vec3) -> Vec3 {
+        return *self + (-*_rhs);
+    }
+}
+
 // Dot product
 impl std::ops::Mul<Vec3> for Vec3 {
     type Output = f32;
@@ -156,10 +164,12 @@ impl Axis {
 pub struct Ray {
     pub origin: Vec3,
     pub direction: Vec3,
+    pub n_inv: Vec3
 }
 
-pub struct Intersection {
-    pub t: f32,
-    pub pos: Vec3,
-    pub normal: Vec3,
+impl Ray {
+    pub fn calc_n_inv(&mut self) {
+        self.n_inv = Vec3::new(1.0 / self.direction.x, 1.0 / self.direction.y, 1.0 / self.direction.z);
+    }
 }
+

src/main.rs

@@ -12,51 +12,8 @@ use crate::random::RNG;
 use crate::mesh::Mesh;
 use crate::geometry::Ray;
 use crate::image::Image;
-use crate::geometry::Intersection;
 use crate::bvh::BVH;
 
-// ╺┳╸┏━┓┏━┓┏━╸┏━╸   ┏━┓┏━┓╻ ╻┏━┓
-//  ┃ ┣┳┛┣━┫┃  ┣╸    ┣┳┛┣━┫┗┳┛┗━┓
-//  ╹ ╹┗╸╹ ╹┗━╸┗━╸   ╹┗╸╹ ╹ ╹ ┗━┛
-
-pub fn intersect_sphere(ray: &Ray, center: &Vec3, radius: f32) -> Option<Intersection> {
-    let v: Vec3 = *center - ray.origin;
-    let cos_alpha = v.normalized() * ray.direction;
-
-    let t = cos_alpha * v.length();
-
-    let closest_point = ray.origin + t * ray.direction;
-    let min_distance = center.distance_to(closest_point);
-
-    if min_distance < radius {
-        let x = (radius * radius - min_distance * min_distance).sqrt();
-        let intersection = closest_point - x * ray.direction;
-        let normal = (intersection - *center).normalized();
-
-        return Some(Intersection {
-            t: t,
-            pos: intersection,
-            normal: normal,
-        });
-    } else {
-        return None;
-    }
-
-}
-
-fn intersect(r: &Ray) -> Option<Intersection> {
-    // Have a simple scene with a singular sphere at (0,0,5) with radius r=1
-    return intersect_sphere(r, &Vec3 {x: 0.0, y: 0.0, z: -5.0}, 1.8);
-}
-
-fn raytrace(r: &Ray) -> Vec3 {
-    let i = intersect(r);
-    
-    match i {
-        None => Vec3 { x: 0.0, y: 0.0, z: 0.0 },
-        Some(intersect) => Vec3 { x: 1.0, y: 1.0, z: 0.0 }
-    }
-}
 
 // ┏┳┓┏━┓╻┏┓╻
 // ┃┃┃┣━┫┃┃┗┫
@@ -66,9 +23,9 @@ fn main() {
 
     let m = Mesh::load("benchmarks/suzanne.obj");
     println!("First vertex: {}", m[0].vertices[m[0].vertices.len() - 1]);
-    let bvh = BVH::from(&m, 1);
-    println!("\n\n{}", bvh);
-    return;
+    println!("Building BVH…");
+    let bvh = BVH::from(&m, 8);
+    println!("BVH created");
     let mut rng = RNG::create();
 
     // Generate camera rays
@@ -85,7 +42,7 @@ fn main() {
             let target_x = (f32::from(x) / f32::from(WIDTH) - 0.5) * (FOV / 2.0).tan() * D;
             let target_y = (f32::from(y) / f32::from(HEIGHT) - 0.5) * (FOV / 2.0).tan() * D;
             let pixel = img.mod_pixel(x.into(), y.into());
-            let r = Ray {
+            let mut r = Ray {
                 origin: Vec3 {
                     x: 0.0,
                     y: 0.0,
@@ -95,34 +52,32 @@ fn main() {
                     x: target_x,
                     y: target_y,
                     z: D
-                }
+                },
+                n_inv: Vec3::new(0.0, 0.0, 0.0)
             };
+            r.calc_n_inv();
+
             //println!("Shooting ray {} {} {}", target_x, target_y, D);
             
             for _ in 0..SPP {
-                let sample : Vec3= 255.0 / f32::from(SPP) * raytrace(&r);
+                let intersection = bvh.intersect(&r);
+
+                if intersection.is_some() {
+                    println!("HIT! {}", r.direction);
+                }
+
+                let color = intersection
+                    .map(|i| Vec3::new(i.t, i.t, i.t))
+                    .unwrap_or(Vec3::new(0.0, 0.0, 0.0));
+
+
+                //let sample : Vec3= 255.0 / f32::from(SPP) * raytrace(&r);
                 //println!("{} --> {}", r.direction, sample);
-                *pixel += sample;
+                *pixel += color;
             }
         }
+        println!("Progress: {}", f32::from(x) / f32::from(WIDTH));
     }
 
     img.save_bmp(&Path::new("Render.bmp")).expect("Failed to save rendered image");
 }
-
-
-// ╺┳╸┏━╸┏━┓╺┳╸┏━┓
-//  ┃ ┣╸ ┗━┓ ┃ ┗━┓
-//  ╹ ┗━╸┗━┛ ╹ ┗━┛
-
-#[test]
-fn sphereIntersect() {
-    let result = intersect_sphere(&Ray {
-        origin: Vec3::new(0.0, 0.0, -5.0),
-        direction: Vec3::new(0.0, 0.0, 1.0)
-
-    }, &Vec3::new(0.0, 0.0, 0.0),
-    1.0);
-
-    assert!(result.is_some(), "Should have found an intersection");
-}

src/mesh.rs

@@ -10,7 +10,6 @@ use std::path::Path;
 
 use crate::geometry::Vec3;
 use crate::geometry::Ray;
-use crate::geometry::Intersection;
 
 
 /** Data structure for a single triangle face. */
@@ -298,12 +297,6 @@ impl Mesh {
         return objects;
     }
 
-    pub fn intersect(&self, ray: &Ray) -> Option<Intersection> {
-
-
-        None
-    }
-
     pub fn face_coords(&self, f: &Face) -> [&Vec3; 3] {
         return [
             &self.vertices[f.v[0]],