rtree.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2017 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.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>
 
#include <boost/geometry/algorithms/distance.hpp>
#include <boost/geometry/index/rtree.hpp>
#include <boost/geometry/strategies/strategies.hpp>
 
#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(Container const &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(
    Translator const & 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()(InternalNode const &node);
 
  inline void
  operator()(Leaf const &);
 
  Translator const &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);
 
 
 
// 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)
{
  std_cxx20::ranges::iota_view<typename ContainerType::size_type,
                               typename ContainerType::size_type>
    indices(0, container.size());
  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();
}
 
 
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif