diff --git a/ray-tracing-core/src/aabb.rs b/ray-tracing-core/src/aabb.rs
index 65b36e9..c1a36cb 100644
--- a/ray-tracing-core/src/aabb.rs
+++ b/ray-tracing-core/src/aabb.rs
@@ -49,6 +49,14 @@ impl AABB {
         self.max
     }
 
+    pub fn center(self) -> Pos3 {
+        Pos3::new(
+            (self.min.x + self.max.x) * 0.5,
+            (self.min.y + self.max.y) * 0.5,
+            (self.min.z + self.max.z) * 0.5,
+        )
+    }
+
     pub fn contains_point(self, pos: Pos3) -> bool {
         self.min.x <= pos.x
             && pos.x <= self.max.x
diff --git a/ray-tracing-pbrt-scene/src/bvh.rs b/ray-tracing-pbrt-scene/src/bvh.rs
new file mode 100644
index 0000000..2f5d0e4
--- /dev/null
+++ b/ray-tracing-pbrt-scene/src/bvh.rs
@@ -0,0 +1,182 @@
+use ray_tracing_core::{aabb::AABB, prelude::*};
+
+pub type Index = u32;
+
+pub trait BvhData {
+    type Intersect: AsRef<Float>;
+
+    // fn get_slice(&mut self) -> &mut [Self::Sortable];
+
+    fn len(&self) -> Index;
+
+    fn sort_range_dim(&mut self, id: std::ops::Range<Index>, dim: u8);
+
+    fn get_aabb(&self, id: Index) -> AABB;
+
+    fn get_aabb_range(&self, id: std::ops::Range<Index>) -> AABB {
+        id.map(|i| self.get_aabb(i))
+            .reduce(|a, b| a.extend(b.min()).extend(b.max()))
+            .unwrap()
+    }
+
+    fn intersect(&self, id: Index, ray: Ray, min: Float, max: Float) -> Option<Self::Intersect>;
+}
+
+#[derive(Debug, Clone, Copy)]
+enum Node {
+    Inner {
+        left: Index,
+        left_aabb: AABB,
+        right: Index,
+        right_aabb: AABB,
+    },
+    Leaf {
+        start: Index,
+        count: Index,
+    },
+}
+
+#[derive(Debug)]
+pub struct Bvh<T> {
+    data: T,
+    bvh: Vec<Node>,
+}
+
+impl<T: BvhData> Bvh<T> {
+    pub fn new(mut data: T, min: Index) -> Self {
+        let len = data.len();
+        let mut bvh = vec![Node::Leaf {
+            start: 0,
+            count: len,
+        }];
+
+        let aabb = data.get_aabb_range(0..len);
+
+        Self::build_bvh(&mut data, &mut bvh, 0, aabb, min);
+
+        Self { data, bvh }
+    }
+
+    fn build_bvh(data: &mut T, bvh: &mut Vec<Node>, node: usize, aabb: AABB, min: Index) {
+        let (start, count) = if let Node::Leaf { start, count } = bvh[node] {
+            (start, count)
+        } else {
+            unreachable!()
+        };
+
+        if count <= min {
+            return;
+        }
+
+        let size = aabb.size();
+
+        let dim = if size.x() > Float::max(size.y(), size.z()) {
+            0
+        } else if size.y() > Float::max(size.x(), size.z()) {
+            1
+        } else {
+            2
+        };
+
+        // sort
+        data.sort_range_dim(start..start + count, dim);
+
+        let mid = count / 2;
+
+        let left_aabb = data.get_aabb_range(start..(start + mid));
+        let right_aabb = data.get_aabb_range((start + mid)..(start + count));
+
+        let i = bvh.len() as Index;
+        bvh[node] = Node::Inner {
+            left: i,
+            left_aabb,
+            right: i + 1,
+            right_aabb,
+        };
+        bvh.push(Node::Leaf { start, count: mid });
+        bvh.push(Node::Leaf {
+            start: start + mid,
+            count: count - mid,
+        });
+
+        Self::build_bvh(data, bvh, i as usize, left_aabb, min);
+        Self::build_bvh(data, bvh, i as usize + 1, right_aabb, min);
+    }
+
+    fn intersect_inner(
+        &self,
+        node: Index,
+        ray: Ray,
+        min: Float,
+        max: Float,
+    ) -> Option<T::Intersect> {
+        match self.bvh[node as usize] {
+            Node::Inner {
+                left,
+                left_aabb,
+                right,
+                right_aabb,
+            } => {
+                let left_intersect = left_aabb.intersect_ray(ray, min, max);
+                let right_intersect = right_aabb.intersect_ray(ray, min, max);
+
+                match (left_intersect, right_intersect) {
+                    (None, None) => None,
+                    (None, Some(_)) => self.intersect_inner(right, ray, min, max),
+                    (Some(_), None) => self.intersect_inner(left, ray, min, max),
+                    (Some(l), Some(r)) => {
+                        let close;
+                        let far;
+                        if l < r {
+                            close = left;
+                            far = right;
+                        } else {
+                            close = right;
+                            far = left;
+                        }
+
+                        if let Some(close_intersect) = self.intersect_inner(close, ray, min, max) {
+                            if let Some(far_intersect) = self
+                                .intersect_inner(
+                                    far,
+                                    ray,
+                                    min,
+                                    Float::min(max, *close_intersect.as_ref()),
+                                )
+                                .filter(|far_intersect| {
+                                    *far_intersect.as_ref() < *close_intersect.as_ref()
+                                })
+                            {
+                                Some(far_intersect)
+                            } else {
+                                Some(close_intersect)
+                            }
+                        } else {
+                            self.intersect_inner(far, ray, min, max)
+                        }
+                    }
+                }
+            }
+            Node::Leaf { start, count } => {
+                let mut intersection = None;
+                for i in start..(start + count) {
+                    if let Some(t) = self.data.intersect(i, ray, min, max)
+                        && min <= *t.as_ref()
+                        && *t.as_ref() <= max
+                        && !intersection.as_ref().is_some_and(
+                            |old_t: &<T as BvhData>::Intersect| *t.as_ref() > *old_t.as_ref(),
+                        )
+                    {
+                        intersection.replace(t);
+                    }
+                }
+
+                intersection
+            }
+        }
+    }
+
+    pub fn intersect(&self, ray: Ray, min: Float, max: Float) -> Option<T::Intersect> {
+        self.intersect_inner(0, ray, min, max)
+    }
+}
diff --git a/ray-tracing-pbrt-scene/src/lib.rs b/ray-tracing-pbrt-scene/src/lib.rs
index a3b3f69..2344d78 100644
--- a/ray-tracing-pbrt-scene/src/lib.rs
+++ b/ray-tracing-pbrt-scene/src/lib.rs
@@ -1,7 +1,8 @@
 use crate::{
+    bvh::Bvh,
     either::Either,
     scene::PbrtScene,
-    shape::{Shape, ShapeAlpha, ShapeType},
+    shape::{Shape, ShapeAlpha, ShapeType, TriangleMesh},
     tokenizer::{Token, Tokenizer},
 };
 use material::PbrtMaterial;
@@ -23,6 +24,7 @@ use texture::PbrtTexture;
 #[macro_use]
 mod tokenizer;
 
+mod bvh;
 mod either;
 mod error;
 mod material;
@@ -151,14 +153,14 @@ fn parse_shape<R>(iter: &mut Tokenizer) -> Result<Statement<R>> {
                 n, Vec<Dir3>, Vec::new();
                 s, Vec<Dir3>, Vec::new();
                 uv, Vec<[Float; 2]>, Vec::new();
-                indices, Vec<usize>, Vec::new();
+                indices, Vec<[usize;3]>, Vec::new();
                 alpha, ShapeAlpha, ShapeAlpha::None
                 =>
                 p, "point3 P", iter.parse_list_3(Pos3::new)?;
                 n, "normal N", iter.parse_list_3(Dir3::new)?;
                 s, "normal S", iter.parse_list_3(Dir3::new)?;
                 uv, "point2 uv", iter.parse_list_2(|u, v| [u, v])?;
-                indices, "integer indices", iter.parse_list()?;
+                indices, "integer indices", iter.parse_list_3(|a, b, c|[a, b, c])?;
                 alpha, "float alpha", ShapeAlpha::Value(iter.parse_parameter()?);
                 alpha, "texture alpha", ShapeAlpha::Texture(iter.parse_parameter()?)
             );
@@ -170,13 +172,6 @@ fn parse_shape<R>(iter: &mut Tokenizer) -> Result<Statement<R>> {
                 bail!("Indices required for trianglemesh with more than 3 points.")
             }
 
-            if t.indices.len() % 3 != 0 {
-                bail!(
-                    "number of indices must be divisible by 3. num indices: {}",
-                    t.indices.len()
-                )
-            }
-
             if !t.n.is_empty() && t.n.len() != t.p.len() {
                 bail!("Number of normals not equal to number of positions.")
             }
@@ -189,13 +184,16 @@ fn parse_shape<R>(iter: &mut Tokenizer) -> Result<Statement<R>> {
             }
 
             Ok(Statement::Shape(
-                ShapeType::TriangleMesh {
-                    indices: t.indices,
-                    p: t.p,
-                    n: t.n,
-                    s: t.s,
-                    uv: t.uv,
-                },
+                ShapeType::TriangleMesh(Bvh::new(
+                    TriangleMesh {
+                        indices: t.indices,
+                        p: t.p,
+                        n: t.n,
+                        s: t.s,
+                        uv: t.uv,
+                    },
+                    8,
+                )),
                 t.alpha,
             ))
         }
