pack_create.hpp 20 KB

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