Add SFC-based bottom-up construction of BVHs (#67)

.github/workflows/CI.yml

@@ -26,7 +26,7 @@ jobs:
       - name: Run clang-format
         uses: jidicula/[email protected]
         with:
-          clang-format-version: '10'
+          clang-format-version: '18'
           check-path: ${{ matrix.directory }}
 
   Linux-GNU:

Docs/Sphinx/source/Concepts.rst

@@ -50,7 +50,7 @@ An example of an implicit function for the same sphere is
    I_{\textrm{sph}}\left(\mathbf{x}\right) = \left|\mathbf{x} - \mathbf{x}_0\right|^2 - R^2.
 
 An important difference between these is the Eikonal property in :eq:`Eikonal`, ensuring that the signed distance function always returns the exact distance to the object.
-Signed distance functions are usually the more useful object, but many operations (e.g. CSG unions) do not preserve the signed distance property.
+Signed distance functions are more useful objects, but many operations (e.g. CSG unions) do not preserve the signed distance property (but still provide *bounds* for the signed distance).
 
 .. _Chap:DCEL:
 
@@ -69,7 +69,7 @@ The DCEL structures consist of the following objects:
 
 As shown in :numref:`Fig:DCEL`, half-edges circulate the inside of the facet, with pointer access to the next half-edge.
 A half-edge also stores a reference to its starting vertex, as well as a reference to its pair-edge.
-From the DCEL structure we can easily obtain all edges or vertices belonging to a single facet, and also jump to a neighboring facet by fetching the pair edge. 
+From the DCEL structure we can obtain all edges or vertices belonging to a single facet by iterating through the half-edges, and also jump to a neighboring facet by fetching the pair edge. 
 
 .. _Fig:DCEL:
 .. figure:: /_static/DCEL.png
@@ -111,7 +111,7 @@ Three cases can be distinguished:
 
    .. tip::
 
-      ``EBGeometry`` uses the crossing number algorithm by default.
+      ``EBGeometry`` uses the crossing number algorithm by default. 
 
    If the point projects to the inside of the face, the signed distance is just :math:`\mathbf{n}_f\cdot\left(\mathbf{x} - \mathbf{x}_f\right)` where :math:`\mathbf{n}_f` is the face normal and :math:`\mathbf{x}_f` is a point on the face plane (e.g., a vertex).
    If the point projects to *outside* the polygon face, the closest feature is either an edge or a vertex.
@@ -163,7 +163,7 @@ Bounding volume hierarchies
 Bounding volume hierarchies (BVHs) are tree structures where the regular nodes are bounding volumes that enclose all geometric primitives (e.g. polygon faces or implicit functions) further down in the hierarchy.
 This means that every node in a BVH is associated with a *bounding volume*.
 The bounding volume can, in principle, be any type of volume. 
-Moreover, there are two types of nodes in a BVH:
+There are two types of nodes in a BVH:
 
 * **Regular/interior nodes.** These do not contain any of the primitives/objects, but store references to subtrees (aka child nodes).
 * **Leaf nodes.** These lie at the bottom of the BVH tree and each of them contains a subset of the geometric primitives.
@@ -186,7 +186,7 @@ Note that the bounding volume for :math:`P` encloses the bounding volumes of :ma
 There is no fundamental limitation to what type of primitives/objects can be enclosed in BVHs, which makes BVHs useful beyond triangulated data sets.
 For example, analytic signed distance functions can also be embedded in BVHs, provided that we can construct bounding volumes that enclose them.
 
-.. note::
+.. important::
 
    ``EBGeometry`` is not limited to binary trees, but supports :math:`k` -ary trees where each regular node has :math:`k` child nodes. 
 
@@ -226,7 +226,10 @@ Top-down construction can thus be illustrated as a recursive procedure:
 	 if(enoughPrimitives(child)):
 	    child.topDownConstruction(child.objects)
 
-In practice, the above procedure is supplemented by more sophisticated criteria for terminating the recursion, as well as routines for creating the bounding volumes around the newly inserted nodes. 
+In practice, the above procedure is supplemented by more sophisticated criteria for terminating the recursion, as well as routines for creating the bounding volumes around the newly inserted nodes.
+
+Bottom-up construction is also possible, in which case one constructs the leaf nodes first, and then merge the nodes upward until one reaches a root node.
+In ``EBGeometry``, bottom-up construction is done by means of space-filling curves (e.g., Morton codes).
 
 Tree traversal
 --------------
@@ -237,9 +240,9 @@ For the traversal algorithm we consider the following steps:
 
 * When descending from node :math:`P` we determine that we first investigate the left subtree (node :math:`A`) since its bounding volume is closer than the bounding volumes for the other subtree.
   The other subtree will is investigated after we have recursed to the bottom of the :math:`A` subtree. 
-* Since :math:`A` is a leaf node, we find the signed distance from :math:`\mathbf{x}` to the primitives in :math:`A`.
+* Since :math:`A` is a leaf node, we compute the signed distance from :math:`\mathbf{x}` to the primitives in :math:`A`.
   This requires us to iterate over all the triangles in :math:`A`. 
-* When moving back to :math:`P`, we find that the distance to the primitives in :math:`A` is shorter than the distance from :math:`\mathbf{x}` to the bounding volume that encloses nodes :math:`B` and :math:`C`.
+* When investigating the other child node of :math:`P`, we find that the distance to the primitives in :math:`A` is shorter than the distance from :math:`\mathbf{x}` to the bounding volume that encloses nodes :math:`B` and :math:`C`.
   This immediately permits us to prune the entire subtree containing :math:`B` and :math:`C`.
 
 .. _Fig:TreePruning:
@@ -251,17 +254,19 @@ For the traversal algorithm we consider the following steps:
 
 .. warning::
 
-   Note that all BVH traversal algorithms have linear complexity when the primitives are all at approximately the same distance from the query point.
-   For example, it is necessary to traverse almost the entire tree when one tries to compute the signed distance at the origin of a tessellated sphere.
+   BVH traversal has :math:`\log N` complexity on average.
+   However in the worst case the traversal algorithm may have linear complexity if the primitives are all at approximately the same distance from the query point.
+   For example, it is necessary to traverse almost the entire tree when one tries to compute the signed distance at the origin of a tessellated sphere since all triangles and their bounding volumes are approximately at the same distance from the center.
 
-Note that types of tree traversal (that do not compute the signed distance) are also possible, e.g. we may want to compute the union :math:`I\left(\mathbf{x}\right) = \min\left(I_1\left(\mathbf{x}\right), I_2\left(\mathbf{x}\right), .\ldots\right)`.
-``EBGeometry`` supports a fairly flexible approach to the tree traversal and update algorithms.
+Other types of tree traversal (that do not compute the signed distance) are also possible.
+``EBGeometry`` supports a fairly flexible approach to the tree traversal and update algorithms such that the user is permitted to use the hierarhical traversal algorithm also for other types of operations (e.g., for finding all facets within a specified distance from a point).
 
 Octree
 ======
 
 Octrees are tree-structures where each interior node has exactly eight children.
-Such trees are usually used for spatial partitioning (and in this case the eight children have no spatial overlap), and the leaf nodes may also contain actual data. 
+Such trees are usually used for spatial partitioning.
+Unlike BVH trees, the eight child nodes have no spatial overlap.
 
 Octree construction can be done in (at least) two ways:
 
@@ -301,4 +306,6 @@ Combining objects
 * Difference.
 
 Some of these CSG operations also have smooth equivalents, i.e. for smoothing the transition between combined objects.
-Fast CSG operations are also supported by ``EBGeometry``, e.g. the BVH-accelerated CSG union where one uses the BVH when searching for the relevant geometric primitive(s). 
+Fast CSG operations are also supported by ``EBGeometry``, e.g. the BVH-accelerated CSG union where one uses the BVH when searching for the relevant geometric primitive(s).
+This functionality is motivated by the fact that a CSG union is normally implemented as :math:`\min\left(I_1, I_2, I_3, \ldots,I_N\right)`, which has :math:`\mathcal{O}\left(N\right)` complexity when there are :math:`N` objects.
+BVH trees can reduce this to :math:`\mathcal{O}\left(\log N\right)` complexity.

Docs/Sphinx/source/ImplemBVH.rst

@@ -14,7 +14,7 @@ The full BVH is encapsulated by a class
 
 .. code-block:: c++
 
-   template <class T, class P, class BV, int K>
+   template <class T, class P, class BV, size_t K>
    class NodeT;
 
 The above template parameters are:
@@ -40,24 +40,31 @@ Bounding volumes
 *  **Axis-aligned bounding box**, which is templated as ``EBGeometry::BoundingVolumes::AABBT<T>``.
 
 For full API details, see `the doxygen API <doxygen/html/namespaceBoundingVolumes.html>`_.
+Other types of bounding volumes can in principle be added, with the only requirement being that they conform to the same interface as the ``AABB`` and ``BoundingSphere`` volumes.
 
 .. _Chap:BVHConstruction:
 
 Construction
 ------------
 
-Constructing a BVH is done by
+Constructing a BVH is done by:
 
-*  Creating a root node and providing it with the geometric primitives and their bounding volumes.
-*  Partitioning the BVH by providing a partitioning function. 
+#.  Creating a root node and providing it with the geometric primitives and their bounding volumes.
+#.  Partitioning the BVH by providing a partitioning function. 
 
 The first step is usually a matter of simply constructing the root node using the full constructor, which takes a list of primitives and their associated bounding volumes. 
-The second step is to recursively build the BVH, which is done through the function ``topDownSortAndPartition()``, see the code below:
+The second step is to recursively build the BVH.
+We currently support top-down and bottom-up construction (using space-filling curves).
 
-.. literalinclude:: ../../../Source/EBGeometry_BVH.hpp
-   :language: c++
-   :lines: 29, 62-94, 217-227, 248-257, 263-268, 274-285, 404, 643,644
-   :caption: Header section of the BVH implementation.
+.. tip::
+
+   The default construction methods performs the hierarchical subdivision by only considering the *bounding volumes*.
+   Consequently, the build process is identical regardless of what type of primitives (e.g., triangles or analytic spheres) are contained in the BVH.
+
+Top-down construction
+_____________________
+
+Top-down construction is done through the function ``topDownSortAndPartition()``, `see the doxygen API for the BVH implementation <doxygen/html/doxygen/html/classBVH_1_1NodeT.html>`_.
 
 The optional input arguments to ``topDownSortAndPartition`` are polymorphic functions of type indicated above, and have the following responsibilities:
 
@@ -68,6 +75,23 @@ The optional input arguments to ``topDownSortAndPartition`` are polymorphic func
 
 Default arguments for these are provided, bubt users are free to partition their BVHs in their own way should they choose.
 
+Bottom-up construction
+______________________
+
+The bottom-up construction uses a space-filling curve (e.g., a Morton curve) for first building the leaf nodes.
+This construction is done such that each leaf node contains approximately the number of primitives, and all leaf nodes exist on the same level.
+To use bottom-up construction, one may use the member function
+
+.. literalinclude:: ../../../Source/EBGeometry_BVH.hpp
+   :language: c++
+   :lines: 298-309
+
+The template argument is the space-filling curve that the user wants to apply.
+Currently, we support Morton codes and nested indices.
+For Morton curves, one would e.g. call ``bottomUpSortAndPartition<SFC::Morton>`` while for nested indices (which are not recommended) the signature is likewise ``bottomUpSortAndPartition<SFC::Nested``.
+
+Build times for SFC-based bottom-up construction are generally speaking faster than top-down construction, but tends to produce worse trees such that traversal becomes slower. 
+
 .. _Chap:LinearBVH:
 
 Compact form
@@ -190,7 +214,7 @@ This function has a signature
   using Visiter = std::function<bool(const NodeType& a_node, const Meta a_meta)>;
 
 where ``NodeType`` is the type of node (which is different for full/flat BVHs), and the ``Meta`` template parameter is discussed below.
-If this function returns true, the node will be visisted and if the function returns false then the node will be pruned from the tree traversal.
+If this function returns true, the node will be visisted and if the function returns false then the node will be pruned from the tree traversal. Typically, the ``Meta`` parameter will contain the necessary information that determines whether or not to visit the subtree.
 
 Traversal pattern
 _________________
@@ -206,6 +230,9 @@ This function has the signature:
   template <class NodeType, class Meta, size_t K>
   using Sorter = std::function<void(std::array<std::pair<std::shared_ptr<const NodeType>, Meta>, K>& a_children)>;
 
+Sorting the child nodes is completely optional.
+The user can leave this function empty if it does not matter which subtrees are visited first. 
+
 Update rule
 ___________
 
@@ -217,6 +244,8 @@ These are done by a user-supplied update-rule:
   template <class P>
   using Updater = std::function<void(const PrimitiveListT<P>& a_primitives)>;
 
+Typically, the ``Updater`` will modify parameters that appear in a local scope outside of the tree traversal (e.g. updating the minimum distance to a DCEL mesh).
+
 Meta-data
 _________
 
@@ -229,8 +258,25 @@ The signature for meta-data construction is
   template <class NodeType, class Meta>
   using MetaUpdater = std::function<Meta(const NodeType& a_node)>;
 
-Traversal example
-_________________
+The biggest difference between ``Updater`` and ``MetaUpdater`` is that ``Updater`` is *only* called on leaf nodes whereas ``MetaUpdater`` is also called for internal nodes.
+One typical example for DCEL meshes is that ``Updater`` computes the distance from an input point to the triangles in a leaf node, whereas ``MetaUpdater`` computes the distance from the input point to the bounding volumes of a child nodes.
+This information is then used in ``Sorter`` in order to determine a preferred child visit pattern when descending along subtrees. 
+
+Traversal algorithm
+___________________
+
+The code-block below shows the implementation of the BVH traversal.
+The implementation uses a non-recursive queue-based formulation when descending along subtrees.
+Observe that each entry in the stack contains both the node itself *and* any meta-data we want to attach to the node.
+If the traversal decides to visit a node, it immediately computes the specified meta-data of the node, and the user can then sort the children based on that data.
+
+.. literalinclude:: ../../../Source/EBGeometry_BVHImplem.hpp
+   :language: c++
+   :lines: 284-322
+   :caption: Tree traversal algorithm for the BVH tree.
+
+Traversal examples
+__________________
 
 The DCEL mesh distance fields use a traversal pattern based on
 
@@ -242,6 +288,6 @@ These rules are given below.
 
 .. literalinclude:: ../../../Source/EBGeometry_MeshDistanceFunctionsImplem.hpp
    :language: c++
-   :lines: 127-161
+   :lines: 97-132
    :caption: Tree traversal criterion for computing the signed distance to a DCEL mesh using the BVH accelerator.
 	     See :file:`Source/EBGeometry_MeshDistanceFunctionsImplem.hpp` for details.

Docs/Sphinx/source/ImplemCSG.rst

@@ -37,7 +37,7 @@ These are also available through functions that automatically cast the resulting
 
 .. literalinclude:: ../../../Source/EBGeometry_Transform.hpp   
    :language: c++
-   :lines: 20-90
+   :lines: 20-111
 
 CSG operations
 --------------

Docs/Sphinx/source/ImplemDCEL.rst

@@ -10,7 +10,7 @@ The DCEL functionality exists under the namespace ``EBGeometry::DCEL`` and conta
 
 .. important::
 
-   The DCEL functionality is *not* restricted to triangles, but supports N-sided polygons, including *meta-data* attached to the vertices, edges, and facets. 
+   The DCEL functionality is *not* restricted to triangles, but supports N-sided polygons, including *meta-data* attached to the vertices, edges, and facets. The latter is particularly useful in case on wants to associate e.g. boundary conditions to specific triangles. 
 
 Main types
 ----------
@@ -77,7 +77,7 @@ The above DCEL classes have member functions of the type:
 
 which can be used to compute the distance to the various features on the mesh.
 
-Meta-data can be attached to the DCEL primitives by selecting an appropriate type for ``Meta`` above (which defaults to ``short``). 
+Meta-data can be attached to the DCEL primitives by selecting an appropriate type for ``Meta`` above.
 
 
 .. _Chap:BVHIntegration:

Docs/Sphinx/source/ImplemOctree.rst

@@ -5,7 +5,7 @@ Octree
 
 The octree functionality is encapsulated in the namespace ``EBGeometry::Octree``.
 For the full API, see `the doxygen API <doxygen/html/namespaceOctree.html>`__.
-Currently, only full octrees are supported (i.e. pointer-based representation).
+Currently, only full octrees are supported (i.e., using a pointer-based representation).
 
 Octrees are encapsulated by a class
 

Docs/Sphinx/source/Introduction.rst

@@ -20,13 +20,13 @@ To obtain ``EBGeometry``, clone the code from `github <https://github.com/rmrsk/
 To compile the ``EBGeometry`` example codes, navigate to the :file:`EBGeometry/Examples` folder.
 Folders that are named :file:`EBGeometry_<something>` are pure ``EBGeometry`` examples and can be compiled without any third-party dependencies.
 
-To run the ``EBGeometry`` examples, navigate to one of the folders and execute
+To run the ``EBGeometry`` examples, navigate to one of the folders in :file:`EBGeometry/Examples/EBGeometry_*` and execute
 
 .. code-block:: bash
 
-   g++ -O3 -std=c++17 main.cpp && ./a.out
+   g++ -O3 main.cpp && ./a.out
 
-All ``EBGeometry`` examples should run using this command.
+All ``EBGeometry`` examples can be run using this command.
 README files present in each folder provide more information regarding the functionality and usage of each example code.
 
 Third-party examples

Docs/Sphinx/source/Parsers.rst

@@ -43,10 +43,12 @@ To read one or multiple STL files and turn it into DCEL meshes, use
    :language: c++
    :lines: 53-67
 
+Note that this will only expose the DCEL mesh, but not include any signed distance functionality.
+
 DCEL mesh SDF
 _____________
 
-To read one or multiple STL files and turn it into signed distance representations, use
+To read one or multiple STL files and also turn it into signed distance representations, use
 
 .. literalinclude:: ../../../Source/EBGeometry_Parser.hpp
    :language: c++
@@ -59,7 +61,7 @@ To read one or multiple STL files and turn it into signed distance representatio
 
 .. literalinclude:: ../../../Source/EBGeometry_Parser.hpp
    :language: c++
-   :lines: 85-99
+   :lines: 85-105
 
 .. _Chap:LinearSTL:
 
@@ -70,7 +72,7 @@ To read one or multiple STL files and turn it into signed distance representatio
 
 .. literalinclude:: ../../../Source/EBGeometry_Parser.hpp
    :language: c++
-   :lines: 101-115
+   :lines: 107-127
 
 .. _Chap:PolySoups:
 
@@ -91,11 +93,12 @@ To turn this into a DCEL mesh, one should compress the triangle soup (get rid of
 
 .. literalinclude:: ../../../Source/EBGeometry_Parser.hpp
    :language: c++
-   :lines: 136-165
+   :lines: 146-165
 
 The ``compress`` function will discard duplicate vertices from the soup, while the ``soupToDCEL`` will simply turn the remaining polygon soup into a DCEL mesh.
+This function will also compute the vertex and edge normal vectors.
 
-.. tip::
+.. warning::
 
    ``soupToDCEL`` will issue plenty of warnings if the polygon soup is not watertight and orientable.