EBGeometry_BVH.hpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2022 Robert Marskar <robert.marskar@sintef.no>
//
// SPDX-License-Identifier: GPL-3.0-or-later
 
#ifndef EBGEOMETRY_BVH_HPP
#define EBGEOMETRY_BVH_HPP
 
// Std includes
#include <algorithm>
#include <array>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <iterator>
#include <limits>
#include <memory>
#include <type_traits>
#include <utility>
#include <vector>
 
#if defined(__SSE4_1__) || defined(__AVX__)
#include <immintrin.h>
#endif
 
// Our includes
#include "EBGeometry_BoundingVolumes.hpp"
#include "EBGeometry_Macros.hpp"
#include "EBGeometry_SFC.hpp"
#include "EBGeometry_Vec.hpp"
 
namespace EBGeometry {
 
namespace BVH {
 
enum class Build
{
  TopDown, 
  Morton,  
  Nested,  
  SAH      
};
 
struct ClusterSpec
{
  size_t maxClusterSize = 8; 
};
 
template <typename T>
[[nodiscard]] constexpr size_t
DefaultBranchingRatio() noexcept
{
  static_assert(std::is_floating_point_v<T>, "BVH::DefaultBranchingRatio requires a floating-point T");
#if defined(__AVX512F__)
  if constexpr (std::is_same_v<T, double>) {
    return 8;
  }
  else {
    return 16;
  }
#elif defined(__AVX__)
  if constexpr (std::is_same_v<T, double>) {
    return 4;
  }
  else {
    return 8;
  }
#else
  return 4;
#endif
}
 
template <class T, class P, class BV, size_t K>
class TreeBVH;
 
template <class P>
using PrimitiveList = std::vector<std::shared_ptr<const P>>;
 
template <class P>
struct SharedPtrStorage
{
  using StorageType = std::shared_ptr<const P>;
 
  [[nodiscard]] static const P&
  get(const StorageType& a_stored) noexcept
  {
    return *a_stored;
  }
 
  static void
  appendTreeLeaf(std::vector<StorageType>& a_dst, const PrimitiveList<P>& a_leafPrims)
  {
    a_dst.insert(a_dst.end(), a_leafPrims.begin(), a_leafPrims.end());
  }
 
  static void
  appendAliased(std::vector<StorageType>& a_dst, const std::shared_ptr<std::vector<P>>& a_block)
  {
    a_dst.reserve(a_dst.size() + a_block->size());
 
    for (size_t i = 0; i < a_block->size(); i++) {
      a_dst.emplace_back(a_block, &(*a_block)[i]);
    }
  }
};
 
template <class P>
struct ValueStorage
{
  using StorageType = P;
 
  [[nodiscard]] static const P&
  get(const StorageType& a_stored) noexcept
  {
    return a_stored;
  }
 
  static void
  appendTreeLeaf(std::vector<StorageType>& a_dst, const PrimitiveList<P>& a_leafPrims)
  {
    // NB: do NOT reserve(a_dst.size() + a_leafPrims.size()) here. appendTreeLeaf is called once per
    // leaf as the whole tree is packed, and reserving to an each-time-slightly-larger *exact* size
    // defeats std::vector's geometric growth -- it forces a fresh reallocation (copying every
    // element already appended) on every leaf, making a whole build O(N^2) in the primitive count.
    // Both the identity pack() constructor and the direct top-down PackedBVH constructor append
    // leaves this way, so the quadratic cost would hit every ValueStorage build over many leaves.
    // Plain push_back grows the buffer geometrically and keeps construction linear.
    for (const auto& p : a_leafPrims) {
      a_dst.push_back(*p);
    }
  }
 
  static void
  appendAliased(std::vector<StorageType>& a_dst, const std::shared_ptr<std::vector<P>>& a_block)
  {
    if (a_dst.empty()) {
      a_dst = std::move(*a_block);
    }
    else {
      a_dst.insert(a_dst.end(), std::make_move_iterator(a_block->begin()), std::make_move_iterator(a_block->end()));
    }
  }
};
 
template <class T, class P, size_t K, class StoragePolicy = SharedPtrStorage<P>>
class PackedBVH;
 
template <class P, class BV>
using PrimAndBV = std::pair<std::shared_ptr<const P>, BV>;
 
template <class P, class BV>
using PrimAndBVList = std::vector<PrimAndBV<P, BV>>;
 
template <class P, class BV, size_t K>
using Partitioner = std::function<std::array<PrimAndBVList<P, BV>, K>(PrimAndBVList<P, BV> a_primsAndBVs)>;
 
template <class T, class P, class BV, size_t K>
using LeafPredicate = std::function<bool(const TreeBVH<T, P, BV, K>& a_node)>;
 
template <class P>
using LeafEvaluator = std::function<void(const PrimitiveList<P>& a_primitives)>;
 
template <class P, class StoragePolicy = SharedPtrStorage<P>>
using PackedLeafEvaluator = std::function<void(
  const std::vector<typename StoragePolicy::StorageType>& a_primitives, size_t a_offset, size_t a_count)>;
 
template <class NodeType, class NodeKey>
using PrunePredicate = std::function<bool(const NodeType& a_node, const NodeKey& a_nodeKey)>;
 
template <class NodeType, class NodeKey, size_t K>
using ChildOrderer =
  std::function<void(std::array<std::pair<std::shared_ptr<const NodeType>, NodeKey>, K>& a_children)>;
 
template <class NodeKey, size_t K>
using PackedChildOrderer = std::function<void(std::array<std::pair<uint32_t, NodeKey>, K>& a_children)>;
 
template <class NodeType, class NodeKey>
using NodeKeyFactory = std::function<NodeKey(const NodeType& a_node)>;
 
template <class X, size_t K>
auto EqualCounts = [](std::vector<X> a_primitives) noexcept -> std::array<std::vector<X>, K> {
  static_assert(K >= 2, "EqualCounts<X, K>: branching factor K must be at least 2");
 
  EBGEOMETRY_EXPECT(!a_primitives.empty());
 
  const int length = a_primitives.size() / K;
  int       remain = a_primitives.size() % K;
 
  int begin = 0;
  int end   = 0;
 
  std::array<std::vector<X>, K> chunks;
 
  for (size_t k = 0; k < K; k++) {
    end += (remain > 0) ? length + 1 : length;
    remain--;
 
    // Move the [begin, end) slice out -- the input is taken by value and not reused, so the elements
    // (e.g. shared_ptr primitive handles) transfer without copying or refcount churn.
    chunks[k] = std::vector<X>(std::make_move_iterator(a_primitives.begin() + begin),
                               std::make_move_iterator(a_primitives.begin() + end));
 
    begin = end;
  }
 
  return chunks;
};
 
template <class T, class P, class BV, size_t K>
auto PrimitiveCentroidPartitioner =
  [](PrimAndBVList<P, BV> a_primsAndBVs) noexcept -> std::array<PrimAndBVList<P, BV>, K> {
  EBGEOMETRY_EXPECT(!a_primsAndBVs.empty());
 
  Vec3T<T> lo = +Vec3T<T>::max();
  Vec3T<T> hi = -Vec3T<T>::max();
 
  for (const auto& pbv : a_primsAndBVs) {
    lo = min(lo, pbv.first->getCentroid());
    hi = max(hi, pbv.first->getCentroid());
  }
 
  const size_t splitDir = (hi - lo).maxDir(true);
 
  // The input is taken by value; sort it in place (no working copy) and move it into EqualCounts.
  std::sort(a_primsAndBVs.begin(),
            a_primsAndBVs.end(),
            [splitDir](const PrimAndBV<P, BV>& pbv1, const PrimAndBV<P, BV>& pbv2) -> bool {
              return pbv1.first->getCentroid(splitDir) < pbv2.first->getCentroid(splitDir);
            });
 
  return BVH::EqualCounts<PrimAndBV<P, BV>, K>(std::move(a_primsAndBVs));
};
 
template <class T, class P, class BV, size_t K>
auto BVCentroidPartitioner = [](PrimAndBVList<P, BV> a_primsAndBVs) -> std::array<PrimAndBVList<P, BV>, K> {
  EBGEOMETRY_EXPECT(!a_primsAndBVs.empty());
 
  Vec3T<T> lo = +Vec3T<T>::max();
  Vec3T<T> hi = -Vec3T<T>::max();
 
  for (const auto& pbv : a_primsAndBVs) {
    lo = min(lo, pbv.second.getCentroid());
    hi = max(hi, pbv.second.getCentroid());
  }
 
  const size_t splitDir = (hi - lo).maxDir(true);
 
  // The input is taken by value; sort it in place (no working copy) and move it into EqualCounts.
  std::sort(a_primsAndBVs.begin(),
            a_primsAndBVs.end(),
            [splitDir](const PrimAndBV<P, BV>& pbv1, const PrimAndBV<P, BV>& pbv2) -> bool {
              return pbv1.second.getCentroid()[splitDir] < pbv2.second.getCentroid()[splitDir];
            });
 
  return BVH::EqualCounts<PrimAndBV<P, BV>, K>(std::move(a_primsAndBVs));
};
 
template <class T, class P, class BV>
inline size_t
SAH2WaySplit(PrimAndBVList<P, BV>& a_list,
             const size_t          a_begin,
             const size_t          a_end,
             const bool            a_longestAxisOnly = false) noexcept
{
  static_assert(std::is_same_v<BV, EBGeometry::BoundingVolumes::AABBT<T>>, "SAH2WaySplit requires BV == AABBT<T>");
 
  EBGEOMETRY_EXPECT(a_begin < a_end);
 
  constexpr int BINS = 32;
 
  const size_t N = a_end - a_begin;
 
  // Centroid bounding box
  Vec3T<T> clo = +Vec3T<T>::max();
  Vec3T<T> chi = -Vec3T<T>::max();
 
  for (size_t i = a_begin; i < a_end; i++) {
    const auto& c = a_list[i].second.getCentroid();
 
    clo = min(clo, c);
    chi = max(chi, c);
  }
 
  T   bestCost  = std::numeric_limits<T>::max();
  T   bestPlane = T(0);
  int bestAxis  = -1;
 
  Vec3T<T> binLo[BINS];
  Vec3T<T> binHi[BINS];
 
  int binCnt[BINS];
 
  T   leftArea[BINS - 1];
  int leftCnt[BINS - 1];
 
  const int firstAxis = a_longestAxisOnly ? (chi - clo).maxDir(true) : 0;
  const int lastAxis  = a_longestAxisOnly ? firstAxis : 2;
 
  for (int axis = firstAxis; axis <= lastAxis; axis++) {
    const T lo  = clo[axis];
    const T hi  = chi[axis];
    const T ext = hi - lo;
 
    if (ext <= T(0)) {
      continue;
    }
 
    const T scale = T(BINS) / ext;
 
    for (int b = 0; b < BINS; b++) {
      binLo[b]  = Vec3T<T>::max();
      binHi[b]  = -Vec3T<T>::max();
      binCnt[b] = 0;
    }
 
    for (size_t i = a_begin; i < a_end; i++) {
      const int b = std::min(BINS - 1, (int)((a_list[i].second.getCentroid()[axis] - lo) * scale));
      binLo[b]    = min(binLo[b], a_list[i].second.getLowCorner());
      binHi[b]    = max(binHi[b], a_list[i].second.getHighCorner());
      binCnt[b]   = binCnt[b] + 1;
    }
 
    // Left prefix: accumulated AABB and count for bins [0..b]
    Vec3T<T> rlo  = Vec3T<T>::max();
    Vec3T<T> rhi  = -Vec3T<T>::max();
    int      rcnt = 0;
 
    for (int b = 0; b < BINS - 1; b++) {
      if (binCnt[b] > 0) {
        rlo = min(rlo, binLo[b]);
        rhi = max(rhi, binHi[b]);
      }
 
      rcnt        = rcnt + binCnt[b];
      leftArea[b] = (rcnt > 0) ? BV(rlo, rhi).getArea() : T(0);
      leftCnt[b]  = rcnt;
    }
 
    // Right suffix sweep; evaluate SAH cost at each candidate split boundary
    Vec3T<T> rrlo  = Vec3T<T>::max();
    Vec3T<T> rrhi  = -Vec3T<T>::max();
    int      rrcnt = 0;
 
    for (int b = BINS - 1; b >= 1; b--) {
      if (binCnt[b] > 0) {
        rrlo = min(rrlo, binLo[b]);
        rrhi = max(rrhi, binHi[b]);
      }
 
      rrcnt += binCnt[b];
 
      if (leftCnt[b - 1] > 0 && rrcnt > 0) {
        const T cost = leftArea[b - 1] * T(leftCnt[b - 1]) + BV(rrlo, rrhi).getArea() * T(rrcnt);
 
        if (cost < bestCost) {
          bestCost  = cost;
          bestAxis  = axis;
          bestPlane = lo + T(b) / scale;
        }
      }
    }
  }
 
  // All axes degenerate — equal split
  if (bestAxis < 0) {
    return a_begin + N / 2;
  }
 
  auto mid = std::partition(
    a_list.begin() + a_begin, a_list.begin() + a_end, [bestAxis, bestPlane](const PrimAndBV<P, BV>& pbv) noexcept {
      return pbv.second.getCentroid()[bestAxis] < bestPlane;
    });
 
  const size_t splitIdx = static_cast<size_t>(std::distance(a_list.begin(), mid));
 
  // Guard: if everything ended up on one side, fall back to equal split
  return (splitIdx == a_begin || splitIdx == a_end) ? a_begin + N / 2 : splitIdx;
}
 
template <class T, class P, class BV>
inline void
SAHKWaySplit(PrimAndBVList<P, BV>&                   a_list,
             const size_t                            a_begin,
             const size_t                            a_end,
             const size_t                            a_K,
             std::vector<std::pair<size_t, size_t>>& a_groups,
             const bool                              a_longestAxisOnly = false) noexcept
{
  static_assert(std::is_same_v<BV, EBGeometry::BoundingVolumes::AABBT<T>>, "SAHKWaySplit requires BV == AABBT<T>");
 
  if (a_K <= 1 || a_begin >= a_end) {
    a_groups.emplace_back(a_begin, a_end);
 
    return;
  }
 
  const size_t K1 = a_K / 2;
  const size_t K2 = a_K - K1;
 
  // Clamp the SAH split so the left half has >= K1 elements and the right has >= K2.
  // This guarantees neither recursive call receives an under-populated range, which
  // would cause empty sub-lists to reach the TreeBVH constructor and crash on
  // AABBT(vector::front()) when the vector is empty.
  // The clamp is valid whenever a_end - a_begin >= a_K = K1 + K2.
  const size_t rawMid = SAH2WaySplit<T, P, BV>(a_list, a_begin, a_end, a_longestAxisOnly);
  const size_t mid    = std::max(a_begin + K1, std::min(a_end - K2, rawMid));
 
  SAHKWaySplit<T, P, BV>(a_list, a_begin, mid, K1, a_groups, a_longestAxisOnly);
  SAHKWaySplit<T, P, BV>(a_list, mid, a_end, K2, a_groups, a_longestAxisOnly);
}
 
template <class T, class P, class BV, size_t K, bool LongestAxisOnly = false>
auto BinnedSAHPartitioner = [](PrimAndBVList<P, BV> a_primsAndBVs) -> std::array<PrimAndBVList<P, BV>, K> {
  EBGEOMETRY_EXPECT(!a_primsAndBVs.empty());
 
  // The input is taken by value; partition it in place (no working copy). SAHKWaySplit reorders it
  // and yields K [begin, end) index ranges, which we move out (disjoint, each moved once).
  std::vector<std::pair<size_t, size_t>> groups;
  groups.reserve(K);
 
  SAHKWaySplit<T, P, BV>(a_primsAndBVs, 0, a_primsAndBVs.size(), K, groups, LongestAxisOnly);
 
  std::array<PrimAndBVList<P, BV>, K> result;
  for (size_t k = 0; k < K; k++) {
    const auto [b, e] = groups[k];
    result[k]         = PrimAndBVList<P, BV>(std::make_move_iterator(a_primsAndBVs.begin() + b),
                                     std::make_move_iterator(a_primsAndBVs.begin() + e));
  }
 
  return result;
};
 
template <class T, class P, class BV>
inline size_t
Midpoint2WaySplit(PrimAndBVList<P, BV>& a_list, const size_t a_begin, const size_t a_end) noexcept
{
  EBGEOMETRY_EXPECT(a_begin < a_end);
 
  const size_t N = a_end - a_begin;
 
  Vec3T<T> lo = +Vec3T<T>::max();
  Vec3T<T> hi = -Vec3T<T>::max();
 
  for (size_t i = a_begin; i < a_end; i++) {
    const auto& c = a_list[i].second.getCentroid();
 
    lo = min(lo, c);
    hi = max(hi, c);
  }
 
  const size_t splitDir = (hi - lo).maxDir(true);
  const T      midpoint = T(0.5) * (lo[splitDir] + hi[splitDir]);
 
  auto mid = std::partition(
    a_list.begin() + a_begin, a_list.begin() + a_end, [splitDir, midpoint](const PrimAndBV<P, BV>& pbv) noexcept {
      return pbv.second.getCentroid()[splitDir] < midpoint;
    });
 
  const size_t splitIdx = static_cast<size_t>(std::distance(a_list.begin(), mid));
 
  // Guard: if every centroid ended up on one side (e.g. all coincide on the split axis, or the
  // distribution is skewed entirely to one side of the midpoint), fall back to an equal-count
  // split so neither resulting group is ever empty.
  return (splitIdx == a_begin || splitIdx == a_end) ? a_begin + N / 2 : splitIdx;
}
 
template <class T, class P, class BV>
inline void
MidpointKWaySplit(PrimAndBVList<P, BV>&                   a_list,
                  const size_t                            a_begin,
                  const size_t                            a_end,
                  const size_t                            a_K,
                  std::vector<std::pair<size_t, size_t>>& a_groups) noexcept
{
  if (a_K <= 1 || a_begin >= a_end) {
    a_groups.emplace_back(a_begin, a_end);
 
    return;
  }
 
  const size_t K1 = a_K / 2;
  const size_t K2 = a_K - K1;
 
  // Clamp so the left half has >= K1 elements and the right has >= K2, guaranteeing neither
  // recursive call receives an under-populated range (see SAHKWaySplit's identical clamp for why:
  // an empty sub-list reaching the TreeBVH constructor crashes on AABBT(vector::front())).
  const size_t rawMid = Midpoint2WaySplit<T, P, BV>(a_list, a_begin, a_end);
  const size_t mid    = std::max(a_begin + K1, std::min(a_end - K2, rawMid));
 
  MidpointKWaySplit<T, P, BV>(a_list, a_begin, mid, K1, a_groups);
  MidpointKWaySplit<T, P, BV>(a_list, mid, a_end, K2, a_groups);
}
 
template <class T, class P, class BV, size_t K>
auto MidpointPartitioner = [](PrimAndBVList<P, BV> a_primsAndBVs) -> std::array<PrimAndBVList<P, BV>, K> {
  EBGEOMETRY_EXPECT(!a_primsAndBVs.empty());
 
  // The input is taken by value; partition it in place (no working copy). MidpointKWaySplit reorders
  // it and yields K [begin, end) index ranges, which we move out (disjoint, each moved once). This
  // is the crux of the partitioner's cost: the midpoint split itself is a couple of std::partition
  // passes, so the plumbing (copies) is what dominated before -- now eliminated.
  std::vector<std::pair<size_t, size_t>> groups;
  groups.reserve(K);
 
  MidpointKWaySplit<T, P, BV>(a_primsAndBVs, 0, a_primsAndBVs.size(), K, groups);
 
  std::array<PrimAndBVList<P, BV>, K> result;
  for (size_t k = 0; k < K; k++) {
    const auto [b, e] = groups[k];
    result[k]         = PrimAndBVList<P, BV>(std::make_move_iterator(a_primsAndBVs.begin() + b),
                                     std::make_move_iterator(a_primsAndBVs.begin() + e));
  }
 
  return result;
};
 
template <class T, class P, class BV, size_t K>
auto DefaultLeafPredicate =
  [](const BVH::TreeBVH<T, P, BV, K>& a_node) noexcept -> bool { return (a_node.getPrimitives()).size() < K; };
 
template <class T, class P, class BV, size_t K>
class TreeBVH : public std::enable_shared_from_this<TreeBVH<T, P, BV, K>>
{
  static_assert(std::is_floating_point_v<T>, "TreeBVH: T must be a floating-point type");
  static_assert(K >= 2, "TreeBVH: branching factor K must be at least 2");
 
public:
  using PrimitiveList = BVH::PrimitiveList<P>;
 
