Use SAH to generate BVH & fix ray triangle intersection

9 files changed, 175 insertions, 115 deletions

diff --git a/include/aabb.hpp b/include/aabb.hpp
index 2630a17..3a4e680 100644
--- a/include/aabb.hpp
+++ b/include/aabb.hpp
@@ -10,7 +10,7 @@ class aabb {
         aabb() = default;

         aabb(vec3 min_point, vec3 max_point): minimum(min_point), maximum(max_point) {}

         aabb(vec3 position, float radius): minimum(position - vec3(radius)), maximum(position + vec3(radius)) {}

-        aabb(point3 v1, point3 v2, point3 v3) {

+        aabb(point3 &v1, point3 &v2, point3 &v3) {

             minimum.x = fmin(fmin(v1.x, v2.x), v3.x);

             minimum.y = fmin(fmin(v1.y, v2.y), v3.y);

             minimum.z = fmin(fmin(v1.z, v2.z), v3.z);

@@ -20,20 +20,53 @@ class aabb {
             maximum.z = fmax(fmax(v1.z, v2.z), v3.z);

         }

 

-        static aabb surrounding_box(aabb box0, aabb box1) {

-            point3 min_point(fmin(box0.minimum.x, box1.minimum.x),

-                         fmin(box0.minimum.y, box1.minimum.y),

-                         fmin(box0.minimum.z, box1.minimum.z));

+        void union_with(const aabb &box) {

+            minimum.x = fmin(minimum.x, box.minimum.x);

+            minimum.y = fmin(minimum.y, box.minimum.y);

+            minimum.z = fmin(minimum.z, box.minimum.z);

 

-            point3 max_point(fmax(box0.maximum.x, box1.maximum.x),

-                       fmax(box0.maximum.y, box1.maximum.y),

-                       fmax(box0.maximum.z, box1.maximum.z));

+            maximum.x = fmax(maximum.x, box.maximum.x);

+            maximum.y = fmax(maximum.y, box.maximum.y);

+            maximum.z = fmax(maximum.z, box.maximum.z);

+        }

+

+        void union_with(const point3 &point) {

+            minimum.x = fmin(minimum.x, point.x);

+            minimum.y = fmin(minimum.y, point.y);

+            minimum.z = fmin(minimum.z, point.z);

+

+            maximum.x = fmax(maximum.x, point.x);

+            maximum.y = fmax(maximum.y, point.y);

+            maximum.z = fmax(maximum.z, point.z);

+        }

+

+        int32_t maximum_axis() {

+            vec3 diag = diagonal();

+

+            if (diag.x > diag.y && diag.x > diag.z) {

+                return 0;

+            } else if (diag.y > diag.z) {

+                return 1;

+            } else {

+                return 2;

+            }

+        }

+

+        vec3 diagonal() {

+            return maximum - minimum;

+        }

+

+        vec3 center() {

+            return diagonal() / 2.f;

+        }

 

-            return aabb(min_point, max_point);

+        float surface_area() {

+            auto diag = diagonal();

+            return 2.f * (diag.x * diag.y + diag.x * diag.z + diag.y * diag.z);

         }

 

         vec3 minimum = vec3(std::numeric_limits<float>::max());

-        vec3 maximum = vec3(std::numeric_limits<float>::min());

+        vec3 maximum = vec3(std::numeric_limits<float>::lowest());

 };

 

 #endif // !__AABB_HPP_

diff --git a/include/bvh-node.hpp b/include/bvh-node.hpp
index 4821fc7..72c4bdf 100644
--- a/include/bvh-node.hpp
+++ b/include/bvh-node.hpp
@@ -90,6 +90,8 @@ public:
             }

         }

 

+        // TODO: Pack nodes

+

         // exporter();

     }

 

@@ -146,36 +148,31 @@ public:
         return -1;

     }

 

+    std::vector<bvh_node>& nodes;

+    std::vector<temp_node> temp_nodes;

+    std::vector<temp_node> leafs;

+

+    // std::vector<sphere>& spheres;

+    std::vector<triangle>& triangles;

+

+    uint32_t next_node_id = 0;

+

 private:

 

+    static constexpr uint32_t buckets_count = 12;

+

     aabb compute_bounds(uint32_t begin, uint32_t end) {

         aabb global_box = aabb();

 

         for (uint32_t box_id = begin; box_id < end; box_id++) {

             auto leaf = leafs[box_id];

             // global_box = aabb::surrounding_box(global_box, spheres[leaf.primitive_id].bounding_box);

-            global_box = aabb::surrounding_box(global_box, triangles[leaf.primitive_id].bounding_box);

+            global_box.union_with(triangles[leaf.primitive_id].bounding_box);

         }

 

         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) {

@@ -192,22 +189,96 @@ private:
             temp_nodes[left_id] = leafs[begin];

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

         } else {

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

+            auto parent_bb = temp_nodes[parent_id].bounding_box;

+            int32_t split_axis = parent_bb.maximum_axis();

 

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

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

-            //         return spheres[a.primitive_id].center[split_axis] < spheres[b.primitive_id].center[split_axis];

-            //     }

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

-                    return triangles[a.primitive_id].center()[split_axis] < triangles[b.primitive_id].center()[split_axis];

+            if (count <= 4) {

+                // Subdivide the node by splitting the set in two equal parts

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

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

+                    //     return spheres[a.primitive_id].center[split_axis] < spheres[b.primitive_id].center[split_axis];

+                    // }

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

+                        return triangles[a.primitive_id].center()[split_axis] < triangles[b.primitive_id].center()[split_axis];

+                    }

+                );

+

+                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);

+            } else {

+                // Subdivide based on SAH

+                struct Bucket {

+                    uint32_t count = 0;

+                    aabb bb;

+                };

+                Bucket buckets[buckets_count];

+

+                // Split the set in N spatial buckets

+                // Simplify the optimal split axis search to just an axis between each bucket

+                aabb centroid_bounds;

+

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

+                    auto primitive = triangles[leafs[i].primitive_id];

+                    centroid_bounds.union_with(primitive.center());

                 }

-            );

 

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

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

+                float split_axis_size = centroid_bounds.maximum[split_axis] - centroid_bounds.minimum[split_axis];

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

+                    auto primitive = triangles[leafs[i].primitive_id];

+                    auto normalized_offset = (primitive.center()[split_axis] - centroid_bounds.minimum[split_axis]) / split_axis_size;

+                    uint32_t bucket_index = normalized_offset * buckets_count;

+                    if (bucket_index == buckets_count) {

+                        bucket_index -= 1;

+                    }

 

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

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

+                    buckets[bucket_index].count++;

+                    buckets[bucket_index].bb.union_with(primitive.bounding_box);

+                }

+

+                // Compute each potential split axis cost

+                float buckets_cost[buckets_count - 1] = { 0.f };

+                for (size_t bucket_index = 0; bucket_index < buckets_count - 1; bucket_index++) {

+                    aabb right_side_bb, left_side_bb;

+                    float right_side_count = 0, left_side_count = 0;

+

+                    for (size_t i = 0; i <= bucket_index; i++) {

+                        left_side_bb.union_with(buckets[i].bb);

+                        left_side_count += buckets[i].count;

+                    }

+

+                    for (size_t i = bucket_index + 1; i < buckets_count; i++) {

+                        right_side_bb.union_with(buckets[i].bb);

+                        right_side_count += buckets[i].count;

+                    }

+

+                    float left = left_side_count * left_side_bb.surface_area();

+                    float right = right_side_count * right_side_bb.surface_area();

+                    buckets_cost[bucket_index] = 0.125f + (left + right) / parent_bb.surface_area();

+                }

+

+                float min_cost = buckets_cost[0];

+                uint32_t min_cost_bucket = 0;

