detail/rtree/pack_create.hpp

0001 // Boost.Geometry Index
0002 //
0003 // R-tree initial packing
0004 //
0005 // Copyright (c) 2011-2023 Adam Wulkiewicz, Lodz, Poland.
0006 // Copyright (c) 2020 Caian Benedicto, Campinas, Brazil.
0007 //
0008 // This file was modified by Oracle on 2019-2023.
0009 // Modifications copyright (c) 2019-2023 Oracle and/or its affiliates.
0010 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0011 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0012 //
0013 // Use, modification and distribution is subject to the Boost Software License,
0014 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0015 // http://www.boost.org/LICENSE_1_0.txt)
0016
0017 #ifndef BOOST_GEOMETRY_INDEX_DETAIL_RTREE_PACK_CREATE_HPP
0018 #define BOOST_GEOMETRY_INDEX_DETAIL_RTREE_PACK_CREATE_HPP
0019
0020 #include <boost/core/ignore_unused.hpp>
0021
0022 #include <boost/geometry/algorithms/centroid.hpp>
0023 #include <boost/geometry/algorithms/detail/expand_by_epsilon.hpp>
0024 #include <boost/geometry/algorithms/expand.hpp>
0025
0026 #include <boost/geometry/index/detail/algorithms/bounds.hpp>
0027 #include <boost/geometry/index/detail/algorithms/content.hpp>
0028 #include <boost/geometry/index/detail/algorithms/is_valid.hpp>
0029 #include <boost/geometry/index/detail/algorithms/nth_element.hpp>
0030 #include <boost/geometry/index/detail/rtree/node/node_elements.hpp>
0031 #include <boost/geometry/index/detail/rtree/node/subtree_destroyer.hpp>
0032 #include <boost/geometry/index/parameters.hpp>
0033
0034 #include <boost/geometry/util/constexpr.hpp>
0035
0036 namespace boost { namespace geometry { namespace index { namespace detail { namespace rtree {
0037
0038 namespace pack_utils {
0039
0040 template <std::size_t Dimension>
0041 struct biggest_edge
0042 {
0043     BOOST_STATIC_ASSERT(0 < Dimension);
0044     template <typename Box>
0045     static inline void apply(Box const& box, typename coordinate_type<Box>::type & length, std::size_t & dim_index)
0046     {
0047         biggest_edge<Dimension-1>::apply(box, length, dim_index);
0048         typename coordinate_type<Box>::type curr
0049             = geometry::get<max_corner, Dimension-1>(box) - geometry::get<min_corner, Dimension-1>(box);
0050         if ( length < curr )
0051         {
0052             dim_index = Dimension - 1;
0053             length = curr;
0054         }
0055     }
0056 };
0057
0058 template <>
0059 struct biggest_edge<1>
0060 {
0061     template <typename Box>
0062     static inline void apply(Box const& box, typename coordinate_type<Box>::type & length, std::size_t & dim_index)
0063     {
0064         dim_index = 0;
0065         length = geometry::get<max_corner, 0>(box) - geometry::get<min_corner, 0>(box);
0066     }
0067 };
0068
0069 template <std::size_t I>
0070 struct point_entries_comparer
0071 {
0072     template <typename PointEntry>
0073     bool operator()(PointEntry const& e1, PointEntry const& e2) const
0074     {
0075         return geometry::get<I>(e1.first) < geometry::get<I>(e2.first);
0076     }
0077 };
0078
0079 template <std::size_t I, std::size_t Dimension>
0080 struct nth_element_and_half_boxes
0081 {
0082     template <typename EIt, typename Box>
0083     static inline void apply(EIt first, EIt median, EIt last, Box const& box,
0084                              Box & left, Box & right, std::size_t dim_index)
0085     {
0086         if (I == dim_index)
0087         {
0088             index::detail::nth_element(first, median, last, point_entries_comparer<I>());
0089
0090             geometry::convert(box, left);
0091             geometry::convert(box, right);
0092             auto const mi = geometry::get<min_corner, I>(box);
0093             auto const ma = geometry::get<max_corner, I>(box);
0094             auto const center = mi + (ma - mi) / 2;
0095             geometry::set<max_corner, I>(left, center);
0096             geometry::set<min_corner, I>(right, center);
0097         }
0098         else
0099         {
0100             nth_element_and_half_boxes
0101                 <
0102                     I + 1, Dimension
0103                 >::apply(first, median, last, box, left, right, dim_index);
0104         }
0105     }
0106 };
0107
0108 template <std::size_t Dimension>
0109 struct nth_element_and_half_boxes<Dimension, Dimension>
0110 {
0111     template <typename EIt, typename Box>
0112     static inline void apply(EIt , EIt , EIt , Box const& , Box & , Box & , std::size_t ) {}
0113 };
0114
0115 } // namespace pack_utils
0116
0117 // STR leafs number are calculated as rcount/max
0118 // and the number of splitting planes for each dimension as (count/max)^(1/dimension)
0119 // <-> for dimension==2 -> sqrt(count/max)
0120 //
0121 // The main flaw of this algorithm is that the resulting tree will have bad structure for:
0122 // 1. non-uniformly distributed elements
0123 //      Statistic check could be performed, e.g. based on variance of lengths of elements edges for each dimension
0124 // 2. elements distributed mainly along one axis
0125 //      Calculate bounding box of all elements and then number of dividing planes for a dimension
0126 //      from the length of BB edge for this dimension (more or less assuming that elements are uniformly-distributed squares)
0127 //
0128 // Another thing is that the last node may have less elements than Max or even Min.
0129 // The number of splitting planes must be chosen more carefully than count/max
0130 //
0131 // This algorithm is something between STR and TGS
0132 // it is more similar to the top-down recursive kd-tree creation algorithm
0133 // using object median split and split axis of greatest BB edge
0134 // BB is only used as a hint (assuming objects are distributed uniformly)
0135 //
0136 // Implemented algorithm guarantees that the number of elements in nodes will be between Min and Max
0137 // and that nodes are packed as tightly as possible
0138 // e.g. for 177 values Max = 5 and Min = 2 it will construct the following tree:
0139 // ROOT                 177
0140 // L1          125               52
0141 // L2  25  25  25  25  25   25  17    10
0142 // L3  5x5 5x5 5x5 5x5 5x5  5x5 3x5+2 2x5
0143
0144 template <typename MembersHolder>
0145 class pack
0146 {
0147     typedef typename MembersHolder::node node;
0148     typedef typename MembersHolder::internal_node internal_node;
0149     typedef typename MembersHolder::leaf leaf;
0150
0151     typedef typename MembersHolder::node_pointer node_pointer;
0152     typedef typename MembersHolder::size_type size_type;
0153     typedef typename MembersHolder::parameters_type parameters_type;
0154     typedef typename MembersHolder::translator_type translator_type;
0155     typedef typename MembersHolder::allocators_type allocators_type;
0156
0157     typedef typename MembersHolder::box_type box_type;
0158     typedef typename geometry::point_type<box_type>::type point_type;
0159     typedef typename geometry::coordinate_type<point_type>::type coordinate_type;
0160     typedef typename detail::default_content_result<box_type>::type content_type;
0161     typedef typename detail::strategy_type<parameters_type>::type strategy_type;
0162     static const std::size_t dimension = geometry::dimension<point_type>::value;
0163
0164     typedef typename rtree::container_from_elements_type<
0165         typename rtree::elements_type<leaf>::type,
0166         size_type
0167     >::type values_counts_container;
0168
0169     typedef typename rtree::elements_type<internal_node>::type internal_elements;
0170     typedef typename internal_elements::value_type internal_element;
0171
0172     typedef rtree::subtree_destroyer<MembersHolder> subtree_destroyer;
0173
0174 public:
0175     // Arbitrary iterators
0176     template <typename InIt> inline static
0177     node_pointer apply(InIt first, InIt last,
0178                        size_type & values_count,
0179                        size_type & leafs_level,
0180                        parameters_type const& parameters,
0181                        translator_type const& translator,
0182                        allocators_type & allocators)
0183     {
0184         return apply(first, last, values_count, leafs_level, parameters, translator,
0185                      allocators, boost::container::new_allocator<void>());
0186     }
0187
0188     template <typename InIt, typename TmpAlloc> inline static
0189     node_pointer apply(InIt first, InIt last,
0190                        size_type & values_count,
0191                        size_type & leafs_level,
0192                        parameters_type const& parameters,
0193                        translator_type const& translator,
0194                        allocators_type & allocators,
0195                        TmpAlloc const& temp_allocator)
0196     {
0197         typedef typename std::iterator_traits<InIt>::difference_type diff_type;
0198
0199         diff_type diff = std::distance(first, last);
0200         if ( diff <= 0 )
0201             return node_pointer(0);
0202
0203         typedef std::pair<point_type, InIt> entry_type;
0204         typedef typename boost::container::allocator_traits<TmpAlloc>::
0205             template rebind_alloc<entry_type> temp_entry_allocator_type;
0206
0207         temp_entry_allocator_type temp_entry_allocator(temp_allocator);
0208         boost::container::vector<entry_type, temp_entry_allocator_type> entries(temp_entry_allocator);
0209
0210         values_count = static_cast<size_type>(diff);
0211         entries.reserve(values_count);
0212
0213         auto const& strategy = index::detail::get_strategy(parameters);
0214
0215         expandable_box<box_type, strategy_type> hint_box(strategy);
0216         for ( ; first != last ; ++first )
0217         {
0218             // NOTE: support for iterators not returning true references adapted
0219             // to Geometry concept and default translator returning true reference
0220             // An alternative would be to dereference the iterator and translate
0221             // in one expression each time the indexable was needed.
0222             typename std::iterator_traits<InIt>::reference in_ref = *first;
0223             typename translator_type::result_type indexable = translator(in_ref);
0224
0225             // NOTE: added for consistency with insert()
0226             // CONSIDER: alternative - ignore invalid indexable or throw an exception
0227             BOOST_GEOMETRY_INDEX_ASSERT(detail::is_valid(indexable), "Indexable is invalid");
0228
0229             hint_box.expand(indexable);
0230
0231             point_type pt;
0232             geometry::centroid(indexable, pt, strategy);
0233             entries.push_back(std::make_pair(pt, first));
0234         }
0235
0236         subtree_elements_counts subtree_counts = calculate_subtree_elements_counts(values_count, parameters, leafs_level);
0237         internal_element el = per_level(entries.begin(), entries.end(), hint_box.get(), values_count, subtree_counts,
0238                                         parameters, translator, allocators);
0239
0240         return el.second;
0241     }
0242
0243 private:
0244     template <typename BoxType, typename Strategy>
0245     class expandable_box
0246     {
0247     public:
0248         explicit expandable_box(Strategy const& strategy)
0249             : m_strategy(strategy), m_initialized(false)
0250         {}
0251
0252         template <typename Indexable>
0253         explicit expandable_box(Indexable const& indexable, Strategy const& strategy)
0254             : m_strategy(strategy), m_initialized(true)
0255         {
0256             detail::bounds(indexable, m_box, m_strategy);
0257         }
0258
0259         template <typename Indexable>
0260         void expand(Indexable const& indexable)
0261         {
0262             if ( !m_initialized )
0263             {
0264                 // it's guaranteed that the Box will be initialized
0265                 // only for Points, Boxes and Segments but that's ok
0266                 // since only those Geometries can be stored
0267                 detail::bounds(indexable, m_box, m_strategy);
0268                 m_initialized = true;
0269             }
0270             else
0271             {
0272                 detail::expand(m_box, indexable, m_strategy);
0273             }
0274         }
0275
0276         void expand_by_epsilon()
0277         {
0278             geometry::detail::expand_by_epsilon(m_box);
0279         }
0280
0281         BoxType const& get() const
0282         {
0283             BOOST_GEOMETRY_INDEX_ASSERT(m_initialized, "uninitialized envelope accessed");
0284             return m_box;
0285         }
0286
0287     private:
0288         BoxType m_box;
0289         Strategy m_strategy;
0290         bool m_initialized;
0291     };
0292
0293     struct subtree_elements_counts
0294     {
0295         subtree_elements_counts(size_type ma, size_type mi) : maxc(ma), minc(mi) {}
0296         size_type maxc;
0297         size_type minc;
0298     };
0299
0300     template <typename EIt> inline static
0301     internal_element per_level(EIt first, EIt last,
0302                                box_type const& hint_box,
0303                                size_type values_count,
0304                                subtree_elements_counts const& subtree_counts,
0305                                parameters_type const& parameters,
0306                                translator_type const& translator,
0307                                allocators_type & allocators)
0308     {
0309         BOOST_GEOMETRY_INDEX_ASSERT(0 < std::distance(first, last) && static_cast<size_type>(std::distance(first, last)) == values_count,
0310                                     "unexpected parameters");
0311
0312         if ( subtree_counts.maxc <= 1 )
0313         {
0314             // ROOT or LEAF
0315             BOOST_GEOMETRY_INDEX_ASSERT(values_count <= parameters.get_max_elements(),
0316                                         "too big number of elements");
0317             // if !root check m_parameters.get_min_elements() <= count
0318
0319             // create new leaf node
0320             node_pointer n = rtree::create_node<allocators_type, leaf>::apply(allocators);                       // MAY THROW (A)
0321             subtree_destroyer auto_remover(n, allocators);
0322             leaf & l = rtree::get<leaf>(*n);
0323
0324             // reserve space for values
0325             rtree::elements(l).reserve(values_count);                                                       // MAY THROW (A)
0326
0327             // calculate values box and copy values
0328             //   initialize the box explicitly to avoid GCC-4.4 uninitialized variable warnings with O2
0329             expandable_box<box_type, strategy_type> elements_box(translator(*(first->second)),
0330                                                                  detail::get_strategy(parameters));
0331             rtree::elements(l).push_back(*(first->second));                                                 // MAY THROW (A?,C)
0332             for ( ++first ; first != last ; ++first )
0333             {
0334                 // NOTE: push_back() must be called at the end in order to support move_iterator.
0335                 //       The iterator is dereferenced 2x (no temporary reference) to support
0336                 //       non-true reference types and move_iterator without std::forward<>.
0337                 elements_box.expand(translator(*(first->second)));
0338                 rtree::elements(l).push_back(*(first->second));                                             // MAY THROW (A?,C)
0339             }
0340
0341 #ifdef BOOST_GEOMETRY_INDEX_EXPERIMENTAL_ENLARGE_BY_EPSILON
0342             // Enlarge bounds of a leaf node.
0343             // It's because Points and Segments are compared WRT machine epsilon
0344             // This ensures that leafs bounds correspond to the stored elements
0345             // NOTE: this is done only if the Indexable is a different kind of Geometry
0346             //   than the bounds (only Box for now). Spatial predicates are checked
0347             //   the same way for Geometry of the same kind.
0348             if BOOST_GEOMETRY_CONSTEXPR (! index::detail::is_bounding_geometry
0349                                             <
0350                                                 typename indexable_type<translator_type>::type
0351                                             >::value)
0352             {
0353                 elements_box.expand_by_epsilon();
0354             }
0355 #endif
0356
0357             auto_remover.release();
0358             return internal_element(elements_box.get(), n);
0359         }
0360
0361         // calculate next max and min subtree counts
0362         subtree_elements_counts next_subtree_counts = subtree_counts;
0363         next_subtree_counts.maxc /= parameters.get_max_elements();
0364         next_subtree_counts.minc /= parameters.get_max_elements();
0365
0366         // create new internal node
0367         node_pointer n = rtree::create_node<allocators_type, internal_node>::apply(allocators);                  // MAY THROW (A)
0368         subtree_destroyer auto_remover(n, allocators);
0369         internal_node & in = rtree::get<internal_node>(*n);
0370
0371         // reserve space for values
0372         size_type nodes_count = calculate_nodes_count(values_count, subtree_counts);
0373         rtree::elements(in).reserve(nodes_count);                                                           // MAY THROW (A)
0374         // calculate values box and copy values
0375         expandable_box<box_type, strategy_type> elements_box(detail::get_strategy(parameters));
0376
0377         per_level_packets(first, last, hint_box, values_count, subtree_counts, next_subtree_counts,
0378                           rtree::elements(in), elements_box,
0379                           parameters, translator, allocators);
0380
0381         auto_remover.release();
0382         return internal_element(elements_box.get(), n);
0383     }
0384
0385     template <typename EIt, typename ExpandableBox> inline static
0386     void per_level_packets(EIt first, EIt last,
0387                            box_type const& hint_box,
0388                            size_type values_count,
0389                            subtree_elements_counts const& subtree_counts,
0390                            subtree_elements_counts const& next_subtree_counts,
0391                            internal_elements & elements,
0392                            ExpandableBox & elements_box,
0393                            parameters_type const& parameters,
0394                            translator_type const& translator,
0395                            allocators_type & allocators)
0396     {
0397         BOOST_GEOMETRY_INDEX_ASSERT(0 < std::distance(first, last) && static_cast<size_type>(std::distance(first, last)) == values_count,
0398                                     "unexpected parameters");
0399
0400         BOOST_GEOMETRY_INDEX_ASSERT(subtree_counts.minc <= values_count,
0401                                     "too small number of elements");
0402
0403         // only one packet
0404         if ( values_count <= subtree_counts.maxc )
0405         {
0406             // the end, move to the next level
0407             internal_element el = per_level(first, last, hint_box, values_count, next_subtree_counts,
0408                                             parameters, translator, allocators);
0409
0410             // in case if push_back() do throw here
0411             // and even if this is not probable (previously reserved memory, nonthrowing pairs copy)
0412             // this case is also tested by exceptions test.
0413             subtree_destroyer auto_remover(el.second, allocators);
0414             // this container should have memory allocated, reserve() called outside
0415             elements.push_back(el);                                                 // MAY THROW (A?,C) - however in normal conditions shouldn't
0416             auto_remover.release();
0417
0418             elements_box.expand(el.first);
0419             return;
0420         }
0421
0422         size_type median_count = calculate_median_count(values_count, subtree_counts);
0423         EIt median = first + median_count;
0424
0425         coordinate_type greatest_length;
0426         std::size_t greatest_dim_index = 0;
0427         pack_utils::biggest_edge<dimension>::apply(hint_box, greatest_length, greatest_dim_index);
0428         box_type left, right;
0429         pack_utils::nth_element_and_half_boxes<0, dimension>
0430             ::apply(first, median, last, hint_box, left, right, greatest_dim_index);
0431
0432         per_level_packets(first, median, left,
0433                           median_count, subtree_counts, next_subtree_counts,
0434                           elements, elements_box,
0435                           parameters, translator, allocators);
0436         per_level_packets(median, last, right,
0437                           values_count - median_count, subtree_counts, next_subtree_counts,
0438                           elements, elements_box,
0439                           parameters, translator, allocators);
0440     }
0441
0442     inline static
0443     subtree_elements_counts calculate_subtree_elements_counts(size_type elements_count, parameters_type const& parameters, size_type & leafs_level)
0444     {
0445         boost::ignore_unused(parameters);
0446
0447         subtree_elements_counts res(1, 1);
0448         leafs_level = 0;
0449
0450         size_type smax = parameters.get_max_elements();
0451         for ( ; smax < elements_count ; smax *= parameters.get_max_elements(), ++leafs_level )
0452             res.maxc = smax;
0453
0454         res.minc = parameters.get_min_elements() * (res.maxc / parameters.get_max_elements());
0455
0456         return res;
0457     }
0458
0459     inline static
0460     size_type calculate_nodes_count(size_type count,
0461                                     subtree_elements_counts const& subtree_counts)
0462     {
0463         size_type n = count / subtree_counts.maxc;
0464         size_type r = count % subtree_counts.maxc;
0465
0466         if ( 0 < r && r < subtree_counts.minc )
0467         {
0468             size_type count_minus_min = count - subtree_counts.minc;
0469             n = count_minus_min / subtree_counts.maxc;
0470             r = count_minus_min % subtree_counts.maxc;
0471             ++n;
0472         }
0473
0474         if ( 0 < r )
0475             ++n;
0476
0477         return n;
0478     }
0479
0480     inline static
0481     size_type calculate_median_count(size_type count,
0482                                      subtree_elements_counts const& subtree_counts)
0483     {
0484         // e.g. for max = 5, min = 2, count = 52, subtree_max = 25, subtree_min = 10
0485
0486         size_type n = count / subtree_counts.maxc; // e.g. 52 / 25 = 2
0487         size_type r = count % subtree_counts.maxc; // e.g. 52 % 25 = 2
0488         size_type median_count = (n / 2) * subtree_counts.maxc; // e.g. 2 / 2 * 25 = 25
0489
0490         if ( 0 != r ) // e.g. 0 != 2
0491         {
0492             if ( subtree_counts.minc <= r ) // e.g. 10 <= 2 == false
0493             {
0494                 //BOOST_GEOMETRY_INDEX_ASSERT(0 < n, "unexpected value");
0495                 median_count = ((n+1)/2) * subtree_counts.maxc; // if calculated ((2+1)/2) * 25 which would be ok, but not in all cases
0496             }
0497             else // r < subtree_counts.second  // e.g. 2 < 10 == true
0498             {
0499                 size_type count_minus_min = count - subtree_counts.minc; // e.g. 52 - 10 = 42
0500                 n = count_minus_min / subtree_counts.maxc; // e.g. 42 / 25 = 1
0501                 r = count_minus_min % subtree_counts.maxc; // e.g. 42 % 25 = 17
0502                 if ( r == 0 )                               // e.g. false
0503                 {
0504                     // n can't be equal to 0 because then there wouldn't be any element in the other node
0505                     //BOOST_GEOMETRY_INDEX_ASSERT(0 < n, "unexpected value");
0506                     median_count = ((n+1)/2) * subtree_counts.maxc;     // if calculated ((1+1)/2) * 25 which would be ok, but not in all cases
0507                 }
0508                 else
0509                 {
0510                     if ( n == 0 )                                        // e.g. false
0511                         median_count = r;                                // if calculated -> 17 which is wrong!
0512                     else
0513                         median_count = ((n+2)/2) * subtree_counts.maxc; // e.g. ((1+2)/2) * 25 = 25
0514                 }
0515             }
0516         }
0517
0518         return median_count;
0519     }
0520 };
0521
0522 }}}}} // namespace boost::geometry::index::detail::rtree
0523
0524 #endif // BOOST_GEOMETRY_INDEX_DETAIL_RTREE_PACK_CREATE_HPP