rtree.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2018 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#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/distance.hpp>
#include <boost/geometry/index/rtree.hpp>
#include <boost/geometry/strategies/strategies.hpp>
DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
 
#include <memory>
 
 
DEAL_II_NAMESPACE_OPEN
 
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 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 Value,
          typename Options,
          typename Translator,
          typename Box,
          typename Allocators>
struct ExtractLevelVisitor
  : public boost::geometry::index::detail::rtree::visitor<
      Value,
      typename Options::parameters_type,
      Box,
      Allocators,
      typename Options::node_tag,
      true>::type
{
  inline ExtractLevelVisitor(
    const Translator  &translator,
    const unsigned int target_level,
    std::vector<BoundingBox<boost::geometry::dimension<Box>::value>> &boxes);
 
  using InternalNode =
    typename boost::geometry::index::detail::rtree::internal_node<
      Value,
      typename Options::parameters_type,
      Box,
      Allocators,
      typename Options::node_tag>::type;
 
  using Leaf = typename boost::geometry::index::detail::rtree::leaf<
    Value,
    typename Options::parameters_type,
    Box,
    Allocators,
    typename Options::node_tag>::type;
 
  inline void
  operator()(const InternalNode &node);
 
  inline void
  operator()(const Leaf &);
 
  const Translator &translator;
 
  size_t level;
 
  const size_t target_level;
 
  std::vector<BoundingBox<boost::geometry::dimension<Box>::value>> &boxes;
};
 
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 Value,
          typename Options,
          typename Translator,
          typename Box,
          typename Allocators>
ExtractLevelVisitor<Value, Options, Translator, Box, Allocators>::
  ExtractLevelVisitor(
    const Translator  &translator,
    const unsigned int target_level,
    std::vector<BoundingBox<boost::geometry::dimension<Box>::value>> &boxes)
  : translator(translator)
  , level(0)
  , target_level(target_level)
  , boxes(boxes)
{}
 
 
 
template <typename Value,
          typename Options,
          typename Translator,
          typename Box,
          typename Allocators>
void
ExtractLevelVisitor<Value, Options, Translator, Box, Allocators>::operator()(
  const ExtractLevelVisitor::InternalNode &node)
{
  using ElementsType =
    typename boost::geometry::index::detail::rtree::elements_type<
      InternalNode>::type;
 
  const auto &elements = boost::geometry::index::detail::rtree::elements(node);
 
  if (level == target_level)
    {
      const auto offset = boxes.size();
      boxes.resize(offset + elements.size());
 
      unsigned int i = offset;
      for (typename ElementsType::const_iterator it = elements.begin();
           it != elements.end();
           ++it)
        {
          boost::geometry::convert(it->first, boxes[i]);
          ++i;
        }
      return;
    }
 
  const size_t level_backup = level;
  ++level;
 
  for (typename ElementsType::const_iterator it = elements.begin();
       it != elements.end();
       ++it)
    {
      boost::geometry::index::detail::rtree::apply_visitor(*this, *it->second);
    }
 
  level = level_backup;
}
 
template <typename Value,
          typename Options,
          typename Translator,
          typename Box,
          typename Allocators>
