pack_create.hpp 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376
  1. // Boost.Geometry Index
  2. //
  3. // R-tree initial packing
  4. //
  5. // Copyright (c) 2011-2013 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. namespace boost { namespace geometry { namespace index { namespace detail { namespace rtree {
  13. namespace pack_utils {
  14. template <std::size_t Dimension>
  15. struct biggest_edge
  16. {
  17. BOOST_STATIC_ASSERT(0 < Dimension);
  18. template <typename Box>
  19. static inline void apply(Box const& box, typename coordinate_type<Box>::type & length, std::size_t & dim_index)
  20. {
  21. biggest_edge<Dimension-1>::apply(box, length, dim_index);
  22. typename coordinate_type<Box>::type curr
  23. = geometry::get<max_corner, Dimension-1>(box) - geometry::get<min_corner, Dimension-1>(box);
  24. if ( length < curr )
  25. {
  26. dim_index = Dimension - 1;
  27. length = curr;
  28. }
  29. }
  30. };
  31. template <>
  32. struct biggest_edge<1>
  33. {
  34. template <typename Box>
  35. static inline void apply(Box const& box, typename coordinate_type<Box>::type & length, std::size_t & dim_index)
  36. {
  37. dim_index = 0;
  38. length = geometry::get<max_corner, 0>(box) - geometry::get<min_corner, 0>(box);
  39. }
  40. };
  41. template <std::size_t I>
  42. struct point_entries_comparer
  43. {
  44. template <typename PointEntry>
  45. bool operator()(PointEntry const& e1, PointEntry const& e2) const
  46. {
  47. return geometry::get<I>(e1.first) < geometry::get<I>(e2.first);
  48. }
  49. };
  50. template <std::size_t I, std::size_t Dimension>
  51. struct partial_sort_and_half_boxes
  52. {
  53. template <typename EIt, typename Box>
  54. static inline void apply(EIt first, EIt median, EIt last, Box const& box, Box & left, Box & right, std::size_t dim_index)
  55. {
  56. if ( I == dim_index )
  57. {
  58. std::partial_sort(first, median, last, point_entries_comparer<I>());
  59. geometry::convert(box, left);
  60. geometry::convert(box, right);
  61. typename coordinate_type<Box>::type edge_len
  62. = geometry::get<max_corner, I>(box) - geometry::get<min_corner, I>(box);
  63. typename coordinate_type<Box>::type median
  64. = geometry::get<min_corner, I>(box) + edge_len / 2;
  65. geometry::set<max_corner, I>(left, median);
  66. geometry::set<min_corner, I>(right, median);
  67. }
  68. else
  69. partial_sort_and_half_boxes<I+1, Dimension>::apply(first, median, last, box, left, right, dim_index);
  70. }
  71. };
  72. template <std::size_t Dimension>
  73. struct partial_sort_and_half_boxes<Dimension, Dimension>
  74. {
  75. template <typename EIt, typename Box>
  76. static inline void apply(EIt , EIt , EIt , Box const& , Box & , Box & , std::size_t ) {}
  77. };
  78. } // namespace pack_utils
  79. // STR leafs number are calculated as rcount/max
  80. // and the number of splitting planes for each dimension as (count/max)^(1/dimension)
  81. // <-> for dimension==2 -> sqrt(count/max)
  82. //
  83. // The main flaw of this algorithm is that the resulting tree will have bad structure for:
  84. // 1. non-uniformly distributed elements
  85. // Statistic check could be performed, e.g. based on variance of lengths of elements edges for each dimension
  86. // 2. elements distributed mainly along one axis
  87. // Calculate bounding box of all elements and then number of dividing planes for a dimension
  88. // from the length of BB edge for this dimension (more or less assuming that elements are uniformly-distributed squares)
  89. //
  90. // Another thing is that the last node may have less elements than Max or even Min.
  91. // The number of splitting planes must be chosen more carefully than count/max
  92. //
  93. // This algorithm is something between STR and TGS
  94. // it is more similar to the top-down recursive kd-tree creation algorithm
  95. // using object median split and split axis of greatest BB edge
  96. // BB is only used as a hint (assuming objects are distributed uniformly)
  97. //
  98. // Implemented algorithm guarantees that the number of elements in nodes will be between Min and Max
  99. // and that nodes are packed as tightly as possible
  100. // e.g. for 177 values Max = 5 and Min = 2 it will construct the following tree:
  101. // ROOT 177
  102. // L1 125 52
  103. // L2 25 25 25 25 25 25 17 10
  104. // L3 5x5 5x5 5x5 5x5 5x5 5x5 3x5+2 2x5
  105. template <typename Value, typename Options, typename Translator, typename Box, typename Allocators>
  106. class pack
  107. {
  108. typedef typename rtree::node<Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag>::type node;
  109. typedef typename rtree::internal_node<Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag>::type internal_node;
  110. typedef typename rtree::leaf<Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag>::type leaf;
  111. typedef typename Allocators::node_pointer node_pointer;
  112. typedef rtree::node_auto_ptr<Value, Options, Translator, Box, Allocators> node_auto_ptr;
  113. typedef typename Allocators::size_type size_type;
  114. typedef typename traits::point_type<Box>::type point_type;
  115. typedef typename traits::coordinate_type<point_type>::type coordinate_type;
  116. typedef typename detail::default_content_result<Box>::type content_type;
  117. typedef typename Options::parameters_type parameters_type;
  118. static const std::size_t dimension = traits::dimension<point_type>::value;
  119. typedef typename rtree::container_from_elements_type<
  120. typename rtree::elements_type<leaf>::type,
  121. std::size_t
  122. >::type values_counts_container;
  123. typedef typename rtree::elements_type<internal_node>::type internal_elements;
  124. typedef typename internal_elements::value_type internal_element;
  125. public:
  126. // Arbitrary iterators
  127. template <typename InIt> inline static
  128. node_pointer apply(InIt first, InIt last, size_type & values_count, size_type & leafs_level,
  129. parameters_type const& parameters, Translator const& translator, Allocators & allocators)
  130. {
  131. typedef typename std::iterator_traits<InIt>::difference_type diff_type;
  132. diff_type diff = std::distance(first, last);
  133. if ( diff <= 0 )
  134. return node_pointer(0);
  135. typedef std::pair<point_type, InIt> entry_type;
  136. std::vector<entry_type> entries;
  137. values_count = static_cast<size_type>(diff);
  138. entries.reserve(values_count);
  139. Box hint_box;
  140. geometry::assign_inverse(hint_box);
  141. for ( ; first != last ; ++first )
  142. {
  143. geometry::expand(hint_box, translator(*first));
  144. point_type pt;
  145. geometry::centroid(translator(*first), pt);
  146. entries.push_back(std::make_pair(pt, first));
  147. }
  148. subtree_elements_counts subtree_counts = calculate_subtree_elements_counts(values_count, parameters, leafs_level);
  149. internal_element el = per_level(entries.begin(), entries.end(), hint_box, values_count, subtree_counts,
  150. parameters, translator, allocators);
  151. return el.second;
  152. }
  153. private:
  154. struct subtree_elements_counts
  155. {
  156. subtree_elements_counts(std::size_t ma, std::size_t mi) : maxc(ma), minc(mi) {}
  157. std::size_t maxc;
  158. std::size_t minc;
  159. };
  160. template <typename EIt> inline static
  161. internal_element per_level(EIt first, EIt last, Box const& hint_box, std::size_t values_count, subtree_elements_counts const& subtree_counts,
  162. parameters_type const& parameters, Translator const& translator, Allocators & allocators)
  163. {
  164. BOOST_ASSERT(0 < std::distance(first, last) && static_cast<std::size_t>(std::distance(first, last)) == values_count);
  165. if ( subtree_counts.maxc <= 1 )
  166. {
  167. // ROOT or LEAF
  168. BOOST_ASSERT(values_count <= parameters.get_max_elements());
  169. // if !root check m_parameters.get_min_elements() <= count
  170. // create new leaf node
  171. node_pointer n = rtree::create_node<Allocators, leaf>::apply(allocators); // MAY THROW (A)
  172. node_auto_ptr auto_remover(n, allocators);
  173. leaf & l = rtree::get<leaf>(*n);
  174. // reserve space for values
  175. rtree::elements(l).reserve(values_count); // MAY THROW (A)
  176. // calculate values box and copy values
  177. Box elements_box;
  178. geometry::assign_inverse(elements_box);
  179. for ( ; first != last ; ++first )
  180. {
  181. rtree::elements(l).push_back(*(first->second)); // MAY THROW (A?,C)
  182. geometry::expand(elements_box, translator(*(first->second)));
  183. }
  184. auto_remover.release();
  185. return internal_element(elements_box, n);
  186. }
  187. // calculate next max and min subtree counts
  188. subtree_elements_counts next_subtree_counts = subtree_counts;
  189. next_subtree_counts.maxc /= parameters.get_max_elements();
  190. next_subtree_counts.minc /= parameters.get_max_elements();
  191. // create new internal node
  192. node_pointer n = rtree::create_node<Allocators, internal_node>::apply(allocators); // MAY THROW (A)
  193. node_auto_ptr auto_remover(n, allocators);
  194. internal_node & in = rtree::get<internal_node>(*n);
  195. // reserve space for values
  196. std::size_t nodes_count = calculate_nodes_count(values_count, subtree_counts);
  197. rtree::elements(in).reserve(nodes_count); // MAY THROW (A)
  198. // calculate values box and copy values
  199. Box elements_box;
  200. geometry::assign_inverse(elements_box);
  201. per_level_packets(first, last, hint_box, values_count, subtree_counts, next_subtree_counts,
  202. rtree::elements(in), elements_box,
  203. parameters, translator, allocators);
  204. auto_remover.release();
  205. return internal_element(elements_box, n);
  206. }
  207. template <typename EIt> inline static
  208. void per_level_packets(EIt first, EIt last, Box const& hint_box,
  209. std::size_t values_count,
  210. subtree_elements_counts const& subtree_counts,
  211. subtree_elements_counts const& next_subtree_counts,
  212. internal_elements & elements, Box & elements_box,
  213. parameters_type const& parameters, Translator const& translator, Allocators & allocators)
  214. {
  215. BOOST_ASSERT(0 < std::distance(first, last) && static_cast<std::size_t>(std::distance(first, last)) == values_count);
  216. BOOST_ASSERT_MSG( subtree_counts.minc <= values_count, "too small number of elements");
  217. // only one packet
  218. if ( values_count <= subtree_counts.maxc )
  219. {
  220. // the end, move to the next level
  221. internal_element el = per_level(first, last, hint_box, values_count, next_subtree_counts,
  222. parameters, translator, allocators);
  223. // in case if push_back() do throw here
  224. // and even if this is not probable (previously reserved memory, nonthrowing pairs copy)
  225. // this case is also tested by exceptions test.
  226. node_auto_ptr auto_remover(el.second, allocators);
  227. // this container should have memory allocated, reserve() called outside
  228. elements.push_back(el); // MAY THROW (A?,C) - however in normal conditions shouldn't
  229. auto_remover.release();
  230. geometry::expand(elements_box, el.first);
  231. return;
  232. }
  233. std::size_t median_count = calculate_median_count(values_count, subtree_counts);
  234. EIt median = first + median_count;
  235. coordinate_type greatest_length;
  236. std::size_t greatest_dim_index = 0;
  237. pack_utils::biggest_edge<dimension>::apply(hint_box, greatest_length, greatest_dim_index);
  238. Box left, right;
  239. pack_utils::partial_sort_and_half_boxes<0, dimension>
  240. ::apply(first, median, last, hint_box, left, right, greatest_dim_index);
  241. per_level_packets(first, median, left,
  242. median_count, subtree_counts, next_subtree_counts,
  243. elements, elements_box,
  244. parameters, translator, allocators);
  245. per_level_packets(median, last, right,
  246. values_count - median_count, subtree_counts, next_subtree_counts,
  247. elements, elements_box,
  248. parameters, translator, allocators);
  249. }
  250. inline static
  251. subtree_elements_counts calculate_subtree_elements_counts(std::size_t elements_count, parameters_type const& parameters, size_type & leafs_level)
  252. {
  253. (void)parameters;
  254. subtree_elements_counts res(1, 1);
  255. leafs_level = 0;
  256. std::size_t smax = parameters.get_max_elements();
  257. for ( ; smax < elements_count ; smax *= parameters.get_max_elements(), ++leafs_level )
  258. res.maxc = smax;
  259. res.minc = parameters.get_min_elements() * (res.maxc / parameters.get_max_elements());
  260. return res;
  261. }
  262. inline static
  263. std::size_t calculate_nodes_count(std::size_t count,
  264. subtree_elements_counts const& subtree_counts)
  265. {
  266. std::size_t n = count / subtree_counts.maxc;
  267. std::size_t r = count % subtree_counts.maxc;
  268. if ( 0 < r && r < subtree_counts.minc )
  269. {
  270. std::size_t count_minus_min = count - subtree_counts.minc;
  271. n = count_minus_min / subtree_counts.maxc;
  272. r = count_minus_min % subtree_counts.maxc;
  273. ++n;
  274. }
  275. if ( 0 < r )
  276. ++n;
  277. return n;
  278. }
  279. inline static
  280. std::size_t calculate_median_count(std::size_t count,
  281. subtree_elements_counts const& subtree_counts)
  282. {
  283. // e.g. for max = 5, min = 2, count = 52, subtree_max = 25, subtree_min = 10
  284. std::size_t n = count / subtree_counts.maxc; // e.g. 52 / 25 = 2
  285. std::size_t r = count % subtree_counts.maxc; // e.g. 52 % 25 = 2
  286. std::size_t median_count = (n / 2) * subtree_counts.maxc; // e.g. 2 / 2 * 25 = 25
  287. if ( 0 != r ) // e.g. 0 != 2
  288. {
  289. if ( subtree_counts.minc <= r ) // e.g. 10 <= 2 == false
  290. {
  291. //BOOST_ASSERT_MSG(0 < n, "unexpected value");
  292. median_count = ((n+1)/2) * subtree_counts.maxc; // if calculated ((2+1)/2) * 25 which would be ok, but not in all cases
  293. }
  294. else // r < subtree_counts.second // e.g. 2 < 10 == true
  295. {
  296. std::size_t count_minus_min = count - subtree_counts.minc; // e.g. 52 - 10 = 42
  297. n = count_minus_min / subtree_counts.maxc; // e.g. 42 / 25 = 1
  298. r = count_minus_min % subtree_counts.maxc; // e.g. 42 % 25 = 17
  299. if ( r == 0 ) // e.g. false
  300. {
  301. // n can't be equal to 0 because then there wouldn't be any element in the other node
  302. //BOOST_ASSERT_MSG(0 < n, "unexpected value");
  303. median_count = ((n+1)/2) * subtree_counts.maxc; // if calculated ((1+1)/2) * 25 which would be ok, but not in all cases
  304. }
  305. else
  306. {
  307. if ( n == 0 ) // e.g. false
  308. median_count = r; // if calculated -> 17 which is wrong!
  309. else
  310. median_count = ((n+2)/2) * subtree_counts.maxc; // e.g. ((1+2)/2) * 25 = 25
  311. }
  312. }
  313. }
  314. return median_count;
  315. }
  316. };
  317. }}}}} // namespace boost::geometry::index::detail::rtree
  318. #endif // BOOST_GEOMETRY_INDEX_DETAIL_RTREE_PACK_CREATE_HPP