diff --git a/Makefile b/Makefile
index b0814d7..afd67b5 100644
--- a/Makefile
+++ b/Makefile
@@ -1,9 +1,9 @@
 CXX = g++
 
-CXXFLAGS  = -std=c++20 -Wall -Wextra -Wno-unused-parameter -march=native -Ofast -flto -fopenmp
+CXXFLAGS  = -std=c++20 -Wall -Wextra -Wno-unused-parameter -march=native -O3 -flto -fopenmp
 
-DEPS = util.h vec3.h color.h ray.h camera.h hittable.h hittable_list.h sphere.h material.h lodepng.h moving_sphere.h
-OBJ = main.o material.o vec3.o lodepng.o moving_sphere.o
+DEPS = util.h vec3.h color.h ray.h camera.h hittable.h hittable_list.h sphere.h material.h lodepng.h moving_sphere.h bvh.h aabb.h
+OBJ = main.o material.o vec3.o lodepng.o moving_sphere.o bvh.o
 
 TARGET = toytracer
 
diff --git a/aabb.h b/aabb.h
new file mode 100644
index 0000000..5779c05
--- /dev/null
+++ b/aabb.h
@@ -0,0 +1,50 @@
+#ifndef AABB_H
+#define AABB_H
+
+#include <cmath>
+
+#include "vec3.h"
+
+class Aabb {
+  public:
+    Aabb() {}
+    Aabb(const Point3& a, const Point3& b) {
+        _min = a;
+        _max = b;
+    }
+
+    Point3 min() const { return _min; }
+    Point3 max() const { return _max; }
+
+    bool hit(const Ray& r, double tmin, double tmax) const {
+        for (int a = 0; a < 3; a++) {
+            auto invD = 1.0f / r.direction()[a];
+            auto t0 = (min()[a] - r.origin()[a]) * invD;
+            auto t1 = (max()[a] - r.origin()[a]) * invD;
+            if (invD < 0.0f)
+                std::swap(t0, t1);
+            tmin = t0 > tmin ? t0 : tmin;
+            tmax = t1 < tmax ? t1 : tmax;
+            if (tmax <= tmin)
+                return false;
+        }
+        return true;
+    }
+
+  private:
+    Point3 _min;
+    Point3 _max;
+};
+
+inline Aabb surrounding_box(Aabb box0, Aabb box1) {
+    Point3 small(std::fmin(box0.min().x(), box1.min().x()),
+                 std::fmin(box0.min().y(), box1.min().y()),
+                 std::fmin(box0.min().z(), box1.min().z()));
+
+    Point3 big(std::fmax(box0.max().x(), box1.max().x()),
+               std::fmax(box0.max().y(), box1.max().y()),
+               std::fmax(box0.max().z(), box1.max().z()));
+    return Aabb(small, big);
+}
+
+#endif // AABB_H
diff --git a/bvh.cpp b/bvh.cpp
new file mode 100644
index 0000000..c1f1838
--- /dev/null
+++ b/bvh.cpp
@@ -0,0 +1,67 @@
+#include "bvh.h"
+
+#include <algorithm>
+#include <functional>
+
+bool Bvh_node::bounding_box(double t0, double t1, Aabb& output_box) const {
+    output_box = box;
+    return true;
+}
+
+bool Bvh_node::hit(const Ray& r, double t_min, double t_max, hit_record& rec) const {
+    if (!box.hit(r, t_min, t_max))
+        return false;
+
+    bool hit_left = left->hit(r, t_min, t_max, rec);
+    bool hit_right = right->hit(r, t_min, hit_left ? rec.t : t_max, rec);
+
+    return hit_left || hit_right;
+}
+
+bool box_compare(const std::shared_ptr<Hittable> a, const std::shared_ptr<Hittable> b, int axis, double t0, double t1) {
+    Aabb box_a;
+    Aabb box_b;
+
+    if (!a->bounding_box(t0, t1, box_a) || !b->bounding_box(t0, t1, box_b)) {
+        std::cerr << "No bounding box in Bvh_node constructor.\n";
+        std::exit(1);
+    }
+
+    return box_a.min()[axis] < box_b.min()[axis];
+}
+
+Bvh_node::Bvh_node(std::vector<std::shared_ptr<Hittable>>& objects, int start, int end, double time0, double time1) {
+    using namespace std::placeholders;
+
+    int axis = random_int(0, 2);
+    auto comparator = std::bind(box_compare, _1, _2, axis, time0, time1);
+
+    int object_span = end - start;
+
+    if (object_span == 1) {
+        left = right = objects[start];
+    } else if (object_span == 2) {
+        if (comparator(objects[start], objects[start+1])) {
+            left = objects[start];
+            right = objects[start+1];
+        } else {
+            left = objects[start+1];
+            right = objects[start];
+        }
+    } else {
+        std::sort(objects.begin() + start, objects.begin() + end, comparator);
+
+        auto mid = start + object_span/2;
+        left = std::make_shared<Bvh_node>(objects, start, mid, time0, time1);
+        right = std::make_shared<Bvh_node>(objects, mid, end, time0, time1);
+    }
+
+    Aabb box_left, box_right;
+
+    if (!left->bounding_box(time0, time1, box_left) || !right->bounding_box(time0, time1, box_right)) {
+        std::cerr << "No bounding box in Bvh_node constructor.\n";
+        std::exit(1);
+    }
+
+    box = surrounding_box(box_left, box_right);
+}
diff --git a/bvh.h b/bvh.h
new file mode 100644
index 0000000..2c5e454
--- /dev/null
+++ b/bvh.h
@@ -0,0 +1,23 @@
+#ifndef BVH_H
+#define BVH_H
+
+#include "hittable.h"
+
+class Bvh_node : public Hittable {
+  public:
+    Bvh_node();
+
+    Bvh_node(
+            std::vector<std::shared_ptr<Hittable>>& objects,
+            int start, int end, double time0, double time1);
+
+    virtual bool hit(const Ray&, double tmin, double tmax, hit_record& rec) const;
+    virtual bool bounding_box(double t0, double t1, Aabb& output_box) const;
+
+  private:
+    std::shared_ptr<Hittable> left;
+    std::shared_ptr<Hittable> right;
+    Aabb box;
+};
+
+#endif // BVH_H
diff --git a/hittable.h b/hittable.h
index 6736158..caa95e1 100644
--- a/hittable.h
+++ b/hittable.h
@@ -5,6 +5,7 @@
 
 #include "ray.h"
 #include "material.h"