void
ExtractLevelVisitor<Value, Options, Translator, Box, Allocators>::operator()(
  const ExtractLevelVisitor::Leaf &)
{}
 
 
 
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;
 
  using RtreeView =
    boost::geometry::index::detail::rtree::utilities::view<Rtree>;
  RtreeView rtv(tree);
 
  std::vector<BoundingBox<dim>> boxes;
 
  if (rtv.depth() == 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.resize(1);
      boost::geometry::convert(tree.bounds(), boxes[0]);
    }
  else
    {
      const unsigned int target_level =
        std::min<unsigned int>(level, rtv.depth() - 1);
 
      ExtractLevelVisitor<typename RtreeView::value_type,
                          typename RtreeView::options_type,
                          typename RtreeView::translator_type,
                          typename RtreeView::box_type,
                          typename RtreeView::allocators_type>
        extract_level_visitor(rtv.translator(), target_level, boxes);
      rtv.apply_visitor(extract_level_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 Value,
          typename Options,
          typename Translator,
          typename Box,
          typename Allocators>
struct NodeVisitor : public boost::geometry::index::detail::rtree::visitor<
                       Value,
                       typename Options::parameters_type,
                       Box,
                       Allocators,
                       typename Options::node_tag,
                       true>::type
{
  inline NodeVisitor(
    const Translator &translator,
    unsigned int      target_level,
    std::vector<
      std::vector<BoundingBox<boost::geometry::dimension<Box>::value>>> &boxes);
 
 
  using InternalNode =
    typename boost::geometry::index::detail::rtree::internal_node<
      Value,
      typename Options::parameters_type,
      Box,
      Allocators,
      typename Options::node_tag>::type;
 
  using Leaf = typename boost::geometry::index::detail::rtree::leaf<
    Value,
    typename Options::parameters_type,
    Box,
    Allocators,
    typename Options::node_tag>::type;
 
  inline void
  operator()(const InternalNode &node);
 
  inline void
  operator()(const Leaf &);
 
  const Translator &translator;
 
  size_t level;
 
  size_t node_counter;
 
  const size_t target_level;
 
  std::vector<std::vector<BoundingBox<boost::geometry::dimension<Box>::value>>>
    &boxes_in_boxes;
};
 
 
 
template <typename Value,
          typename Options,
          typename Translator,
          typename Box,
          typename Allocators>
NodeVisitor<Value, Options, Translator, Box, Allocators>::NodeVisitor(
  const Translator  &translator,
  const unsigned int target_level,
  std::vector<std::vector<BoundingBox<boost::geometry::dimension<Box>::value>>>
    &bb_in_boxes)
  : translator(translator)
  , level(0)
  , node_counter(0)
  , target_level(target_level)
  , boxes_in_boxes(bb_in_boxes)
{}
 
 
 
template <typename Value,
          typename Options,
          typename Translator,
          typename Box,
          typename Allocators>
void
NodeVisitor<Value, Options, Translator, Box, Allocators>::operator()(
  const NodeVisitor::InternalNode &node)
{
  using elements_type =
    typename boost::geometry::index::detail::rtree::elements_type<
      InternalNode>::type; //  pairs of bounding box and pointer to child node
  const elements_type &elements =
    boost::geometry::index::detail::rtree::elements(node);
 
  if (level == target_level)
    {
      const unsigned int n_children = elements.size();
      const auto         offset     = boxes_in_boxes.size();
      boxes_in_boxes.resize(offset + n_children);
    }
 
  if (level == target_level + 1)
    {
      // I have now access to children of level target_level
      boxes_in_boxes[node_counter].resize(elements.size()); // number of bboxes
      unsigned int i = 0;
      for (typename elements_type::const_iterator it = elements.begin();
           it != elements.end();
           ++it)
        {
          boost::geometry::convert(it->first, boxes_in_boxes[node_counter][i]);
          ++i;
        }
      // Children have been stored, go to the next parent.
      ++node_counter;
      return;
    }
 
  size_t level_backup = level;
  ++level;
 
  for (typename elements_type::const_iterator it = elements.begin();
       it != elements.end();
       ++it)
    {
      boost::geometry::index::detail::rtree::apply_visitor(*this, *it->second);
    }
 
  level = level_backup;
}
 
 
 
template <typename Value,
          typename Options,
          typename Translator,
          typename Box,
          typename Allocators>
void
NodeVisitor<Value, Options, Translator, Box, Allocators>::operator()(
  const NodeVisitor::Leaf &)
{
  // No children for leaf nodes.
  boxes_in_boxes.clear();
}
 
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;
 
  using RtreeView =
    boost::geometry::index::detail::rtree::utilities::view<Rtree>;
  RtreeView rtv(tree);
 
  std::vector<std::vector<BoundingBox<dim>>> boxes_in_boxes;
 
  if (rtv.depth() == 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, rtv.depth() - 1);
 
      NodeVisitor<typename RtreeView::value_type,
                  typename RtreeView::options_type,
                  typename RtreeView::translator_type,
                  typename RtreeView::box_type,
                  typename RtreeView::allocators_type>
        node_visitor(rtv.translator(), target_level, boxes_in_boxes);
      rtv.apply_visitor(node_visitor);
    }
 
  return boxes_in_boxes;
}
 
 
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif