boost/geometry/index/detail/rtree/pack_create.hpp

   1 // Boost.Geometry Index
   2 //
   3 // R-tree initial packing
   4 //
   5 // Copyright (c) 2011-2017 Adam Wulkiewicz, Lodz, Poland.
   6 //
   7 // Use, modification and distribution is subject to the Boost Software License,
   8 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
   9 // http://www.boost.org/LICENSE_1_0.txt)
  10
  11 #ifndef BOOST_GEOMETRY_INDEX_DETAIL_RTREE_PACK_CREATE_HPP
  12 #define BOOST_GEOMETRY_INDEX_DETAIL_RTREE_PACK_CREATE_HPP
  13
  14 #include <boost/geometry/algorithms/expand.hpp>
  15 #include <boost/geometry/index/detail/algorithms/bounds.hpp>
  16 #include <boost/geometry/index/detail/algorithms/nth_element.hpp>
  17
  18 #include <boost/geometry/algorithms/detail/expand_by_epsilon.hpp>
  19
  20 namespace boost { namespace geometry { namespace index { namespace detail { namespace rtree {
  21
  22 namespace pack_utils {
  23
  24 template <std::size_t Dimension>
  25 struct biggest_edge
  26 {
  27     BOOST_STATIC_ASSERT(0 < Dimension);
  28     template <typename Box>
  29     static inline void apply(Box const& box, typename coordinate_type<Box>::type & length, std::size_t & dim_index)
  30     {
  31         biggest_edge<Dimension-1>::apply(box, length, dim_index);
  32         typename coordinate_type<Box>::type curr
  33             = geometry::get<max_corner, Dimension-1>(box) - geometry::get<min_corner, Dimension-1>(box);
  34         if ( length < curr )
  35         {
  36             dim_index = Dimension - 1;
  37             length = curr;
  38         }
  39     }
  40 };
  41
  42 template <>
  43 struct biggest_edge<1>
  44 {
  45     template <typename Box>
  46     static inline void apply(Box const& box, typename coordinate_type<Box>::type & length, std::size_t & dim_index)
  47     {
  48         dim_index = 0;
  49         length = geometry::get<max_corner, 0>(box) - geometry::get<min_corner, 0>(box);
  50     }
  51 };
  52
  53 template <std::size_t I>
  54 struct point_entries_comparer
  55 {
  56     template <typename PointEntry>
  57     bool operator()(PointEntry const& e1, PointEntry const& e2) const
  58     {
  59         return geometry::get<I>(e1.first) < geometry::get<I>(e2.first);
  60     }
  61 };
  62
  63 template <std::size_t I, std::size_t Dimension>
  64 struct nth_element_and_half_boxes
  65 {
  66     template <typename EIt, typename Box>
  67     static inline void apply(EIt first, EIt median, EIt last, Box const& box, Box & left, Box & right, std::size_t dim_index)
  68     {
  69         if ( I == dim_index )
  70         {
  71             index::detail::nth_element(first, median, last, point_entries_comparer<I>());
  72
  73             geometry::convert(box, left);
  74             geometry::convert(box, right);
  75             typename coordinate_type<Box>::type edge_len
  76                 = geometry::get<max_corner, I>(box) - geometry::get<min_corner, I>(box);
  77             typename coordinate_type<Box>::type median
  78                 = geometry::get<min_corner, I>(box) + edge_len / 2;
  79             geometry::set<max_corner, I>(left, median);
  80             geometry::set<min_corner, I>(right, median);
  81         }
  82         else
  83             nth_element_and_half_boxes<I+1, Dimension>::apply(first, median, last, box, left, right, dim_index);
  84     }
  85 };
  86
  87 template <std::size_t Dimension>
  88 struct nth_element_and_half_boxes<Dimension, Dimension>
  89 {
  90     template <typename EIt, typename Box>
  91     static inline void apply(EIt , EIt , EIt , Box const& , Box & , Box & , std::size_t ) {}
  92 };
  93
  94 } // namespace pack_utils
  95
  96 // STR leafs number are calculated as rcount/max
  97 // and the number of splitting planes for each dimension as (count/max)^(1/dimension)
  98 // <-> for dimension==2 -> sqrt(count/max)
  99 //
 100 // The main flaw of this algorithm is that the resulting tree will have bad structure for:
 101 // 1. non-uniformly distributed elements
 102 //      Statistic check could be performed, e.g. based on variance of lengths of elements edges for each dimension
 103 // 2. elements distributed mainly along one axis
 104 //      Calculate bounding box of all elements and then number of dividing planes for a dimension
 105 //      from the length of BB edge for this dimension (more or less assuming that elements are uniformly-distributed squares)
 106 //
 107 // Another thing is that the last node may have less elements than Max or even Min.
 108 // The number of splitting planes must be chosen more carefully than count/max
 109 //
 110 // This algorithm is something between STR and TGS
 111 // it is more similar to the top-down recursive kd-tree creation algorithm
 112 // using object median split and split axis of greatest BB edge
 113 // BB is only used as a hint (assuming objects are distributed uniformly)
 114 //
 115 // Implemented algorithm guarantees that the number of elements in nodes will be between Min and Max
 116 // and that nodes are packed as tightly as possible
 117 // e.g. for 177 values Max = 5 and Min = 2 it will construct the following tree:
 118 // ROOT                 177
 119 // L1          125               52
 120 // L2  25  25  25  25  25   25  17    10
 121 // L3  5x5 5x5 5x5 5x5 5x5  5x5 3x5+2 2x5
 122
 123 template <typename Value, typename Options, typename Translator, typename Box, typename Allocators>
 124 class pack
 125 {
 126     typedef typename rtree::node<Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag>::type node;
 127     typedef typename rtree::internal_node<Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag>::type internal_node;
 128     typedef typename rtree::leaf<Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag>::type leaf;
 129
 130     typedef typename Allocators::node_pointer node_pointer;
 131     typedef rtree::subtree_destroyer<Value, Options, Translator, Box, Allocators> subtree_destroyer;
 132     typedef typename Allocators::size_type size_type;
 133
 134     typedef typename geometry::point_type<Box>::type point_type;
 135     typedef typename geometry::coordinate_type<point_type>::type coordinate_type;
 136     typedef typename detail::default_content_result<Box>::type content_type;
 137     typedef typename Options::parameters_type parameters_type;
 138     static const std::size_t dimension = geometry::dimension<point_type>::value;
 139
 140     typedef typename rtree::container_from_elements_type<
 141         typename rtree::elements_type<leaf>::type,
 142         std::size_t
 143     >::type values_counts_container;
 144
 145     typedef typename rtree::elements_type<internal_node>::type internal_elements;
 146     typedef typename internal_elements::value_type internal_element;
 147
 148 public:
 149     // Arbitrary iterators
 150     template <typename InIt> inline static
 151     node_pointer apply(InIt first, InIt last, size_type & values_count, size_type & leafs_level,
 152                        parameters_type const& parameters, Translator const& translator, Allocators & allocators)
 153     {
 154         typedef typename std::iterator_traits<InIt>::difference_type diff_type;
 155
 156         diff_type diff = std::distance(first, last);
 157         if ( diff <= 0 )
 158             return node_pointer(0);
 159
 160         typedef std::pair<point_type, InIt> entry_type;
 161         std::vector<entry_type> entries;
 162
 163         values_count = static_cast<size_type>(diff);
 164         entries.reserve(values_count);
 165
 166         expandable_box<Box> hint_box;
 167         for ( ; first != last ; ++first )
 168         {
 169             // NOTE: support for iterators not returning true references adapted
 170             // to Geometry concept and default translator returning true reference
 171             // An alternative would be to dereference the iterator and translate
 172             // in one expression each time the indexable was needed.
 173             typename std::iterator_traits<InIt>::reference in_ref = *first;
 174             typename Translator::result_type indexable = translator(in_ref);
 175
 176             // NOTE: added for consistency with insert()
 177             // CONSIDER: alternative - ignore invalid indexable or throw an exception
 178             BOOST_GEOMETRY_INDEX_ASSERT(detail::is_valid(indexable), "Indexable is invalid");
 179
 180             hint_box.expand(indexable);
 181
 182             point_type pt;
 183             geometry::centroid(indexable, pt);
 184             entries.push_back(std::make_pair(pt, first));
 185         }
 186
 187         subtree_elements_counts subtree_counts = calculate_subtree_elements_counts(values_count, parameters, leafs_level);
 188         internal_element el = per_level(entries.begin(), entries.end(), hint_box.get(), values_count, subtree_counts,
 189                                         parameters, translator, allocators);
 190
 191         return el.second;
 192     }
 193
 194 private:
 195     template <typename BoxType>
 196     class expandable_box
 197     {
 198     public:
 199         expandable_box()
 200             : m_initialized(false)
 201         {}
 202
 203         template <typename Indexable>
 204         explicit expandable_box(Indexable const& indexable)
 205             : m_initialized(true)
 206         {
 207             detail::bounds(indexable, m_box);
 208         }
 209
 210         template <typename Indexable>
 211         void expand(Indexable const& indexable)
 212         {
 213             if ( !m_initialized )
 214             {
 215                 // it's guaranteed that the Box will be initialized
 216                 // only for Points, Boxes and Segments but that's ok
 217                 // since only those Geometries can be stored
 218                 detail::bounds(indexable, m_box);
 219                 m_initialized = true;
 220             }
 221             else
 222             {
 223                 geometry::expand(m_box, indexable);
 224             }
 225         }
 226
 227         void expand_by_epsilon()
 228         {
 229             geometry::detail::expand_by_epsilon(m_box);
 230         }
 231
 232         BoxType const& get() const
 233         {
 234             BOOST_GEOMETRY_INDEX_ASSERT(m_initialized, "uninitialized envelope accessed");
 235             return m_box;
 236         }
 237
 238     private:
 239         bool m_initialized;
 240         BoxType m_box;
 241     };
 242
 243     struct subtree_elements_counts
 244     {
 245         subtree_elements_counts(std::size_t ma, std::size_t mi) : maxc(ma), minc(mi) {}
 246         std::size_t maxc;
 247         std::size_t minc;
 248     };
 249
 250     template <typename EIt> inline static
 251     internal_element per_level(EIt first, EIt last, Box const& hint_box, std::size_t values_count, subtree_elements_counts const& subtree_counts,
 252                                parameters_type const& parameters, Translator const& translator, Allocators & allocators)
 253     {
 254         BOOST_GEOMETRY_INDEX_ASSERT(0 < std::distance(first, last) && static_cast<std::size_t>(std::distance(first, last)) == values_count,
 255                                     "unexpected parameters");
 256
 257         if ( subtree_counts.maxc <= 1 )
 258         {
 259             // ROOT or LEAF
 260             BOOST_GEOMETRY_INDEX_ASSERT(values_count <= parameters.get_max_elements(),
 261                                         "too big number of elements");
 262             // if !root check m_parameters.get_min_elements() <= count
 263
 264             // create new leaf node
 265             node_pointer n = rtree::create_node<Allocators, leaf>::apply(allocators);                       // MAY THROW (A)
 266             subtree_destroyer auto_remover(n, allocators);
 267             leaf & l = rtree::get<leaf>(*n);
 268
 269             // reserve space for values
 270             rtree::elements(l).reserve(values_count);                                                       // MAY THROW (A)
 271
 272             // calculate values box and copy values
 273             //   initialize the box explicitly to avoid GCC-4.4 uninitialized variable warnings with O2
 274             expandable_box<Box> elements_box(translator(*(first->second)));
 275             rtree::elements(l).push_back(*(first->second));                                                 // MAY THROW (A?,C)
 276             for ( ++first ; first != last ; ++first )
 277             {
 278                 // NOTE: push_back() must be called at the end in order to support move_iterator.
 279                 //       The iterator is dereferenced 2x (no temporary reference) to support
 280                 //       non-true reference types and move_iterator without boost::forward<>.
 281                 elements_box.expand(translator(*(first->second)));
 282                 rtree::elements(l).push_back(*(first->second));                                             // MAY THROW (A?,C)
 283             }
 284
 285 #ifdef BOOST_GEOMETRY_INDEX_EXPERIMENTAL_ENLARGE_BY_EPSILON
 286             // Enlarge bounds of a leaf node.
 287             // It's because Points and Segments are compared WRT machine epsilon
 288             // This ensures that leafs bounds correspond to the stored elements
 289             // NOTE: this is done only if the Indexable is a different kind of Geometry
 290             //   than the bounds (only Box for now). Spatial predicates are checked
 291             //   the same way for Geometry of the same kind.
 292             if ( BOOST_GEOMETRY_CONDITION((
 293                     ! index::detail::is_bounding_geometry
 294                         <
 295                             typename indexable_type<Translator>::type
 296                         >::value )) )
 297             {
 298                 elements_box.expand_by_epsilon();
 299             }
 300 #endif
 301
 302             auto_remover.release();
 303             return internal_element(elements_box.get(), n);
 304         }
 305
 306         // calculate next max and min subtree counts
 307         subtree_elements_counts next_subtree_counts = subtree_counts;
 308         next_subtree_counts.maxc /= parameters.get_max_elements();
 309         next_subtree_counts.minc /= parameters.get_max_elements();
 310
 311         // create new internal node
 312         node_pointer n = rtree::create_node<Allocators, internal_node>::apply(allocators);                  // MAY THROW (A)
 313         subtree_destroyer auto_remover(n, allocators);
 314         internal_node & in = rtree::get<internal_node>(*n);
 315
 316         // reserve space for values
 317         std::size_t nodes_count = calculate_nodes_count(values_count, subtree_counts);
 318         rtree::elements(in).reserve(nodes_count);                                                           // MAY THROW (A)
 319         // calculate values box and copy values
 320         expandable_box<Box> elements_box;
 321
 322         per_level_packets(first, last, hint_box, values_count, subtree_counts, next_subtree_counts,
 323                           rtree::elements(in), elements_box,
 324                           parameters, translator, allocators);
 325
 326         auto_remover.release();
 327         return internal_element(elements_box.get(), n);
 328     }
 329
 330     template <typename EIt, typename ExpandableBox> inline static
 331     void per_level_packets(EIt first, EIt last, Box const& hint_box,
 332                            std::size_t values_count,
 333                            subtree_elements_counts const& subtree_counts,
 334                            subtree_elements_counts const& next_subtree_counts,
 335                            internal_elements & elements, ExpandableBox & elements_box,
 336                            parameters_type const& parameters, Translator const& translator, Allocators & allocators)
 337     {
 338         BOOST_GEOMETRY_INDEX_ASSERT(0 < std::distance(first, last) && static_cast<std::size_t>(std::distance(first, last)) == values_count,
 339                                     "unexpected parameters");
 340
 341         BOOST_GEOMETRY_INDEX_ASSERT(subtree_counts.minc <= values_count,
 342                                     "too small number of elements");
 343
 344         // only one packet
 345         if ( values_count <= subtree_counts.maxc )
 346         {
 347             // the end, move to the next level
 348             internal_element el = per_level(first, last, hint_box, values_count, next_subtree_counts,
 349                                             parameters, translator, allocators);
 350
 351             // in case if push_back() do throw here
 352             // and even if this is not probable (previously reserved memory, nonthrowing pairs copy)
 353             // this case is also tested by exceptions test.
 354             subtree_destroyer auto_remover(el.second, allocators);
 355             // this container should have memory allocated, reserve() called outside
 356             elements.push_back(el);                                                 // MAY THROW (A?,C) - however in normal conditions shouldn't
 357             auto_remover.release();
 358
 359             elements_box.expand(el.first);
 360             return;
 361         }
 362
 363         std::size_t median_count = calculate_median_count(values_count, subtree_counts);
 364         EIt median = first + median_count;
 365
 366         coordinate_type greatest_length;
 367         std::size_t greatest_dim_index = 0;
 368         pack_utils::biggest_edge<dimension>::apply(hint_box, greatest_length, greatest_dim_index);
 369         Box left, right;
 370         pack_utils::nth_element_and_half_boxes<0, dimension>
 371             ::apply(first, median, last, hint_box, left, right, greatest_dim_index);
 372
 373         per_level_packets(first, median, left,
 374                           median_count, subtree_counts, next_subtree_counts,
 375                           elements, elements_box,
 376                           parameters, translator, allocators);
 377         per_level_packets(median, last, right,
 378                           values_count - median_count, subtree_counts, next_subtree_counts,
 379                           elements, elements_box,
 380                           parameters, translator, allocators);
 381     }
 382
 383     inline static
 384     subtree_elements_counts calculate_subtree_elements_counts(std::size_t elements_count, parameters_type const& parameters, size_type & leafs_level)
 385     {
 386         boost::ignore_unused_variable_warning(parameters);
 387
 388         subtree_elements_counts res(1, 1);
 389         leafs_level = 0;
 390
 391         std::size_t smax = parameters.get_max_elements();
 392         for ( ; smax < elements_count ; smax *= parameters.get_max_elements(), ++leafs_level )
 393             res.maxc = smax;
 394
 395         res.minc = parameters.get_min_elements() * (res.maxc / parameters.get_max_elements());
 396
 397         return res;
 398     }
 399
 400     inline static
 401     std::size_t calculate_nodes_count(std::size_t count,
 402                                       subtree_elements_counts const& subtree_counts)
 403     {
 404         std::size_t n = count / subtree_counts.maxc;
 405         std::size_t r = count % subtree_counts.maxc;
 406
 407         if ( 0 < r && r < subtree_counts.minc )
 408         {
 409             std::size_t count_minus_min = count - subtree_counts.minc;
 410             n = count_minus_min / subtree_counts.maxc;
 411             r = count_minus_min % subtree_counts.maxc;
 412             ++n;
 413         }
 414
 415         if ( 0 < r )
 416             ++n;
 417
 418         return n;
 419     }
 420
 421     inline static
 422     std::size_t calculate_median_count(std::size_t count,
 423                                        subtree_elements_counts const& subtree_counts)
 424     {
 425         // e.g. for max = 5, min = 2, count = 52, subtree_max = 25, subtree_min = 10
 426
 427         std::size_t n = count / subtree_counts.maxc; // e.g. 52 / 25 = 2
 428         std::size_t r = count % subtree_counts.maxc; // e.g. 52 % 25 = 2
 429         std::size_t median_count = (n / 2) * subtree_counts.maxc; // e.g. 2 / 2 * 25 = 25
 430
 431         if ( 0 != r ) // e.g. 0 != 2
 432         {
 433             if ( subtree_counts.minc <= r ) // e.g. 10 <= 2 == false
 434             {
 435                 //BOOST_GEOMETRY_INDEX_ASSERT(0 < n, "unexpected value");
 436                 median_count = ((n+1)/2) * subtree_counts.maxc; // if calculated ((2+1)/2) * 25 which would be ok, but not in all cases
 437             }
 438             else // r < subtree_counts.second  // e.g. 2 < 10 == true
 439             {
 440                 std::size_t count_minus_min = count - subtree_counts.minc; // e.g. 52 - 10 = 42
 441                 n = count_minus_min / subtree_counts.maxc; // e.g. 42 / 25 = 1
 442                 r = count_minus_min % subtree_counts.maxc; // e.g. 42 % 25 = 17
 443                 if ( r == 0 )                               // e.g. false
 444                 {
 445                     // n can't be equal to 0 because then there wouldn't be any element in the other node
 446                     //BOOST_GEOMETRY_INDEX_ASSERT(0 < n, "unexpected value");
 447                     median_count = ((n+1)/2) * subtree_counts.maxc;     // if calculated ((1+1)/2) * 25 which would be ok, but not in all cases
 448                 }
 449                 else
 450                 {
 451                     if ( n == 0 )                                        // e.g. false
 452                         median_count = r;                                // if calculated -> 17 which is wrong!
 453                     else
 454                         median_count = ((n+2)/2) * subtree_counts.maxc; // e.g. ((1+2)/2) * 25 = 25
 455                 }
 456             }
 457         }
 458
 459         return median_count;
 460     }
 461 };
 462
 463 }}}}} // namespace boost::geometry::index::detail::rtree
 464
 465 #endif // BOOST_GEOMETRY_INDEX_DETAIL_RTREE_PACK_CREATE_HPP