+                for (size_t bucket_index = 1; bucket_index < buckets_count - 1; bucket_index++) {

+                    if (buckets_cost[bucket_index] < min_cost) {

+                        min_cost = buckets_cost[bucket_index];

+                        min_cost_bucket = bucket_index;

+                    }

+                }

+

+                uint32_t middle_primitive_index = buckets[0].count;

+                for (size_t bucket_index = 1; bucket_index <= min_cost_bucket; bucket_index++) {

+                    middle_primitive_index += buckets[bucket_index].count;

+                }

+

+                temp_nodes[left_id].bounding_box = compute_bounds(begin, begin + middle_primitive_index);

+                temp_nodes[right_id].bounding_box = compute_bounds(begin + middle_primitive_index, end);

+

+                subdivide(left_id, begin, begin + middle_primitive_index);

+                subdivide(right_id, begin + middle_primitive_index, end);

+            }

         }

     }

 

@@ -227,16 +298,6 @@ private:
             depth_first_order(node.left_id + 1, next_id, new_id);

         }

     }

-

-public:

-    std::vector<bvh_node>& nodes;

-    std::vector<temp_node> temp_nodes;

-    std::vector<temp_node> leafs;

-

-    // std::vector<sphere>& spheres;

-    std::vector<triangle>& triangles;

-

-    uint32_t next_node_id = 0;

 };

 

 #endif // !__BVH_NODE_HPP_