  using Vec3 = Vec3T<T>;
 
  using Node = TreeBVH<T, P, BV, K>;
 
  using NodePtr = std::shared_ptr<Node>;
 
  using Partitioner = BVH::Partitioner<P, BV, K>;
 
  using LeafPredicate = BVH::LeafPredicate<T, P, BV, K>;
 
  TreeBVH() noexcept;
 
  TreeBVH(const std::vector<PrimAndBV<P, BV>>& a_primsAndBVs);
 
  TreeBVH(std::vector<PrimAndBV<P, BV>>&& a_primsAndBVs);
 
  ~TreeBVH() noexcept = default;
 
  TreeBVH(const TreeBVH& a_other) = delete;
 
  TreeBVH&
  operator=(const TreeBVH& a_other) = delete;
 
  TreeBVH(TreeBVH&& a_other) noexcept = default;
 
  TreeBVH&
  operator=(TreeBVH&& a_other) noexcept = default;
 
  inline void
  topDownSortAndPartition(const Partitioner&   a_partitioner = BVCentroidPartitioner<T, P, BV, K>,
                          const LeafPredicate& a_stopCrit    = DefaultLeafPredicate<T, P, BV, K>);
 
  template <typename S>
  inline void
  bottomUpSortAndPartition();
 
  [[nodiscard]] inline bool
  isLeaf() const noexcept;
 
  [[nodiscard]] inline bool
  isPartitioned() const noexcept;
 
  [[nodiscard]] inline const BV&
  getBoundingVolume() const noexcept;
 
  [[nodiscard]] inline const PrimitiveList&
  getPrimitives() const noexcept;
 
  [[nodiscard]] inline const std::vector<BV>&
  getBoundingVolumes() const noexcept;
 
  [[nodiscard]] inline T
  getDistanceToBoundingVolume(const Vec3& a_point) const noexcept;
 
  [[nodiscard]] inline const std::array<std::shared_ptr<TreeBVH<T, P, BV, K>>, K>&
  getChildren() const noexcept;
 
  [[nodiscard]] inline std::shared_ptr<TreeBVH<T, P, BV, K>>
  deepCopy() const;
 
  template <class NodeKey>
  inline void
  traverse(const BVH::LeafEvaluator<P>&               a_leafEvaluator,
           const BVH::PrunePredicate<Node, NodeKey>&  a_prunePredicate,
           const BVH::ChildOrderer<Node, NodeKey, K>& a_childOrderer,
           const BVH::NodeKeyFactory<Node, NodeKey>&  a_nodeKeyFactory) const noexcept;
 
  template <class StoragePolicy = BVH::SharedPtrStorage<P>>
  [[nodiscard]] inline std::shared_ptr<PackedBVH<T, P, K, StoragePolicy>>
  pack() const;
 
