diff options
Diffstat (limited to 'core/aabbtree.cpp')
| -rw-r--r-- | core/aabbtree.cpp | 795 |
1 files changed, 795 insertions, 0 deletions
diff --git a/core/aabbtree.cpp b/core/aabbtree.cpp new file mode 100644 index 0000000..f7cb0c7 --- /dev/null +++ b/core/aabbtree.cpp @@ -0,0 +1,795 @@ +// This code contains NVIDIA Confidential Information and is disclosed to you +// under a form of NVIDIA software license agreement provided separately to you. +// +// Notice +// NVIDIA Corporation and its licensors retain all intellectual property and +// proprietary rights in and to this software and related documentation and +// any modifications thereto. Any use, reproduction, disclosure, or +// distribution of this software and related documentation without an express +// license agreement from NVIDIA Corporation is strictly prohibited. +// +// ALL NVIDIA DESIGN SPECIFICATIONS, CODE ARE PROVIDED "AS IS.". NVIDIA MAKES +// NO WARRANTIES, EXPRESSED, IMPLIED, STATUTORY, OR OTHERWISE WITH RESPECT TO +// THE MATERIALS, AND EXPRESSLY DISCLAIMS ALL IMPLIED WARRANTIES OF NONINFRINGEMENT, +// MERCHANTABILITY, AND FITNESS FOR A PARTICULAR PURPOSE. +// +// Information and code furnished is believed to be accurate and reliable. +// However, NVIDIA Corporation assumes no responsibility for the consequences of use of such +// information or for any infringement of patents or other rights of third parties that may +// result from its use. No license is granted by implication or otherwise under any patent +// or patent rights of NVIDIA Corporation. Details are subject to change without notice. +// This code supersedes and replaces all information previously supplied. +// NVIDIA Corporation products are not authorized for use as critical +// components in life support devices or systems without express written approval of +// NVIDIA Corporation. +// +// Copyright (c) 2013-2016 NVIDIA Corporation. All rights reserved. + +#include "aabbtree.h" + +#include "maths.h" +#include "platform.h" + +#include <algorithm> +#include <iostream> + +using namespace std; + +#if _WIN32 +_declspec (thread) uint32_t AABBTree::s_traceDepth; +#endif + +AABBTree::AABBTree(const Vec3* vertices, uint32_t numVerts, const uint32_t* indices, uint32_t numFaces) + : m_vertices(vertices) + , m_numVerts(numVerts) + , m_indices(indices) + , m_numFaces(numFaces) +{ + // build stats + m_treeDepth = 0; + m_innerNodes = 0; + m_leafNodes = 0; + + Build(); +} + +namespace +{ + + struct FaceSorter + { + FaceSorter(const Vec3* positions, const uint32_t* indices, uint32_t n, uint32_t axis) + : m_vertices(positions) + , m_indices(indices) + , m_numIndices(n) + , m_axis(axis) + { + } + + inline bool operator()(uint32_t lhs, uint32_t rhs) const + { + float a = GetCentroid(lhs); + float b = GetCentroid(rhs); + + if (a == b) + return lhs < rhs; + else + return a < b; + } + + inline float GetCentroid(uint32_t face) const + { + const Vec3& a = m_vertices[m_indices[face*3+0]]; + const Vec3& b = m_vertices[m_indices[face*3+1]]; + const Vec3& c = m_vertices[m_indices[face*3+2]]; + + return (a[m_axis] + b[m_axis] + c[m_axis])/3.0f; + } + + const Vec3* m_vertices; + const uint32_t* m_indices; + uint32_t m_numIndices; + uint32_t m_axis; + }; + + inline uint32_t LongestAxis(const Vector3& v) + { + if (v.x > v.y && v.x > v.z) + return 0; + else + return (v.y > v.z) ? 1 : 2; + } + +} // anonymous namespace + +void AABBTree::CalculateFaceBounds(uint32_t* faces, uint32_t numFaces, Vector3& outMinExtents, Vector3& outMaxExtents) +{ + Vector3 minExtents(FLT_MAX); + Vector3 maxExtents(-FLT_MAX); + + // calculate face bounds + for (uint32_t i=0; i < numFaces; ++i) + { + Vector3 a = Vector3(m_vertices[m_indices[faces[i]*3+0]]); + Vector3 b = Vector3(m_vertices[m_indices[faces[i]*3+1]]); + Vector3 c = Vector3(m_vertices[m_indices[faces[i]*3+2]]); + + minExtents = Min(a, minExtents); + maxExtents = Max(a, maxExtents); + + minExtents = Min(b, minExtents); + maxExtents = Max(b, maxExtents); + + minExtents = Min(c, minExtents); + maxExtents = Max(c, maxExtents); + } + + outMinExtents = minExtents; + outMaxExtents = maxExtents; +} + +// track current tree depth +static uint32_t s_depth = 0; + +void AABBTree::Build() +{ + assert(m_numFaces*3); + + //const double startTime = GetSeconds(); + + const uint32_t numFaces = m_numFaces; + + // build initial list of faces + m_faces.reserve(numFaces); + /* + for (uint32_t i=0; i < numFaces; ++i) + { + m_faces[i] = i; + } + */ + + // calculate bounds of each face and store + m_faceBounds.reserve(numFaces); + + std::vector<Bounds> stack; + for (uint32_t i=0; i < numFaces; ++i) + { + Bounds top; + CalculateFaceBounds(&i, 1, top.m_min, top.m_max); + + m_faces.push_back(i); + m_faceBounds.push_back(top); + /* + stack.push_back(top); + + while (!stack.empty()) + { + Bounds b = stack.back(); + stack.pop_back(); + + const float kAreaThreshold = 200.0f; + + if (b.GetSurfaceArea() < kAreaThreshold) + { + // node is good, append to our face list + m_faces.push_back(i); + m_faceBounds.push_back(b); + } + else + { + // split along longest axis + uint32_t a = LongestAxis(b.m_max-b.m_min); + + float splitPos = (b.m_min[a] + b.m_max[a])*0.5f; + Bounds left(b); + left.m_max[a] = splitPos; + + assert(left.GetSurfaceArea() < b.GetSurfaceArea()); + + + Bounds right(b); + right.m_min[a] = splitPos; + + assert(right.GetSurfaceArea() < b.GetSurfaceArea()); + + stack.push_back(left); + stack.push_back(right); + } + } + */ + } + + m_nodes.reserve(uint32_t(numFaces*1.5f)); + + // allocate space for all the nodes + m_freeNode = 1; + + // start building + BuildRecursive(0, &m_faces[0], numFaces); + + assert(s_depth == 0); + + + /* + const double buildTime = (GetSeconds()-startTime); + cout << "AAABTree Build Stats:" << endl; + cout << "Node size: " << sizeof(Node) << endl; + cout << "Build time: " << buildTime << "s" << endl; + cout << "Inner nodes: " << m_innerNodes << endl; + cout << "Leaf nodes: " << m_leafNodes << endl; + cout << "Alloc nodes: " << m_nodes.size() << endl; + cout << "Avg. tris/leaf: " << m_faces.size() / float(m_leafNodes) << endl; + cout << "Max depth: " << m_treeDepth << endl; + */ + // free some memory + FaceBoundsArray f; + m_faceBounds.swap(f); +} + +// partion faces around the median face +uint32_t AABBTree::PartitionMedian(Node& n, uint32_t* faces, uint32_t numFaces) +{ + FaceSorter predicate(&m_vertices[0], &m_indices[0], m_numFaces*3, LongestAxis(n.m_maxExtents-n.m_minExtents)); + std::nth_element(faces, faces+numFaces/2, faces+numFaces, predicate); + + return numFaces/2; +} + +// partion faces based on the surface area heuristic +uint32_t AABBTree::PartitionSAH(Node& n, uint32_t* faces, uint32_t numFaces) +{ + + /* + Vector3 mean(0.0f); + Vector3 variance(0.0f); + + // calculate best axis based on variance + for (uint32_t i=0; i < numFaces; ++i) + { + mean += 0.5f*(m_faceBounds[faces[i]].m_min + m_faceBounds[faces[i]].m_max); + } + + mean /= float(numFaces); + + for (uint32_t i=0; i < numFaces; ++i) + { + Vector3 v = 0.5f*(m_faceBounds[faces[i]].m_min + m_faceBounds[faces[i]].m_max) - mean; + v *= v; + variance += v; + } + + uint32_t bestAxis = LongestAxis(variance); + */ + + uint32_t bestAxis = 0; + uint32_t bestIndex = 0; + float bestCost = FLT_MAX; + + for (uint32_t a=0; a < 3; ++a) + //uint32_t a = bestAxis; + { + // sort faces by centroids + FaceSorter predicate(&m_vertices[0], &m_indices[0], m_numFaces*3, a); + std::sort(faces, faces+numFaces, predicate); + + // two passes over data to calculate upper and lower bounds + vector<float> cumulativeLower(numFaces); + vector<float> cumulativeUpper(numFaces); + + Bounds lower; + Bounds upper; + + for (uint32_t i=0; i < numFaces; ++i) + { + lower.Union(m_faceBounds[faces[i]]); + upper.Union(m_faceBounds[faces[numFaces-i-1]]); + + cumulativeLower[i] = lower.GetSurfaceArea(); + cumulativeUpper[numFaces-i-1] = upper.GetSurfaceArea(); + } + + float invTotalSA = 1.0f / cumulativeUpper[0]; + + // test all split positions + for (uint32_t i=0; i < numFaces-1; ++i) + { + float pBelow = cumulativeLower[i] * invTotalSA; + float pAbove = cumulativeUpper[i] * invTotalSA; + + float cost = 0.125f + (pBelow*i + pAbove*(numFaces-i)); + if (cost <= bestCost) + { + bestCost = cost; + bestIndex = i; + bestAxis = a; + } + } + } + + // re-sort by best axis + FaceSorter predicate(&m_vertices[0], &m_indices[0], m_numFaces*3, bestAxis); + std::sort(faces, faces+numFaces, predicate); + + return bestIndex+1; +} + +void AABBTree::BuildRecursive(uint32_t nodeIndex, uint32_t* faces, uint32_t numFaces) +{ + const uint32_t kMaxFacesPerLeaf = 6; + + // if we've run out of nodes allocate some more + if (nodeIndex >= m_nodes.size()) + { + uint32_t s = std::max(uint32_t(1.5f*m_nodes.size()), 512U); + + //cout << "Resizing tree, current size: " << m_nodes.size()*sizeof(Node) << " new size: " << s*sizeof(Node) << endl; + + m_nodes.resize(s); + } + + // a reference to the current node, need to be careful here as this reference may become invalid if array is resized + Node& n = m_nodes[nodeIndex]; + + // track max tree depth + ++s_depth; + m_treeDepth = max(m_treeDepth, s_depth); + + CalculateFaceBounds(faces, numFaces, n.m_minExtents, n.m_maxExtents); + + // calculate bounds of faces and add node + if (numFaces <= kMaxFacesPerLeaf) + { + n.m_faces = faces; + n.m_numFaces = numFaces; + + ++m_leafNodes; + } + else + { + ++m_innerNodes; + + // face counts for each branch + //const uint32_t leftCount = PartitionMedian(n, faces, numFaces); + const uint32_t leftCount = PartitionSAH(n, faces, numFaces); + const uint32_t rightCount = numFaces-leftCount; + + // alloc 2 nodes + m_nodes[nodeIndex].m_children = m_freeNode; + + // allocate two nodes + m_freeNode += 2; + + // split faces in half and build each side recursively + BuildRecursive(m_nodes[nodeIndex].m_children+0, faces, leftCount); + BuildRecursive(m_nodes[nodeIndex].m_children+1, faces+leftCount, rightCount); + } + + --s_depth; +} + +struct StackEntry +{ + uint32_t m_node; + float m_dist; +}; + + +#define TRACE_STATS 0 + +/* +bool AABBTree::TraceRay(const Vec3& start, const Vector3& dir, float& outT, float& outU, float& outV, float& outW, float& outSign, uint32_t& outIndex) const +{ +#if _WIN32 + // reset stats + s_traceDepth = 0; +#endif + + const Vector3 rcp_dir(1.0f/dir.x, 1.0f/dir.y, 1.0f/dir.z); + + // some temp variables + Vector3 normal; + float t, u, v, w, s; + + float minT, minU, minV, minW, minSign; + minU = minV = minW = minSign = minT = FLT_MAX; + + uint32_t minIndex = 0; + Vector3 minNormal; + + const uint32_t kStackDepth = 50; + + StackEntry stack[kStackDepth]; + stack[0].m_node = 0; + stack[0].m_dist = 0.0f; + + uint32_t stackCount = 1; + + while (stackCount) + { + // pop node from back + StackEntry& e = stack[--stackCount]; + + // ignore if another node has already come closer + if (e.m_dist >= minT) + { + continue; + } + + const Node* node = &m_nodes[e.m_node]; + +filth: + + if (node->m_faces == NULL) + { + +#if TRACE_STATS + extern uint32_t g_nodesChecked; + ++g_nodesChecked; +#endif + +#if _WIN32 + ++s_traceDepth; +#endif + // find closest node + const Node& leftChild = m_nodes[node->m_children+0]; + const Node& rightChild = m_nodes[node->m_children+1]; + + float dist[2] = {FLT_MAX, FLT_MAX}; + + IntersectRayAABBOmpf(start, rcp_dir, leftChild.m_minExtents, leftChild.m_maxExtents, dist[0]); + IntersectRayAABBOmpf(start, rcp_dir, rightChild.m_minExtents, rightChild.m_maxExtents, dist[1]); + + const uint32_t closest = dist[1] < dist[0]; // 0 or 1 + const uint32_t furthest = closest ^ 1; + + if (dist[furthest] < minT) + { + StackEntry& e = stack[stackCount++]; + e.m_node = node->m_children+furthest; + e.m_dist = dist[furthest]; + } + + // lifo + if (dist[closest] < minT) + { + node = &m_nodes[node->m_children+closest]; + goto filth; + } + } + else + { + for (uint32_t i=0; i < node->m_numFaces; ++i) + { + const uint32_t faceIndex = node->m_faces[i]; + const uint32_t indexStart = faceIndex*3; + + const Vec3& a = m_vertices[m_indices[indexStart+0]]; + const Vec3& b = m_vertices[m_indices[indexStart+1]]; + const Vec3& c = m_vertices[m_indices[indexStart+2]]; +#if TRACE_STATS + extern uint32_t g_trisChecked; + ++g_trisChecked; +#endif + + if (IntersectRayTriTwoSided(start, dir, a, b, c, t, u, v, w, s)) + { + if (t < minT && t > 0.01f) + { + minT = t; + minU = u; + minV = v; + minW = w; + minSign = s; + minIndex = faceIndex; + } + } + } + } + } + + // copy to outputs + outT = minT; + outU = minU; + outV = minV; + outW = minW; + outSign = minSign; + outIndex = minIndex; + + return (outT != FLT_MAX); +} +*/ + +bool AABBTree::TraceRay(const Vec3& start, const Vector3& dir, float& outT, float& u, float& v, float& w, float& faceSign, uint32_t& faceIndex) const +{ + //s_traceDepth = 0; + + Vector3 rcp_dir(1.0f/dir.x, 1.0f/dir.y, 1.0f/dir.z); + + outT = FLT_MAX; + TraceRecursive(0, start, dir, outT, u, v, w, faceSign, faceIndex); + + return (outT != FLT_MAX); +} + + +void AABBTree::TraceRecursive(uint32_t nodeIndex, const Vec3& start, const Vector3& dir, float& outT, float& outU, float& outV, float& outW, float& faceSign, uint32_t& faceIndex) const +{ + const Node& node = m_nodes[nodeIndex]; + + if (node.m_faces == NULL) + { +#if _WIN32 + ++s_traceDepth; +#endif + +#if TRACE_STATS + extern uint32_t g_nodesChecked; + ++g_nodesChecked; +#endif + + // find closest node + const Node& leftChild = m_nodes[node.m_children+0]; + const Node& rightChild = m_nodes[node.m_children+1]; + + float dist[2] = {FLT_MAX, FLT_MAX}; + + IntersectRayAABB(start, dir, leftChild.m_minExtents, leftChild.m_maxExtents, dist[0], NULL); + IntersectRayAABB(start, dir, rightChild.m_minExtents, rightChild.m_maxExtents, dist[1], NULL); + + uint32_t closest = 0; + uint32_t furthest = 1; + + if (dist[1] < dist[0]) + { + closest = 1; + furthest = 0; + } + + if (dist[closest] < outT) + TraceRecursive(node.m_children+closest, start, dir, outT, outU, outV, outW, faceSign, faceIndex); + + if (dist[furthest] < outT) + TraceRecursive(node.m_children+furthest, start, dir, outT, outU, outV, outW, faceSign, faceIndex); + + } + else + { + Vector3 normal; + float t, u, v, w, s; + + for (uint32_t i=0; i < node.m_numFaces; ++i) + { + uint32_t indexStart = node.m_faces[i]*3; + + const Vec3& a = m_vertices[m_indices[indexStart+0]]; + const Vec3& b = m_vertices[m_indices[indexStart+1]]; + const Vec3& c = m_vertices[m_indices[indexStart+2]]; +#if TRACE_STATS + extern uint32_t g_trisChecked; + ++g_trisChecked; +#endif + + if (IntersectRayTriTwoSided(start, dir, a, b, c, t, u, v, w, s)) + { + if (t < outT) + { + outT = t; + outU = u; + outV = v; + outW = w; + faceSign = s; + faceIndex = node.m_faces[i]; + } + } + } + } +} + +/* +bool AABBTree::TraceRay(const Vec3& start, const Vector3& dir, float& outT, Vector3* outNormal) const +{ + outT = FLT_MAX; + TraceRecursive(0, start, dir, outT, outNormal); + + return (outT != FLT_MAX); +} + +void AABBTree::TraceRecursive(uint32_t n, const Vec3& start, const Vector3& dir, float& outT, Vector3* outNormal) const +{ + const Node& node = m_nodes[n]; + + if (node.m_numFaces == 0) + { + extern _declspec(thread) uint32_t g_traceDepth; + ++g_traceDepth; +#if _DEBUG + extern uint32_t g_nodesChecked; + ++g_nodesChecked; +#endif + float t; + if (IntersectRayAABB(start, dir, node.m_minExtents, node.m_maxExtents, t, NULL)) + { + if (t <= outT) + { + TraceRecursive(n*2+1, start, dir, outT, outNormal); + TraceRecursive(n*2+2, start, dir, outT, outNormal); + } + } + } + else + { + Vector3 normal; + float t, u, v, w; + + for (uint32_t i=0; i < node.m_numFaces; ++i) + { + uint32_t indexStart = node.m_faces[i]*3; + + const Vec3& a = m_vertices[m_indices[indexStart+0]]; + const Vec3& b = m_vertices[m_indices[indexStart+1]]; + const Vec3& c = m_vertices[m_indices[indexStart+2]]; +#if _DEBUG + extern uint32_t g_trisChecked; + ++g_trisChecked; +#endif + + if (IntersectRayTri(start, dir, a, b, c, t, u, v, w, &normal)) + { + if (t < outT) + { + outT = t; + + if (outNormal) + *outNormal = normal; + } + } + } + } +} +*/ +bool AABBTree::TraceRaySlow(const Vec3& start, const Vector3& dir, float& outT, float& outU, float& outV, float& outW, float& faceSign, uint32_t& faceIndex) const +{ + const uint32_t numFaces = GetNumFaces(); + + float minT, minU, minV, minW, minS; + minT = minU = minV = minW = minS = FLT_MAX; + + Vector3 minNormal(0.0f, 1.0f, 0.0f); + + Vector3 n(0.0f, 1.0f, 0.0f); + float t, u, v, w, s; + bool hit = false; + uint32_t minIndex = 0; + + for (uint32_t i=0; i < numFaces; ++i) + { + const Vec3& a = m_vertices[m_indices[i*3+0]]; + const Vec3& b = m_vertices[m_indices[i*3+1]]; + const Vec3& c = m_vertices[m_indices[i*3+2]]; + + if (IntersectRayTriTwoSided(start, dir, a, b, c, t, u, v, w, s)) + { + if (t < minT) + { + minT = t; + minU = u; + minV = v; + minW = w; + minS = s; + minNormal = n; + minIndex = i; + hit = true; + } + } + } + + outT = minT; + outU = minU; + outV = minV; + outW = minW; + faceSign = minS; + faceIndex = minIndex; + + return hit; +} + +void AABBTree::DebugDraw() +{ + /* + glPolygonMode( GL_FRONT_AND_BACK, GL_LINE ); + + DebugDrawRecursive(0, 0); + + glPolygonMode( GL_FRONT_AND_BACK, GL_FILL ); + */ + +} + +void AABBTree::DebugDrawRecursive(uint32_t nodeIndex, uint32_t depth) +{ + static uint32_t kMaxDepth = 3; + + if (depth > kMaxDepth) + return; + + + /* + Node& n = m_nodes[nodeIndex]; + + Vector3 minExtents = FLT_MAX; + Vector3 maxExtents = -FLT_MAX; + + // calculate face bounds + for (uint32_t i=0; i < m_vertices.size(); ++i) + { + Vector3 a = m_vertices[i]; + + minExtents = Min(a, minExtents); + maxExtents = Max(a, maxExtents); + } + + + glBegin(GL_QUADS); + glVertex3f(minExtents.x, maxExtents.y, 0.0f); + glVertex3f(maxExtents.x, maxExtents.y, 0.0f); + glVertex3f(maxExtents.x, minExtents.y, 0.0f); + glVertex3f(minExtents.x, minExtents.y, 0.0f); + glEnd(); + + + n.m_center = Vec3(minExtents+maxExtents)/2; + n.m_extents = (maxExtents-minExtents)/2; + */ +/* + if (n.m_minEtextents != Vector3(0.0f)) + { + Vec3 corners[8]; + corners[0] = n.m_center + Vector3(-n.m_extents.x, n.m_extents.y, n.m_extents.z); + corners[1] = n.m_center + Vector3(n.m_extents.x, n.m_extents.y, n.m_extents.z); + corners[2] = n.m_center + Vector3(n.m_extents.x, -n.m_extents.y, n.m_extents.z); + corners[3] = n.m_center + Vector3(-n.m_extents.x, -n.m_extents.y, n.m_extents.z); + + corners[4] = n.m_center + Vector3(-n.m_extents.x, n.m_extents.y, -n.m_extents.z); + corners[5] = n.m_center + Vector3(n.m_extents.x, n.m_extents.y, -n.m_extents.z); + corners[6] = n.m_center + Vector3(n.m_extents.x, -n.m_extents.y, -n.m_extents.z); + corners[7] = n.m_center + Vector3(-n.m_extents.x, -n.m_extents.y, -n.m_extents.z); + + glBegin(GL_QUADS); + glColor3f(0.0f, 1.0f, 0.0f); + glVertex3fv(corners[0]); + glVertex3fv(corners[1]); + glVertex3fv(corners[2]); + glVertex3fv(corners[3]); + + glVertex3fv(corners[1]); + glVertex3fv(corners[5]); + glVertex3fv(corners[6]); + glVertex3fv(corners[2]); + + glVertex3fv(corners[0]); + glVertex3fv(corners[4]); + glVertex3fv(corners[5]); + glVertex3fv(corners[1]); + + glVertex3fv(corners[4]); + glVertex3fv(corners[5]); + glVertex3fv(corners[6]); + glVertex3fv(corners[7]); + + glVertex3fv(corners[0]); + glVertex3fv(corners[4]); + glVertex3fv(corners[7]); + glVertex3fv(corners[3]); + + glVertex3fv(corners[3]); + glVertex3fv(corners[7]); + glVertex3fv(corners[6]); + glVertex3fv(corners[2]); + + glEnd(); + + DebugDrawRecursive(nodeIndex*2+1, depth+1); + DebugDrawRecursive(nodeIndex*2+2, depth+1); + } + */ +} |