diff --git a/Cargo.toml b/Cargo.toml
index 4a319a5..6d39549 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -5,3 +5,5 @@ edition = "2021"
 
 [dependencies]
 anyhow = "*"
+vecmath = "*"
+itertools = "*"
diff --git a/src/bvh.rs b/src/bvh.rs
new file mode 100644
index 0000000..8e1d557
--- /dev/null
+++ b/src/bvh.rs
@@ -0,0 +1,76 @@
+use vecmath::vec3_add;
+
+use crate::{primitives::Box3, stl::parser::StlSolid};
+
+/// A node representing a node single in the bounding volume hierarchy
+#[derive(Default)]
+struct BvhNode {
+    aabb: Box3,
+    left_child: usize,
+    right_child: usize,
+    prim_idx: usize,
+    prim_count: usize,
+}
+
+/// Struct representing the total bounding volume hierarchy
+pub struct Bvh {
+    nodes: Vec<BvhNode>,
+}
+
+/// Return type for get_subdivision_position_and_axis
+struct SplitPlane {
+    axis: usize,
+    pos: f64,
+}
+
+impl Bvh {
+    /// Create a new Bounding volume hierarchy from an STL solid
+    pub fn new(solid: &StlSolid) -> Self {
+        // First we need to get all of the triangle centroids
+        let centroids: Vec<[f32; 3]> = solid
+            .triangles
+            .iter()
+            .map(|tri| {
+                tri.iter()
+                    .fold([0.0; 3], |acc, idx| vec3_add(solid.vertices[*idx], acc))
+            })
+            .collect();
+
+        // Initialize the root node
+        let mut nodes = vec![BvhNode {
+            prim_idx: 0,
+            prim_count: solid.triangles.len(),
+            ..Default::default()
+        }];
+        nodes[0]
+            .aabb
+            .shrink_wrap_primitives(&solid.vertices, &solid.triangles);
+
+        let mut bvh = Bvh { nodes };
+
+        bvh.subdivide(0);
+        bvh
+    }
+
+    /// Subdivide a node, subdivision doesn't occur if the node has only 2 primitives
+    pub fn subdivide(&mut self, node_idx: usize) {
+        let split_plane = self.get_split_plane(node_idx);
+    }
+
+    /// Calculate the optimal split plane axis and position.
+    /// This function should use the surface area heuristic to select the optimal split
+    fn get_split_plane(&self, node_idx: usize) -> SplitPlane {
+        // Initial approach just splits halfway along the longest axis
+        // TODO: Improve this to use the surface area heuristic
+        let axis = (0..3).fold(0, |acc, a| {
+            if self.nodes[node_idx].aabb.extents[acc] < self.nodes[node_idx].aabb.extents[a] {
+                a
+            } else {
+                acc
+            }
+        });
+
+        let pos = self.nodes[node_idx].aabb.extents[axis] / 2.0;
+        SplitPlane { axis, pos }
+    }
+}
diff --git a/src/lib.rs b/src/lib.rs
index fa3c805..d7364e1 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,2 +1,5 @@
+pub mod bvh;
+pub mod octree;
+pub mod primitives;
 pub mod stl;
 pub mod utilities;