  template <class Q, class Converter, class StoragePolicy = BVH::SharedPtrStorage<Q>>
  [[nodiscard]] inline std::shared_ptr<PackedBVH<T, Q, K, StoragePolicy>>
  packWith(Converter&& a_converter) const;
 
protected:
  BV m_boundingVolume;
 
  bool m_partitioned;
 
  std::vector<std::shared_ptr<const P>> m_primitives;
 
  std::vector<BV> m_boundingVolumes;
 
  std::array<std::shared_ptr<TreeBVH<T, P, BV, K>>, K> m_children;
 
  inline PrimitiveList&
  getPrimitives() noexcept;
 
  inline std::vector<BV>&
  getBoundingVolumes() noexcept;
 
  inline void
  setChildren(const std::array<std::shared_ptr<TreeBVH<T, P, BV, K>>, K>& a_children) noexcept;
};
 
template <class T, class P, size_t K, class StoragePolicy>
class PackedBVH
{
  static_assert(std::is_floating_point_v<T>, "PackedBVH: T must be a floating-point type");
  static_assert(K >= 2, "PackedBVH: branching factor K must be at least 2");
 
public:
  using BV = EBGeometry::BoundingVolumes::AABBT<T>;
 
  using StorageType = typename StoragePolicy::StorageType;
 
  struct Node
  {
    BV m_bv{};
 
    uint32_t m_primOff{};
 
    uint32_t m_numPrims{};
 
    std::array<uint32_t, K> m_childOff{};
 
    inline void
    setBoundingVolume(const BV& a_bv) noexcept
    {
      m_bv = a_bv;
    }
 
    inline void
    setPrimitivesOffset(uint32_t a_off) noexcept
    {
      m_primOff = a_off;
    }
 
    inline void
    setNumPrimitives(uint32_t a_n) noexcept
    {
      m_numPrims = a_n;
    }
 
    inline void
    setChildOffset(uint32_t a_off, size_t a_k) noexcept
    {
      EBGEOMETRY_EXPECT(a_k < K);
      m_childOff[a_k] = a_off;
    }
 
    [[nodiscard]] inline const BV&
    getBoundingVolume() const noexcept
    {
      return m_bv;
    }
 
    [[nodiscard]] inline uint32_t
    getPrimitivesOffset() const noexcept
    {
      return m_primOff;
    }
 
    [[nodiscard]] inline uint32_t
    getNumPrimitives() const noexcept
    {
      return m_numPrims;
    }
 
    [[nodiscard]] inline const std::array<uint32_t, K>&
    getChildOffsets() const noexcept
    {
      return m_childOff;
    }
 
