rtree.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2018 - 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_numerics_rtree_h
#define dealii_numerics_rtree_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/point.h>
#include <deal.II/base/std_cxx20/iota_view.h>
 
#include <deal.II/boost_adaptors/bounding_box.h>
#include <deal.II/boost_adaptors/point.h>
#include <deal.II/boost_adaptors/segment.h>
 
DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
#include <boost/geometry/algorithms/buffer.hpp>
#include <boost/geometry/algorithms/distance.hpp>
#include <boost/geometry/index/rtree.hpp>
#include <boost/geometry/strategies/strategies.hpp>
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
 
#include <memory>
#include <type_traits>
#include <utility>
 
 
DEAL_II_NAMESPACE_OPEN
 
namespace internal
{
  namespace RTreeImplementation
  {
    struct EmptyVisitor
    {};
  } // namespace RTreeImplementation
} // namespace internal
 
template <typename LeafType,
          typename IndexType = boost::geometry::index::linear<16>,
          typename IndexableGetter =
            boost::geometry::index::indexable<LeafType>>
using RTree =
  boost::geometry::index::rtree<LeafType, IndexType, IndexableGetter>;
 
template <typename RTreeType>
using RTreeView =
  boost::geometry::index::detail::rtree::utilities::view<RTreeType>;
 
template <typename IndexType = boost::geometry::index::linear<16>,
          typename LeafTypeIterator,
          typename IndexableGetter = boost::geometry::index::indexable<
            typename LeafTypeIterator::value_type>>
RTree<typename LeafTypeIterator::value_type, IndexType, IndexableGetter>
pack_rtree(const LeafTypeIterator &begin, const LeafTypeIterator &end);
 
template <typename IndexType = boost::geometry::index::linear<16>,
          typename ContainerType,
          typename IndexableGetter = boost::geometry::index::indexable<
            typename ContainerType::value_type>>
RTree<typename ContainerType::value_type, IndexType, IndexableGetter>
pack_rtree(const ContainerType &container);
 
template <typename Container>
class IndexableGetterFromIndices
{
public:
  using IndexableGetter =
    typename boost::geometry::index::indexable<typename Container::value_type>;
 
  using result_type = typename IndexableGetter::result_type;
 
  using size_t = typename Container::size_type;
 
  explicit IndexableGetterFromIndices(const Container &c)
    : container(c)
  {}
 
  result_type
  operator()(size_t i) const
  {
    return getter(container[i]);
  }
 
private:
  const Container &container;
 
  IndexableGetter getter;
};
 
template <typename IndexType = boost::geometry::index::linear<16>,
          typename ContainerType>
RTree<typename ContainerType::size_type,
      IndexType,
      IndexableGetterFromIndices<ContainerType>>
pack_rtree_of_indices(const ContainerType &container);
 
template <
  typename RTreeType,
  typename InternalVisitor  = internal::RTreeImplementation::EmptyVisitor,
  typename LeafVisitor      = internal::RTreeImplementation::EmptyVisitor,
  typename ElementVisitor   = internal::RTreeImplementation::EmptyVisitor,
  typename IndexableVisitor = internal::RTreeImplementation::EmptyVisitor>
