Depth first ordered bvh

3 files changed, 120 insertions, 100 deletions

diff --git a/include/bvh-node.hpp b/include/bvh-node.hpp
index 2affa38..16c3cb3 100644
--- a/include/bvh-node.hpp
+++ b/include/bvh-node.hpp
@@ -14,16 +14,22 @@ class bvh_node {
 

         aabb bounding_box;

 

-        int32_t left_id    = -1;

+        int32_t next_id    = -1;

         int32_t sphere_id   = -1;

         uint32_t padding[2];

 };

 

+struct temp_node : public bvh_node {

+    int32_t left_id;

+    int32_t df_id;

+};

+

 class bvh {

 public: 

     void build(sphere* spheres, uint32_t spheres_count) {

         nodes = std::vector<bvh_node>(2 * spheres_count - 1);

-        leafs = std::vector<bvh_node>(spheres_count);

+        temp_nodes = std::vector<temp_node>(2 * spheres_count - 1);

+        leafs = std::vector<temp_node>(spheres_count);

         this->spheres = spheres;

 

         for (uint32_t id = 0; id < spheres_count; id++) {

@@ -31,10 +37,26 @@ public:
             leafs[id].sphere_id = id;

         }

 

-        nodes[0] = bvh_node();

+        temp_nodes[0] = temp_node();

 

-        nodes[0].bounding_box = compute_bounds(0, spheres_count);

+        temp_nodes[0].bounding_box = compute_bounds(0, spheres_count);

         subdivide(0, 0, spheres_count);

+

+        set_depth_first_order();

+

+        for (uint32_t id = 0; id < temp_nodes.size(); id++) {

+            auto old_node = temp_nodes[id];

+            auto &node = nodes[old_node.df_id];

+

+            node.sphere_id = old_node.sphere_id;

+            node.bounding_box = old_node.bounding_box;

+

+            if (old_node.next_id != -1) {

+                node.next_id = temp_nodes[old_node.next_id].df_id;

+            }

+        }

+

+        // exporter();

     }

 

     //-------------------------

@@ -43,30 +65,53 @@ public:
     void exporter() {

         std::cout << "digraph bvh {" << std::endl;

 

-        std::cout << "\t 0";

-        traverse(nodes[0].left_id);

+        traverse_depth_first(temp_nodes, 0);

 

-        std::cout << "\t 0";

-        traverse(nodes[0].left_id + 1);

+        // traverse_depth_first_ordered(-1, 0);

 

         std::cout << "}" << std::endl;

     }

 

-    bool traverse(int32_t id) {

+    bool traverse_depth_first(std::vector<temp_node>& nodes, int32_t id) {

         if (id <= 0) {

             std::cout << ";" << std::endl;

             return false;

         }

 

         std::cout << " -> " << id;

-        if (traverse(nodes[id].left_id)) {

+        if (traverse_depth_first(nodes, nodes[id].left_id)) {

             std::cout << "\t " << id;

-            traverse(nodes[id].left_id + 1);

+            traverse_depth_first(nodes, nodes[id].left_id + 1);

         }

 

         return true;

     }

 

+    int32_t traverse_depth_first_ordered(std::vector<bvh_node>& nodes, int32_t parent_id, int32_t id) {

+        if (id == -1) {

+            return -1;

+        }

+

+        if (nodes[id].sphere_id == -1) {

+            int32_t left = id + 1;

+            int32_t next_id;

+

+            std::cout << "\t " << parent_id << " -> " << id << ";\n";

+            next_id = traverse_depth_first_ordered(nodes, id, left);

+

+            if (next_id != -1) {

+                return traverse_depth_first_ordered(nodes, id, next_id);

+            }

+        } else if (nodes[id].sphere_id != -1 && nodes[id].next_id != -1) {

+            std::cout << "\t " << parent_id << " -> " << id << ";\n";

+            return nodes[id].next_id;

+        } else {

+            std::cout << "\t " << parent_id << " -> " << id << ";\n";

+        }

+

+        return -1;

+    }

+

 private:

 

     aabb compute_bounds(uint32_t begin, uint32_t end) {

@@ -80,34 +125,39 @@ private:
         return global_box;

     }

 

+    int32_t find_split_axis(aabb bounding_box) {

+        int32_t split_axis = 0;

+        float size = 0;

+

+        // Find largest axis

+        for (size_t id = 0; id < 3; id++) {

+            auto current_size = bounding_box.maximum[id] - bounding_box.minimum[id];

+            if (current_size > size) {

+                size = current_size;

+                split_axis = id;

+            }

+        }

+

+        return split_axis;

+    }

+

     void subdivide(uint32_t parent_id, uint32_t begin, uint32_t end) {

         size_t count = end - begin;

          if (count == 1) {

-            nodes[parent_id] = leafs[begin];

+            temp_nodes[parent_id] = leafs[begin];

             return;

         }

 

         auto left_id = ++next_node_id;

         auto right_id = ++next_node_id;

 

-        nodes[parent_id].left_id = left_id;

+        temp_nodes[parent_id].left_id = left_id;

 

         if (count == 2) {

-            nodes[left_id] = leafs[begin];

-            nodes[right_id] = leafs[begin + 1];

+            temp_nodes[left_id] = leafs[begin];

+            temp_nodes[right_id] = leafs[begin + 1];

         } else {

-            auto bound = nodes[parent_id].bounding_box;

-            int split_axis = 0;

-            float size = 0;

-

-            // Find largest axis

-            for (size_t id = 0; id < 3; id++) {

-                auto current_size = bound.maximum[id] - bound.minimum[id];

-                if (current_size > size) {

-                    size = current_size;

-                    split_axis = id;

-                }

-            }

+            int32_t split_axis = find_split_axis(temp_nodes[parent_id].bounding_box);

 

             std::sort(leafs.begin() + begin, leafs.begin() + end,

                 [&](const bvh_node& a, const bvh_node& b) {

@@ -115,28 +165,35 @@ private:
                 }

             );

 

-            float split_axis_pos = (bound.maximum.x + bound.minimum.x) * 0.5f;

-            int mid = end - 1;

+            temp_nodes[left_id].bounding_box = compute_bounds(begin, begin + (count / 2));

+            temp_nodes[right_id].bounding_box = compute_bounds(begin + (count / 2), end);

+

+            subdivide(left_id, begin, begin + (count / 2));

+            subdivide(right_id, begin + (count / 2), end);

+        }

+    }

+

+    void set_depth_first_order() {

+        int32_t df_index = 0;

+        depth_first_order(0, -1, df_index);

+    }

 

-            for (size_t i = begin + 1; i < end; i++) {

-                auto leaf = leafs[i];

-                if (spheres[leaf.sphere_id].center[split_axis] >= split_axis_pos) {

-                    mid = i;

-                    break;

-                }

-            }

 

-            nodes[left_id].bounding_box = compute_bounds(begin, mid);

-            nodes[right_id].bounding_box = compute_bounds(mid, end);

+    void depth_first_order(int32_t id, int32_t next_id, int32_t &new_id) {

+        auto &node = temp_nodes[id];

+        node.next_id = next_id;

+        node.df_id = new_id++;

 

-            subdivide(left_id, begin, mid);

-            subdivide(right_id, mid, end);

+        if (node.left_id > 0) {

+            depth_first_order(node.left_id, node.left_id + 1, new_id);

+            depth_first_order(node.left_id + 1, next_id, new_id);

         }

     }

 

 public:

     std::vector<bvh_node> nodes;

-    std::vector<bvh_node> leafs;

+    std::vector<temp_node> temp_nodes;

+    std::vector<temp_node> leafs;

     sphere* spheres;

     uint32_t next_node_id = 0;

 };

diff --git a/shaders/compute.comp b/shaders/compute.comp
index 2c56fbf..2f6c5f5 100644
--- a/shaders/compute.comp
+++ b/shaders/compute.comp
@@ -60,7 +60,7 @@ struct camera {
 
 struct bvh_node {
     aabb bb;
-    int left;
+    int next_id;
     int sphere_id;
     uint[2] padding;
 };
@@ -121,26 +121,17 @@ vec2 disk_vec(inout uint seed) {
 }
 
 bool hit_aabb(aabb b, ray r, float min_t, float max_t) {
-    for (int axis_index = 0; axis_index < 3; axis_index++) {
-        float inv_dir = 1.f / r.direction[axis_index];
-        float t0 = (b.min[axis_index] - r.origin[axis_index]) * inv_dir;
-        float t1 = (b.max[axis_index] - r.origin[axis_index]) * inv_dir;
+    const vec3 invdir = 1.0 / r.direction;
+    const vec3 f = (b.max.xyz - r.origin.xyz) * invdir;
+    const vec3 n = (b.min.xyz - r.origin.xyz) * invdir;
 
-        if (inv_dir < 0.f) {
-            float temp = t0;
-            t0 = t1;
-            t1 = temp;
-        }
-
-        min_t = t0 > min_t ? t0 : min_t;
-        max_t = t1 < max_t ? t1 : max_t;
+    const vec3 tmax = max(f, n);
+    const vec3 tmin = min(f, n);
 
-        if (max_t <= min_t) {
-            return false;
-        }
-    }
+    const float t1 = min(tmax.x, min(tmax.y, tmax.z));
+    const float t0 = max(max(tmin.x, max(tmin.y, tmin.z)), 0.0f);
 
-    return true;
+    return t1 >= t0;
 }
 
 uint hit_aabbs(ray r, float min_t, float max_t) {
@@ -192,56 +183,28 @@ bool hit_sphere(sphere s, ray r, float min_t, float max_t, out hit_info info) {
 bool hit_node(ray r, float min_t, float max_t, out hit_info info) {
     hit_info temp_info;
     float closest_t = max_t;
-
-    if (!hit_aabb(nodes[0].bb, r, min_t, max_t)) {
-        return false;
-    }
-
-    int stack_id = -1;
-    int stack[512];
-    bvh_node current_node = nodes[0];
     bool hit = false;
+    int id = 0;
 
-    while(true) {
-        // Reached a leaf
-        if (current_node.sphere_id != -1) {
-            if (hit_sphere(spheres[current_node.sphere_id], r, min_t, closest_t, temp_info)) {
+    while(id != -1) {
+        // Leafs
+        if (nodes[id].sphere_id != -1) {
+            if (hit_sphere(spheres[nodes[id].sphere_id], r, min_t, closest_t, temp_info)) {
                 info = temp_info;
                 closest_t = temp_info.t;
                 hit = true;
             }
-
-            if (stack_id >= 0) {
-                current_node = nodes[stack[stack_id]];
-                stack_id--;
-                continue;
-            }
-
-            return hit;
-        }
-
-        bool hit_left = hit_aabb(nodes[current_node.left].bb, r, min_t, max_t);
-        bool hit_right = hit_aabb(nodes[current_node.left + 1].bb, r, min_t, max_t);
-
-        if (hit_left && hit_right) {
-            stack_id++;
-            stack[stack_id] = current_node.left;
-
-            current_node = nodes[current_node.left + 1];
-        } else if (hit_left) {
-            current_node = nodes[current_node.left];
-        } else if (hit_right) {
-            current_node = nodes[current_node.left + 1];
-        } else {
-            if (stack_id >= 0) {
-                current_node = nodes[stack[stack_id]];
-                stack_id--;
+            id = nodes[id].next_id;
+        } else { // Internal nodes
+            if (hit_aabb(nodes[id].bb, r, min_t, closest_t)) {
+                id++;
             } else {
-                return hit;
+                id = nodes[id].next_id;
             }
         }
     }
 
+    return hit;
 }
 
 bool hit_spheres(ray r, float min_t, float max_t, out hit_info info) {
@@ -351,7 +314,7 @@ void main() {
     vec4 out_color = ray_color(r, seed);
 
     if (debug_bvh) {
-        out_color -= hit_aabbs(r, 0.001, 1e15) * vec4(0.0, 0.05, 0.05, 0.0);
+        out_color -= hit_aabbs(r, 0.001, 1e15) * vec4(0.0, 0.025, 0.025, 0.0);
     }
 
     float spp_scale = (1.0 / sample_index);
diff --git a/src/utils.cpp b/src/utils.cpp
index a754b1d..a46cd6b 100644
--- a/src/utils.cpp
+++ b/src/utils.cpp
@@ -54,7 +54,7 @@ color sky_color { 0.5, 0.7, 1.0 };
 uint32_t random_scene(sphere *world) {
     material ground_material = { { 0.5, 0.5, 0.5 } };
     ground_material.type = MATERIAL_TYPE::LAMBERTIAN;
-    world[0] = sphere{ { 0,-1000,0 }, ground_material, 1000 };
+    world[0] = sphere{ { 0, -1000, 0 }, ground_material, 1000 };
 
     uint32_t world_sphere_index = 1;