    [[nodiscard]] inline bool
    isLeaf() const noexcept
    {
      return m_numPrims > 0;
    }
 
    [[nodiscard]] inline T
    getDistanceToBoundingVolume(const Vec3T<T>& a_point) const noexcept
    {
      return m_bv.getDistance(a_point);
    }
 
    [[nodiscard]] inline T
    getDistanceToBoundingVolume2(const Vec3T<T>& a_point) const noexcept
    {
      return m_bv.getDistance2(a_point);
    }
  };
 
  PackedBVH() = delete;
 
  inline PackedBVH(const TreeBVH<T, P, BV, K>& a_tree);
 
  template <class Q, class Converter>
  inline PackedBVH(const TreeBVH<T, Q, BV, K>& a_tree, Converter&& a_converter);
 
  template <class S = SFC::Morton>
  inline PackedBVH(std::vector<std::pair<P, BV>> a_primsAndBVs, size_t a_targetLeafSize, S a_sfc = S{});
 
  inline PackedBVH(std::vector<std::pair<P, BV>>          a_primsAndBVs,
                   const BVH::Partitioner<P, BV, K>&      a_partitioner = BVCentroidPartitioner<T, P, BV, K>,
                   const BVH::LeafPredicate<T, P, BV, K>& a_stopCrit    = DefaultLeafPredicate<T, P, BV, K>);
 
  inline PackedBVH(std::vector<std::pair<P, BV>> a_primsAndBVs, BVH::ClusterSpec a_spec);
 
  inline ~PackedBVH() = default;
 
  PackedBVH(const PackedBVH& a_other) = default;
 
  PackedBVH&
  operator=(const PackedBVH& a_other) = default;
 
  PackedBVH(PackedBVH&& a_other) noexcept = default;
 
  PackedBVH&
  operator=(PackedBVH&& a_other) noexcept = default;
 
  [[nodiscard]] inline const std::vector<StorageType>&
  getPrimitives() const noexcept;
 
  [[nodiscard]] inline const BV&
  getBoundingVolume() const noexcept;
 
  [[nodiscard]] inline BV
  computeBoundingVolume() const noexcept;
 
  template <class NodeKey>
  inline void
  traverse(const BVH::PackedLeafEvaluator<P, StoragePolicy>& a_leafEvaluator,
           const BVH::PrunePredicate<Node, NodeKey>&         a_prunePredicate,
           const BVH::PackedChildOrderer<NodeKey, K>&        a_childOrderer,
           const BVH::NodeKeyFactory<Node, NodeKey>&         a_nodeKeyFactory) const noexcept;
 
  template <class State, class LeafEvaluator, class PruneDistSquared>
  inline void
  pruneTraverse(const Vec3T<T>&    a_point,
                State&             a_state,
                LeafEvaluator&&    a_evalLeaf,
                PruneDistSquared&& a_pruneDist2) const noexcept;
 
protected:
  inline PackedBVH(std::vector<Node>&& a_linearNodes, std::vector<StorageType>&& a_primitives) noexcept
    : m_linearNodes(std::move(a_linearNodes)), m_primitives(std::move(a_primitives))
  {
    this->buildSoA();
  }
 
  std::vector<Node> m_linearNodes;
 
  std::vector<StorageType> m_primitives;
 
  struct ChildAABBSoA
  {
    alignas(sizeof(T) * K) T m_lo[3][K];
 
    alignas(sizeof(T) * K) T m_hi[3][K];
  };
 
  std::vector<ChildAABBSoA> m_childAabbSoA;
 
  inline void
  buildSoA();
};
 
} // namespace BVH
 
} // namespace EBGeometry
 
#include "EBGeometry_BVHImplem.hpp"
 
#endif