Implemented bounding spheres for objects

Cargo.toml

@@ -8,3 +8,4 @@ edition = "2018"
 
 [dependencies]
 nalgebra = "0.18"
+rand = "0.7.3"

README.md

@@ -18,14 +18,14 @@ This list may be changed or extended in the future.
   - [x] Plane struct
   - [x] Plane intersection test
   - [x] Color mapping on planes
-- [ ] Triangle objects
+- [x] Triangle objects
   - [x] Triangle struct
   - [x] Triangle intersection test
   - [x] Triangle normal generation
   - [x] Color mapping on triangles
   - [x] Triangle mesh struct
   - [x] Triangle mesh intersection test
-- [ ] Bounding boxes
+- [x] Bounding spheres
 - [ ] Direct lighting
   - [ ] Point light sources
     - [ ] Point source struct

src/main.rs

@@ -55,7 +55,7 @@ fn main() -> std::io::Result<()> {
     let camera = Camera::new(Point3::new(0.0,0.0,0.0), Vector3::new(0.0,0.0,1.0), 1.0, 16.0 / 9.0, 2.0, 480);
 
     let scene = vec![
-        Object::new(TriangleMesh::singleton(Point3::new(-1.0, -1.0, 3.0), Point3::new(0.0, 1.0, 3.0), Point3::new(1.0, -1.0, 3.0), |t, u, v| Color::new(t, u, v)))
+        Object::new_boundless(TriangleMesh::singleton(Point3::new(-1.0, -1.0, 2.0), Point3::new(0.0, 1.0, 2.0), Point3::new(1.0, -1.0, 2.0), |t, u, v| Color::new(t, u, v)))
     ];
 
     render(&camera, &scene, "out.ppm")

src/object.rs

@@ -2,6 +2,7 @@
 mod sphere; pub use sphere::*;
 mod plane; pub use plane::*;
 mod triangle; pub use triangle::*;
+mod bound; pub use bound::*;
 
 use na::*;
 
@@ -22,18 +23,49 @@ pub trait Surface {
     // Takes in a point (assumed to be on the object's surface)
     // and returns the color information on that point.
     fn getcolor(&self, point: Point3<f32>) -> Color;
+
+    // Creates a bounding sphere around the object.
+    fn bound(&self) -> Bound;
 }
 
 pub struct Object {
-    pub surface: Box<dyn Surface>
+    pub surface: Box<dyn Surface>,
+    bound: Bound
 }
 
+#[allow(dead_code)]
 impl Object {
-    pub fn new(surface: impl 'static + Surface) -> Self {
-        Object { surface: Box::new(surface) }
+    // Creates a new object with a custom bounding sphere.
+    pub fn new_(surface: impl 'static + Surface, center: Point3<f32>, radius: f32) -> Self {
+        Object {
+            surface: Box::new(surface),
+            bound: Bound { center: center, radius: radius, bypass: false }
+        }
     }
 
-    pub fn intersect(&self, ray: Ray) -> Option<f32> { self.surface.intersect(ray) }
+    // Creates a new object with no bounding sphere.
+    pub fn new_boundless(surface: impl 'static + Surface) -> Self {
+        Object {
+            surface: Box::new(surface),
+            bound: Bound::bypass()
+        }
+    }
+
+    // Creates a new object with the default bounding sphere.
+    pub fn new(surface: impl 'static + Surface) -> Self {
+        let bound = surface.bound();
+        Object {
+            surface: Box::new(surface),
+            bound: bound
+        }
+    }
+
+
+    pub fn intersect(&self, ray: Ray) -> Option<f32> {
+        if self.bound.is_intersected(ray) {
+            self.surface.intersect(ray)
+        } else { None }
+    }
     pub fn normal(&self, point: Point3<f32>) -> Unit<Vector3<f32>> { self.surface.normal(point) }
     pub fn getcolor(&self, point: Point3<f32>) -> Color { self.surface.getcolor(point) }
 }

src/object/bound.rs

@@ -0,0 +1,37 @@
+extern crate nalgebra as na;
+
+use na::distance;
+use na::geometry::Point3;
+
+use crate::types::Ray;
+
+// A bounding sphere, used for
+// intersection test optimization.
+#[derive(Debug)]
+pub struct Bound {
+    pub center: Point3<f32>,
+    pub radius: f32,
+
+    // If true, then the bounding sphere is disabled.
+    pub bypass: bool
+}
+
+impl Bound {
+    pub fn is_intersected(&self, ray: Ray) -> bool {
+
+        if self.bypass { return true; }
+
+        let l = ray.origin - self.center;
+
+        let b_2 = ray.direction.dot(&l);
+        let c = l.norm_squared() - self.radius * self.radius;
+
+        let discr = b_2 * b_2 * c;
+
+        discr >= 0.0
+    }
+
+    pub fn contains(&self, point: &Point3<f32>) -> bool { distance(&self.center, point) < self.radius }
+
+    pub fn bypass() -> Self { Bound { center: Point3::origin(), radius: 0.0, bypass: true } }
+}

src/object/plane.rs

@@ -4,7 +4,7 @@ use na::*;
 use na::geometry::Point3;
 
 use crate::types::*;
-use super::Surface;
+use super::{Surface, bound::*};
 
 pub struct Plane {
     pub center: Point3<f32>,        // Plane origin (used for texture mapping).
@@ -86,6 +86,10 @@ impl Surface for Plane {
 
         (*self.texture)(x, y)
     }
+
+    // Planes are infinite, so no finite
+    // bounding sphere could possibly contain one.
+    fn bound(&self) -> Bound { Bound::bypass() }
 }
 
 #[cfg(test)]

src/object/sphere.rs

@@ -6,7 +6,7 @@ use na::*;
 use na::geometry::Point3;
 
 use crate::types::*;
-use super::Surface;
+use super::{Surface, bound::*};
 
 pub struct Sphere {
     pub center: Point3<f32>, // Center point of the sphere.
@@ -36,25 +36,24 @@ impl Sphere {
 
 impl Surface for Sphere {
     fn intersect(&self, ray: Ray) -> Option<f32> {
-        fn solve_quadratic(a: f32, b: f32, c: f32) -> Option<(f32, f32)> {
-            let discr = b * b - 4.0 * a * c;
+        fn solve_quadratic(b: f32, c: f32) -> Option<(f32, f32)> {
+            let discr = b * b - 4.0 * c;
 
             if discr < 0.0 { None }
             else if discr == 0.0 {
-                let x = -0.5 * b / a;
+                let x = -0.5 * b;
                 Some((x, x))
             } else {
                 let q = if b > 0.0 { -0.5 * (b + discr.sqrt()) } else { -0.5 * (b - discr.sqrt()) };
-                Some((q / a, c / q))
+                Some((q, c / q))
             }
         }
 
         let l = ray.origin - self.center;
-        let a = ray.direction.dot(&ray.direction);
         let b = 2.0 * ray.direction.dot(&l);
-        let c = l.dot(&l) - self.radius * self.radius;
+        let c = l.normsquared() - self.radius * self.radius;
 
-        let (mut t0, mut t1) = solve_quadratic(a, b, c)?;
+        let (mut t0, mut t1) = solve_quadratic(b, c)?;
 
         if t0 > t1 { std::mem::swap(&mut t0, &mut t1); }
 
@@ -78,4 +77,6 @@ impl Surface for Sphere {
 
         (*self.texture)(x, y)
     }
+
+    fn bound(&self) -> Bound { Bound::bypass() }
 }

src/object/triangle.rs

@@ -6,7 +6,7 @@ use na::*;
 use na::geometry::Point3;
 
 use crate::types::*;
-use super::Surface;
+use super::{Surface, bound::*};
 
 pub struct Triangle {
     pub v1: usize, // Handles to 3 vertices.
@@ -162,6 +162,85 @@ impl Surface for TriangleMesh {
     fn getcolor(&self, point: Point3<f32>) -> Color {
         self.closest_tri(point).getcolor(&self.vertices, point)
     }
+
+    // Uses Welzl's algorithm to solve the bounding sphere problem
+    fn bound(&self) -> Bound {
+        fn triangle_sphere(point1: &Point3<f32>, point2: &Point3<f32>, point3: &Point3<f32>) -> (Point3<f32>, f32) {
+            let a = point3 - point1;
+            let b = point2 - point1;
+
+            let crs = b.cross(&a);
+
+            let to_center = (crs.cross(&b) * a.norm_squared() + a.cross(&crs) * b.norm_squared())
+                            / (2.0 * crs.norm_squared());
+
+            let radius = to_center.norm();
+
+            (point1 + to_center, radius)
+        }
+
+        fn tetrahedron_sphere(point1: &Point3<f32>, point2: &Point3<f32>, point3: &Point3<f32>, point4: &Point3<f32>) -> (Point3<f32>, f32) {
+            let matrix = Matrix4::from_rows(&[point1.to_homogeneous().transpose(),
+                                            point2.to_homogeneous().transpose(),
+                                            point3.to_homogeneous().transpose(),
+                                            point4.to_homogeneous().transpose()]);
+
+            let a = matrix.determinant() * 2.0;
+            let mut matrix_mut = matrix.clone();
+
+            let squares = Vector4::new(point1.coords.norm_squared(), point2.coords.norm_squared(), point3.coords.norm_squared(), point4.coords.norm_squared());
+            matrix_mut.set_column(0, &squares);
+            let center_x = matrix_mut.determinant();
+
+            matrix_mut.set_column(1, &matrix.index((.., 0)));
+            let center_y = -matrix_mut.determinant();
+
+            matrix_mut.set_column(2, &matrix.index((.., 1)));
+            let center_z = matrix_mut.determinant();
+
+            let center = Point3::new(center_x / a, center_y / a, center_z / a);
+            let radius = distance(point1, &center);
+
+            (center, radius)
+        }
+
+        fn smallest_sphere(points: Vec<&Point3<f32>>, boundary: Vec<&Point3<f32>>) -> (Point3<f32>, f32) {
+            println!("{:?}\n{:?}\n", points, boundary);
+            if points.len() == 0 || boundary.len() == 4 {
+                match boundary.len() {
+                    0 => (Point3::new(0.0, 0.0, 0.0), 0.0),
+                    1 => (*boundary[0], 0.0),
+                    2 => { let half_span = 0.5 * (boundary[1] - boundary[0]);
+                            (*boundary[0] + half_span, half_span.norm()) },
+                    3 => triangle_sphere(boundary[0], boundary[1], boundary[2]),
+                    4 => tetrahedron_sphere(boundary[0], boundary[1], boundary[2], boundary[3]),
+                    _ => unreachable!()
+                }
+            } else {
+                let removed = points[0];
+                let points = Vec::from(&points[1..]);
+
+                let bound = smallest_sphere(points.clone(), boundary.clone());
+                if distance(&bound.0, removed) < bound.1 { return bound; }
+
+                let mut boundary = boundary.clone();
+                boundary.push(removed);
+
+                smallest_sphere(points, boundary)
+            }
+        }
+
+        extern crate rand;
+        use rand::thread_rng;
+        use rand::seq::SliceRandom;
+
+        let mut points: Vec<&Point3<f32>> = self.vertices.iter().collect();
+        points.shuffle(&mut thread_rng());
+
+        let (center, radius) = smallest_sphere(points, Vec::new());
+
+        Bound { center: center, radius: radius + 0.01, bypass: false }
+    }
 }
 
 #[cfg(test)]
@@ -199,4 +278,66 @@ mod tests {
 
         assert_eq!(roundcolor(triangle.getcolor(point)), roundcolor(Color::new(t, u, v)));
     }
+
+    #[test]
+    fn triangle_bounds() {
+        let point1 = Point3::new(0.0, 0.0, 0.0);
+        let point2 = Point3::new(1.0, 0.0, 0.0);
+        let point3 = Point3::new(0.0, 1.0, 0.0);
+
+        let triangle = TriangleMesh::singleton_solid(point1, point2, point3, Color::black());
+
+        let bound = triangle.bound();
+
+        println!("{:?}", bound);
+
+        assert!(bound.contains(&point1));
+        assert!(bound.contains(&point2));
+        assert!(bound.contains(&point3));
+    }
+    /*
+    #[test]
+    fn triangle_tobound() {
+        let point1 = Point3::new(-3.0, 4.0, -6.0);
+        let point2 = Point3::new(5.0, -2.0, -7.0);
+        let point3 = Point3::new(9.0, -7.0, 3.0);
+
+        let (center, radius) = triangle_sphere(&point1, &point2, &point3);
+        let bound = Bound { center: center, radius: radius + 0.01, bypass: false };
+
+        println!("{:?}", bound);
+
+        println!("{}\n{}\n{}", distance(&bound.center, &point1),
+                                distance(&bound.center, &point2),
+                                distance(&bound.center, &point3));
+
+        assert!(bound.contains(&point1));
+        assert!(bound.contains(&point2));
+        assert!(bound.contains(&point3));
+    }
+
+    #[test]
+    fn triangle_tetrabound() {
+        let point1 = Point3::new(8.0, -2.0, -5.0);
+        let point2 = Point3::new(-3.0, 4.0, -6.0);
+        let point3 = Point3::new(-3.0, -9.0, 3.0);
+        let point4 = Point3::new(-6.0, 5.0, -9.0);
+
+        let (center, radius) = tetrahedron_sphere(&point1, &point2, &point3, &point4);
+
+        let bound = Bound { center: center, radius: radius + 0.01, bypass: false };
+
+        println!("{:?}", bound);
+
+        println!("{}\n{}\n{}\n{}", distance(&bound.center, &point1),
+                                distance(&bound.center, &point2),
+                                distance(&bound.center, &point3),
+                                distance(&bound.center, &point4));
+
+        assert!(bound.contains(&point1));
+        assert!(bound.contains(&point2));
+        assert!(bound.contains(&point3));
+        assert!(bound.contains(&point4));
+    }
+    */
 }

src/types.rs

@@ -30,6 +30,7 @@ pub struct Color {
     _private: () // Private field prevents direct construction
 }
 
+#[allow(dead_code)]
 impl Color {
     pub fn new(red: f32, green: f32, blue: f32) -> Self {
         Color {