diff --git a/src/octree/base.rs b/src/octree/base.rs
new file mode 100644
index 0000000..0f8ba62
--- /dev/null
+++ b/src/octree/base.rs
@@ -0,0 +1,9 @@
+pub struct Octree {
+    _octants: Vec<Octant>,
+    _min_size: f64,
+}
+
+pub struct Octant {
+    _morton_id: usize,
+    _level: usize,
+}
diff --git a/src/octree/mod.rs b/src/octree/mod.rs
new file mode 100644
index 0000000..2859ed3
--- /dev/null
+++ b/src/octree/mod.rs
@@ -0,0 +1,2 @@
+pub mod base;
+pub mod morton_ids;
diff --git a/src/octree/morton_ids.rs b/src/octree/morton_ids.rs
new file mode 100644
index 0000000..65f0969
--- /dev/null
+++ b/src/octree/morton_ids.rs
@@ -0,0 +1,17 @@
+/// Calculate the morton id in 3d space based on integer position
+pub fn encode(mut x: usize, mut y: usize, mut z: usize) {
+    let mut id = 0;
+    let mut level = 0;
+    while x > 0 || y > 0 || z > 0 {
+        let xbit = x % 2;
+        let ybit = x % 2;
+        let zbit = x % 2;
+        let curr_code = zbit + 2 * ybit + 4 * xbit;
+
+        id += curr_code << (3 * level);
+        x /= 2;
+        y /= 2;
+        z /= 2;
+        level += 1;
+    }
+}
diff --git a/src/primitives.rs b/src/primitives.rs
new file mode 100644
index 0000000..51779a5
--- /dev/null
+++ b/src/primitives.rs
@@ -0,0 +1,73 @@
+use std::collections::BTreeSet;
+
+use itertools::Itertools;
+use vecmath::{vec3_add, vec3_scale};
+
+#[derive(Default)]
+pub struct Box3 {
+    pub center: [f64; 3],
+    pub extents: [f64; 3],
+}
+
+impl Box3 {
+    /// Create a new box from the bounding values
+    pub fn new_from_bounds(min_bound: [f32; 3], max_bound: [f32; 3]) -> Self {
+        let mut center = [0.0; 3];
+        let mut extents = [0.0; 3];
+        for (i, (min_val, max_val)) in min_bound.iter().zip(max_bound).enumerate() {
+            center[i] = ((min_val + max_val) / 2.0) as f64;
+            extents[i] = (max_val - min_val) as f64;
+        }
+        Box3 { center, extents }
+    }
+
+    /// This function returns the 8 vertices of an axis-aligned bounding box.
+    pub fn get_vertices(&self) -> Vec<[f64; 3]> {
+        let half_extent = vec3_scale(self.extents, 0.5);
+
+        // Define all eight corners relative to the center.
+        [-half_extent[0], half_extent[0]]
+            .into_iter()
+            .cartesian_product([-half_extent[1], half_extent[1]])
+            .cartesian_product([-half_extent[2], half_extent[2]])
+            .map(|((x, y), z)| vec3_add(self.center, [x, y, z]))
+            .collect()
+    }
+
+    /// Expand a box to fit all primitive inside of it.
+    /// NOTE: vertices does not need to include only vertices being considered. Instead only
+    /// vertices included as indices in the primitives list are considered
+    pub fn shrink_wrap_primitives<T>(&mut self, vertices: &[[f32; 3]], primitives: &[T])
+    where
+        for<'a> &'a T: IntoIterator<Item = &'a usize>,
+    {
+        // If the current primitive list is empty, zero bounding box and exit
+        if primitives.is_empty() {
+            self.center = [0.0; 3];
+            self.extents = [0.0; 3]
+        }
+        // Get all of the unique vertex indices so we don't have to check them multiple times
+        let unique_vertex_idxs: BTreeSet<&usize> = BTreeSet::from_iter(primitives.iter().flatten());
+
+        // Now expand the box
+        let mut min_bd = [f32::INFINITY; 3];
+        let mut max_bd = [f32::NEG_INFINITY; 3];
+
+        for idx in unique_vertex_idxs.into_iter() {
+            for ((vertex_pos, min_pos), max_pos) in vertices[*idx]
+                .iter()
+                .zip(min_bd.iter_mut())
+                .zip(max_bd.iter_mut())
+            {
+                *min_pos = vertex_pos.min(*min_pos);
+                *max_pos = vertex_pos.max(*max_pos);
+            }
+        }
+    }
+}
+
+pub struct Triangle {
+    pub v1: [f64; 3],
+    pub v2: [f64; 3],
+    pub v3: [f64; 3],
+}
diff --git a/src/stl/parser.rs b/src/stl/parser.rs
index 09bc43f..063d4a7 100644
--- a/src/stl/parser.rs
+++ b/src/stl/parser.rs
@@ -10,9 +10,9 @@ use std::path::Path;
 /// [f32; 3]. Triangles are represented as a [usize; 3] where each entry is an index in the
 /// vertices Vec
 pub struct StlSolid {
-    vertices: Vec<[f32; 3]>,
-    triangles: Vec<[usize; 3]>,
-    normals: Vec<[f32; 3]>,
+    pub vertices: Vec<[f32; 3]>,
+    pub triangles: Vec<[usize; 3]>,
+    pub normals: Vec<[f32; 3]>,
 }
 
 pub fn parse_ascii_stl(file_path: &impl AsRef<Path>) -> Result<StlSolid> {
diff --git a/src/utilities/intersection_checks.rs b/src/utilities/intersection_checks.rs
new file mode 100644
index 0000000..e43dccf
--- /dev/null
+++ b/src/utilities/intersection_checks.rs
@@ -0,0 +1,50 @@
+use core::f64;
+
+use vecmath::{vec3_cross, vec3_dot, vec3_sub};
+
+use crate::primitives::{Box3, Triangle};
+
+/// Determine whether two shapes interesect using the separating axis theorem
+pub fn tri_aabb_intersect(bx: &Box3, tri: &Triangle) -> bool {
+    // Get edge vectors for triangle
+    let e1 = vec3_sub(tri.v2, tri.v1);
+    let e2 = vec3_sub(tri.v3, tri.v2);
+    let e3 = vec3_sub(tri.v1, tri.v3);
+
+    // Get face and edge normals
+    let tri_face = vec3_cross(e1, e2);
+    let n1 = vec3_cross(e1, e1);
+    let n2 = vec3_cross(e2, e1);
+    let n3 = vec3_cross(e3, e1);
+
+    // Face normals AABB
+    let u0 = [1.0, 0.0, 0.0];
+    let u1 = [0.0, 1.0, 0.0];
+    let u2 = [0.0, 0.0, 1.0];
+
+    // Now check all of the axes
+    let box_vertices = bx.get_vertices();
+    for axis in [u0, u1, u2, tri_face, n1, n2, n3] {
+        let box_projection = project_pnts(&box_vertices, axis);
+        let tri_projection = project_pnts(&[tri.v1, tri.v2, tri.v3], axis);
+
+        let start = f64::max(box_projection.0, tri_projection.0);
+        let end = f64::min(box_projection.1, tri_projection.1);
+        if start < end {
+            return false;
+        }
+    }
+    return true;
+}
+
+/// Project points unto an Axis and return the max and min
+fn project_pnts(pnts: &[[f64; 3]], axis: [f64; 3]) -> (f64, f64) {
+    let mut min_proj = f64::INFINITY;
+    let mut max_proj = f64::NEG_INFINITY;
+    for pnt in pnts {
+        let projection = vec3_dot(axis, *pnt);
+        min_proj = min_proj.min(projection);
+        max_proj = max_proj.max(projection);
+    }
+    (min_proj, max_proj)
+}
diff --git a/src/utilities/mod.rs b/src/utilities/mod.rs
index 6755121..8ee97d8 100644
--- a/src/utilities/mod.rs
+++ b/src/utilities/mod.rs
@@ -1,2 +1,3 @@
+pub mod intersection_checks;
 #[cfg(test)]
 pub mod test_macros;