struct RTreeFunctionalVisitor
  : public boost::geometry::index::detail::rtree::visitor<
      typename RTreeView<RTreeType>::value_type,
      typename RTreeView<RTreeType>::options_type::parameters_type,
      typename RTreeView<RTreeType>::box_type,
      typename RTreeView<RTreeType>::allocators_type,
      typename RTreeView<RTreeType>::options_type::node_tag,
      true>
{
  static constexpr unsigned int dim =
    boost::geometry::dimension<typename RTreeType::indexable_type>::value;
 
  using internal_node_type =
    typename boost::geometry::index::detail::rtree::internal_node<
      typename RTreeView<RTreeType>::value_type,
      typename RTreeView<RTreeType>::options_type::parameters_type,
      typename RTreeView<RTreeType>::box_type,
      typename RTreeView<RTreeType>::allocators_type,
      typename RTreeView<RTreeType>::options_type::node_tag>::type;
 
  using leaf_type = typename boost::geometry::index::detail::rtree::leaf<
    typename RTreeView<RTreeType>::value_type,
    typename RTreeView<RTreeType>::options_type::parameters_type,
    typename RTreeView<RTreeType>::box_type,
    typename RTreeView<RTreeType>::allocators_type,
    typename RTreeView<RTreeType>::options_type::node_tag>::type;
 
  using value_type = typename RTreeType::value_type;
 
  using indexable_type = typename RTreeType::indexable_type;
 
  using translator_type = typename RTreeView<RTreeType>::translator_type;
 
  using empty_visitor = internal::RTreeImplementation::EmptyVisitor;
 
  static constexpr bool has_internal_visitor =
    !std::is_same_v<InternalVisitor, empty_visitor>;
 
  static constexpr bool has_leaf_visitor =
    !std::is_same_v<LeafVisitor, empty_visitor>;
 
  static constexpr bool has_element_visitor =
    !std::is_same_v<ElementVisitor, empty_visitor>;
 
  static constexpr bool has_indexable_visitor =
    !std::is_same_v<IndexableVisitor, empty_visitor>;
 
  InternalVisitor internal_node_visitor;
 
  LeafVisitor leaf_visitor;
 
  ElementVisitor element_visitor;
 
  IndexableVisitor indexable_visitor;
 
  void
  operator()(internal_node_type &node);
 
  void
  operator()(leaf_type &leaf);
 
  RTreeFunctionalVisitor(const RTreeType &tree,
                         InternalVisitor  internal_node_visitor = {},
                         LeafVisitor      leaf_visitor          = {},
                         ElementVisitor   element_visitor       = {},
                         IndexableVisitor indexable_visitor     = {});
 
  std::size_t     current_level;
  translator_type translator;
};
 
template <
  typename RTreeType,
  typename InternalVisitor  = internal::RTreeImplementation::EmptyVisitor,
  typename LeafVisitor      = internal::RTreeImplementation::EmptyVisitor,
  typename ElementVisitor   = internal::RTreeImplementation::EmptyVisitor,
  typename IndexableVisitor = internal::RTreeImplementation::EmptyVisitor>
void
visit_rtree(const RTreeType &tree,
            InternalVisitor  internal_node_visitor = {},
            LeafVisitor      leaf_visitor          = {},
            ElementVisitor   element_visitor       = {},
            IndexableVisitor indexable_visitor     = {});
 
 
template <typename Rtree>
inline std::vector<BoundingBox<
  boost::geometry::dimension<typename Rtree::indexable_type>::value>>
extract_rtree_level(const Rtree &tree, const unsigned int level);
 
 
 
template <typename Rtree>
std::vector<std::vector<BoundingBox<
  boost::geometry::dimension<typename Rtree::indexable_type>::value>>>
extract_children_of_level(const Rtree &tree, const unsigned int level);
 
 
 
// Inline and template functions
#ifndef DOXYGEN
 
template <typename IndexType,
          typename LeafTypeIterator,
          typename IndexableGetter>
RTree<typename LeafTypeIterator::value_type, IndexType, IndexableGetter>
pack_rtree(const LeafTypeIterator &begin, const LeafTypeIterator &end)
{
  return RTree<typename LeafTypeIterator::value_type,
               IndexType,
               IndexableGetter>(begin, end);
}
 
 
 
template <typename IndexType, typename ContainerType, typename IndexableGetter>
RTree<typename ContainerType::value_type, IndexType, IndexableGetter>
pack_rtree(const ContainerType &container)
{
  return pack_rtree<IndexType, decltype(container.begin()), IndexableGetter>(
    container.begin(), container.end());
}
 
 
 
template <typename IndexType, typename ContainerType>
RTree<typename ContainerType::size_type,
      IndexType,
      IndexableGetterFromIndices<ContainerType>>
pack_rtree_of_indices(const ContainerType &container)
{
  // We need an array that holds the indices we want to pack. The rtree
  // implementation in BOOST, for reasons not entirely clear, insists
  // on using a reference to the elements of the range. This is fine if
  // the indices are stored in a container, so that's what we do.
  // (It would be nice if we could just pass a std::ranges::iota_view
  // instead, but that has no nested 'reference' type, and this then
  // trips up BOOST rtree.)
  std::vector<typename ContainerType::size_type> indices(container.size());
  for (typename ContainerType::size_type i = 0; i < container.size(); ++i)
    indices[i] = i;
 
  return RTree<typename ContainerType::size_type,
               IndexType,
               IndexableGetterFromIndices<ContainerType>>(
    indices.begin(),
    indices.end(),
    IndexType(),
    IndexableGetterFromIndices<ContainerType>(container));
}
 
 
 
template <typename RTreeType,
          typename InternalVisitor,
          typename LeafVisitor,
          typename ElementVisitor,
          typename IndexableVisitor>
