From 0dafba17af044532fa2f14c61176e7994336f56c Mon Sep 17 00:00:00 2001 From: ladevieq Date: Sun, 1 Aug 2021 16:48:27 +0200 Subject: Use SAH to generate BVH & fix ray triangle intersection --- include/aabb.hpp | 53 +++++++++++++++---- include/bvh-node.hpp | 141 ++++++++++++++++++++++++++++++++++++--------------- include/gltf.hpp | 2 +- include/triangle.hpp | 3 +- shaders/compute.comp | 21 ++++---- src/gltf.cpp | 31 +++++------ src/main.cpp | 2 +- src/scene.cpp | 10 ++-- src/vec3.cpp | 29 ----------- 9 files changed, 176 insertions(+), 116 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::max()); - vec3 maximum = vec3(std::numeric_limits::min()); + vec3 maximum = vec3(std::numeric_limits::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& nodes; + std::vector temp_nodes; + std::vector leafs; + + // std::vector& spheres; + std::vector& 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); - - 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]; + auto parent_bb = temp_nodes[parent_id].bounding_box; + int32_t split_axis = parent_bb.maximum_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; + } + + buckets[bucket_index].count++; + buckets[bucket_index].bb.union_with(primitive.bounding_box); + } - subdivide(left_id, begin, begin + (count / 2)); - subdivide(right_id, begin + (count / 2), end); + // 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& nodes; - std::vector temp_nodes; - std::vector leafs; - - // std::vector& spheres; - std::vector& 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 indices; std::vector positions; - std::vector normals; + // std::vector 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(); 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(); - auto normals_view_index = accessors[normals_accessor_index]["bufferView"].get(); - auto normals_view = gltf_json["bufferViews"][normals_view_index]; + // { + // auto normals_accessor_index = primitive["attributes"]["NORMAL"].get(); + // auto normals_view_index = accessors[normals_accessor_index]["bufferView"].get(); + // auto normals_view = gltf_json["bufferViews"][normals_view_index]; - auto normals_buffer_index = normals_view["buffer"].get(); - auto normals_offset = normals_view["byteOffset"].get(); - auto normals_length = normals_view["byteLength"].get(); + // auto normals_buffer_index = normals_view["buffer"].get(); + // auto normals_offset = normals_view["byteOffset"].get(); + // auto normals_length = normals_view["byteLength"].get(); - 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; - } -} -- cgit v1.3