@@ -863,7 +861,7 @@ fn inner_parse_pbrt<R: Rng + std::fmt::Debug>(
 
     // dbg!(context);
 
-    Ok(dbg!(pbrt))
+    Ok(pbrt)
 }
 
 pub fn parse_pbrt_v4<R: Rng + std::fmt::Debug>(
diff --git a/ray-tracing-pbrt-scene/src/shape.rs b/ray-tracing-pbrt-scene/src/shape.rs
index 36f676b..c2aff45 100644
--- a/ray-tracing-pbrt-scene/src/shape.rs
+++ b/ray-tracing-pbrt-scene/src/shape.rs
@@ -2,7 +2,13 @@ use std::sync::Arc;
 
 use ray_tracing_core::{affine_transform::AffineTransform, prelude::*, scene::Intersection};
 
-use crate::{AreaLight, either::Either, material::PbrtMaterial, scene::UVMaterial};
+use crate::{
+    AreaLight,
+    bvh::{Bvh, BvhData},
+    either::Either,
+    material::PbrtMaterial,
+    scene::UVMaterial,
+};
 
 #[derive(Debug)]
 #[allow(dead_code)]
@@ -21,6 +27,122 @@ pub(crate) enum ShapeAlpha {
     Texture(String),
 }
 
+pub(crate) struct InnerIntersect {
+    t: Float,
+    n: Dir3,
+    u: Float,
+    v: Float,
+}
+
+impl AsRef<Float> for InnerIntersect {
+    fn as_ref(&self) -> &Float {
+        &self.t
+    }
+}
+#[derive(Debug)]
+pub(crate) struct TriangleMesh {
+    pub(crate) indices: Vec<[usize; 3]>,
+    pub(crate) p: Vec<Pos3>,
+    pub(crate) n: Vec<Dir3>,
+    #[allow(dead_code)]
+    pub(crate) s: Vec<Dir3>,
+    pub(crate) uv: Vec<[Float; 2]>,
+}
+
+impl BvhData for TriangleMesh {
+    type Intersect = InnerIntersect;
+
+    fn len(&self) -> crate::bvh::Index {
+        self.indices.len() as crate::bvh::Index
+    }
+
+    fn sort_range_dim(&mut self, id: std::ops::Range<crate::bvh::Index>, dim: u8) {
+        let get_key =
+            |t: &[usize; 3]| t.iter().map(|&v| self.p[v][dim as usize]).sum::<Float>() / 3.0;
+        self.indices[id.start as usize..id.end as usize]
+            .sort_by(|a, b| get_key(a).partial_cmp(&get_key(b)).unwrap());
+    }
+
+    fn get_aabb(&self, id: crate::bvh::Index) -> AABB {
+        let v = self.indices[id as usize];
+
+        AABB::new(self.p[v[0]], self.p[v[1]]).extend(self.p[v[2]])
+    }
+
+    fn intersect(
+        &self,
+        id: crate::bvh::Index,
+        ray: Ray,
+        min: Float,
+        max: Float,
+    ) -> Option<Self::Intersect> {
+        let _ = max;
+        let _ = min;
+        let v = [
+            self.p[self.indices[id as usize][0]],
+            self.p[self.indices[id as usize][1]],
+            self.p[self.indices[id as usize][2]],
+        ];
+
+        let e1 = v[1] - v[0];
+        let e2 = v[2] - v[0];
+
+        let ray_cross_e2 = Dir3::cross(ray.dir(), e2);
+        let det = e1.dot(ray_cross_e2);
+
+        if det > -f32::EPSILON && det < f32::EPSILON {
+            return None; // This ray is parallel to this triangle.
+        }
+
+        let inv_det = 1.0 / det;
+        let s = ray.start() - v[0];
+        let u = inv_det * s.dot(ray_cross_e2);
+        if !(0.0..=1.0).contains(&u) {
+            return None;
+        }
+
+        let s_cross_e1 = s.cross(e1);
+        let v = inv_det * Dir3::dot(ray.dir(), s_cross_e1);
+        if v < 0.0 || u + v > 1.0 {
+            return None;
+        }
+        // At this stage we can compute t to find out where the intersection point is on the line.
+        let t = inv_det * e2.dot(s_cross_e1);
+
+        if t > Float::EPSILON {
+            // ray intersection
+            let w = 1.0 - u - v;
+            Some(InnerIntersect {
+                t,
+                n: if self.n.is_empty() {
+                    Dir3::cross(e1, e2)
+                } else {
+                    w * self.n[self.indices[id as usize][0]]
+                        + u * self.n[self.indices[id as usize][1]]
+                        + v * self.n[self.indices[id as usize][2]]
+                },
+                u: if self.uv.is_empty() {
+                    0.0
+                } else {
+                    w * self.uv[self.indices[id as usize][0]][0]
+                        + u * self.uv[self.indices[id as usize][1]][0]
+                        + v * self.uv[self.indices[id as usize][2]][0]
+                },
+                v: if self.uv.is_empty() {
+                    0.0
+                } else {
+                    w * self.uv[self.indices[id as usize][0]][1]
+                        + u * self.uv[self.indices[id as usize][1]][1]
+                        + v * self.uv[self.indices[id as usize][2]][1]
+                },
+            })
+        } else {
+            // This means that there is a line intersection but not a ray intersection.
+            None
+        }
+    }
+}
+
 #[derive(Debug)]
 #[allow(dead_code)]
 pub(crate) enum ShapeType {
@@ -30,13 +152,7 @@ pub(crate) enum ShapeType {
         zmax: Float,
         phimax: Float,
     },
-    TriangleMesh {
-        indices: Vec<usize>,
-        p: Vec<Pos3>,
-        n: Vec<Dir3>,
-        s: Vec<Dir3>,
-        uv: Vec<[Float; 2]>,
-    },
+    TriangleMesh(Bvh<TriangleMesh>),
     BilinearMesh {
         indices: Vec<usize>,
         p: Vec<Pos3>,
@@ -156,12 +272,7 @@ fn bilinear_intersection(
 }
 
 impl<R: Rng + std::fmt::Debug> Shape<R> {
-    fn inner_intersect(
-        &self,
-        ray: Ray,
-        min: Float,
-        max: Float,
-    ) -> Option<(Float, Dir3, Float, Float)> {
+    fn inner_intersect(&self, ray: Ray, min: Float, max: Float) -> Option<InnerIntersect> {
         match &self.obj {
             &ShapeType::Sphere {
                 radius,
@@ -213,11 +324,20 @@ impl<R: Rng + std::fmt::Debug> Shape<R> {
                         if phi < 0.0 {
                             phi += 2.0 * FloatConsts::PI;
                         }
-                        return Some((t, p.as_dir(), phi, Float::acos(p.z() / radius)));
+                        Some(InnerIntersect {
+                            t,
+                            n: p.as_dir(),
+                            u: phi,
+                            v: Float::acos(p.z() / radius),
+                        })
+                    } else {
+                        None
                     }
+                } else {
+                    None
                 }
             }
-            ShapeType::TriangleMesh { .. } => (),
+            ShapeType::TriangleMesh(bvh) => bvh.intersect(ray, min, max),
             ShapeType::BilinearMesh {
                 indices,
                 p,
@@ -228,15 +348,18 @@ impl<R: Rng + std::fmt::Debug> Shape<R> {
                     && let Some((t, u, v)) =
                         bilinear_intersection(ray, min, max, [p[0], p[1], p[2], p[3]])
                 {
-                    return Some((t, Dir3::new(0.0, 0.0, 1.0), u, v));
+                    Some(InnerIntersect {
+                        t,
+                        n: Dir3::new(0.0, 0.0, 1.0),
+                        u,
+                        v,
+                    })
+                } else {
+                    None
                 }
             }
-            ShapeType::LoopSubDiv { .. } => todo!(),
-            ShapeType::Disk { .. } => (),
-            ShapeType::PlyMesh { .. } => (),
+            _ => None,
         }
-
-        None
     }
     pub(crate) fn intersect(
         &'_ self,
@@ -247,10 +370,10 @@ impl<R: Rng + std::fmt::Debug> Shape<R> {
         let ray = self.ctm.transform_ray(ray);
 
         self.inner_intersect(ray, min, max)
-            .map(|(t, normal, u, v)| {
+            .map(|InnerIntersect { t, n, u, v }| {
                 Intersection::new(
                     t,
-                    self.ctm.inverse_transform_normal(normal).normalize(),
+                    self.ctm.inverse_transform_normal(n).normalize(),
                     self.material.get_a().map(|m| UVMaterial {
                         u,
                         v,
diff --git a/ray-tracing-tev/src/main.rs b/ray-tracing-tev/src/main.rs
index b9aa5fe..93b69eb 100644
--- a/ray-tracing-tev/src/main.rs
+++ b/ray-tracing-tev/src/main.rs
@@ -222,7 +222,13 @@ fn main() {
     }
 
     if let Some(pbrt_filename) = &args.pbrt_filename {
-        let pbrt = parse_pbrt_v4(pbrt_filename).unwrap();
+        let pbrt = match parse_pbrt_v4(pbrt_filename) {
+            Ok(pbrt) => pbrt,
+            Err(e) => {
+                println!("{e:?}");
+                std::process::exit(1);
+            }
+        };
 
         let c = BasicCamera::from_look_at(
             args.width,