pack_create.hpp 20 KB

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