BVH

The BVH class implements a bounding volume hierarchy. This data structure recursively subdivides a rectilinear region of interest into a “tree” of subregions.

BVH is well-suited for particle-mesh or ray-mesh intersection tests. It is also well-suited to data sets where the contents are unevenly distributed, since the bins are subdivided based on their contents.

The following code example shows how a BVH can be used to accelerate a point-mesh intersection algorithm. First, we generate bounding boxes for all triangles in the mesh, and call BVH::initialize() with the bounding boxes.

#include "axom/spin/BVH.hpp"
// the BVH is in 2D, storing an index to 2D triangles
using BVH2DType = axom::spin::BVH<in2D>;
// supporting classes
using BoundingBox2DType = axom::primal::BoundingBox<double, in2D>;
using Point2DType = axom::primal::Point<double, in2D>;
using Triangle2DType = axom::primal::Triangle<double, in2D>;
BoundingBox2DType findBbox(Triangle2DType& tri);

BVH2DType* buildBVHTree(std::vector<Triangle2DType>& tris)
{
  // Create a BVH
  BVH2DType* tree = new BVH2DType;

  std::vector<BoundingBox2DType> bboxes(tris.size());

  // Get bounding boxes of each element
  for(size_t i = 0; i < tris.size(); ++i)
  {
    bboxes[i] = findBbox(tris[i]);
  }

  // Build bounding volume hierarchy from bounding boxes
  tree->initialize(bboxes.data(), bboxes.size());

  return tree;
}

After the structure is built, we can use the BVH to generate a list of element IDs that are candidate neighbors to the query points. Call BVH::findPoints() to get the list of element IDs that the query point intersects. The key idea of findPoints() is that testing for probe intersection with a child node (bounding box) is cheap. If a node intersection test fails (misses), the child node, and all of its corresponding elements, can be skipped during the BVH traversal. If the probe does intersect a child node, the node’s children are also tested for probe intersection. Without the acceleration data structure, each probe point must be tested against each triangle.

void findCandidateBVHTreeBins(BVH2DType* tree,
                              Point2DType ppoint,
                              std::vector<int>& candidates)
{
  axom::IndexType offsets;
  axom::IndexType counts;
  // Get the candidates for a given probe point:
  // BVH::findPoints takes an array of points, and allocates and fills an array
  // for all the candidate intersections with the points in a packed manner.
  axom::Array<axom::IndexType> candidatesArray;
  tree->findPoints(axom::ArrayView<axom::IndexType>(&offsets, 1),
                   axom::ArrayView<axom::IndexType>(&counts, 1),
                   candidatesArray,
                   1,
                   &ppoint);

  // Since we are only querying one point, offsets == 0 and
  // len(candidatesPtr) == counts
  candidates =
    std::vector<int>(candidatesArray.data(), candidatesArray.data() + counts);
}

Note that the returned packed candidate intersection array (candidatesPtr above) needs to be deallocated by the caller.

Finally, test the point against all candidate neighbor triangles.

void findIntersectionsWithCandidates(std::vector<Triangle2DType>& tris,
                                     std::vector<int>& candidates,
                                     Point2DType ppoint,
                                     std::vector<int>& intersections)
{
  // Test if ppoint lands in any of its neighbor triangles.
  int csize = static_cast<int>(candidates.size());
  for(int i = 0; i < csize; ++i)
  {
    Triangle2DType& t = tris[candidates[i]];
    if(t.checkInTriangle(ppoint))
    {
      intersections.push_back(candidates[i]);
    }
  }
}

Device Traversal API

The BVH class contains a getTraverser() method, which returns an object that can be used to traverse a BVH with user-defined actions.

The returned traverser type has one function, traverse_tree(), which takes the following arguments:

  • const QueryObject& p: the object to traverse the BVH with. This is passed into each invocation of the traversal predicate.

  • LeafAction&& lf: a function or lambda which is executed on each leaf node of the BVH that is reached during traversal. It should take in two arguments, the index of the leaf node in the BVH, and a pointer to an array mapping leaf node indices to the original index of elements.

  • Predicate&& predicate: a function which determines whether to traverse down to a given internal node. It should take in two arguments: the query object, and the tentative node’s bounding box.

This object may be used within a CUDA kernel, so long as the execution space parameter of BVH is set correctly.

