Add octree for faster mesh rendering

CMakeLists.txt

@@ -16,6 +16,8 @@ add_executable(raytracer
     include/image.h
     include/tiny_obj_loader.h
     include/scene.h
+    include/aabb.h
+    include/octree.h
     src/main.cpp
     src/scene.cpp)
 target_include_directories(raytracer PUBLIC include PRIVATE ${GLM_INCLUDE_DIR})

include/aabb.h

@@ -0,0 +1,51 @@
+#pragma once
+
+struct AABB {
+    glm::vec3 min, max;
+    
+    glm::vec3 center() const {
+        return 0.5f * (max + min);
+    }
+    
+    glm::vec3 extent() const {
+        return max - center();
+    }
+    
+    bool contains(const glm::vec3 point) const {
+        return glm::all(glm::lessThan(point, max)) && glm::all(glm::greaterThan(point, min));
+    }
+    
+    bool inside(const AABB extent) const {
+        return(max.x > extent.min.x &&
+               min.x < extent.max.x &&
+               max.y > extent.min.y &&
+               min.y < extent.max.y &&
+               max.z > extent.min.z &&
+               min.z < extent.max.z);
+    }
+    
+    bool contains(const Ray& ray) const {
+        const float t1 = (min.x - ray.origin.x) / ray.direction.x;
+        const float t2 = (max.x - ray.origin.x) / ray.direction.x;
+        const float t3 = (min.y - ray.origin.y) / ray.direction.y;
+        const float t4 = (max.y - ray.origin.y) / ray.direction.y;
+        const float t5 = (min.z - ray.origin.z) / ray.direction.z;
+        const float t6 = (max.z - ray.origin.z) / ray.direction.z;
+        
+        const float tmin = std::min(std::max(std::min(t1, t2), std::min(t3, t4)), std::min(t5, t6));
+        const float tmax = std::min(std::min(std::max(t1, t2), std::max(t3, t4)), std::max(t5, t6));
+        
+        // if tmax < 0, ray (line) is intersecting AABB, but whole AABB is behing us
+        if(tmax < 0)
+            return false;
+        
+        // if tmin > tmax, ray doesn't intersect AABB
+        if(tmin > tmax)
+            return false;
+        
+        if(tmin < 0.0f)
+            return true;
+
+        return true;
+    }
+};

include/octree.h

@@ -0,0 +1,89 @@
+#pragma once
+
+#include "aabb.h"
+
+constexpr int max_contained_types = 16;
+constexpr int max_octree_depth = 2;
+
+constexpr auto child_pattern = {
+    glm::vec3(+1, -1, -1),
+    glm::vec3(+1, -1, +1),
+    glm::vec3(+1, +1, -1),
+    glm::vec3(+1, +1, +1),
+    glm::vec3(-1, -1, -1),
+    glm::vec3(-1, -1, +1),
+    glm::vec3(-1, +1, -1),
+    glm::vec3(-1, +1, +1),
+};
+
+template<typename UnderlyingType>
+struct Node {
+    using NodeType = Node<UnderlyingType>;
+    
+    AABB extent;
+    int depth = 0;
+    
+    std::vector<UnderlyingType> contained_objects;
+    
+    std::array<std::unique_ptr<NodeType>, 8> children;
+    bool is_split = false;
+    
+    void add(UnderlyingType& object, const AABB extent) {
+        if(is_split) {
+            for(auto& child : children) {
+                if(extent.inside(child->extent))
+                    child->add(object, extent);
+            }
+        } else {
+            contained_objects.push_back(object);
+            
+            if(contained_objects.size() >= max_contained_types && depth < max_octree_depth)
+                split();
+        }
+    }
+    
+    void split() {
+        is_split = true;
+        
+        const auto center = extent.center();
+        const auto split_size = extent.extent().x / 2.0f;
+        
+        int i = 0;
+        for(auto& point : child_pattern) {
+            auto child = std::make_unique<NodeType>();
+            child->depth = depth + 1;
+            
+            const auto position = center + point * split_size;
+            child->extent.min = glm::vec3(position.x - split_size, position.y - split_size, position.z - split_size);
+            child->extent.max = glm::vec3(position.x + split_size, position.y + split_size, position.z + split_size);
+            
+            children[i++] = std::move(child);
+        }
+        
+        for(auto& object : contained_objects) {
+            for(auto& child : children) {
+                if(object.extent.inside(child->extent))
+                    child->add(object, extent);
+            }
+        }
+        
+        contained_objects.clear();
+    }
+};
+
+template<typename UnderlyingType>
+struct Octree {
+    using NodeType = Node<UnderlyingType>;
+    
+    NodeType root;
+    
+    Octree(const glm::vec3 min, const glm::vec3 max) {
+        root = NodeType();
+        root.extent.min = min;
+        root.extent.max = max;
+    }
+    
+    void add(UnderlyingType& object, const AABB extent) {
+        root.add(object, extent);
+    }
+};

include/scene.h

@@ -2,18 +2,32 @@
 
 #include <optional>
 #include <glm/glm.hpp>
+#include <array>
+#include <iostream>
 
 #include <tiny_obj_loader.h>
 
 #include "ray.h"
 #include "intersections.h"
 #include "lighting.h"
+#include "octree.h"
 
 constexpr glm::vec3 light_position = glm::vec3(5);
 constexpr float light_bias = 0.01f;
 constexpr int max_depth = 2;
 constexpr int num_indirect_samples = 4;
 
+struct TriangleBox {
+    uint64_t vertice_index = 0;
+    tinyobj::mesh_t* mesh;
+    AABB extent;
+};
+
+struct Object;
+
+glm::vec3 fetch_position(const Object& object, const tinyobj::mesh_t& mesh, const int32_t index, const int32_t vertex);
+glm::vec3 fetch_normal(const Object& object, const tinyobj::mesh_t& mesh, const int32_t index, const int32_t vertex);
+
 struct Object {
     glm::vec3 position = glm::vec3(0);
     glm::vec3 color = glm::vec3(1);
@@ -21,30 +35,62 @@ struct Object {
     tinyobj::attrib_t attrib;
     std::vector<tinyobj::shape_t> shapes;
     std::vector<tinyobj::material_t> materials;
-};
     
-struct Scene {
+    std::unique_ptr<Octree<TriangleBox>> octree;
-    std::vector<Object> objects;
 
-    Object& load_from_file(const std::string_view path) {
+    void create_octree() {
-        Object o;
+        octree = std::make_unique<Octree<TriangleBox>>(glm::vec3(-1), glm::vec3(1));
         
-        tinyobj::LoadObj(&o.attrib, &o.shapes, &o.materials, nullptr, path.data());
+        for(auto& shape : shapes) {
+            for(size_t i = 0; i < shape.mesh.num_face_vertices.size(); i++) {
+                const glm::vec3 v0 = fetch_position(*this, shape.mesh, i, 0);
+                const glm::vec3 v1 = fetch_position(*this, shape.mesh, i, 1);
+                const glm::vec3 v2 = fetch_position(*this, shape.mesh, i, 2);
                 
-        return objects.emplace_back(o);
+                AABB extent;
+                extent.min = glm::min(v0, v1);
+                extent.min = glm::min(extent.min, v2);
+                
+                extent.max = glm::max(v0, v1);
+                extent.max = glm::max(extent.max, v2);
+                
+                TriangleBox box = {};
+                box.vertice_index = i;
+                box.extent = extent;
+                box.mesh = &shape.mesh;
+                
+                octree->add(box, box.extent);
+            }
+        }
     }
 };
 
-glm::vec3 fetch_position(const Object& object, const tinyobj::mesh_t& mesh, const int32_t index, const int32_t vertex);
+struct Scene {
-glm::vec3 fetch_normal(const Object& object, const tinyobj::mesh_t& mesh, const int32_t index, const int32_t vertex);
+    std::vector<std::unique_ptr<Object>> objects;
+    
+    Object& load_from_file(const std::string_view path) {
+        auto o = std::make_unique<Object>();
+        
+        tinyobj::LoadObj(&o->attrib, &o->shapes, &o->materials, nullptr, path.data());
+      
+        return *objects.emplace_back(std::move(o));
+    }
+    
+    void generate_acceleration() {
+        for(auto& object : objects) {
+            object->create_octree();
+        }
+    }
+};
 
 struct HitResult {
     glm::vec3 position, normal;
-    Object object;
+    const Object* object = nullptr;
 };
 
 std::optional<HitResult> test_mesh(const Ray ray, const Object& object, const tinyobj::mesh_t& mesh, float& tClosest);
 std::optional<HitResult> test_scene(const Ray ray, const Scene& scene, float tClosest = std::numeric_limits<float>::infinity());
+std::optional<HitResult> test_scene_octree(const Ray ray, const Scene& scene, float tClosest = std::numeric_limits<float>::infinity());
 
 struct SceneResult {
     HitResult hit;
@@ -52,3 +98,4 @@ struct SceneResult {
 };
 
 std::optional<SceneResult> cast_scene(const Ray ray, const Scene& scene, const int depth = 0);
+

src/main.cpp

@@ -23,7 +23,7 @@
 #include <tiny_obj_loader.h>
 
 // scene information
-constexpr int32_t width = 128, height = 128;
+constexpr int32_t width = 256, height = 256;
 
 const Camera camera = [] {
     Camera camera;
@@ -178,6 +178,25 @@ void dump_to_file() {
     stbi_write_png("output.png", width, height, 3, pixels, width * 3);
 }
 
+template<typename UnderlyingType>
+void walk_node(Node<UnderlyingType>& node) {
+    if(ImGui::TreeNode(&node, "min: (%f %f %f)\n max: (%f %f %f)", node.extent.min.x, node.extent.min.y, node.extent.min.z, node.extent.max.x, node.extent.max.y, node.extent.max.z)) {
+        ImGui::Text("Is split: %i", node.is_split);
+        ImGui::Text("Contained triangles: %lu", node.contained_objects.size());
+
+        if(node.is_split) {
+            for(auto& child : node.children)
+                walk_node(*child);
+        }
+        
+        ImGui::TreePop();
+    }
+}
+
+void walk_object(Object& object) {
+    walk_node(object.octree->root);
+}
+
 int main(int, char*[]) {
     SDL_Init(SDL_INIT_VIDEO | SDL_INIT_TIMER);
     
@@ -224,6 +243,8 @@ int main(int, char*[]) {
                     auto& plane = scene.load_from_file("plane.obj");
                     plane.position.y = -1;
                     plane.color = {1, 0, 0};
+                    
+                    scene.generate_acceleration();
                 }
                 
                 ImGui::EndMenu();
@@ -243,6 +264,14 @@ int main(int, char*[]) {
             image_dirty = false;
         }
         
+        for(auto& object : scene.objects) {
+            if(ImGui::TreeNode("Object")) {
+                walk_object(*object);
+                
+                ImGui::TreePop();
+            }
+        }
+        
         ImGui::Render();
         
         auto& io = ImGui::GetIO();

src/scene.cpp

@@ -52,18 +52,91 @@ std::optional<HitResult> test_mesh(const Ray ray, const Object& object, const ti
         return {};
 }
 
+std::optional<HitResult> test_triangle(const Ray ray, const Object& object, const tinyobj::mesh_t& mesh, const size_t i, float& tClosest) {
+    bool intersection = false;
+    HitResult result = {};
+    
+    const glm::vec3 v0 = fetch_position(object, mesh, i, 0) + object.position;
+    const glm::vec3 v1 = fetch_position(object, mesh, i, 1) + object.position;
+    const glm::vec3 v2 = fetch_position(object, mesh, i, 2) + object.position;
+    
+    float t = std::numeric_limits<float>::infinity(), u, v;
+    if(intersections::ray_triangle(ray, v0, v1, v2, t, u, v)) {
+        if(t < tClosest && t > epsilon) {
+            const glm::vec3 n0 = fetch_normal(object, mesh, i, 0);
+            const glm::vec3 n1 = fetch_normal(object, mesh, i, 1);
+            const glm::vec3 n2 = fetch_normal(object, mesh, i, 2);
+            
+            result.normal = (1 - u - v) * n0 + u * n1 + v * n2;
+            result.position = ray.origin + ray.direction * t;
+            
+            tClosest = t;
+            
+            intersection = true;
+        }
+    }
+    
+    if(intersection)
+        return result;
+    else
+        return {};
+}
+
 std::optional<HitResult> test_scene(const Ray ray, const Scene& scene, float tClosest) {
     bool intersection = false;
     HitResult result = {};
     
     for(auto& object : scene.objects) {
-        for(uint32_t i = 0; i < object.shapes.size(); i++) {
+        for(uint32_t i = 0; i < object->shapes.size(); i++) {
-            auto mesh = object.shapes[i].mesh;
+            auto mesh = object->shapes[i].mesh;
             
-            if(const auto hit = test_mesh(ray, object, mesh, tClosest)) {
+            if(const auto hit = test_mesh(ray, *object, mesh, tClosest)) {
                 intersection = true;
                 result = hit.value();
-                result.object = object;
+                result.object = object.get();
+            }
+        }
+    }
+    
+    if(intersection)
+        return result;
+    else
+        return {};
+}
+
+template<typename T>
+Node<T>* find_hit_ray_search(Node<T>& node, const Ray ray, std::vector<Node<T>*>* out) {
+    if(node.extent.contains(ray)) {
+        if(node.is_split) {
+            for(auto& child : node.children)
+                find_hit_ray_search(*child, ray, out);
+        } else {
+            out->push_back(&node);
+        }
+    }
+}
+
+template<typename T>
+std::vector<Node<T>*> find_hit_ray(Node<T>& node, const Ray ray) {
+    std::vector<Node<T>*> vec;
+    
+    find_hit_ray_search(node, ray, &vec);
+    
+    return vec;
+}
+
+std::optional<HitResult> test_scene_octree(const Ray ray, const Scene& scene, float tClosest) {
+    bool intersection = false;
+    HitResult result = {};
+    
+    for(auto& object : scene.objects) {
+        for(auto& node : find_hit_ray(object->octree->root, ray)) {
+            for(auto& triangle_object : node->contained_objects) {
+                if(const auto hit = test_triangle(ray, *object, *triangle_object.mesh, triangle_object.vertice_index, tClosest)) {
+                    intersection = true;
+                    result = hit.value();
+                    result.object = object.get();
+                }
             }
         }
     }
@@ -75,7 +148,7 @@ std::optional<HitResult> test_scene(const Ray ray, const Scene& scene, float tCl
 }
 
 // methods adapated from https://users.cg.tuwien.ac.at/zsolnai/gfx/smallpaint/
-inline std::tuple<glm::vec3, glm::vec3> orthogonal_system(const glm::vec3& v1) {
+std::tuple<glm::vec3, glm::vec3> orthogonal_system(const glm::vec3& v1) {
     glm::vec3 v2;
     if(glm::abs(v1.x) > glm::abs(v1.y)) {
         // project to the y = 0 plane and construct a normalized orthogonal vector in this plane
@@ -101,7 +174,7 @@ std::optional<SceneResult> cast_scene(const Ray ray, const Scene& scene, const i
     if(depth > max_depth)
         return {};
     
-    if(auto hit = test_scene(ray, scene)) {
+    if(auto hit = test_scene_octree(ray, scene)) {
         const float diffuse = lighting::point_light(hit->position, light_position, hit->normal);
         
         // direct lighting calculation
@@ -112,7 +185,7 @@ std::optional<SceneResult> cast_scene(const Ray ray, const Scene& scene, const i
             
             const Ray shadow_ray(hit->position + (hit->normal * light_bias), light_dir);
             
-            const float shadow = test_scene(shadow_ray, scene) ? 0.0f : 1.0f;
+            const float shadow = test_scene_octree(shadow_ray, scene) ? 0.0f : 1.0f;
             
             direct = diffuse * shadow * glm::vec3(1);
         }
@@ -143,7 +216,7 @@ std::optional<SceneResult> cast_scene(const Ray ray, const Scene& scene, const i
         
         SceneResult result = {};
         result.hit = *hit;
-        result.color = (indirect + direct) * hit->object.color;
+        result.color = (indirect + direct) * hit->object->color;
         result.indirect = indirect;
         
         return result;