diff --git a/rtiow/renderer/src/bvh_triangles.rs b/rtiow/renderer/src/bvh_triangles.rs
new file mode 100644
index 0000000..1fe2dc9
--- /dev/null
+++ b/rtiow/renderer/src/bvh_triangles.rs
@@ -0,0 +1,294 @@
+/// Implementation based on blog post @
+/// https://jacco.ompf2.com/2022/04/13/how-to-build-a-bvh-part-1-basics/
+use std::f32::EPSILON;
+
+use stl::STL;
+
+use crate::{
+    aabb::AABB,
+    hitable::{Hit, HitRecord},
+    material::Material,
+    ray::Ray,
+    vec3::{cross, dot, Vec3},
+};
+
+#[derive(Debug)]
+struct BVHNode {
+    aabb: AABB,
+    left_child: usize,
+    right_child: usize,
+    first_prim: usize,
+    prim_count: usize,
+}
+
+impl BVHNode {
+    fn is_leaf(&self) -> bool {
+        self.prim_count > 0
+    }
+}
+
+#[derive(Debug)]
+pub struct Triangle {
+    centroid: Vec3,
+    verts: [Vec3; 3],
+}
+
+#[derive(Debug)]
+pub struct BVHTriangles<M>
+where
+    M: Material,
+{
+    pub triangles: Vec<Triangle>,
+    material: M,
+    bvh_nodes: Vec<BVHNode>,
+}
+
+const ROOT_NODE_IDX: usize = 0;
+impl<M> BVHTriangles<M>
+where
+    M: Material,
+{
+    pub fn new(stl: &STL, material: M) -> BVHTriangles<M> {
+        let triangles: Vec<_> = stl
+            .triangles
+            .iter()
+            .map(|t| {
+                let v0 = t.verts[0];
+                let v1 = t.verts[1];
+                let v2 = t.verts[2];
+                let centroid = (v0 + v1 + v2) * 0.3333;
+                Triangle {
+                    centroid,
+                    verts: [v0, v1, v2],
+                }
+            })
+            .collect();
+
+        let n = 2 * triangles.len() - 2;
+        let bvh_nodes = Vec::with_capacity(n);
+        let mut bvh = BVHTriangles {
+            triangles,
+            bvh_nodes,
+            material,
+        };
+        bvh.build_bvh();
+        bvh
+    }
+
+    fn build_bvh(&mut self) {
+        // assign all triangles to root node
+        let root = BVHNode {
+            aabb: AABB::new(0f32.into(), 0f32.into()),
+            left_child: 0,
+            right_child: 0,
+            first_prim: 0,
+            prim_count: self.triangles.len() - 1,
+        };
+        self.bvh_nodes.push(root);
+        self.update_node_bounds(ROOT_NODE_IDX);
+        // subdivide recursively
+        self.subdivide(ROOT_NODE_IDX);
+    }
+    fn update_node_bounds(&mut self, node_idx: usize) {
+        let node = &mut self.bvh_nodes[node_idx];
+        let mut aabb_min: Vec3 = f32::MAX.into();
+        let mut aabb_max: Vec3 = f32::MIN.into();
+        for i in node.first_prim..node.prim_count {
+            let leaf_tri = &self.triangles[i];
+            aabb_min = vec3::min(aabb_min, leaf_tri.verts[0]);
+            aabb_min = vec3::min(aabb_min, leaf_tri.verts[1]);
+            aabb_min = vec3::min(aabb_min, leaf_tri.verts[2]);
+            aabb_max = vec3::max(aabb_max, leaf_tri.verts[0]);
+            aabb_max = vec3::max(aabb_max, leaf_tri.verts[1]);
+            aabb_max = vec3::max(aabb_max, leaf_tri.verts[2]);
+        }
+        node.aabb = AABB::new(aabb_min, aabb_max);
+    }
+    fn subdivide(&mut self, idx: usize) {
+        // Early out if we're down to just 2 or less triangles
+        if self.bvh_nodes[idx].prim_count <= 2 {
+            return;
+        }
+        let (first_prim, prim_count, left_count, i) = {
+            let node = &self.bvh_nodes[idx];
+
+            // Compute split plane and position.
+            let extent = node.aabb.min() - node.aabb.max();
+            let axis = node.aabb.longest_axis();
+            let split_pos = node.aabb.min()[axis] + extent[axis] * 0.5;
+
+            // Split the group in two halves.
+            let mut i = node.first_prim as isize;
+            let mut j = i + node.prim_count as isize - 1;
+            while i <= j {
+                if self.triangles[i as usize].centroid[axis] < split_pos {
+                    i += 1;
+                } else {
+                    self.triangles.swap(i as usize, j as usize);
+                    j -= 1;
+                }
+            }
+
+            // Create child nodes for each half.
+            let left_count = i as usize - node.first_prim;
+            if left_count == 0 || left_count == node.prim_count {
+                return;
+            }
+            (node.first_prim, node.prim_count, left_count, i)
+        };
+
+        // create child nodes
+        let left_child_idx = self.bvh_nodes.len();
+        let right_child_idx = self.bvh_nodes.len() + 1;
+        let left = BVHNode {
+            aabb: AABB::new(0f32.into(), 0f32.into()),
+            left_child: 0,
+            right_child: 0,
+            first_prim: first_prim,
+            prim_count: left_count,
+        };
+        let right = BVHNode {
+            aabb: AABB::new(0f32.into(), 0f32.into()),
+            left_child: 0,
+            right_child: 0,
+            first_prim: i as usize,
+            prim_count: prim_count - left_count,
+        };
+        self.bvh_nodes.push(left);
+        self.bvh_nodes.push(right);
+        let node = &mut self.bvh_nodes[idx];
+        node.left_child = left_child_idx;
+        node.right_child = right_child_idx;
+        node.prim_count = 0;
+
+        // Recurse
+        self.update_node_bounds(left_child_idx);
+        self.update_node_bounds(right_child_idx);
+        self.subdivide(left_child_idx);
+        self.subdivide(right_child_idx);
+    }
+
+    fn intersect_bvh(&self, r: Ray, node_idx: usize, t_min: f32, t_max: f32) -> Option<HitRecord> {
+        let node = &self.bvh_nodes[node_idx];
+        if !node.aabb.hit(r, t_min, t_max) {
+            return None;
+        }
+
+        if node.is_leaf() {
+            for tri in &self.triangles {
+                if let Some(RayTriangleResult { t, p }) =
+                    ray_triangle_intersect_moller_trumbore(r, tri)
+                {
+                    //if let Some(RayTriangleResult { t, p }) = ray_triangle_intersect_geometric(r, tri) {
+                    // We don't support UV (yet?).
+                    let uv = (0.5, 0.5);
+                    let v0 = tri.verts[0];
+                    let v1 = tri.verts[1];
+                    let v2 = tri.verts[2];
+
+                    let v0v1 = v1 - v0;
+                    let v0v2 = v2 - v0;
+                    let normal = cross(v0v1, v0v2);
+                    return Some(HitRecord {
+                        t,
+                        uv,
+                        p,
+                        normal,
+                        material: &self.material,
+                    });
+                }
+            }
+        } else {
+            let hr = self.intersect_bvh(r, node.left_child, t_min, t_max);
+            if hr.is_some() {
+                return hr;
+            }
+            let hr = self.intersect_bvh(r, node.left_child + 1, t_min, t_max);
+            if hr.is_some() {
+                return hr;
+            }
+        }
+        None
+    }
+}
+
+impl<M> Hit for BVHTriangles<M>
+where
+    M: Material,
+{
+    fn hit(&self, r: Ray, t_min: f32, t_max: f32) -> Option<HitRecord> {
+        self.intersect_bvh(r, 0, t_min, t_max)
+    }
+
+    fn bounding_box(&self, _t_min: f32, _t_max: f32) -> Option<AABB> {
+        Some(self.bvh_nodes[0].aabb)
+    }
+}
+
+struct RayTriangleResult {
+    t: f32,
+    p: Vec3,
+}
+///
+/// Based on https://www.scratchapixel.com/lessons/3d-basic-rendering/ray-tracing-rendering-a-triangle/moller-trumbore-ray-triangle-intersection.html
+fn ray_triangle_intersect_moller_trumbore(r: Ray, tri: &Triangle) -> Option<RayTriangleResult> {
+    // #ifdef MOLLER_TRUMBORE
+    //     Vec3f v0v1 = v1 - v0;
+    //     Vec3f v0v2 = v2 - v0;
+    //     Vec3f pvec = dir.crossProduct(v0v2);
+    //     float det = v0v1.dotProduct(pvec);
+    // #ifdef CULLING
+    //     // if the determinant is negative, the triangle is 'back facing'
+    //     // if the determinant is close to 0, the ray misses the triangle
+    //     if (det < kEpsilon) return false;
+    // #else
+    //     // ray and triangle are parallel if det is close to 0
+    //     if (fabs(det) < kEpsilon) return false;
+    // #endif
+    //     float invDet = 1 / det;
+    //
+    //     Vec3f tvec = orig - v0;
+    //     u = tvec.dotProduct(pvec) * invDet;
+    //     if (u < 0 || u > 1) return false;
+    //
+    //     Vec3f qvec = tvec.crossProduct(v0v1);
+    //     v = dir.dotProduct(qvec) * invDet;
+    //     if (v < 0 || u + v > 1) return false;
+    //
+    //     t = v0v2.dotProduct(qvec) * invDet;
+    //
+    let v0 = tri.verts[0];
+    let v1 = tri.verts[1];
+    let v2 = tri.verts[2];
+
+    let v0v1 = v1 - v0;
+    let v0v2 = v2 - v0;
+    let p = cross(r.direction, v0v2);
+    let det = dot(v0v1, p);
+    if det < EPSILON {
+        return None;
+    }
+
+    let inv_det = 1. / det;
+
+    let t = r.origin - v0;
+    let u = dot(t, p) * inv_det;
+    if u < 0. || u > 1. {
+        return None;
+    }
+
+    let q = cross(t, v0v1);
+    let v = dot(r.direction, q) * inv_det;
+    if v < 0. || u + v > 1. {
+        return None;
+    }
+    let t = dot(v0v2, q) * inv_det;
+
+    if t > EPSILON {
+        return Some(RayTriangleResult {
+            t,
+            p: r.origin + r.direction * t,
+        });
+    }
+    None
+}
diff --git a/rtiow/renderer/src/lib.rs b/rtiow/renderer/src/lib.rs
index d88add8..0f40cbc 100644
--- a/rtiow/renderer/src/lib.rs
+++ b/rtiow/renderer/src/lib.rs
@@ -1,5 +1,6 @@
 pub mod aabb;
 pub mod bvh;
+pub mod bvh_triangles;
 pub mod camera;
 pub mod constant_medium;
 pub mod cuboid;
diff --git a/rtiow/renderer/src/scenes/stltest.rs b/rtiow/renderer/src/scenes/stltest.rs
index 92cbee6..a11803a 100644
--- a/rtiow/renderer/src/scenes/stltest.rs
+++ b/rtiow/renderer/src/scenes/stltest.rs
@@ -2,6 +2,7 @@ use std::io::{BufReader, Cursor};
 use stl::STL;
 
 use crate::{
+    bvh_triangles::BVHTriangles,
     camera::Camera,
     hitable::Hit,
     hitable_list::HitableList,
@@ -12,7 +13,6 @@ use crate::{
     scale::Scale,
     sphere::Sphere,
     texture::ConstantTexture,
-    triangles::Triangles,
     vec3::Vec3,
 };
 
@@ -102,12 +102,11 @@ pub fn new(opt: &Opt) -> Scene {
         )),
         // STL Mesh
         Box::new(Scale::new(
-            Triangles::new(
+            BVHTriangles::new(
                 &stl_cube,
                 Lambertian::new(ConstantTexture::new(Vec3::new(0.6, 0.6, 0.6))),
-                1.,
             ),
-            0.5,
+            200.,
         )),
     ];
     let world: Box<dyn Hit> = if opt.use_accel {