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> applyboost::geometry::index::detail::rtree::pack_utils::biggest_edge29 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> applyboost::geometry::index::detail::rtree::pack_utils::biggest_edge46 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> operator ()boost::geometry::index::detail::rtree::pack_utils::point_entries_comparer57 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> applyboost::geometry::index::detail::rtree::pack_utils::nth_element_and_half_boxes67 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> applyboost::geometry::index::detail::rtree::pack_utils::nth_element_and_half_boxes91 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 apply(InIt first,InIt last,size_type & values_count,size_type & leafs_level,parameters_type const & parameters,Translator const & translator,Allocators & allocators)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: expandable_box()199 expandable_box() 200 : m_initialized(false) 201 {} 202 203 template <typename Indexable> expandable_box(Indexable const & 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> expand(Indexable const & 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 expand_by_epsilon()227 void expand_by_epsilon() 228 { 229 geometry::detail::expand_by_epsilon(m_box); 230 } 231 get() const232 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 { subtree_elements_countsboost::geometry::index::detail::rtree::pack::subtree_elements_counts245 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 per_level(EIt first,EIt last,Box const & hint_box,std::size_t values_count,subtree_elements_counts const & subtree_counts,parameters_type const & parameters,Translator const & translator,Allocators & allocators)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 per_level_packets(EIt first,EIt last,Box const & hint_box,std::size_t values_count,subtree_elements_counts const & subtree_counts,subtree_elements_counts const & next_subtree_counts,internal_elements & elements,ExpandableBox & elements_box,parameters_type const & parameters,Translator const & translator,Allocators & allocators)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 calculate_subtree_elements_counts(std::size_t elements_count,parameters_type const & parameters,size_type & leafs_level)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 calculate_nodes_count(std::size_t count,subtree_elements_counts const & subtree_counts)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 calculate_median_count(std::size_t count,subtree_elements_counts const & subtree_counts)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 466