Implement BVH

src/bvh.rs

@@ -4,6 +4,7 @@
 
 
 use crate::geometry::Vec3;
+use crate::geometry::Axis;
 use crate::mesh::Face;
 use crate::mesh::Mesh;
 
@@ -16,21 +17,6 @@ pub struct AABB {
     pub max: Vec3
 }
 
-/**
- * Bounding Volume Hierarchy
- */
-pub struct FaceReference<'a> {
-    m: &'a Mesh,
-    f: &'a Face
-}
-
-pub enum BVH<'a> {
-    Node(&'a BVH<'a>, &'a BVH<'a>),
-    Leaf(&'a BVH<'a>, Vec<FaceReference<'a>>)
-
-
-}
-
 fn f32_min(a: f32, b: f32) -> f32 {
     return if a < b { a } else { b };
 }
@@ -52,16 +38,103 @@ impl AABB {
         self.max.z = f32::max(self.max.z, coords.iter().map(|c| c.z).reduce(f32_max).unwrap());
     }
 
-    pub fn from(m: &Mesh, f: &Face) -> AABB {
+    pub fn from(fr: &[FaceReference]) -> AABB {
         let mut aabb = AABB {
             min: Vec3::new(f32::INFINITY, f32::INFINITY, f32::INFINITY),
             max: Vec3::new(-f32::INFINITY, -f32::INFINITY, -f32::INFINITY),
         };
 
-        aabb.extend(m, f);
-        
+        for face_reference in fr {
+            aabb.extend(face_reference.m, face_reference.f);
+        }
+
         return aabb;
     }
 }
 
 
+
+/**
+ * Bounding Volume Hierarchy
+ */
+pub enum BVH<'a> {
+    Node(AABB, Box<BVH<'a>>, Box<BVH<'a>>),
+    Leaf(AABB, Vec<FaceReference<'a>>)
+}
+
+
+#[derive(Clone)]
+pub struct FaceReference<'a> {
+    m: &'a Mesh,
+    f: &'a Face
+}
+
+impl BVH<'_> {
+    fn build_recursive<'a, 'b>(frs: &'b mut [FaceReference<'a>], recursion_limit: u16, split_axis: Axis) -> BVH<'a> {
+        let aabb = AABB::from(frs);
+        let frs_count = frs.len();
+
+        if recursion_limit == 0 || frs_count <= 3 {
+            // We have reached the maximum of recursions necessary.
+            return BVH::Leaf(aabb, Vec::from(frs));
+        }
+
+        // Sort along the current axis
+        frs.sort_by(|a,b| a.f.center.select(split_axis).partial_cmp(
+                              &b.f.center.select(split_axis)).unwrap());
+
+        // Select using surface area heuristic
+        let total_surface: f32= frs.iter().map(|fr| fr.f.area).sum();
+
+        let mut split_index : usize = 0;
+        let mut surface_area = 0.0;
+        while surface_area < total_surface / 2.0 {
+            surface_area += frs[split_index].f.area;
+            split_index += 1;
+        }
+
+        // Recurse
+        let (frs1, frs2) = frs.split_at_mut(split_index);
+        let left = Self::build_recursive(frs1, recursion_limit - 1, split_axis.next());
+        let right = Self::build_recursive(frs2, recursion_limit - 1, split_axis.next());
+
+        return BVH::Node(aabb, Box::new(left), Box::new(right));
+    }
+
+    pub fn from<'a>(meshes: &'a Vec<Mesh>, recursion_depth: u16) -> BVH {
+        // Create references for each triangle
+        let mut references = vec!();
+
+        for m in meshes {
+            for f in &m.faces {
+                references.push(FaceReference {
+                    m: &m,
+                    f: &f
+                });
+            }
+        }
+
+
+        let bvh = Self::build_recursive(&mut references, recursion_depth, Axis::x);
+
+        return bvh;
+    }
+}
+
+impl std::fmt::Display for AABB {
+    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
+        return write!(f, "{} -> {}", self.min, self.max);
+    }
+}
+impl std::fmt::Display for BVH<'_> {
+    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
+        match self {
+            BVH::Leaf(aabb,_) => {
+                return write!(f, "(Leaf {})", aabb);
+            }
+            BVH::Node(aabb, left, right) => {
+                return write!(f, "(Node {}; {} {})", aabb, left, right);
+            }
+        }
+    }
+}

src/geometry.rs

