diff --git a/EBGeometry.hpp b/EBGeometry.hpp index 8392ed56..da0c3817 100644 --- a/EBGeometry.hpp +++ b/EBGeometry.hpp @@ -9,8 +9,6 @@ #include "Source/EBGeometry_DcelIterator.hpp" #include "Source/EBGeometry_DcelParser.hpp" #include "Source/EBGeometry_SignedDistanceFunction.hpp" -#include "Source/EBGeometry_SignedDistanceDcel.hpp" -#include "Source/EBGeometry_SignedDistanceBVH.hpp" #include "Source/EBGeometry_Union.hpp" #include "Source/EBGeometry_UnionBVH.hpp" diff --git a/Examples/AMReX/main.cpp b/Examples/AMReX/main.cpp index a02a1487..97f5134e 100644 --- a/Examples/AMReX/main.cpp +++ b/Examples/AMReX/main.cpp @@ -15,133 +15,175 @@ using namespace amrex; -// BVH tree degree -constexpr int K = 2; - namespace amrex { namespace EB2 { - using prec = float; - using Vec3 = EBGeometry::Vec3T; - using BV = EBGeometry::BoundingVolumes::AABBT; - using Node = EBGeometry::BVH::NodeT, BV, K>; - - - class SignedDistanceBVH : public EBGeometry::SignedDistanceBVH { + /*! + @brief This is an AMReX-capable version of the EBGeometry BVH accelerator. It is templated as T, BV, K which indicate + the EBGeometry precision, bounding volume, and tree degree. + */ + template + class SignedDistanceBVH { public: - + /*! + @brief Alias for builder node, for encapsulating a "standard" BVH node + */ + using BuilderNode = EBGeometry::BVH::NodeT, BV, K>; + + /*! + @brief Alias for linearized BVH node + */ + using LinearNode = EBGeometry::BVH::LinearBVH, BV, K>; + + /*! + @brief Alias for always-3D vector + */ + using Vec3 = EBGeometry::Vec3T; + + /*! + @brief Full constructor. + @param[in] a_filename File name. Must be a PLY file and will be parser by the PLY parser. + @param[in] a_flipSign Hook for swapping inside/outside. + */ SignedDistanceBVH(const std::string a_filename, const bool a_flipSign) { - auto mesh = EBGeometry::Dcel::Parser::PLY::readASCII(a_filename); + // 1. Read mesh from file. + auto mesh = EBGeometry::Dcel::Parser::PLY::readASCII(a_filename); - m_rootNode = std::make_shared (mesh->getFaces()); + // 2. Create a standard BVH hierarchy. This is not a compact ree. + auto root = std::make_shared(mesh->getFaces()); + root->topDownSortAndPartitionPrimitives(EBGeometry::Dcel::defaultBVConstructor, + EBGeometry::Dcel::spatialSplitPartitioner, + EBGeometry::Dcel::defaultStopFunction); - m_rootNode->topDownSortAndPartitionPrimitives(EBGeometry::Dcel::defaultBVConstructor, - EBGeometry::Dcel::spatialSplitPartitioner, - EBGeometry::Dcel::defaultStopFunction); + // 3. Flatten the tree onto a tighter memory representation. + m_rootNode = root->flattenTree(); } + /*! + @brief Copy constructor. + @param[in] a_other Other SDF. + */ SignedDistanceBVH(const SignedDistanceBVH& a_other) { this->m_rootNode = a_other.m_rootNode; this->m_flipSign = a_other.m_flipSign; - this->m_transformOps = a_other.m_transformOps; } + /*! + @brief AMReX's implicit function definition. + */ Real operator() (AMREX_D_DECL(Real x, Real y, Real z)) const noexcept { - const Real sign = (m_flipSign) ? -1.0 : 1.0; - return sign * m_rootNode->signedDistance(this->transformPoint(Vec3(x,y,z))); + return sign * m_rootNode->signedDistance(Vec3(x,y,z)); }; + /*! + @brief Also an AMReX implicit function implementation + */ inline Real operator() (const RealArray& p) const noexcept { return this->operator()(AMREX_D_DECL(p[0],p[1],p[2])); - } + } + + protected: + + /*! + @brief Root node of the linearized BVH hierarchy. + */ + std::shared_ptr m_rootNode; + + /*! + @brief Hook for flipping the sign + */ + bool m_flipSign; }; } } -int main (int argc, char* argv[]) -{ - amrex::Initialize(argc, argv); +int main (int argc, char* argv[]) { + amrex::Initialize(argc, argv); - { - int n_cell = 128; - int max_grid_size = 32; - int which_geom = 0; + int n_cell = 128; + int max_grid_size = 32; + int which_geom = 0; - std::string filename; + std::string filename; - // read parameters - ParmParse pp; - pp.query("n_cell", n_cell); - pp.query("max_grid_size", max_grid_size); - pp.query("which_geom", which_geom); + // read parameters + ParmParse pp; + pp.query("n_cell", n_cell); + pp.query("max_grid_size", max_grid_size); + pp.query("which_geom", which_geom); - Geometry geom; - { - RealBox rb; + Geometry geom; + { + RealBox rb; - if(which_geom == 0){// Airfoil case - rb = RealBox({-100,-100,-75}, {400,100,125}); - filename = "../PLY/airfoil.ply"; - } - else if (which_geom == 1){ // Sphere case - rb = RealBox({-400,-400,-400}, {400,400,400}); - filename = "../PLY/sphere.ply"; - } - else if (which_geom == 2){ // Dodecahedron - rb = RealBox({-2.,-2.,-2.}, {2.,2.,2.}); - filename = "../PLY/dodecahedron.ply"; - } - else if (which_geom == 3){ // Horse - rb = RealBox({-0.12,-0.12,-0.12}, {0.12,0.12,0.12}); - filename = "../PLY/horse.ply"; - } - else if (which_geom == 4){ // Car - // rb = RealBox({-20,-20,-20}, {20,20,20}); // Doesn't work. - rb = RealBox({-10,-5,-5}, {10,5,5}); // Works. - filename = "../PLY/porsche.ply"; - } - else if (which_geom == 5){ // Orion - rb = RealBox({-10,-5,-10}, {10,10,10}); - filename = "../PLY/orion.ply"; - } - else if (which_geom == 6){ // Armadillo - rb = RealBox({-100,-75,-100}, {100,125,100}); - filename = "../PLY/armadillo.ply"; - } - - Array is_periodic{false, false, false}; - Geometry::Setup(&rb, 0, is_periodic.data()); - Box domain(IntVect(0), IntVect(n_cell-1)); - geom.define(domain); - } - - // Create the signed distance function. - EB2::SignedDistanceBVH sdf(filename, false); - - auto gshop = EB2::makeShop(sdf); - EB2::Build(gshop, geom, 0, 0); - - // Put some data - MultiFab mf; - { - BoxArray ba(geom.Domain()); - ba.maxSize(max_grid_size); - DistributionMapping dm{ba}; - - std::unique_ptr factory - = amrex::makeEBFabFactory(geom, ba, dm, {2,2,2}, EBSupport::full); - - mf.define(ba, dm, 1, 0, MFInfo(), *factory); - mf.setVal(1.0); - } - - EB_WriteSingleLevelPlotfile("plt", mf, {"rho"}, geom, 0.0, 0); + if(which_geom == 0){// Airfoil case + rb = RealBox({-100,-100,-75}, {400,100,125}); + filename = "../PLY/airfoil.ply"; + } + else if (which_geom == 1){ // Sphere case + rb = RealBox({-400,-400,-400}, {400,400,400}); + filename = "../PLY/sphere.ply"; + } + else if (which_geom == 2){ // Dodecahedron + rb = RealBox({-2.,-2.,-2.}, {2.,2.,2.}); + filename = "../PLY/dodecahedron.ply"; + } + else if (which_geom == 3){ // Horse + rb = RealBox({-0.12,-0.12,-0.12}, {0.12,0.12,0.12}); + filename = "../PLY/horse.ply"; + } + else if (which_geom == 4){ // Car + // rb = RealBox({-20,-20,-20}, {20,20,20}); // Doesn't work. + rb = RealBox({-10,-5,-5}, {10,5,5}); // Works. + filename = "../PLY/porsche.ply"; } + else if (which_geom == 5){ // Orion + rb = RealBox({-10,-5,-10}, {10,10,10}); + filename = "../PLY/orion.ply"; + } + else if (which_geom == 6){ // Armadillo + rb = RealBox({-100,-75,-100}, {100,125,100}); + filename = "../PLY/armadillo.ply"; + } + + Array is_periodic{false, false, false}; + Geometry::Setup(&rb, 0, is_periodic.data()); + Box domain(IntVect(0), IntVect(n_cell-1)); + geom.define(domain); + } + + // Create our signed distance function. K is the tree degree while T is the EBGeometry precision. + constexpr int K = 2; + + using T = float; + using Vec3 = EBGeometry::Vec3T; + using BV = EBGeometry::BoundingVolumes::AABBT; + + EB2::SignedDistanceBVH sdf(filename, false); + + auto gshop = EB2::makeShop(sdf); + EB2::Build(gshop, geom, 0, 0); + + // Put some data + MultiFab mf; + { + BoxArray ba(geom.Domain()); + ba.maxSize(max_grid_size); + DistributionMapping dm{ba}; + + std::unique_ptr factory + = amrex::makeEBFabFactory(geom, ba, dm, {2,2,2}, EBSupport::full); + + mf.define(ba, dm, 1, 0, MFInfo(), *factory); + mf.setVal(1.0); + } + + EB_WriteSingleLevelPlotfile("plt", mf, {"rho"}, geom, 0.0, 0); - amrex::Finalize(); + amrex::Finalize(); } diff --git a/Examples/Basic/README.md b/Examples/Basic/README.md index a015a6be..18a87c37 100644 --- a/Examples/Basic/README.md +++ b/Examples/Basic/README.md @@ -1,4 +1,11 @@ -This folder contains a basic example of using EBGeometry. +This folder contains a basic example of using EBGeometry, with three different representations of the signed distance field. + +* A naive approach which iterates through all facets and computed the signed distance. +* Using a conventional bounding volume hierarchy with (references to) primitives stored directly in the nodes. +* A flatten, more compact bounding volume hierarchy. + +Note that SDF queries have different complexity for different geometries and input points. +For example, a tessellated sphere has a "blind spot" in it's center where even BVHs will visit most, if not all, primitives. Compiling --------- diff --git a/Examples/Basic/main.cpp b/Examples/Basic/main.cpp index 29e91005..c140793f 100644 --- a/Examples/Basic/main.cpp +++ b/Examples/Basic/main.cpp @@ -5,14 +5,16 @@ #include #include +#include #include "../../EBGeometry.hpp" using namespace EBGeometry; using namespace EBGeometry::Dcel; -// Degree of bounding volume hierarchies. -constexpr int K = 2; +// Degree of bounding volume hierarchies. We use a 4-ary tree here, where each +// regular node has four children. +constexpr int K = 4; int main(int argc, char *argv[]) { @@ -20,53 +22,67 @@ int main(int argc, char *argv[]) { std::vector all_args; std::string file; - + + // Read the input file. if (argc == 2) { file = "../PLY/" + std::string(argv[1]); } else{ std::cerr << "Missing file name. Use ./a.out 'filename' where 'filename' is one of the files in ../PLY\n"; - - return 1; } - // Declare precision. - using precision = float; + // Declare the precision T as float. + using T = float; - // Aliases for cutting down on typing. - using BV = BoundingVolumes::AABBT; - using Vec3 = Vec3T; - using Face = FaceT; - using SDF = SignedDistanceFunction; - using slowSDF = SignedDistanceDcel; - using fastSDF = SignedDistanceBVH; + // Aliases for cutting down on typing. + using BV = BoundingVolumes::AABBT; + using Vec3 = Vec3T; - std::string inputFile; - - // Create an empty DCEL mesh + // Parse the mesh from file. One can call the signed distance function directly on the mesh, but it will + // iterate through every facet. std::cout << "Parsing input file\n"; - auto mesh = EBGeometry::Dcel::Parser::PLY::readASCII(file); - - // Create a signed distance function from the mesh. This is the object - // that will iterate through each and every facet in the input mesh. - auto slow = std::make_shared(mesh, false); + std::shared_ptr > directSDF = EBGeometry::Dcel::Parser::PLY::readASCII(file); // Create a bounding-volume hierarchy of the same mesh type. We begin by create the root node and supplying all the mesh faces to it. Here, // our bounding volume hierarchy bounds the facets in a binary tree. - auto root = std::make_shared > (mesh->getFaces()); - - std::cout << "Partitioning BVH\n"; - root->topDownSortAndPartitionPrimitives(EBGeometry::Dcel::defaultBVConstructor, - EBGeometry::Dcel::spatialSplitPartitioner, - EBGeometry::Dcel::defaultStopFunction); - - - auto fast = std::make_shared(root, false); - - - // Query the distance to a point. - std::cout << "Distance to point using direct method = " << (*slow)(Vec3::one()) << std::endl; - std::cout << "Distance to point using bounding volumes = " << (*fast)(Vec3::one()) << std::endl; + std::cout << "Partitioning BVH\n"; + auto bvhSDF = std::make_shared, BV, K> > (directSDF->getFaces()); + bvhSDF->topDownSortAndPartitionPrimitives(EBGeometry::Dcel::defaultBVConstructor, + EBGeometry::Dcel::spatialSplitPartitioner, + EBGeometry::Dcel::defaultStopFunction); + + // Create the linear representation of the conventional BVH SDF above. + std::cout << "Flattening BVH tree\n"; + auto linSDF = bvhSDF->flattenTree(); + + // Compute signed distance for this position and time all SDF representations. + const Vec3 point = Vec3::one(); + + const auto t1 = std::chrono::high_resolution_clock::now(); + const T directDist = directSDF->signedDistance(point); + const auto t2 = std::chrono::high_resolution_clock::now(); + + const auto t3 = std::chrono::high_resolution_clock::now(); + const T bvhDist = bvhSDF->signedDistance(point); + const auto t4 = std::chrono::high_resolution_clock::now(); + + const auto t5 = std::chrono::high_resolution_clock::now(); + const T linDist = linSDF->signedDistance(point); + const auto t6 = std::chrono::high_resolution_clock::now(); + + // Kill all the SDF representations. + directSDF = nullptr; + bvhSDF = nullptr; + linSDF = nullptr; + + // Get the timings. + const std::chrono::duration directTime = t2-t1; + const std::chrono::duration bvhTime = t4-t3; + const std::chrono::duration linTime = t6-t5; + + std::cout << "Distance and time using direct query = " << directDist << ", which took " << directTime.count() << " us\n"; + std::cout << "Distance and time using bvh query = " << bvhDist << ", which took " << bvhTime. count() << " us\n"; + std::cout << "Distance and time using linear bvh query = " << linDist << ", which took " << linTime. count() << " us\n"; return 0; } diff --git a/Examples/Union/README.md b/Examples/Union/README.md index c4778805..9b95f297 100644 --- a/Examples/Union/README.md +++ b/Examples/Union/README.md @@ -6,7 +6,7 @@ Two unions are defined: : * A union that uses BVHs for finding the closest object(s). The scene is defined as an array of analytic spheres bound by axis-aligned bounding boxes. -To shift the computational load to the signed distance function itself, computing the distance function has been given a delay of about 1 millisecond. +To shift the computational load to the signed distance function itself, computing the distance function has been given a delay of about 1 microsecond. Compiling --------- diff --git a/README.md b/README.md index bd10d99f..51120969 100644 --- a/README.md +++ b/README.md @@ -4,22 +4,29 @@ EBGeometry A compact code for computing signed distance functions to watertight and orientable surface grids. Can be used with embedded-boundary (EB) codes like Chombo or AMReX. -The tesselations must consist of planar polygons, but these polygons are not necessarily restricted to triangles. +The tesselations must consist of planar polygons (not necessarily triangles). Internally, the surface mesh is stored in a doubly-connected edge list (DCEL), i.e. a half-edge data structure. -On watertight and orientable grids, the distance to any feature (facet, edge, vertex) is well defined, and can naively be computed by computing the distance to every facet. +On watertight and orientable grids, the distance to any feature (facet, edge, vertex) is well defined, and can naively be computed in various ways: -EBGeometry provides bounding volume hierarchies (BVHs) for accelerating the signed distance computation. -In the DCEL context the BVHs are used for bounding the facets on the surface mesh, but there are no fundamental limitations on which objects that can be bounded. -Thus, multiple objects (e.g., surface grids or analytic functions) can also be bound in the BVHs. -Querying the distance to the mesh through the BVH is much faster than the naive approach. -On average, if the mesh consists of N facets then a BVH has O(log(N)) complexity while a direct search has O(N) complexity. +* Directly, by iterating through all facets. +* With conventional bounding volume hierarchies (BVHs). +* With compact (linearized) BVHs. + +The BVHs in EBGeometry are not limited to facets. +Users can also embed entire objects (e.g., analytic functions) in the BVHs, e.g. the BVH accelerator can be used for a packed sphere geometry. +BVHs can also be nested so that the BVH accelerator is used to embed objects that are themselves described by a BVH. +For example, a scene consisting of many objects described by surface grids can be embedded as a BVH-of-BVH type of scene. + + + +In addition, EBGeometry provides standard operators for signed distance fields like rotations, translations, and scalings. +Multi-object scenes can be constructed with conventional unions, or with BVH-enabled unions (which can be orders of magnitudes faster). Requirements ------------ * A C++ compiler which supports C++14. -* EBGeometry takes a watertight and orientable surface as input (only PLY files currently supported). - Although EBGeometry does it's best at processing grids that contain self-intersections, holes, and hanging vertices the signed distance function is not well-defined. +* Watertight and orientable surfaces (only PLY files currently supported). Basic usage ----------- diff --git a/Source/EBGeometry_BVH.hpp b/Source/EBGeometry_BVH.hpp index 44f2010f..3afc8745 100644 --- a/Source/EBGeometry_BVH.hpp +++ b/Source/EBGeometry_BVH.hpp @@ -34,6 +34,18 @@ namespace BVH { template class NodeT; + /*! + @brief Forward declare linear node class + */ + template + class LinearNodeT; + + /*! + @brief Forward declare linear BVH class. + */ + template + class LinearBVH; + /*! @brief Alias to cut down on typing. */ @@ -67,7 +79,7 @@ namespace BVH { /*! @brief Enum for determining if a BVH node is a leaf or a regular node (only leaf nodes contain data) */ - enum class NodeType { + enum class NodeType : bool { Regular, Leaf, }; @@ -259,6 +271,13 @@ namespace BVH { inline T pruneUnordered2(const Vec3& a_point) const noexcept; + /*! + @brief Flatten everything beneath this node into a depth-first sorted BVH hierarchy. + @details This will compute the flattening of the standard BVH tree and return a pointer to the root node. + */ + inline + std::shared_ptr > flattenTree(); + protected: /*! @@ -411,6 +430,243 @@ namespace BVH { */ inline void pruneUnordered2(T& a_minDist2, std::shared_ptr& a_closest, const Vec3& a_point) const noexcept; + + /*! + @brief Flatten tree method. + @details This function will flatten everything beneath the current node and linearize all the nodes and primitives beneath it to + a_linearNodes and a_sortedPrimitives. This function is called recursively. + @param[inout] a_linearNodes BVH nodes, linearized onto a vector. + @param[inout] a_sortedPrimitives Sorted primitives (in leaf node order). + @param[inout] a_offset Supporting integer for figuring out where in the tree we are. + @note When called from the root node, a_linearNodes and a_sortedPrimitives should be empty and a_offset=0UL. + */ + inline + unsigned long flattenTree(std::vector >& a_linearNodes, + std::vector >& a_sortedPrimitives, + unsigned long& a_offset) const noexcept; + }; + + /*! + @brief Node type for linearized (flattened) BVH. This will be constructed from the other (conventional) BVH type. + + @details T is the precision for Vec3, P is the primitive type you want to enclose, BV is the bounding volume you use for it. + + @note P MUST supply function signedDistance(...) and unsignedDistance2(Vec3). BV must supply a + function getDistance (had this been C++20, we would have use concepts to enforce this). Note that LinearNode is the result + of a flattnened BVH hierarchy where nodes are stored with depth-first ordering for improved cache-location in the downward + traversal. + + @note This class exists so that we can fit the nodes with a smaller memory footprint. The standard BVH node (NodeT) is very useful + when building the tree but less useful when traversing it since it stores references to the primitives in the node itself. It will span + multiple cache lines. This node exists so that we can fit all the BVH info onto fewer cache lines. The number of cache lines will depend + on the tree degree, precision, and bounding volume that is chosen. + + @todo There's a minor optimization that can be made to the memory alignment, which is as follows: For a leaf node we never really need + the m_childOffsets array, and for a regular node we never really need the m_primitivesOffset member. Moreover, m_childOffsets could be + made into a K-1 sized array because we happen to know that the linearized hierarchy will store the first child node immediately after + the regular node. We could shave off 16 bytes of storage, which would mean that a double-precision binary tree only takes up one word + of CPU memory. + */ + template + class LinearNodeT { + public: + + /*! + @brief Alias for cutting down on typing. + */ + using Vec3 = Vec3T; + + /*! + @brief Constructor. + */ + inline + LinearNodeT(); + + /*! + @brief Destructor. + */ + inline + virtual ~LinearNodeT(); + + /*! + @brief Set the bounding volume + @param[in] a_bv Bounding volume for this node. + */ + inline + void setBoundingVolume(const BV& a_boundingVolume) noexcept; + + /*! + @brief Set the offset into the primitives array. + */ + inline + void setPrimitivesOffset(const unsigned long a_primitivesOffset) noexcept; + + /*! + @brief Set number of primitives. + @param[in] a_numPrimitives Number of primitives. + */ + inline + void setNumPrimitives(const int a_numPrimitives) noexcept; + + /*! + @brief Set the child offsets. + @param[in] a_childOffset Offset in node array. + @param[in] a_whichChild Child index in m_childrenOffsets. Must be [0,K-1] + */ + inline + void setChildOffset(const unsigned long a_childOffset, const int a_whichChild) noexcept; + + /*! + @brief Get the node bounding volume. + return m_boundingVolume + */ + inline + const BV& getBoundingVolume() const noexcept; + + /*! + @brief Get the primitives offset + @return Returns m_primitivesOffset + */ + inline + const unsigned long& getPrimitivesOffset() const noexcept; + + /*! + @brief Get the number of primitives. + @return Returns m_numPrimitives + */ + inline + const unsigned long& getNumPrimitives() const noexcept; + + /*! + @brief Get the child offsets + @return Returns m_childOffsets + */ + inline + const std::array& getChildOffsets() const noexcept; + + /*! + @brief Is leaf or not + */ + inline + bool isLeaf() const noexcept; + + /*! + @brief Get the distance from a 3D point to the bounding volume + @param[in] a_point 3D point + @return Returns distance to bounding volume. A zero distance implies that the input point is inside the bounding volume. + */ + inline + T getDistanceToBoundingVolume(const Vec3& a_point) const noexcept; + + /*! + @brief Get the unsigned square from a 3D point to the bounding volume + @param[in] a_point 3D point + @return Returns squared distance to bounding volume. A zero distance implies that the input point is inside the bounding volume. + */ + inline + T getDistanceToBoundingVolume2(const Vec3& a_point) const noexcept; + + /*! + @brief Compute signed distance to primitives. + @note Only call if this is a leaf node. + */ + inline + T getDistanceToPrimitives(const Vec3& a_point, const std::vector >& a_primitives) const noexcept; + + /*! + @brief Pruning algorithm. This is the same algorithm as NodeT::pruneOrdered2, except that the nodes and primitives come in as arguments (and the node + has been collapsed onto one cache line). + */ + inline + void pruneOrdered2(T& a_shortestSquareDistanceSoFar, + unsigned long& a_closestPrimitiveSoFar, + const Vec3& a_point, + const std::vector >& a_linearNodes, + const std::vector >& a_primitives) const noexcept; + + protected: + + /*! + @brief Bounding volume. + */ + BV m_boundingVolume; + + /*! + @brief Offset into primitives array + */ + unsigned long m_primitivesOffset; + + /*! + @brief Number of primitives + */ + int m_numPrimitives; + + /*! + @brief Offset to child nodes. + */ + std::array m_childOffsets; + }; + + /*! + @brief Linear root node for BVH hierarchy + */ + template + class LinearBVH { + public: + + /*! + @brief Cut down on typing + */ + using Vec3 = Vec3T; + + /*! + @brief Alias for cutting down on typing + */ + using LinearNode = LinearNodeT; + + /*! + @brief List of primitives + */ + using PrimitiveList = std::vector >; + + /*! + @brief Disallowed. Use the full constructor please. + */ + LinearBVH() = delete; + + /*! + @brief Full constructor. Associates the nodes and primitives. + @param[in] a_linearNodes Linearized BVH nodes. + @param[in] a_primitives Primitives. + */ + inline + LinearBVH(const std::vector& a_linearNodes, + const PrimitiveList& a_primitives); + + /*! + @brief Destructor. Does nothing + */ + inline + virtual ~LinearBVH(); + + /*! + @brief Function which computes the signed distance + @param[in] a_point 3D point in space + */ + inline + T signedDistance(const Vec3& a_point, const Prune a_pruning = Prune::Ordered2) const noexcept; + + protected: + + /*! + @brief List of linearly stored nodes + */ + std::vector > m_linearNodes; + + /*! + @brief Global list of primitives. Note that this is ALL primitives, sorted so that LinearNodeT can interface into it. + */ + PrimitiveList m_primitives; }; } diff --git a/Source/EBGeometry_BVHImplem.hpp b/Source/EBGeometry_BVHImplem.hpp index 332b88da..5aa8a135 100644 --- a/Source/EBGeometry_BVHImplem.hpp +++ b/Source/EBGeometry_BVHImplem.hpp @@ -489,6 +489,274 @@ namespace BVH { } } } + + template + inline + std::shared_ptr > NodeT::flattenTree() { + + // Create a list of sorted primitives and nodes. + std::vector > sortedPrimitives; + std::vector > linearNodes; + + // Track the offset into the linearized node array. + unsigned long offset = 0; + + // Flatten recursively. + this->flattenTree(linearNodes, sortedPrimitives, offset); + + // Return the root node. + return std::make_shared >(linearNodes, sortedPrimitives); + } + + template + inline + unsigned long NodeT::flattenTree(std::vector >& a_linearNodes, + std::vector >& a_sortedPrimitives, + unsigned long& a_offset) const noexcept { + + // TLDR: This is the main routine for flattening the hierarchy beneath the current node. When this is called we insert + // this node into a_linearNodes and associate the array offsets so that we can find the children in the linearized array. + + // Current node we are dealing with. + const auto curNode = a_offset; + + // Insert a new node corresponding to this node and provide it with the current bounding volume. + a_linearNodes.emplace_back(LinearNodeT()); + a_linearNodes[curNode].setBoundingVolume(m_boundingVolume); + + a_offset++; + + switch(m_nodeType){ + case NodeType::Leaf: + { + // Insert primitives and offsets. + a_linearNodes[curNode].setNumPrimitives (m_primitives. size()); + a_linearNodes[curNode].setPrimitivesOffset(a_sortedPrimitives.size()); + + a_sortedPrimitives.insert(a_sortedPrimitives.end(), m_primitives.begin(), m_primitives.end()); + + break; + } + case NodeType::Regular: + { + a_linearNodes[curNode].setNumPrimitives (0 ); + a_linearNodes[curNode].setPrimitivesOffset(0UL); + + // Go through the children nodes and + for (int k = 0; k < K; k++){ + const int offset = m_children[k]->flattenTree(a_linearNodes, a_sortedPrimitives, a_offset); + + a_linearNodes[curNode].setChildOffset(offset, k); + } + + break; + } + } + + return curNode; + } + + template + inline + LinearNodeT::LinearNodeT() { + + // Initialize everything. + m_boundingVolume = BV(); + m_primitivesOffset = 0UL; + m_numPrimitives = 0; + + for (auto& offset : m_childOffsets){ + offset = 0UL; + } + } + + template + inline + LinearNodeT::~LinearNodeT() { + } + + template + inline + void LinearNodeT::setBoundingVolume(const BV& a_boundingVolume) noexcept { + m_boundingVolume = a_boundingVolume; + } + + template + inline + void LinearNodeT::setPrimitivesOffset(const unsigned long a_primitivesOffset) noexcept { + m_primitivesOffset = a_primitivesOffset; + } + + template + inline + void LinearNodeT::setNumPrimitives(const int a_numPrimitives) noexcept { + m_numPrimitives = a_numPrimitives; + } + + template + inline + void LinearNodeT::setChildOffset(const unsigned long a_childOffset, const int a_whichChild) noexcept { + m_childOffsets[a_whichChild] = a_childOffset; + } + + template + inline + const BV& LinearNodeT::getBoundingVolume() const noexcept { + return m_boundingVolume; + } + + template + inline + const unsigned long& LinearNodeT::getPrimitivesOffset() const noexcept { + return m_primitivesOffset; + } + + template + inline + const unsigned long& LinearNodeT::getNumPrimitives() const noexcept { + return m_numPrimitives; + } + + template + inline + const std::array& LinearNodeT::getChildOffsets() const noexcept { + return m_childOffsets; + } + + template + inline + bool LinearNodeT::isLeaf() const noexcept { + return m_numPrimitives > 0; + } + + template + inline + T LinearNodeT::getDistanceToBoundingVolume(const Vec3& a_point) const noexcept{ + return m_boundingVolume.getDistance(a_point); + } + + template + inline + T LinearNodeT::getDistanceToBoundingVolume2(const Vec3& a_point) const noexcept{ + return m_boundingVolume.getDistance2(a_point); + } + + template + inline + T LinearNodeT::getDistanceToPrimitives(const Vec3T& a_point, const std::vector >& a_primitives) const noexcept { + T minDist = std::numeric_limits::infinity(); + + for (unsigned int i = 0; i < m_numPrimitives; i++){ + const T curDist = a_primitives[m_primitivesOffset + i]->signedDistance(a_point); + + if(std::abs(curDist) < std::abs(minDist)){ + minDist = curDist; + } + } + + return minDist; + } + + template + inline + void LinearNodeT::pruneOrdered2(T& a_shortestSquareDistanceSoFar, + unsigned long& a_closestPrimitiveSoFar, + const Vec3& a_point, + const std::vector >& a_linearNodes, + const std::vector >& a_primitives) const noexcept { + + // If this is a leaf node, see if one of the primitives in this node is closer than the closest distance this far. If it is, + // update the closest primitive. + + if(m_numPrimitives > 0){ + + // See if this node has primitives that are closer than the closest one we've found so far. + for (unsigned int i = 0; i < m_numPrimitives; i++){ + const std::shared_ptr& curPrim = a_primitives[m_primitivesOffset + i]; + + const T curDist2 = curPrim->unsignedDistance2(a_point); + + if(curDist2 < a_shortestSquareDistanceSoFar){ + a_shortestSquareDistanceSoFar = curDist2; + a_closestPrimitiveSoFar = m_primitivesOffset + i; + } + } + } + else{ + // In this case we are at a regular node and we need to determine which branch to descend first. One we've selected a branch we won't visit + // the other branch until we've drilled all the way down into the bottom of the tree, so this choice is critical. We assume that the closest + // bounding volume is more likely to also contain the closest primitive, so we sort the distance to the bounding volumes and take the closest + // bounding volume first. + + // Compute the distance to the child node bounding volumes. + std::array, K > distancesAndNodes; + for (int k = 0; k < K; k++){ + const unsigned long& childOffset = m_childOffsets[k]; + const LinearNodeT& childNode = a_linearNodes[childOffset]; + + distancesAndNodes[k] = std::make_pair(childNode.getDistanceToBoundingVolume2(a_point), childOffset); + } + + // Sort, closest node goes first. + std::sort(distancesAndNodes.begin(), + distancesAndNodes.end(), + [](const std::pair& node1, const std::pair& node2){ + return node1.first < node2.first; + }); + + + // Go through the child nodes and call this function again. + for (int k = 0; k < K; k++){ + const unsigned long& childOffset = distancesAndNodes[k].second; + const LinearNodeT& childNode = a_linearNodes[childOffset]; + + if(distancesAndNodes[k].first < a_shortestSquareDistanceSoFar){ + childNode.pruneOrdered2(a_shortestSquareDistanceSoFar, a_closestPrimitiveSoFar, a_point, a_linearNodes, a_primitives); + } + else{ // Prune the other subtree. + break; + } + } + } + } + + template + inline + LinearBVH::LinearBVH(const std::vector& a_linearNodes, + const PrimitiveList& a_primitives) { + m_linearNodes = a_linearNodes; + m_primitives = a_primitives ; + } + + template + inline + LinearBVH::~LinearBVH() { + + } + + template + T LinearBVH::signedDistance(const Vec3& a_point, const Prune a_pruning) const noexcept { + T minDist = std::numeric_limits::infinity(); + + switch(a_pruning){ + case Prune::Ordered2: + { + unsigned long closestPrimitiveSoFar = 0UL; + + m_linearNodes[0].pruneOrdered2(minDist, closestPrimitiveSoFar, a_point, m_linearNodes, m_primitives); + + minDist = m_primitives[closestPrimitiveSoFar]->signedDistance(a_point); + + break; + } + default: + std::cerr << "In file EBGeometry_BVHImplem.hpp function LinearBVH::signedDistance(Vec3, Prune) -- bad input enum for 'Prune'\n"; + }; + + return minDist; + } + + } #include "EBGeometry_NamespaceFooter.hpp" diff --git a/Source/EBGeometry_DcelMeshImplem.hpp b/Source/EBGeometry_DcelMeshImplem.hpp index ddf95f38..efcc7b48 100644 --- a/Source/EBGeometry_DcelMeshImplem.hpp +++ b/Source/EBGeometry_DcelMeshImplem.hpp @@ -289,7 +289,7 @@ namespace Dcel { template inline T MeshT::signedDistance(const Vec3& a_point, SearchAlgorithm a_algorithm) const noexcept { - T minDist; + T minDist = std::numeric_limits::max(); switch(a_algorithm){ case SearchAlgorithm::Direct: diff --git a/Source/EBGeometry_SignedDistanceBVH.hpp b/Source/EBGeometry_SignedDistanceBVH.hpp deleted file mode 100644 index d4524d59..00000000 --- a/Source/EBGeometry_SignedDistanceBVH.hpp +++ /dev/null @@ -1,98 +0,0 @@ -/* EBGeometry - * Copyright © 2022 Robert Marskar - * Please refer to Copyright.txt and LICENSE in the EBGeometry root directory. - */ - -/*! - @file EBGeometry_SignedDistanceBVH.hpp - @brief Declaration of a signed-distance object that uses bounding volume hiearchices for accelerating signed queries from DCEL tesselations. - @author Robert Marskar -*/ - -#ifndef EBGeometry_SignedDistanceBVH -#define EBGeometry_SignedDistanceBVH - -// Std includes -#include - -// Our includes -#include "EBGeometry_SignedDistanceFunction.hpp" -#include "EBGeometry_DcelMesh.hpp" -#include "EBGeometry_NamespaceHeader.hpp" - -/*! - @brief Signed distance function from a DCEL mesh, using BVHs for accelerating the query. - @details This class stores a reference to the node in the BVH tree. For that reason, the class is templated using the precision (T), the bounding volume type - that encloses the primitives (BV), and the tree degree (K. -*/ -template -class SignedDistanceBVH : public SignedDistanceFunction { -public: - - /*! - @brief Alias for DCEL faces - */ - using Face = EBGeometry::Dcel::FaceT; - - /*! - @brief Alias for BVH node - */ - using Node = EBGeometry::BVH::NodeT; - - /*! - @brief Alias for DCEl mesh with precision T - */ - using Mesh = Dcel::MeshT; - - /*! - @brief Weak constructor. - */ - SignedDistanceBVH() = default; - - /*! - @brief Full constructor - @param[in] a_rootNode Reference to root node in BVH. - @param[in] a_flipSign Hook for turning inside to outside - */ - SignedDistanceBVH(const std::shared_ptr& a_rootNode, const bool a_flipSign); - - /*! - @brief Copy constructor - @param[in] a_object Other distance function - */ - SignedDistanceBVH(const SignedDistanceBVH& a_object); - - /*! - @brief Destructor (does nothing) - */ - virtual ~SignedDistanceBVH(); - - /*! - @brief Return the root node bounding volume - */ - const BV& getBoundingVolume() const noexcept; - - /*! - @brief Value function - @param[in] a_point 3D point. - */ - T signedDistance(const Vec3T& a_point) const noexcept override; - -protected: - - /*! - @brief Reference to root node in DCEL mesh. - */ - std::shared_ptr m_rootNode; - - /*! - @brief Hook for turning inside to outside - */ - bool m_flipSign; -}; - -#include "EBGeometry_NamespaceFooter.hpp" - -#include "EBGeometry_SignedDistanceBVHImplem.hpp" - -#endif diff --git a/Source/EBGeometry_SignedDistanceBVHImplem.hpp b/Source/EBGeometry_SignedDistanceBVHImplem.hpp deleted file mode 100644 index 96da204d..00000000 --- a/Source/EBGeometry_SignedDistanceBVHImplem.hpp +++ /dev/null @@ -1,54 +0,0 @@ -/* EBGeometry - * Copyright © 2022 Robert Marskar - * Please refer to Copyright.txt and LICENSE in the EBGeometry root directory. - */ - -/*! - @file EBGeometry_SignedDistanceBVHImplem.hpp - @brief Implementation of EBGeometry_SignedDistanceBVH.hpp - @author Robert Marskar -*/ - -#ifndef EBGeometry_SignedDistanceBVHImplem -#define EBGeometry_SignedDistanceBVHImplem - -// Our includes -#include "EBGeometry_SignedDistanceBVH.hpp" -#include "EBGeometry_NamespaceHeader.hpp" - -using namespace Dcel; -using namespace BVH; - -template -SignedDistanceBVH::SignedDistanceBVH(const std::shared_ptr& a_rootNode, const bool a_flipSign) { - this->m_rootNode = a_rootNode; - this->m_flipSign = a_flipSign; -} - -template -SignedDistanceBVH::SignedDistanceBVH(const SignedDistanceBVH& a_other) { - this->m_mesh = a_other.m_mesh; - this->m_flipSign = a_other.m_flipSign; - this->m_transformOps = a_other.m_transformOps; -} - -template -SignedDistanceBVH::~SignedDistanceBVH() { - -} - -template -const BV& SignedDistanceBVH::getBoundingVolume() const noexcept { - return m_rootNode->getBoundingVolume(); -} - -template -T SignedDistanceBVH::signedDistance(const Vec3T& a_point) const noexcept { - const T sign = (m_flipSign) ? -1.0 : 1.0; - - return sign * m_rootNode->signedDistance(this->transformPoint(a_point)); -} - -#include "EBGeometry_NamespaceFooter.hpp" - -#endif diff --git a/Source/EBGeometry_SignedDistanceDcel.hpp b/Source/EBGeometry_SignedDistanceDcel.hpp deleted file mode 100644 index dfb5d223..00000000 --- a/Source/EBGeometry_SignedDistanceDcel.hpp +++ /dev/null @@ -1,82 +0,0 @@ -/* EBGeometry - * Copyright © 2022 Robert Marskar - * Please refer to Copyright.txt and LICENSE in the EBGeometry root directory. - */ - -/*! - @file EBGeometry_SignedDistanceDcel.hpp - @brief Declaration of a signed-distance object that gets its value function from a DCEL surface tesselation - @author Robert Marskar -*/ - -#ifndef EBGeometry_SignedDistanceDcel -#define EBGeometry_SignedDistanceDcel - -// Std includes -#include - -// Our includes -#include "EBGeometry_SignedDistanceFunction.hpp" -#include "EBGeometry_DcelMesh.hpp" -#include "EBGeometry_NamespaceHeader.hpp" - -/*! - @brief Signed distance function from a DCEL mesh. - @details This does not use a bounding volume hierarchy for performance. The template parameter is the floating point precision. -*/ -template -class SignedDistanceDcel : public SignedDistanceFunction { -public: - - /*! - @brief Alias for DCEl mesh with precision T - */ - using Mesh = Dcel::MeshT; - - /*! - @brief Weak constructor. - */ - SignedDistanceDcel() = default; - - /*! - @brief Full constructor - @param[in] a_mesh DCEL mesh - @param[in] a_flipSign Hook for turning inside to outside - */ - SignedDistanceDcel(const std::shared_ptr& a_mesh, const bool a_flipSign); - - /*! - @brief Copy constructor - @param[in] a_object Other distance function - */ - SignedDistanceDcel(const SignedDistanceDcel& a_object); - - /*! - @brief Destructor (does nothing) - */ - virtual ~SignedDistanceDcel(); - - /*! - @brief Value function - @param[in] a_point 3D point. - */ - T signedDistance(const Vec3T& a_point) const noexcept override; - -protected: - - /*! - @brief Half-edge mesh structure. - */ - std::shared_ptr m_mesh; - - /*! - @brief Hook for turning inside to outside - */ - bool m_flipSign; -}; - -#include "EBGeometry_NamespaceFooter.hpp" - -#include "EBGeometry_SignedDistanceDcelImplem.hpp" - -#endif diff --git a/Source/EBGeometry_SignedDistanceDcelImplem.hpp b/Source/EBGeometry_SignedDistanceDcelImplem.hpp deleted file mode 100644 index 64abafac..00000000 --- a/Source/EBGeometry_SignedDistanceDcelImplem.hpp +++ /dev/null @@ -1,48 +0,0 @@ -/* EBGeometry - * Copyright © 2022 Robert Marskar - * Please refer to Copyright.txt and LICENSE in the EBGeometry root directory. - */ - -/*! - @file EBGeometry_SignedDistanceDcelImplem.H - @brief Implementation of EBGeometry_SignedDistanceDcel.H - @author Robert Marskar -*/ - -#ifndef EBGeometry_SignedDistanceDcelImplem_H -#define EBGeometry_SignedDistanceDcelImplem_H - -// Our includes -#include "EBGeometry_Vec.hpp" -#include "EBGeometry_SignedDistanceDcel.hpp" -#include "EBGeometry_NamespaceHeader.hpp" - -template -SignedDistanceDcel::SignedDistanceDcel(const std::shared_ptr& a_mesh, const bool a_flipSign) { - this->m_mesh = a_mesh; - this->m_flipSign = a_flipSign; -} - -template -SignedDistanceDcel::SignedDistanceDcel(const SignedDistanceDcel& a_object) { - this->m_mesh = a_object.m_mesh; - this->m_flipSign = a_object.m_flipSign; - this->m_transformOps = a_object.m_transformOps; -} - -template -SignedDistanceDcel::~SignedDistanceDcel(){ - -} - -template -T SignedDistanceDcel::signedDistance(const Vec3T& a_point) const noexcept { - - T sign = (m_flipSign) ? -1.0 : 1.0; - - return sign * m_mesh->signedDistance(this->transformPoint(a_point)); -} - -#include "EBGeometry_NamespaceFooter.hpp" - -#endif diff --git a/Source/EBGeometry_SignedDistanceFunction.hpp b/Source/EBGeometry_SignedDistanceFunction.hpp index 853e9ade..63470d8a 100644 --- a/Source/EBGeometry_SignedDistanceFunction.hpp +++ b/Source/EBGeometry_SignedDistanceFunction.hpp @@ -21,6 +21,9 @@ /*! @brief Abstract representation of a signed distance function. + @details Users can put whatever they like in here, e.g. analytic functions, DCEL meshes, or DCEL meshes stored in full or compact BVH trees. The signedDistance + function must be implemented by the user. When computing it, the user can apply transformation operators (rotations, scaling, translations) by calling transformPoint + on the input coordinate. */ template class SignedDistanceFunction { diff --git a/Source/EBGeometry_UnionBVH.hpp b/Source/EBGeometry_UnionBVH.hpp index c888d30d..ba8cd9a4 100644 --- a/Source/EBGeometry_UnionBVH.hpp +++ b/Source/EBGeometry_UnionBVH.hpp @@ -82,9 +82,14 @@ class UnionBVH : public SignedDistanceFunction { using SDFList = std::vector >; /*! - @brief Node type in BVH tree + @brief Builder node type in BVH tree. Tree is constructed in "full". */ - using Node = EBGeometry::BVH::NodeT; + using BuilderNode = EBGeometry::BVH::NodeT; + + /*! + @brief Node type in BVH tree. We use a flattened tree. + */ + using LinearNode = EBGeometry::BVH::LinearBVH; /*! @brief List of distance functions @@ -94,7 +99,7 @@ class UnionBVH : public SignedDistanceFunction { /*! @brief Root node for BVH tree */ - std::shared_ptr m_rootNode; + std::shared_ptr m_rootNode; /*! @brief Is good or not diff --git a/Source/EBGeometry_UnionBVHImplem.hpp b/Source/EBGeometry_UnionBVHImplem.hpp index 7840ba31..fd4a4396 100644 --- a/Source/EBGeometry_UnionBVHImplem.hpp +++ b/Source/EBGeometry_UnionBVHImplem.hpp @@ -124,18 +124,20 @@ void UnionBVH::buildTree(const BVConstructor& a_bvConstructor) { }; // Stop function. Exists subdivision if there are not enough primitives left to keep subdividing. We set the limit at 10 primitives. - EBGeometry::BVH::StopFunctionT stopFunc = [] (const Node& a_node) -> bool { + EBGeometry::BVH::StopFunctionT stopFunc = [] (const BuilderNode& a_node) -> bool { const int numPrimsInNode = (a_node.getPrimitives()).size(); return numPrimsInNode < K; }; // Init the root node and partition the primitives. - m_rootNode = std::make_shared(m_distanceFunctions); + auto root = std::make_shared(m_distanceFunctions); - m_rootNode->topDownSortAndPartitionPrimitives(a_bvConstructor, - partitioner, - stopFunc); + root->topDownSortAndPartitionPrimitives(a_bvConstructor, + partitioner, + stopFunc); + m_rootNode = root->flattenTree(); + m_isGood = true; } diff --git a/example.png b/example.png new file mode 100644 index 00000000..c3d681b2 Binary files /dev/null and b/example.png differ