This method can be used to avoid the extra memory allocation needed for holding an array of candidate intersections. For example, if we only wish to count the number of intersections for a given query point, our leaf action could get the underlying mesh element based on the element index, check whether the query point intersects it and then increment a per-point counter.

This method also allows uses of the BVH beyond intersection testing. quest::SignedDistance uses the BVH traversal object to search for the closest surface elements to a query point. The leaf action that is used checks each candidate leaf against a current-minimum candidate; if closer, the current-minimum candidate is set to the new surface element. The predicate used for traversal also utilizes the current-minimum candidate data to avoid traversing internal nodes that are farther than the current minimum squared distance.

Example: Broad-phase collision detection

The following example tests each element of a surface mesh for intersection with each other element. It provides an example of how the device traversal object might be used in a broad-phase collision detection problem. First, we initialize the BVH with the bounding boxes of all the query objects (mesh elements), and create a traverser object:

  // Initialize BVH
  spin::BVH<3, ExecSpace, double> bvh;
  bvh.setAllocatorID(allocatorId);
  bvh.initialize(v_aabbs, v_aabbs.size());

  // Create a traverser object from the BVH
  const auto bvh_device = bvh.getTraverser();

Next, we define a traversal predicate. In this case, since we are testing objects in the mesh against each other, our query objects are bounding boxes. Thus, the traversal predicate tests if the query bounding box intersects with the bounding boxes of nodes in the BVH:

  auto bbIsect = [] AXOM_HOST_DEVICE(const BoxType& queryBbox,
                                     const BoxType& bvhBbox) -> bool {
    return queryBbox.intersectsWith(bvhBbox);
  };

Since we do not know the total number of candidate intersections yet, we must traverse the BVH twice; for the first traversal we count the number of candidate intersections for each query object, allowing us to compute offset indices and total storage requirements for the collision pairs:

  // Allocate arrays for offsets and counts
  axom::Array<IndexType> offsets(ncells, ncells, allocatorId);
  axom::Array<IndexType> counts(ncells, ncells, allocatorId);
  const auto v_counts = counts.view();
  const auto v_offsets = offsets.view();

  // First pass: get number of bounding box collisions for each surface element
  axom::for_all<ExecSpace>(
    ncells,
    AXOM_LAMBDA(IndexType icell) {
      IndexType count = 0;

      // Define a function that is called on every leaf node reached during
      // traversal. The function below simply counts the number of candidate
      // collisions with the given query object.
      auto countCollisions = [&](std::int32_t currentNode,
                                 const std::int32_t* leafNodes) {
        AXOM_UNUSED_VAR(leafNodes);
        if(currentNode > icell)
        {
          count++;
        }
      };

      // Call traverse_tree() to run the counting query.
      bvh_device.traverse_tree(v_aabbs[icell], countCollisions, bbIsect);

      // Afterwards, we can store the number of collisions for each surface
      // element, as well as an overall count of intersections.
      v_counts[icell] = count;
      total_count_reduce += count;
    });

  // Generate offsets
  RAJA::exclusive_scan<exec_pol>(RAJA::make_span(counts.data(), ncells),
                                 RAJA::make_span(offsets.data(), ncells),
                                 RAJA::operators::plus<IndexType> {});

After computing offset indices and allocating output arrays, we can then perform a second traversal through the BVH. This time, we will store candidate collision pairs when we reach a leaf node:

  // Allocate arrays for intersection pairs
  firstPair = axom::Array<IndexType>(ncollisions, ncollisions, allocatorId);
  secondPair = axom::Array<IndexType>(ncollisions, ncollisions, allocatorId);
  const auto v_first_pair = firstPair.view();
  const auto v_second_pair = secondPair.view();

  // Second pass: fill broad-phase collisions array
  axom::for_all<ExecSpace>(
    ncells,
    AXOM_LAMBDA(IndexType icell) {
      IndexType offset = v_offsets[icell];

      // Define a leaf node function that stores the intersection candidate.
      auto fillCollisions = [&](std::int32_t currentNode,
                                const std::int32_t* leafs) {
        if(currentNode > icell)
        {
          v_first_pair[offset] = icell;
          v_second_pair[offset] = leafs[currentNode];
          offset++;
        }
      };

      // Call traverse_tree() a second time to run the counting query.
      bvh_device.traverse_tree(v_aabbs[icell], fillCollisions, bbIsect);
    });

The result of the two-pass query is a list of candidate collision pairs. A code could then do further operations on these pairs of elements, such as test them for actual intersection.