diff --git a/include/gltf.hpp b/include/gltf.hpp
index 80f2749..0a34888 100644
--- a/include/gltf.hpp
+++ b/include/gltf.hpp
@@ -13,7 +13,7 @@ using json = nlohmann::json;
 struct mesh_part {

     std::vector<uint16_t> indices;

     std::vector<vec3> positions;

-    std::vector<vec3> normals;

+    // std::vector<vec3> normals;

 };

 

 struct mesh {

diff --git a/include/triangle.hpp b/include/triangle.hpp
index 57b42ca..f1af241 100644
--- a/include/triangle.hpp
+++ b/include/triangle.hpp
@@ -6,7 +6,6 @@
 

 struct vertex {

     vec3 position;

-    vec3 normal;

 };

 

 struct triangle {

@@ -16,7 +15,7 @@ struct triangle {
     {}

 

     vec3 center() {

-        return (v2.position + v1.position + v3.position) / 3.f;

+        return (v1.position + v2.position + v3.position) / 3.f;

     }

 

     vertex v1;

diff --git a/shaders/compute.comp b/shaders/compute.comp
index 060d213..dc82b91 100644
--- a/shaders/compute.comp
+++ b/shaders/compute.comp
@@ -43,7 +43,6 @@ struct sphere {
 
 struct vertex {
     vec4 position;
-    vec4 normal;
 };
 struct triangle {
     vertex v1;
@@ -218,18 +217,17 @@ bool hit_triangle(triangle tri, ray r, float min_t, float max_t, out hit_info in
     float v = d * dot(q, v2v1);
     float t = d * dot(-normal, originv1);
 
-    if (u < 0.0 || u > 1.0 || v < 0.0 || u + v > 1.0) {
+    if (u < 0.0 || u > 1.0 || v < 0.0 || u + v > 1.0 || t < min_t || t > max_t) {
         return false;
     }
 
     info.t = t;
-    // info.point = at(r, info.t);
+    info.point = at(r, info.t);
 
     vec3 outward_normal = normal;
-    info.front_face = dot(outward_normal, r.direction) < 0.0;
-    // info.normal = info.front_face ? outward_normal : -outward_normal;
-    info.normal = vec3(u, v, 0.0);
-    // info.mat = material(vec4(u, v, 0.0, 0.0), vec4(0.0), 1.0, 1.5, LAMBERTIAN_MATERIAL, 0);
+    // info.front_face = dot(outward_normal, r.direction) < 0.0;
+    info.normal = normal;
+    info.mat = material(vec4(u, v, 0.0, 0.0), vec4(0.0), 1.0, 1.5, LAMBERTIAN_MATERIAL, 0);
 
     return true;
 }
@@ -328,10 +326,10 @@ vec4 ray_color(ray r, uint seed) {
             // final_color += info.mat.emissive * throughput;
             throughput *= attenuation;
 
-            // float NdotR = dot(info.normal, scattered.direction);
-            // float PDF = NdotR / PI;
-            // vec4 BRDF = attenuation / PI;
-            // throughput = (throughput * NdotR / PDF) * BRDF;
+            float NdotR = dot(info.normal, scattered.direction);
+            float PDF = NdotR / PI;
+            vec4 BRDF = attenuation / PI;
+            throughput = (throughput * NdotR / PDF) * BRDF;
             current_ray = scattered;
 
             // Russian Roulette
@@ -372,7 +370,6 @@ void main() {
     vec4 proj_plane_pos = cam.first_pixel + jittered_uvs.x * cam.horizontal + jittered_uvs.y * cam.vertical;
 
     ray r = ray(cam.position.xyz + offset, proj_plane_pos.xyz - cam.position.xyz - offset);
-
     vec4 out_color = ray_color(r, seed);
 
     if (debug_bvh) {
diff --git a/src/gltf.cpp b/src/gltf.cpp
index 72219d6..055f2a6 100644
--- a/src/gltf.cpp
+++ b/src/gltf.cpp
@@ -5,6 +5,8 @@
 

 #include "utils.hpp"

 

+const vec3 scale (0.00800000037997961);

+

 mesh gltf::load(std::filesystem::path &path, std::string &filename) {

     std::fstream f(path / filename);

     json gltf_json;

@@ -61,7 +63,7 @@ void gltf::load_primitive(json &gltf_json, json &primitive, std::vector<std::vec
         auto indices_length = indices_view["byteLength"].get<size_t>();

 

         auto indices_buffer = buffers_content[indices_buffer_index];

-        void* indices_start = indices_buffer.data() + indices_offset;

+        uint16_t* indices_start = (uint16_t*)(indices_buffer.data() + indices_offset);

 

         mesh_part.indices.resize(indices_length / sizeof(uint16_t));

 

@@ -83,25 +85,24 @@ void gltf::load_primitive(json &gltf_json, json &primitive, std::vector<std::vec
         mesh_part.positions.resize(positions_length / (sizeof(float) * 3));

 

         for (size_t i = 0; i < positions_length / (sizeof(float) * 3); i++) {

-            mesh_part.positions[i] = vec3(positions_start[i * 3], positions_start[i * 3 + 1], positions_start[i * 3 + 2]);

+            mesh_part.positions[i] = vec3(positions_start[i * 3], positions_start[i * 3 + 1], positions_start[i * 3 + 2]) * scale;

         }

-        // std::memcpy(mesh_part.positions.data(), positions_start, positions_length);

     }

 

-    {

-        auto normals_accessor_index = primitive["attributes"]["NORMAL"].get<uint32_t>();

-        auto normals_view_index = accessors[normals_accessor_index]["bufferView"].get<uint32_t>();

-        auto normals_view = gltf_json["bufferViews"][normals_view_index];

+    // {

+    //     auto normals_accessor_index = primitive["attributes"]["NORMAL"].get<uint32_t>();

+    //     auto normals_view_index = accessors[normals_accessor_index]["bufferView"].get<uint32_t>();

+    //     auto normals_view = gltf_json["bufferViews"][normals_view_index];

 

-        auto normals_buffer_index = normals_view["buffer"].get<uint32_t>();

-        auto normals_offset = normals_view["byteOffset"].get<size_t>();

-        auto normals_length = normals_view["byteLength"].get<size_t>();

+    //     auto normals_buffer_index = normals_view["buffer"].get<uint32_t>();

+    //     auto normals_offset = normals_view["byteOffset"].get<size_t>();

+    //     auto normals_length = normals_view["byteLength"].get<size_t>();

 

-        auto normals_buffer = buffers_content[normals_buffer_index];

-        void* normals_start = normals_buffer.data() + normals_offset;

+    //     auto normals_buffer = buffers_content[normals_buffer_index];

+    //     void* normals_start = normals_buffer.data() + normals_offset;

 

-        mesh_part.normals.resize(normals_length / (sizeof(float) * 3));

+    //     mesh_part.normals.resize(normals_length / (sizeof(float) * 3));

 

-        std::memcpy(mesh_part.normals.data(), normals_start, normals_length);

-    }

+    //     std::memcpy(mesh_part.normals.data(), normals_start, normals_length);

+    // }

 }

diff --git a/src/main.cpp b/src/main.cpp
index 2634371..d0500e7 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -57,7 +57,7 @@ int main() {
 
     point3 position = { 13.0, 2.0, -3.0 };
     point3 target = { 0.0, 0.0, 0.0 };
-    auto v_fov = 20;
+    auto v_fov = 90;
     auto aperture = 0.1;
     auto focus_distance = 10;
     auto cam = camera(position, target, v_fov, (float)width / (float)height, aperture, focus_distance);
diff --git a/src/scene.cpp b/src/scene.cpp
index 8e25be7..b44abd5 100644
--- a/src/scene.cpp
+++ b/src/scene.cpp
@@ -80,9 +80,9 @@ scene::scene(camera cam, uint32_t width, uint32_t height)
             auto i2 = mesh_part.indices[i + 1];

             auto i3 = mesh_part.indices[i + 2];

 

-            auto v1 = vertex { .position = mesh_part.positions[i1], .normal = mesh_part.normals[i1] };

-            auto v2 = vertex { .position = mesh_part.positions[i2], .normal = mesh_part.normals[i2] };

-            auto v3 = vertex { .position = mesh_part.positions[i3], .normal = mesh_part.normals[i3] };

+            auto v1 = vertex { .position = mesh_part.positions[i1] };

+            auto v2 = vertex { .position = mesh_part.positions[i2] };

+            auto v3 = vertex { .position = mesh_part.positions[i3] };

 

             triangles.emplace_back(v1, v2, v3);

         }

@@ -94,7 +94,5 @@ scene::scene(camera cam, uint32_t width, uint32_t height)
     // bvh builder(spheres, nodes);

     // builder.exporter();

 

-    for (uint32_t idx = 0; idx < builder.nodes.size(); idx++) {

-        nodes[idx] = builder.nodes[idx];

-    }

+    std::memcpy(nodes.data(), builder.nodes.data(), sizeof(nodes[0]) * nodes.size());

 }

diff --git a/src/vec3.cpp b/src/vec3.cpp
index 8f7a720..166f681 100644
--- a/src/vec3.cpp
+++ b/src/vec3.cpp
@@ -139,32 +139,3 @@ vec3 refract(const vec3& uv, const vec3& n, float etai_over_etat) {
 std::ostream& operator<<(std::ostream& out, vec3 vec) {
     return out << vec.x << ' ' << vec.y << ' ' << vec.z << std::endl;
 }
-
-// TODO: Use something else to find random point onto a sphere
-vec3 random_in_unit_sphere() {
-    while (true) {
-        auto p = vec3::random(-1,1);
-        if (p.length_sq() >= 1) continue;
-        return p;
-    }
-}
-
-vec3 random_unit_vector() {
-    return random_in_unit_sphere().unit();
-}
-
-vec3 random_in_hemisphere(const vec3& normal) {
-    vec3 in_unit_sphere = random_in_unit_sphere();
-    if (in_unit_sphere.dot(normal) > 0.0) // In the same hemisphere as the normal
-        return in_unit_sphere;
-    else
-        return -in_unit_sphere;
-}
-
-vec3 random_in_unit_disk() {
-    while (true) {
-        auto p = vec3(randd(-1, 1), randd(-1, 1), 0);
-        if (p.length_sq() >= 1) continue;
-        return p;
-    }
-}