+#include "aabb.h"
 
 class Material;
 
@@ -24,6 +25,7 @@ struct hit_record {
 class Hittable {
 public:
   virtual bool hit(const Ray &r, double tmin, double tmax, hit_record &rec) const = 0;
+  virtual bool bounding_box(double t0, double t1, Aabb& output_box) const = 0;
 };
 
 #endif // HITTABLE_H
diff --git a/hittable_list.h b/hittable_list.h
index 12e4c4f..4ca29f8 100644
--- a/hittable_list.h
+++ b/hittable_list.h
@@ -5,6 +5,7 @@
 #include <vector>
 
 #include "hittable.h"
+#include "bvh.h"
 
 class Hittable_list : public Hittable {
 public:
@@ -15,9 +16,13 @@ public:
   void add(std::shared_ptr<Hittable> object) { objects.push_back(object); }
 
   virtual bool hit(const Ray &r, double tmin, double tmax, hit_record &rec) const;
+  bool bounding_box(double t0, double t1, Aabb& output_box) const;
+
+  Bvh_node generate_bvh(double t0, double t1);
 
 private:
   std::vector<std::shared_ptr<Hittable>> objects;
+
 };
 
 bool Hittable_list::hit(const Ray &r, double tmin, double tmax, hit_record &rec) const {
@@ -36,4 +41,22 @@ bool Hittable_list::hit(const Ray &r, double tmin, double tmax, hit_record &rec)
     return hit_anything;
 }
 
+bool Hittable_list::bounding_box(double t0, double t1, Aabb& output_box) const {
+    if (objects.empty()) return false;
+
+    Aabb temp_box;
+    bool first_box = true;
+
+    for (const auto& object : objects) {
+        if (!object->bounding_box(t0, t1, temp_box)) return false;
+        output_box = first_box ? temp_box : surrounding_box(output_box, temp_box);
+        first_box = false;
+    }
+    return true;
+}
+
+Bvh_node Hittable_list::generate_bvh(double t0, double t1) {
+    return Bvh_node(objects, 0, objects.size(), t0, t1);
+}
+
 #endif // HITTABLE_LIST_H
diff --git a/main.cpp b/main.cpp
index 47b5cac..aea607f 100644
--- a/main.cpp
+++ b/main.cpp
@@ -36,7 +36,7 @@ Color ray_color(const Ray& r, const Hittable& world, int depth) {
     return (1.0 - t) * Color(1.0, 1.0, 1.0) + t * Color(0.5, 0.7, 1.0);
 }
 
-Hittable_list setup_random_scene(const int sph_i) {
+Bvh_node setup_random_scene(const int sph_i) {
     Hittable_list world;
 
     auto ground_material = std::make_shared<Lambertian>(Color(0.5, 0.5, 0.5));
@@ -79,7 +79,7 @@ Hittable_list setup_random_scene(const int sph_i) {
     auto material3 = std::make_shared<Metal>(Color(0.7, 0.6, 0.5), 0.0);
     world.add(std::make_shared<Sphere>(Point3(4, 1, 0), 1.0, material3));
 
-    return world;
+    return world.generate_bvh(0, 1);
 }
 
 struct render_tile {
@@ -92,8 +92,8 @@ struct render_tile {
 int main() {
     const auto aspect_ratio = 16.0 / 9.0;
     //const int image_width = 1920;
-    const int image_width = 1280;
-    //const int image_width = 768;
+    //const int image_width = 1280;
+    const int image_width = 768;
     //const int image_width = 384;
     const int image_height = static_cast<int>(image_width / aspect_ratio);
     //const int samples_per_pixel = 2000;
diff --git a/moving_sphere.cpp b/moving_sphere.cpp
index 603f205..c947431 100644
--- a/moving_sphere.cpp
+++ b/moving_sphere.cpp
@@ -33,3 +33,14 @@ bool Moving_sphere::hit(const Ray& r, double tmin, double tmax, hit_record& rec)
     }
     return false;
 }
+
+bool Moving_sphere::bounding_box(double t0, double t1, Aabb& output_box) const {
+    Aabb box0(
+            center(t0) - Vec3(radius, radius, radius),
+            center(t0) + Vec3(radius, radius, radius));
+    Aabb box1(
+            center(t1) - Vec3(radius, radius, radius),
+            center(t1) + Vec3(radius, radius, radius));
+    output_box = surrounding_box(box0, box1);
+    return true;
+}
diff --git a/moving_sphere.h b/moving_sphere.h
index b45c089..86214e6 100644
--- a/moving_sphere.h
+++ b/moving_sphere.h
@@ -13,6 +13,7 @@ class Moving_sphere : public Hittable {
         : center0(cen0), center1(cen1), time0(t0), time1(t1), radius(r), mat_ptr(m) {}
 
     virtual bool hit(const Ray& r, double tmin, double tmax, hit_record& rec) const;
+    virtual bool bounding_box(double t0, double t1, Aabb& output_box) const;
 
     Point3 center(double time) const;
 
diff --git a/sphere.h b/sphere.h
index 2b22fc4..cfa846f 100644
--- a/sphere.h
+++ b/sphere.h
@@ -14,6 +14,7 @@ public:
   Sphere(Point3 cen, double r, std::shared_ptr<Material> m) : center(cen), radius(r), mat_ptr(m) {}
 
   virtual bool hit(const Ray& r, double tmin, double tmax, hit_record& rec) const;
+  virtual bool bounding_box(double t0, double t1, Aabb& output_box) const;
 
 private:
   Point3 center;
@@ -51,5 +52,12 @@ bool Sphere::hit(const Ray& r, double tmin, double tmax, hit_record& rec) const
     return false;
 }
 
+bool Sphere::bounding_box(double t0, double t1, Aabb& output_box) const {
+    output_box = Aabb(
+            center - Vec3(radius, radius, radius),
+            center + Vec3(radius, radius, radius));
+    return true;
+}
+
 
 #endif // SPHERE_H
diff --git a/util.h b/util.h
index e8ba046..ad69adc 100644
--- a/util.h
+++ b/util.h
@@ -22,6 +22,10 @@ inline double random_double(double min, double max) {
     return min + (max-min) * random_double();
 }
 
+inline int random_int(int min, int max) {
+    return static_cast<int>(random_double(min, max+1));
+}
+
 inline double clamp(double x, double min, double max) {
     if (x < min) return min;
     if (x > max) return max;