Chris@16: // Boost.Geometry Index
Chris@16: //
Chris@16: // R-tree initial packing
Chris@16: //
Chris@101: // Copyright (c) 2011-2015 Adam Wulkiewicz, Lodz, Poland.
Chris@16: //
Chris@16: // Use, modification and distribution is subject to the Boost Software License,
Chris@16: // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
Chris@16: // http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_GEOMETRY_INDEX_DETAIL_RTREE_PACK_CREATE_HPP
Chris@16: #define BOOST_GEOMETRY_INDEX_DETAIL_RTREE_PACK_CREATE_HPP
Chris@16: 
Chris@16: namespace boost { namespace geometry { namespace index { namespace detail { namespace rtree {
Chris@16: 
Chris@16: namespace pack_utils {
Chris@16: 
Chris@16: template <std::size_t Dimension>
Chris@16: struct biggest_edge
Chris@16: {
Chris@16:     BOOST_STATIC_ASSERT(0 < Dimension);
Chris@16:     template <typename Box>
Chris@16:     static inline void apply(Box const& box, typename coordinate_type<Box>::type & length, std::size_t & dim_index)
Chris@16:     {
Chris@16:         biggest_edge<Dimension-1>::apply(box, length, dim_index);
Chris@16:         typename coordinate_type<Box>::type curr
Chris@16:             = geometry::get<max_corner, Dimension-1>(box) - geometry::get<min_corner, Dimension-1>(box);
Chris@16:         if ( length < curr )
Chris@16:         {
Chris@16:             dim_index = Dimension - 1;
Chris@16:             length = curr;
Chris@16:         }
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: template <>
Chris@16: struct biggest_edge<1>
Chris@16: {
Chris@16:     template <typename Box>
Chris@16:     static inline void apply(Box const& box, typename coordinate_type<Box>::type & length, std::size_t & dim_index)
Chris@16:     {
Chris@16:         dim_index = 0;
Chris@16:         length = geometry::get<max_corner, 0>(box) - geometry::get<min_corner, 0>(box);
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: template <std::size_t I>
Chris@16: struct point_entries_comparer
Chris@16: {
Chris@16:     template <typename PointEntry>
Chris@16:     bool operator()(PointEntry const& e1, PointEntry const& e2) const
Chris@16:     {
Chris@16:         return geometry::get<I>(e1.first) < geometry::get<I>(e2.first);
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: template <std::size_t I, std::size_t Dimension>
Chris@101: struct nth_element_and_half_boxes
Chris@16: {
Chris@16:     template <typename EIt, typename Box>
Chris@16:     static inline void apply(EIt first, EIt median, EIt last, Box const& box, Box & left, Box & right, std::size_t dim_index)
Chris@16:     {
Chris@16:         if ( I == dim_index )
Chris@16:         {
Chris@101:             std::nth_element(first, median, last, point_entries_comparer<I>());
Chris@16: 
Chris@16:             geometry::convert(box, left);
Chris@16:             geometry::convert(box, right);
Chris@16:             typename coordinate_type<Box>::type edge_len
Chris@16:                 = geometry::get<max_corner, I>(box) - geometry::get<min_corner, I>(box);
Chris@16:             typename coordinate_type<Box>::type median
Chris@16:                 = geometry::get<min_corner, I>(box) + edge_len / 2;
Chris@16:             geometry::set<max_corner, I>(left, median);
Chris@16:             geometry::set<min_corner, I>(right, median);
Chris@16:         }
Chris@16:         else
Chris@101:             nth_element_and_half_boxes<I+1, Dimension>::apply(first, median, last, box, left, right, dim_index);
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: template <std::size_t Dimension>
Chris@101: struct nth_element_and_half_boxes<Dimension, Dimension>
Chris@16: {
Chris@16:     template <typename EIt, typename Box>
Chris@16:     static inline void apply(EIt , EIt , EIt , Box const& , Box & , Box & , std::size_t ) {}
Chris@16: };
Chris@16: 
Chris@16: } // namespace pack_utils
Chris@16: 
Chris@16: // STR leafs number are calculated as rcount/max
Chris@16: // and the number of splitting planes for each dimension as (count/max)^(1/dimension)
Chris@16: // <-> for dimension==2 -> sqrt(count/max)
Chris@16: //
Chris@16: // The main flaw of this algorithm is that the resulting tree will have bad structure for:
Chris@16: // 1. non-uniformly distributed elements
Chris@16: //      Statistic check could be performed, e.g. based on variance of lengths of elements edges for each dimension
Chris@16: // 2. elements distributed mainly along one axis
Chris@16: //      Calculate bounding box of all elements and then number of dividing planes for a dimension
Chris@16: //      from the length of BB edge for this dimension (more or less assuming that elements are uniformly-distributed squares)
Chris@16: //
Chris@16: // Another thing is that the last node may have less elements than Max or even Min.
Chris@16: // The number of splitting planes must be chosen more carefully than count/max
Chris@16: //
Chris@16: // This algorithm is something between STR and TGS
Chris@16: // it is more similar to the top-down recursive kd-tree creation algorithm
Chris@16: // using object median split and split axis of greatest BB edge
Chris@16: // BB is only used as a hint (assuming objects are distributed uniformly)
Chris@16: //
Chris@16: // Implemented algorithm guarantees that the number of elements in nodes will be between Min and Max
Chris@16: // and that nodes are packed as tightly as possible
Chris@16: // e.g. for 177 values Max = 5 and Min = 2 it will construct the following tree:
Chris@16: // ROOT                 177
Chris@16: // L1          125               52
Chris@16: // L2  25  25  25  25  25   25  17    10
Chris@16: // L3  5x5 5x5 5x5 5x5 5x5  5x5 3x5+2 2x5
Chris@16: 
Chris@16: template <typename Value, typename Options, typename Translator, typename Box, typename Allocators>
Chris@16: class pack
Chris@16: {
Chris@16:     typedef typename rtree::node<Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag>::type node;
Chris@16:     typedef typename rtree::internal_node<Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag>::type internal_node;
Chris@16:     typedef typename rtree::leaf<Value, typename Options::parameters_type, Box, Allocators, typename Options::node_tag>::type leaf;
Chris@16: 
Chris@16:     typedef typename Allocators::node_pointer node_pointer;
Chris@101:     typedef rtree::subtree_destroyer<Value, Options, Translator, Box, Allocators> subtree_destroyer;
Chris@16:     typedef typename Allocators::size_type size_type;
Chris@16: 
Chris@101:     typedef typename geometry::point_type<Box>::type point_type;
Chris@101:     typedef typename geometry::coordinate_type<point_type>::type coordinate_type;
Chris@16:     typedef typename detail::default_content_result<Box>::type content_type;
Chris@16:     typedef typename Options::parameters_type parameters_type;
Chris@101:     static const std::size_t dimension = geometry::dimension<point_type>::value;
Chris@16: 
Chris@16:     typedef typename rtree::container_from_elements_type<
Chris@16:         typename rtree::elements_type<leaf>::type,
Chris@16:         std::size_t
Chris@16:     >::type values_counts_container;
Chris@16: 
Chris@16:     typedef typename rtree::elements_type<internal_node>::type internal_elements;
Chris@16:     typedef typename internal_elements::value_type internal_element;
Chris@16: 
Chris@16: public:
Chris@16:     // Arbitrary iterators
Chris@16:     template <typename InIt> inline static
Chris@16:     node_pointer apply(InIt first, InIt last, size_type & values_count, size_type & leafs_level,
Chris@16:                        parameters_type const& parameters, Translator const& translator, Allocators & allocators)
Chris@16:     {
Chris@16:         typedef typename std::iterator_traits<InIt>::difference_type diff_type;
Chris@16:             
Chris@16:         diff_type diff = std::distance(first, last);
Chris@16:         if ( diff <= 0 )
Chris@16:             return node_pointer(0);
Chris@16: 
Chris@16:         typedef std::pair<point_type, InIt> entry_type;
Chris@16:         std::vector<entry_type> entries;
Chris@16: 
Chris@16:         values_count = static_cast<size_type>(diff);
Chris@16:         entries.reserve(values_count);
Chris@16:         
Chris@16:         Box hint_box;
Chris@16:         geometry::assign_inverse(hint_box);
Chris@16:         for ( ; first != last ; ++first )
Chris@16:         {
Chris@101:             // NOTE: support for iterators not returning true references adapted
Chris@101:             // to Geometry concept and default translator returning true reference
Chris@101:             // An alternative would be to dereference the iterator and translate
Chris@101:             // in one expression each time the indexable was needed.
Chris@101:             typename std::iterator_traits<InIt>::reference in_ref = *first;
Chris@101:             typename Translator::result_type indexable = translator(in_ref);
Chris@101: 
Chris@101:             // NOTE: added for consistency with insert()
Chris@101:             // CONSIDER: alternative - ignore invalid indexable or throw an exception
Chris@101:             BOOST_GEOMETRY_INDEX_ASSERT(detail::is_valid(indexable), "Indexable is invalid");
Chris@101: 
Chris@101:             geometry::expand(hint_box, indexable);
Chris@16: 
Chris@16:             point_type pt;
Chris@101:             geometry::centroid(indexable, pt);
Chris@16:             entries.push_back(std::make_pair(pt, first));
Chris@16:         }
Chris@16: 
Chris@16:         subtree_elements_counts subtree_counts = calculate_subtree_elements_counts(values_count, parameters, leafs_level);
Chris@16:         internal_element el = per_level(entries.begin(), entries.end(), hint_box, values_count, subtree_counts,
Chris@16:                                         parameters, translator, allocators);
Chris@16: 
Chris@16:         return el.second;
Chris@16:     }
Chris@16: 
Chris@16: private:
Chris@16:     struct subtree_elements_counts
Chris@16:     {
Chris@16:         subtree_elements_counts(std::size_t ma, std::size_t mi) : maxc(ma), minc(mi) {}
Chris@16:         std::size_t maxc;
Chris@16:         std::size_t minc;
Chris@16:     };
Chris@16: 
Chris@16:     template <typename EIt> inline static
Chris@16:     internal_element per_level(EIt first, EIt last, Box const& hint_box, std::size_t values_count, subtree_elements_counts const& subtree_counts,
Chris@16:                                parameters_type const& parameters, Translator const& translator, Allocators & allocators)
Chris@16:     {
Chris@101:         BOOST_GEOMETRY_INDEX_ASSERT(0 < std::distance(first, last) && static_cast<std::size_t>(std::distance(first, last)) == values_count,
Chris@101:                                     "unexpected parameters");
Chris@16: 
Chris@16:         if ( subtree_counts.maxc <= 1 )
Chris@16:         {
Chris@16:             // ROOT or LEAF
Chris@101:             BOOST_GEOMETRY_INDEX_ASSERT(values_count <= parameters.get_max_elements(),
Chris@101:                                         "too big number of elements");
Chris@16:             // if !root check m_parameters.get_min_elements() <= count
Chris@16: 
Chris@16:             // create new leaf node
Chris@16:             node_pointer n = rtree::create_node<Allocators, leaf>::apply(allocators);                       // MAY THROW (A)
Chris@101:             subtree_destroyer auto_remover(n, allocators);
Chris@16:             leaf & l = rtree::get<leaf>(*n);
Chris@16: 
Chris@16:             // reserve space for values
Chris@16:             rtree::elements(l).reserve(values_count);                                                       // MAY THROW (A)
Chris@16:             // calculate values box and copy values
Chris@16:             Box elements_box;
Chris@16:             geometry::assign_inverse(elements_box);
Chris@16:             for ( ; first != last ; ++first )
Chris@16:             {
Chris@101:                 // NOTE: push_back() must be called at the end in order to support move_iterator.
Chris@101:                 //       The iterator is dereferenced 2x (no temporary reference) to support
Chris@101:                 //       non-true reference types and move_iterator without boost::forward<>.
Chris@101:                 geometry::expand(elements_box, translator(*(first->second)));
Chris@16:                 rtree::elements(l).push_back(*(first->second));                                             // MAY THROW (A?,C)
Chris@16:             }
Chris@16: 
Chris@16:             auto_remover.release();
Chris@16:             return internal_element(elements_box, n);
Chris@16:         }
Chris@16: 
Chris@16:         // calculate next max and min subtree counts
Chris@16:         subtree_elements_counts next_subtree_counts = subtree_counts;
Chris@16:         next_subtree_counts.maxc /= parameters.get_max_elements();
Chris@16:         next_subtree_counts.minc /= parameters.get_max_elements();
Chris@16: 
Chris@16:         // create new internal node
Chris@16:         node_pointer n = rtree::create_node<Allocators, internal_node>::apply(allocators);                  // MAY THROW (A)
Chris@101:         subtree_destroyer auto_remover(n, allocators);
Chris@16:         internal_node & in = rtree::get<internal_node>(*n);
Chris@16: 
Chris@16:         // reserve space for values
Chris@16:         std::size_t nodes_count = calculate_nodes_count(values_count, subtree_counts);
Chris@16:         rtree::elements(in).reserve(nodes_count);                                                           // MAY THROW (A)
Chris@16:         // calculate values box and copy values
Chris@16:         Box elements_box;
Chris@16:         geometry::assign_inverse(elements_box);
Chris@16: 
Chris@16:         per_level_packets(first, last, hint_box, values_count, subtree_counts, next_subtree_counts,
Chris@16:                           rtree::elements(in), elements_box,
Chris@16:                           parameters, translator, allocators);
Chris@16: 
Chris@16:         auto_remover.release();
Chris@16:         return internal_element(elements_box, n);
Chris@16:     }
Chris@16: 
Chris@16:     template <typename EIt> inline static
Chris@16:     void per_level_packets(EIt first, EIt last, Box const& hint_box,
Chris@16:                            std::size_t values_count,
Chris@16:                            subtree_elements_counts const& subtree_counts,
Chris@16:                            subtree_elements_counts const& next_subtree_counts,
Chris@16:                            internal_elements & elements, Box & elements_box,
Chris@16:                            parameters_type const& parameters, Translator const& translator, Allocators & allocators)
Chris@16:     {
Chris@101:         BOOST_GEOMETRY_INDEX_ASSERT(0 < std::distance(first, last) && static_cast<std::size_t>(std::distance(first, last)) == values_count,
Chris@101:                                     "unexpected parameters");
Chris@16: 
Chris@101:         BOOST_GEOMETRY_INDEX_ASSERT(subtree_counts.minc <= values_count,
Chris@101:                                     "too small number of elements");
Chris@16: 
Chris@16:         // only one packet
Chris@16:         if ( values_count <= subtree_counts.maxc )
Chris@16:         {
Chris@16:             // the end, move to the next level
Chris@16:             internal_element el = per_level(first, last, hint_box, values_count, next_subtree_counts,
Chris@16:                                             parameters, translator, allocators);
Chris@16: 
Chris@16:             // in case if push_back() do throw here
Chris@16:             // and even if this is not probable (previously reserved memory, nonthrowing pairs copy)
Chris@16:             // this case is also tested by exceptions test.
Chris@101:             subtree_destroyer auto_remover(el.second, allocators);
Chris@16:             // this container should have memory allocated, reserve() called outside
Chris@16:             elements.push_back(el);                                                 // MAY THROW (A?,C) - however in normal conditions shouldn't
Chris@16:             auto_remover.release();
Chris@16: 
Chris@16:             geometry::expand(elements_box, el.first);
Chris@16:             return;
Chris@16:         }
Chris@16:         
Chris@16:         std::size_t median_count = calculate_median_count(values_count, subtree_counts);
Chris@16:         EIt median = first + median_count;
Chris@16: 
Chris@16:         coordinate_type greatest_length;
Chris@16:         std::size_t greatest_dim_index = 0;
Chris@16:         pack_utils::biggest_edge<dimension>::apply(hint_box, greatest_length, greatest_dim_index);
Chris@16:         Box left, right;
Chris@101:         pack_utils::nth_element_and_half_boxes<0, dimension>
Chris@16:             ::apply(first, median, last, hint_box, left, right, greatest_dim_index);
Chris@16:         
Chris@16:         per_level_packets(first, median, left,
Chris@16:                           median_count, subtree_counts, next_subtree_counts,
Chris@16:                           elements, elements_box,
Chris@16:                           parameters, translator, allocators);
Chris@16:         per_level_packets(median, last, right,
Chris@16:                           values_count - median_count, subtree_counts, next_subtree_counts,
Chris@16:                           elements, elements_box,
Chris@16:                           parameters, translator, allocators);
Chris@16:     }
Chris@16: 
Chris@16:     inline static
Chris@16:     subtree_elements_counts calculate_subtree_elements_counts(std::size_t elements_count, parameters_type const& parameters, size_type & leafs_level)
Chris@16:     {
Chris@101:         boost::ignore_unused_variable_warning(parameters);
Chris@16: 
Chris@16:         subtree_elements_counts res(1, 1);
Chris@16:         leafs_level = 0;
Chris@16: 
Chris@16:         std::size_t smax = parameters.get_max_elements();
Chris@16:         for ( ; smax < elements_count ; smax *= parameters.get_max_elements(), ++leafs_level )
Chris@16:             res.maxc = smax;
Chris@16: 
Chris@16:         res.minc = parameters.get_min_elements() * (res.maxc / parameters.get_max_elements());
Chris@16: 
Chris@16:         return res;
Chris@16:     }
Chris@16: 
Chris@16:     inline static
Chris@16:     std::size_t calculate_nodes_count(std::size_t count,
Chris@16:                                       subtree_elements_counts const& subtree_counts)
Chris@16:     {
Chris@16:         std::size_t n = count / subtree_counts.maxc;
Chris@16:         std::size_t r = count % subtree_counts.maxc;
Chris@16: 
Chris@16:         if ( 0 < r && r < subtree_counts.minc )
Chris@16:         {
Chris@16:             std::size_t count_minus_min = count - subtree_counts.minc;
Chris@16:             n = count_minus_min / subtree_counts.maxc;
Chris@16:             r = count_minus_min % subtree_counts.maxc;
Chris@16:             ++n;
Chris@16:         }
Chris@16: 
Chris@16:         if ( 0 < r )
Chris@16:             ++n;
Chris@16: 
Chris@16:         return n;
Chris@16:     }
Chris@16: 
Chris@16:     inline static
Chris@16:     std::size_t calculate_median_count(std::size_t count,
Chris@16:                                        subtree_elements_counts const& subtree_counts)
Chris@16:     {
Chris@16:         // e.g. for max = 5, min = 2, count = 52, subtree_max = 25, subtree_min = 10
Chris@16: 
Chris@16:         std::size_t n = count / subtree_counts.maxc; // e.g. 52 / 25 = 2
Chris@16:         std::size_t r = count % subtree_counts.maxc; // e.g. 52 % 25 = 2
Chris@16:         std::size_t median_count = (n / 2) * subtree_counts.maxc; // e.g. 2 / 2 * 25 = 25
Chris@16: 
Chris@16:         if ( 0 != r ) // e.g. 0 != 2
Chris@16:         {
Chris@16:             if ( subtree_counts.minc <= r ) // e.g. 10 <= 2 == false
Chris@16:             {
Chris@101:                 //BOOST_GEOMETRY_INDEX_ASSERT(0 < n, "unexpected value");
Chris@16:                 median_count = ((n+1)/2) * subtree_counts.maxc; // if calculated ((2+1)/2) * 25 which would be ok, but not in all cases
Chris@16:             }
Chris@16:             else // r < subtree_counts.second  // e.g. 2 < 10 == true
Chris@16:             {
Chris@16:                 std::size_t count_minus_min = count - subtree_counts.minc; // e.g. 52 - 10 = 42
Chris@16:                 n = count_minus_min / subtree_counts.maxc; // e.g. 42 / 25 = 1
Chris@16:                 r = count_minus_min % subtree_counts.maxc; // e.g. 42 % 25 = 17
Chris@16:                 if ( r == 0 )                               // e.g. false
Chris@16:                 {
Chris@16:                     // n can't be equal to 0 because then there wouldn't be any element in the other node
Chris@101:                     //BOOST_GEOMETRY_INDEX_ASSERT(0 < n, "unexpected value");
Chris@16:                     median_count = ((n+1)/2) * subtree_counts.maxc;     // if calculated ((1+1)/2) * 25 which would be ok, but not in all cases
Chris@16:                 }
Chris@16:                 else
Chris@16:                 {
Chris@16:                     if ( n == 0 )                                        // e.g. false
Chris@16:                         median_count = r;                                // if calculated -> 17 which is wrong!
Chris@16:                     else
Chris@16:                         median_count = ((n+2)/2) * subtree_counts.maxc; // e.g. ((1+2)/2) * 25 = 25
Chris@16:                 }
Chris@16:             }
Chris@16:         }
Chris@16: 
Chris@16:         return median_count;
Chris@16:     }
Chris@16: };
Chris@16: 
Chris@16: }}}}} // namespace boost::geometry::index::detail::rtree
Chris@16: 
Chris@16: #endif // BOOST_GEOMETRY_INDEX_DETAIL_RTREE_PACK_CREATE_HPP