void
RTreeFunctionalVisitor<RTreeType,
                       InternalVisitor,
                       LeafVisitor,
                       ElementVisitor,
                       IndexableVisitor>::operator()(internal_node_type &node)
{
  const std::size_t level_backup = current_level;
  for (const auto &element :
       boost::geometry::index::detail::rtree::elements(node))
    {
      BoundingBox<dim> box;
      boost::geometry::convert(element.first, box);
 
      const std::size_t child_level = level_backup;
 
      if (auto *leaf_ptr = boost::get<leaf_type>(&*element.second))
        {
          if constexpr (has_leaf_visitor)
            leaf_visitor(child_level, node, box, *leaf_ptr);
        }
      else if (auto *internal_ptr =
                 boost::get<internal_node_type>(&*element.second))
        {
          if constexpr (has_internal_visitor)
            internal_node_visitor(child_level, node, box, *internal_ptr);
        }
 
      current_level = child_level + 1;
      boost::geometry::index::detail::rtree::apply_visitor(*this,
                                                           *element.second);
    }
  current_level = level_backup;
}
 
 
 
template <typename RTreeType,
          typename InternalVisitor,
          typename LeafVisitor,
          typename ElementVisitor,
          typename IndexableVisitor>
void
RTreeFunctionalVisitor<RTreeType,
                       InternalVisitor,
                       LeafVisitor,
                       ElementVisitor,
                       IndexableVisitor>::operator()(leaf_type &leaf)
{
  for (const auto &element :
       boost::geometry::index::detail::rtree::elements(leaf))
    {
      if constexpr (has_element_visitor)
        element_visitor(leaf, element);
      if constexpr (has_indexable_visitor)
        {
          auto idx = translator(element);
          indexable_visitor(leaf, idx);
        }
    }
}
 
 
 
template <typename RTreeType,
          typename InternalVisitor,
          typename LeafVisitor,
          typename ElementVisitor,
          typename IndexableVisitor>
RTreeFunctionalVisitor<RTreeType,
                       InternalVisitor,
                       LeafVisitor,
                       ElementVisitor,
                       IndexableVisitor>::
  RTreeFunctionalVisitor(const RTreeType &tree,
                         InternalVisitor  internal_node_visitor,
                         LeafVisitor      leaf_visitor,
                         ElementVisitor   element_visitor,
                         IndexableVisitor indexable_visitor)
  : internal_node_visitor(std::move(internal_node_visitor))
  , leaf_visitor(std::move(leaf_visitor))
  , element_visitor(std::move(element_visitor))
  , indexable_visitor(std::move(indexable_visitor))
  , current_level(0)
  , translator(RTreeView<RTreeType>(tree).translator())
{
  static_assert(!has_internal_visitor ||
                  std::is_invocable_r_v<void,
                                        InternalVisitor,
                                        const std::size_t &,
                                        const internal_node_type &,
                                        const BoundingBox<dim> &,
                                        const internal_node_type &>,
                "internal_node_visitor must be callable with (std::size_t, "
                "internal_node_type, BoundingBox<dim>, internal_node_type)");
 
  static_assert(
    !has_leaf_visitor || std::is_invocable_r_v<void,
                                               LeafVisitor,
                                               const std::size_t &,
                                               const internal_node_type &,
                                               const BoundingBox<dim> &,
                                               const leaf_type &>,
    "leaf_visitor must be callable with (std::size_t, internal_node_type, "
    "BoundingBox<dim>, leaf_type)");
 
  static_assert(
    !has_element_visitor || std::is_invocable_r_v<void,
                                                  ElementVisitor,
                                                  const leaf_type &,
                                                  const value_type &>,
    "element_visitor must be callable with (leaf_type, value_type)");
 
  static_assert(
    !has_indexable_visitor || std::is_invocable_r_v<void,
                                                    IndexableVisitor,
                                                    const leaf_type &,
                                                    const indexable_type &>,
    "indexable_visitor must be callable with (leaf_type, indexable_type)");
 
  RTreeView<RTreeType> rtv(tree);
  rtv.apply_visitor(*this);
}
 
 
 
template <typename RTreeType,
          typename InternalVisitor,
          typename LeafVisitor,
          typename ElementVisitor,
          typename IndexableVisitor>
