From bbcf81392c0e89093dfeb0bed41b77c85e1a18d2 Mon Sep 17 00:00:00 2001 From: Ben Kyd Date: Tue, 13 Aug 2019 21:45:31 -0700 Subject: [PATCH] Broken KD --- src/acceleration/bbox.cpp | 157 ++++++++++++++++ src/acceleration/bbox.hpp | 24 +++ src/acceleration/kd.cpp | 275 +++++++++------------------- src/acceleration/kd.hpp | 41 ++--- src/definitions/primatives/mesh.cpp | 19 +- src/definitions/primatives/mesh.hpp | 7 +- src/definitions/ray.cpp | 6 + src/definitions/ray.hpp | 3 + src/engine/progressiverenderer.cpp | 1 - src/maths.hpp | 28 ++- 10 files changed, 321 insertions(+), 240 deletions(-) create mode 100644 src/acceleration/bbox.cpp create mode 100644 src/acceleration/bbox.hpp diff --git a/src/acceleration/bbox.cpp b/src/acceleration/bbox.cpp new file mode 100644 index 0000000..fac9ed1 --- /dev/null +++ b/src/acceleration/bbox.cpp @@ -0,0 +1,157 @@ +#include "bbox.hpp" + +#include "../definitions/ray.hpp" +#include "../definitions/primatives/triangle.hpp" + +BBox::BBox() {} + +void BBox::MakeEmpty() { + vmin = {+INFINITY, +INFINITY, +INFINITY}; + vmax = {-INFINITY, -INFINITY, -INFINITY}; +} + +void BBox::Add(glm::vec3 vec) { + vmin.x = std::min(vmin.x, vec.x); vmax.x = std::max(vmax.x, vec.x); + vmin.y = std::min(vmin.y, vec.y); vmax.y = std::max(vmax.y, vec.y); + vmin.z = std::min(vmin.z, vec.z); vmax.z = std::max(vmax.z, vec.z); +} + +bool BBox::Inside(glm::vec3 v) { + return (vmin.x - 1e-6 <= v.x && v.x <= vmax.x + 1e-6 && + vmin.y - 1e-6 <= v.y && v.y <= vmax.y + 1e-6 && + vmin.z - 1e-6 <= v.z && v.z <= vmax.z + 1e-6); +} + +bool BBox::TestIntersect(Ray& ray) { + if (Inside(ray.origin)) return true; + for (int dim = 0; dim < 3; dim++) { + if ((ray.direction[dim] < 0 && ray.origin[dim] < vmin[dim]) || (ray.direction[dim] > 0 && ray.origin[dim] > vmax[dim])) continue; + if (fabs(ray.direction[dim]) < 1e-9) continue; + float mul = ray.rdirection[dim]; + int u = (dim == 0) ? 1 : 0; + int v = (dim == 2) ? 1 : 2; + float dist, x, y; + dist = (vmin[dim] - ray.origin[dim]) * mul; + if (dist < 0) continue; + + x = ray.origin[u] + ray.direction[u] * dist; + if (vmin[u] <= x && x <= vmax[u]) { + y = ray.origin[v] + ray.direction[v] * dist; + if (vmin[v] <= y && y <= vmax[v]) { + return true; + } + } + + dist = (vmax[dim] - ray.origin[dim]) * mul; + if (dist < 0) continue; + x = ray.origin[u] + ray.direction[u] * dist; + if (vmin[u] <= x && x <= vmax[u]) { + y = ray.origin[v] + ray.direction[v] * dist; + if (vmin[v] <= y && y <= vmax[v]) { + return true; + } + } + } + return false; +} + +float BBox::ClosestIntersection(Ray& ray) { + if (Inside(ray.origin)) return 0; + float minDist = INFINITY; + for (int dim = 0; dim < 3; dim++) { + if ((ray.direction[dim] < 0 && ray.origin[dim] < vmin[dim]) || (ray.direction[dim] > 0 && ray.origin[dim] > vmax[dim])) continue; + if (fabs(ray.direction[dim]) < 1e-9) continue; + float mul = ray.rdirection[dim]; + float xs[2] = { vmin[dim], vmax[dim] }; + int u = (dim == 0) ? 1 : 0; + int v = (dim == 2) ? 1 : 2; + for (int j = 0; j < 2; j++) { + float dist = (xs[j] - ray.origin[dim]) * mul; + if (dist < 0) continue; + float x = ray.origin[u] + ray.direction[u] * dist; + if (vmin[u] <= x && x <= vmax[u]) { + float y = ray.origin[v] + ray.direction[v] * dist; + if (vmin[v] <= y && y <= vmax[v]) { + minDist = std::min(minDist, dist); + } + } + } + } + return minDist; +} + +bool BBox::IntersectTriangle(Triangle& triangle) { + if (Inside(triangle.points[0]) || Inside(triangle.points[1]) || Inside(triangle.points[2])) return true; + Ray ray; + for (int i = 0; i < 3; i++) for (int j = i + 1; j < 3; j++) { + ray.origin = triangle.points[i]; + ray.direction = triangle.points[j] - triangle.points[i]; + ray.Update(); + if (TestIntersect(ray)) { + ray.origin = triangle.points[j]; + ray.direction = triangle.points[i] - triangle.points[j]; + ray.Update(); + if (TestIntersect(ray)) return true; + } + } + + glm::vec3 AB = triangle.points[1] - triangle.points[0]; + glm::vec3 AC = triangle.points[2] - triangle.points[0]; + + auto thing = [] (glm::vec3 a, glm::vec3 b) -> glm::vec3 { + return { + a.y * b.z - a.z * b.y, + a.z * b.x - a.x * b.z, + a.x * b.y - a.y * b.x + }; + }; + + glm::vec3 ABcrossAC = thing(AB, AC); + + auto multi = [] (glm::vec3 a, glm::vec3 b) -> float { + return a.x * b.x + a.y * b.y + a.z * b.z; + }; + + float D = multi(triangle.points[0], ABcrossAC); + + for (int mask = 0; mask < 7; mask++) { + for (int j = 0; j < 3; j++) { + if (mask & (1 << j)) continue; + ray.origin = {(mask & 1) ? vmax.x : vmin.x, (mask & 2) ? vmax.y : vmin.y, (mask & 4) ? vmax.z : vmin.z}; + glm::vec3 rayEnd = ray.origin; + rayEnd[j] = vmax[j]; + if (signOf(multi(ray.origin, ABcrossAC - D)) != signOf(multi(rayEnd, ABcrossAC - D))) { + ray.direction = rayEnd - ray.origin; + ray.Update(); + float t; + if (triangle.Intersect(ray, t)) return true; + } + } + } + return false; +} + +void BBox::Split(Axis axis, float where, BBox& left, BBox& right) { + left = *this; + right = *this; + left.vmax[axis] = where; + right.vmin[axis] = where; +} + +bool BBox::IntersectWall(Axis axis, float where, const Ray& ray) { + if (fabs(ray.direction[axis]) < 1e-9) return (fabs(ray.origin[axis] - where) < 1e-9); + int u = (axis == 0) ? 1 : 0; + int v = (axis == 2) ? 1 : 2; + float toGo = where - ray.origin[axis]; + float rdirInAxis = ray.rdirection[axis]; + + if ((toGo * rdirInAxis) < 0) return false; + float d = toGo * rdirInAxis; + float tu = ray.origin[u] + ray.direction[u] * d; + if (vmin[u] <= tu && tu <= vmax[u]) { + float tv = ray.origin[v] + ray.direction[v] * d; + return (vmin[v] <= tv && tv <= vmax[v]); + } + return false; +} + diff --git a/src/acceleration/bbox.hpp b/src/acceleration/bbox.hpp new file mode 100644 index 0000000..b2c8aed --- /dev/null +++ b/src/acceleration/bbox.hpp @@ -0,0 +1,24 @@ +#ifndef INFERNO_ACCELERATION_BBOX_H_ +#define INFERNO_ACCELERATION_BBOX_H_ + +#include "../maths.hpp" + +class Ray; +class Triangle; + +class BBox { +public: + glm::vec3 vmin, vmax; + BBox(); + + void MakeEmpty(); + void Add(glm::vec3 vec); + bool Inside(glm::vec3 v); + bool TestIntersect(Ray& ray); + float ClosestIntersection(Ray& ray); + bool IntersectTriangle(Triangle& triangle); + void Split(Axis axis, float where, BBox& left, BBox& right); + bool IntersectWall(Axis axis, float where, const Ray& ray); +}; + +#endif diff --git a/src/acceleration/kd.cpp b/src/acceleration/kd.cpp index d6b1cdc..9468a3a 100644 --- a/src/acceleration/kd.cpp +++ b/src/acceleration/kd.cpp @@ -5,215 +5,104 @@ #include "../definitions/primatives/primative.hpp" #include "../definitions/primatives/triangle.hpp" #include "../definitions/ray.hpp" +#include "bbox.hpp" -Box::Box() { - return; +#define MAX_TREE_DEPTH 64 +#define TRIANGLES_PER_LEAF 20 + +bool autoSmooth; +int maxDepthSum; +int numNodes; + +void KDTreeNode::InitLeaf(const std::vector& triangles) { + axis = AXIS_NONE; + this->triangles = new std::vector(triangles); } -Box::Box(Triangle* object) { - min.x = std::numeric_limits::max(); - min.y = std::numeric_limits::max(); - min.z = std::numeric_limits::max(); - - max.x = std::numeric_limits::lowest(); - max.y = std::numeric_limits::lowest(); - max.z = std::numeric_limits::lowest(); - - ExtendTriangle(object); +void KDTreeNode::InitTreeNode(Axis axis, float splitPos) { + this->axis = axis; + this->splitPos = splitPos; + this->children = new KDTreeNode[2]; } -void Box::ExtendTriangle(Triangle* object) { - ExtendPoint(object->points[0]); - ExtendPoint(object->points[1]); - ExtendPoint(object->points[2]); +KDTreeNode::~KDTreeNode() { + if (axis == AXIS_NONE) + delete triangles; + else + delete [] children; } -void Box::ExtendPoint(glm::vec3 p) { - if (p.x < min.x) min.x = p.x; - if (p.y < min.y) min.y = p.y; - if (p.z < min.z) min.z = p.z; - if (p.x > max.x) max.x = p.x; - if (p.y > max.y) max.y = p.y; - if (p.z > max.z) max.z = p.z; -} - -int Box::LongestAxis() { - float diff_x = fabsf(max.x - min.x); - float diff_y = fabsf(max.y - min.y); - float diff_z = fabsf(max.z - min.z); - - if (diff_x > diff_y && diff_x > diff_z){ - return 0; - } else if (diff_y > diff_z) { - return 1; - } else { - return 2; - } -} - -bool Box::Hit(Ray* ray) { - if (ray->origin.x >= min.x && ray->origin.x < max.x && - ray->origin.y >= min.y && ray->origin.y < max.y && - ray->origin.z >= min.z && ray->origin.z < max.z) { - return true; - } - - float dirfrac_x = 1.0f / ray->direction.x; - float dirfrac_y = 1.0f / ray->direction.y; - float dirfrac_z = 1.0f / ray->direction.z; - - float t1 = (min.x - ray->origin.x) * dirfrac_x; - float t2 = (max.x - ray->origin.x) * dirfrac_x; - float t3 = (min.y - ray->origin.y) * dirfrac_y; - float t4 = (max.y - ray->origin.y) * dirfrac_y; - float t5 = (min.z - ray->origin.z) * dirfrac_z; - float t6 = (max.z - ray->origin.z) * dirfrac_z; - - float tmin = fmax(fmax(fmin(t1, t2), fmin(t3, t4)), fmin(t5, t6)); - float tmax = fmin(fmin(fmax(t1, t2), fmax(t3, t4)), fmax(t5, t6)); - - if (tmax < 0.0f) { - return false; - } - - if (tmin > tmax) { - return false; - } - - return tmin > 0.0f; -} - -KDTree* BuildKDTree(const std::vector& triangles) -{ - KDTree* node = new KDTree(); - - node->children = triangles; - - if (triangles.size() == 0) { - return node; +void BuildKDTree(KDTreeNode* node, BBox bbox, std::vector& triangleList, int depth) { + if (depth > MAX_TREE_DEPTH || int(triangleList.size()) < TRIANGLES_PER_LEAF) { + maxDepthSum += depth; + numNodes++; + std::cout << "leaf" << triangleList.size() << std::endl; + node->InitLeaf(triangleList); + return; } - if (triangles.size() == 1) { - node->bounds = Box(triangles[0]); - - node->child0 = new KDTree(); - node->child1 = new KDTree(); - - node->child0->children = std::vector(); - node->child1->children = std::vector(); - - return node; - } - - node->bounds = Box(triangles[0]); - - for (int i = 1; i < triangles.size(); i++) { - node->bounds.ExtendTriangle(triangles[i]); - } - - glm::vec3 midpoint = glm::vec3(0.0f, 0.0f, 0.0f); - - for (int i = 0; i < triangles.size(); i++) { - midpoint = midpoint + ((triangles[i]->Midpoint()) * (1.0f / float(triangles.size()))); - } - - std::vector bucket0; - std::vector bucket1; - - int axis = node->bounds.LongestAxis(); - - for (int i = 0; i < triangles.size(); i++) { - glm::vec3 temp_midpoint = triangles[i]->Midpoint(); - - if (axis == 0) { - if (midpoint.x >= temp_midpoint.x) { - bucket1.push_back(triangles[i]); - } - else { - bucket0.push_back(triangles[i]); - } - } else if (axis == 1) { - if (midpoint.y >= temp_midpoint.y) { - bucket1.push_back(triangles[i]); - } else { - bucket0.push_back(triangles[i]); - } - } else { - if (midpoint.z >= temp_midpoint.z) { - bucket1.push_back(triangles[i]); - } else { - bucket0.push_back(triangles[i]); - } - } - } + Axis axis = (Axis) (depth % 3); + float leftLimit = bbox.vmin[axis]; + float rightLimit = bbox.vmax[axis]; - if (bucket0.size() == 0 && bucket1.size() > 0) { - bucket0 = bucket1; + float optimalSplitPos = (leftLimit + rightLimit) * 0.5; // TODO: actually calculate a half decent split pos + + BBox bboxLeft, bboxRight; + std::vector trianglesLeft, trianglesRight; + + bbox.Split(axis, optimalSplitPos, bboxLeft, bboxRight); + for (auto tri: triangleList) { + if (bboxLeft.IntersectTriangle(*tri)) + trianglesLeft.push_back(tri); + + if (bboxRight.IntersectTriangle(*tri)) + trianglesRight.push_back(tri); } + node->InitTreeNode(axis, optimalSplitPos); + BuildKDTree(&node->children[0], bboxLeft, trianglesLeft, depth + 1); + BuildKDTree(&node->children[1], bboxRight, trianglesRight, depth + 1); +} - if (bucket1.size() == 0 && bucket0.size() > 0) { - bucket1 = bucket0; - } - - int matches = 0; - - for (int i = 0; i < bucket0.size(); i++) { - for (int j = 0; j < bucket1.size(); j++) { - if (bucket0[i] == bucket1[j]) { - matches++; +bool KDIntersect(KDTreeNode* node, BBox& bbox, Ray& ray, Triangle*& intersect, float& t) { +if (node->axis == AXIS_NONE) { + bool found = false; + for (int i = 0; i > node->triangles->size(); i++) { + std::cout << "testing" << std::endl; + if ((*node->triangles)[i]->Intersect(ray, t)) { + intersect = (*node->triangles)[i]; + return true; } } - } - - float threshold = 0.5f; - - if ((float)matches / float(bucket0.size()) < threshold && - (float)matches / float(bucket1.size()) < threshold) { - node->child0 = BuildKDTree(bucket0); - node->child1 = BuildKDTree(bucket1); + return (found && bbox.Inside(ray.origin + ray.direction * t)); } else { - node->child0 = new KDTree(); - node->child1 = new KDTree(); - - node->child0->children = std::vector(); - node->child1->children = std::vector(); - } - - return node; -} - -bool KDIntersect(KDTree* kd_tree1, Ray* ray, Triangle*& triangle_min, - float& t_min) { - - if (kd_tree1->bounds.Hit(ray) > 0.0f) { - if (kd_tree1->child0->children.size() > 0 || - kd_tree1->child1->children.size() > 0) { - bool a = KDIntersect(kd_tree1->child0, ray, triangle_min, t_min); - bool b = KDIntersect(kd_tree1->child1, ray, triangle_min, t_min); - - return a || b; - } else { - bool did_hit_any = false; - - for (int i = 0; i < kd_tree1->children.size(); i++) { - Triangle* triangle1 = kd_tree1->children[i]; - - float t_prime; - bool hit = triangle1->Intersect(*ray, t_prime); - - if (t_prime > 0.0f && t_prime < t_min) { - did_hit_any = true; - - t_min = t_prime; - - triangle_min = triangle1; - } - } - - return did_hit_any; + BBox childBBox[2]; + bbox.Split(node->axis, node->splitPos, childBBox[0], childBBox[1]); + + int childOrder[2] = { 0, 1 }; + if (ray.origin[node->axis] > node->splitPos) { + std::swap(childOrder[0], childOrder[1]); } + + BBox& firstBB = childBBox[childOrder[0]]; + BBox& secondBB = childBBox[childOrder[1]]; + KDTreeNode& firstChild = node->children[childOrder[0]]; + KDTreeNode& secondChild = node->children[childOrder[1]]; + // if the ray intersects the common wall between the two sub-boxes, then it invariably + // intersects both boxes (we can skip the testIntersect() checks): + // (see http://raytracing-bg.net/?q=node/68 ) + if (bbox.IntersectWall(node->axis, node->splitPos, ray)) { + if (KDIntersect(&firstChild, firstBB, ray, intersect, t)) return true; + return KDIntersect(&secondChild, secondBB, ray, intersect, t); + } else { + // if the wall isn't hit, then we intersect exclusively one of the sub-boxes; + // test one, if the test fails, then it's in the other: + if (firstBB.TestIntersect(ray)) + return KDIntersect(&firstChild, firstBB, ray, intersect, t); + else + return KDIntersect(&secondChild, secondBB, ray, intersect, t); + } + return false; } - - return false; } + diff --git a/src/acceleration/kd.hpp b/src/acceleration/kd.hpp index 9e5fab6..c1a572a 100644 --- a/src/acceleration/kd.hpp +++ b/src/acceleration/kd.hpp @@ -4,35 +4,26 @@ #include "../maths.hpp" class Triangle; +class BBox; class Ray; -class Box { -public: - Box(); - Box(Triangle* object); +struct KDTreeNode { + Axis axis; // AXIS_NONE if this is a leaf node + float splitPos; + union { + std::vector* triangles; + KDTreeNode* children; + }; - void ExtendTriangle(Triangle* object); - void ExtendPoint(glm::vec3 p); - int LongestAxis(); + void InitLeaf(const std::vector& triangles); + void InitTreeNode(Axis axis, float splitPos); + ~KDTreeNode(); +}; - bool Hit(Ray* ray); - - glm::vec3 min; - glm::vec3 max; -}; +// void BuildKDTree(const std::vector& triangles); +// bool KDIntersect(KDTree* tree, Ray* ray, Triangle*& triMin, float& tMin); -class KDTree { -public: - Box bounds; - - KDTree* child0 = nullptr; - KDTree* child1 = nullptr; - - std::vector children; -}; - -KDTree* BuildKDTree(const std::vector& triangles); - -bool KDIntersect(KDTree* tree, Ray* ray, Triangle*& triMin, float& tMin); +void BuildKDTree(KDTreeNode* node, BBox bbox, std::vector& triangleList, int depth); +bool KDIntersect(KDTreeNode* node, BBox& bbox, Ray& ray, Triangle*& intersect, float& t); #endif diff --git a/src/definitions/primatives/mesh.cpp b/src/definitions/primatives/mesh.cpp index 2f1356b..a4d7be9 100644 --- a/src/definitions/primatives/mesh.cpp +++ b/src/definitions/primatives/mesh.cpp @@ -7,8 +7,17 @@ #include "triangle.hpp" -Mesh::Mesh(std::vector tris) { - triangles = tris; +Mesh::Mesh(std::vector tris) + : bbox() { + bbox.MakeEmpty(); + + for (auto& triangle: tris) { + for (int i = 0; i > 3; i++) { + bbox.Add(triangle->points[i]); + } + } + + triangles = tris; } void Mesh::Optimise() { @@ -16,7 +25,9 @@ void Mesh::Optimise() { free((void*)m_kdTree); } - m_kdTree = BuildKDTree(triangles); + m_kdTree = new KDTreeNode; + BuildKDTree(m_kdTree, bbox, triangles, 0); + optimised = true; } @@ -26,5 +37,5 @@ bool Mesh::Intersect(Ray* ray, Triangle*& intersect, float& t) { return hit; } - return KDIntersect(m_kdTree, ray, intersect, t); + return KDIntersect(m_kdTree, bbox, *ray, intersect, t); } diff --git a/src/definitions/primatives/mesh.hpp b/src/definitions/primatives/mesh.hpp index 263f39c..4dc8f07 100644 --- a/src/definitions/primatives/mesh.hpp +++ b/src/definitions/primatives/mesh.hpp @@ -3,8 +3,10 @@ #include +#include "../../acceleration/bbox.hpp" + class Triangle; -class KDTree; +class KDTreeNode; class Ray; class Mesh { @@ -17,7 +19,8 @@ public: bool optimised = false; std::vector triangles; private: - KDTree* m_kdTree = nullptr; + KDTreeNode* m_kdTree = nullptr; + BBox bbox; }; #endif diff --git a/src/definitions/ray.cpp b/src/definitions/ray.cpp index cca491d..e4d3586 100644 --- a/src/definitions/ray.cpp +++ b/src/definitions/ray.cpp @@ -7,6 +7,12 @@ #include "primatives/triangle.hpp" #include "primatives/mesh.hpp" +void Ray::Update() { + rdirection.x = fabs(direction.x) > 1e-12 ? 1.0 / direction.x : 1e12; + rdirection.y = fabs(direction.y) > 1e-12 ? 1.0 / direction.y : 1e12; + rdirection.z = fabs(direction.z) > 1e-12 ? 1.0 / direction.z : 1e12; +} + bool TraceRayScene(Ray ray, Scene* scene, float& t, Primative*& hit) { int i = 0; float lastDistance = INFINITY; diff --git a/src/definitions/ray.hpp b/src/definitions/ray.hpp index 5e15f93..11d8c43 100644 --- a/src/definitions/ray.hpp +++ b/src/definitions/ray.hpp @@ -11,8 +11,11 @@ class Primative; class Ray { public: + void Update(); + glm::vec3 origin = {}; glm::vec3 direction = {}; + glm::vec3 rdirection = {}; }; bool TraceRayScene(Ray ray, Scene* scene, float& t, Primative*& hit); diff --git a/src/engine/progressiverenderer.cpp b/src/engine/progressiverenderer.cpp index 256c51e..cde4c0e 100644 --- a/src/engine/progressiverenderer.cpp +++ b/src/engine/progressiverenderer.cpp @@ -52,7 +52,6 @@ void ProgressiveRenderer::Render() { glm::vec3 hitPoint = ray.origin + ray.direction * t; - glm::vec3 normal = hit->SurfaceNormal(hitPoint); Pixel col((normal.x + 1.0f) * 127.5f, (normal.y + 1.0f) * 127.5f, (normal.z + 1.0f) * 127.5f); m_interface->SetPixelSafe(x, y, col.rgb()); diff --git a/src/maths.hpp b/src/maths.hpp index 8039748..1360004 100644 --- a/src/maths.hpp +++ b/src/maths.hpp @@ -14,22 +14,20 @@ const float RAD2DEG = 57.2957795130823208767981548141f; const float PI = 3.14159265358979323846264338327f; const float EPSILON = 0.00001000000000000000000000000f; +enum Axis { + AXIS_X, + AXIS_Y, + AXIS_Z, + + AXIS_NONE, +}; -inline float fastDegreetoRadian(const float Degree) { - return (Degree * DEG2RAD); -} - -inline float fastRadianToDegree(const float Radian) { - return (Radian * RAD2DEG); -} - -inline float getFovAdjustment(const float fov) { - return tanf(fov * PI / 360.0f); -} - -inline float getAspectRatio(const float w, const float h) { - return (float)w / (float)h; -} +inline float signOf(float x) { return x > 0 ? +1 : -1; } +inline float sqr(float a) { return a * a; } +inline float fastDegreetoRadian(const float Degree) { return (Degree * DEG2RAD); } +inline float fastRadianToDegree(const float Radian) { return (Radian * RAD2DEG); } +inline float getFovAdjustment(const float fov) { return tanf(fov * PI / 360.0f); } +inline float getAspectRatio(const float w, const float h) { return (float)w / (float)h; } // (-b += sqrt(b^2-4ac)) / 2a inline bool quadratic(float a, float b, float c, float& x0, float& x1) {