@@ -47,6 +47,14 @@ impl Vec3 {
             z: self.x * _rhs.y - self.y * _rhs.x,
         };
     }
+
+    pub fn select(&self, a: Axis) -> f32 {
+        match a {
+            x => self.x,
+            y => self.y,
+            z => self.z
+        }
+    }
 }
 
 impl std::fmt::Display for Vec3 {
@@ -124,6 +132,27 @@ impl std::ops::Mul<Vec3> for f32 {
     }
 }
 
+/**
+ * Coordinate space axis
+ */
+#[derive(Clone,Copy)]
+pub enum Axis {
+    x,
+    y,
+    z
+}
+
+impl Axis {
+    pub fn next(&self) -> Axis {
+        match self {
+            Axis::x => Axis::y,
+            Axis::y => Axis::z,
+            Axis::z => Axis::x
+        }
+    }
+
+}
+
 pub struct Ray {
     pub origin: Vec3,
     pub direction: Vec3,

src/main.rs

@@ -13,6 +13,7 @@ use crate::mesh::Mesh;
 use crate::geometry::Ray;
 use crate::image::Image;
 use crate::geometry::Intersection;
+use crate::bvh::BVH;
 
 // ╺┳╸┏━┓┏━┓┏━╸┏━╸   ┏━┓┏━┓╻ ╻┏━┓
 //  ┃ ┣┳┛┣━┫┃  ┣╸    ┣┳┛┣━┫┗┳┛┗━┓
@@ -65,8 +66,8 @@ fn main() {
 
     let m = Mesh::load("benchmarks/suzanne.obj");
     println!("First vertex: {}", m[0].vertices[m[0].vertices.len() - 1]);
-    let aabb = crate::bvh::AABB::from(&m[0], &m[0].faces[0]);
-    println!("AABB: {} -> {}", aabb.min, aabb.max);
+    let bvh = BVH::from(&m, 1);
+    println!("\n\n{}", bvh);
     return;
     let mut rng = RNG::create();
 

src/mesh.rs

@@ -23,6 +23,39 @@ pub struct Face {
     
     /** Texture coordinates */
     tex: [usize; 3],
+
+    /** Center point. Needed for BVH construction */
+    pub center: Vec3,
+
+    /** Surface Area */
+    pub area: f32
+
+
+}
+
+impl Face {
+    pub fn from(m: &Mesh, t: &Vec<[usize; 3]>) -> Face {
+        let mut face = Face {
+            v:   [t[0][0] - 1, t[1][0] - 1, t[2][0] - 1],
+            n:   [t[0][1] - 1, t[1][1] - 1, t[2][1] - 1],
+            tex: [t[0][2] - 1, t[1][2] - 1, t[2][2] - 1],
+            center: Vec3::new(0.0, 0.0, 0.0),
+            area: 0.0
+        };
+
+
+        // Calculate centroid
+        face.center.x = face.v.iter().map(|i| m.vertices[*i].x).sum::<f32>() / 3.0;
+        face.center.y = face.v.iter().map(|i| m.vertices[*i].y).sum::<f32>() / 3.0;
+        face.center.z = face.v.iter().map(|i| m.vertices[*i].z).sum::<f32>() / 3.0;
+
+        // Calculate surface area
+        let ab = m.vertices[face.v[1]] - m.vertices[face.v[0]];
+        let ac = m.vertices[face.v[2]] - m.vertices[face.v[0]];
+        face.area = ab.cross(ac).length() / 2.0;
+
+        return face;
+    }
 }
 
 #[derive(Copy,Clone,Debug)]
@@ -229,15 +262,8 @@ impl Mesh {
                         panic!("{}", error_msg);
                     }
 
-                    // We want the vertices ordered by type, so we have to transpose the triple
-                    // vector. Furthermore, Wavefront OBJ files are 1-based indices, so we subtract
-                    // one to convert them to 0-based indices.
-                    o.faces.push(Face {
-                        v:   [t[0][0] - 1, t[1][0] - 1, t[2][0] - 1],
-                        n:   [t[0][1] - 1, t[1][1] - 1, t[2][1] - 1],
-                        tex: [t[0][2] - 1, t[1][2] - 1, t[2][2] - 1],
-
-                    });
+                    let face = Face::from(&o, &t);
+                    o.faces.push(face);
                 }
                 "mtllib" => {
                     let filename = arguments.next().expect(format!("Missing filename in line {}", line_counter).as_str());