void
visit_rtree(const RTreeType &tree,
            InternalVisitor  internal_node_visitor,
            LeafVisitor      leaf_visitor,
            ElementVisitor   element_visitor,
            IndexableVisitor indexable_visitor)
{
  RTreeFunctionalVisitor<RTreeType,
                         std::decay_t<InternalVisitor>,
                         std::decay_t<LeafVisitor>,
                         std::decay_t<ElementVisitor>,
                         std::decay_t<IndexableVisitor>>
    visitor(tree,
            std::move(internal_node_visitor),
            std::move(leaf_visitor),
            std::move(element_visitor),
            std::move(indexable_visitor));
}
 
 
 
template <typename Rtree>
inline std::vector<BoundingBox<
  boost::geometry::dimension<typename Rtree::indexable_type>::value>>
extract_rtree_level(const Rtree &tree, const unsigned int level)
{
  constexpr unsigned int dim =
    boost::geometry::dimension<typename Rtree::indexable_type>::value;
 
  std::vector<BoundingBox<dim>> boxes;
 
  unsigned int n_levels_in_tree = n_levels(tree);
 
  if (n_levels_in_tree == 0)
    {
      // The algorithm below does not work for `n_levels_in_tree==0`, which
      // might happen if the number entries in the tree is too small. In this
      // case, simply return a single bounding box.
      boxes.resize(1);
      boost::geometry::convert(tree.bounds(), boxes[0]);
    }
  else
    {
      const unsigned int target_level =
        std::min<unsigned int>(level, n_levels_in_tree - 1);
 
      using Visitor           = RTreeFunctionalVisitor<Rtree>;
      using InternalNode      = typename Visitor::internal_node_type;
      using LeafNode          = typename Visitor::leaf_type;
      const auto node_visitor = [&](const std::size_t &current_level,
                                    const InternalNode &,
                                    const BoundingBox<dim> &box,
                                    const InternalNode &) {
        if (current_level == target_level)
          boxes.push_back(box);
      };
      const auto leaf_visitor = [&](const std::size_t &current_level,
                                    const InternalNode &,
                                    const BoundingBox<dim> &box,
                                    const LeafNode &) {
        if (current_level == target_level)
          boxes.push_back(box);
      };
      visit_rtree(tree, node_visitor, leaf_visitor);
    }
 
  return boxes;
}
 
 
 
template <class Rtree>
unsigned int
n_levels(const Rtree &tree)
{
  boost::geometry::index::detail::rtree::utilities::view<Rtree> rtv(tree);
  return rtv.depth();
}
 
 
template <typename Rtree>
inline std::vector<std::vector<BoundingBox<
  boost::geometry::dimension<typename Rtree::indexable_type>::value>>>
extract_children_of_level(const Rtree &tree, const unsigned int level)
{
  constexpr unsigned int dim =
    boost::geometry::dimension<typename Rtree::indexable_type>::value;
 
  unsigned int n_levels_in_tree = n_levels(tree);
 
  std::vector<std::vector<BoundingBox<dim>>> boxes_in_boxes;
 
  if (n_levels_in_tree == 0)
    {
      // The below algorithm does not work for `rtv.depth()==0`, which might
      // happen if the number entries in the tree is too small.
      // In this case, simply return a single bounding box.
      boxes_in_boxes.resize(1);
      boxes_in_boxes[0].resize(1);
      boost::geometry::convert(tree.bounds(), boxes_in_boxes[0][0]);
    }
  else
    {
      const unsigned int target_level =
        std::min<unsigned int>(level, n_levels_in_tree - 1);
 
      // Track all nodes on target_level and collect their children (level
      // target_level+1) grouped per parent. Depth-first visitation guarantees
      // children immediately follow their parent, so a simple counter suffices.
      std::size_t current_parent_index = 0;
      std::size_t node_counter         = 0;
 
      const auto internal = [&](const std::size_t current_level,
                                const auto &,
                                const BoundingBox<dim> &child_box,
                                const auto &) {
        // Register nodes that live on the target level so their children can
        // be grouped.
        if (current_level == target_level)
          {
            current_parent_index = node_counter++;
            boxes_in_boxes.emplace_back();
            return;
          }
 
        // Record children of nodes on the target level.
        if (current_level == target_level + 1)
          boxes_in_boxes[current_parent_index].push_back(child_box);
      };
 
      const auto leaf = [&](const std::size_t current_level,
                            const auto &,
                            const BoundingBox<dim> &child_box,
                            const auto &) {
        if (current_level == target_level + 1)
          boxes_in_boxes[current_parent_index].push_back(child_box);
      };
 
      visit_rtree(tree, internal, leaf);
    }
 
  return boxes_in_boxes;
}
 
 
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif