pack_create.hpp 20 KB

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