Skip to content
This repository
Browse code

+ add boost/geometry/extensions/index/rtree

  (TODO: remove when it's part of boost release)
  • Loading branch information...
commit 457afbdbc44fcc08194b6d1a6d9ab1eb73d96c7e 1 parent b42e498
Artem Pavlenko artemp authored
68 boost/geometry/extensions/index/rtree/helpers.hpp
... ... @@ -0,0 +1,68 @@
  1 +// Boost.Geometry (aka GGL, Generic Geometry Library)
  2 +
  3 +// Boost.SpatialIndex - geometry helper functions
  4 +//
  5 +// Copyright 2008 Federico J. Fernandez.
  6 +// Use, modification and distribution is subject to the Boost Software License,
  7 +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  8 +// http://www.boost.org/LICENSE_1_0.txt)
  9 +
  10 +#ifndef BOOST_GEOMETRY_GGL_INDEX_RTREE_HELPERS_HPP
  11 +#define BOOST_GEOMETRY_GGL_INDEX_RTREE_HELPERS_HPP
  12 +
  13 +#include <boost/geometry/algorithms/area.hpp>
  14 +#include <boost/geometry/algorithms/disjoint.hpp>
  15 +#include <boost/geometry/core/point_type.hpp>
  16 +
  17 +namespace boost { namespace geometry { namespace index {
  18 +
  19 +/**
  20 + * \brief Given two boxes, returns the minimal box that contains them
  21 + */
  22 +// TODO: use geometry::expand
  23 +template <typename Box>
  24 +inline Box enlarge_box(Box const& b1, Box const& b2)
  25 +{
  26 + // TODO: mloskot - Refactor to readable form. Fix VC++8.0 min/max warnings:
  27 + // warning C4002: too many actual parameters for macro 'min
  28 +
  29 + typedef typename geometry::point_type<Box>::type point_type;
  30 +
  31 + point_type pmin(
  32 + geometry::get<min_corner, 0>(b1) < geometry::get<min_corner, 0>(b2)
  33 + ? geometry::get<min_corner, 0>(b1) : geometry::get<min_corner, 0>(b2),
  34 + geometry::get<min_corner, 1>(b1) < geometry::get<min_corner, 1>(b2)
  35 + ? geometry::get<min_corner, 1>(b1) : geometry::get<min_corner, 1>(b2));
  36 +
  37 + point_type pmax(
  38 + geometry::get<max_corner, 0>(b1) > geometry::get<max_corner, 0>(b2)
  39 + ? geometry::get<max_corner, 0>(b1) : geometry::get<max_corner, 0>(b2),
  40 + geometry::get<max_corner, 1>(b1) > geometry::get<max_corner, 1>(b2)
  41 + ? geometry::get<max_corner, 1>(b1) : geometry::get<max_corner, 1>(b2));
  42 +
  43 + return Box(pmin, pmax);
  44 +}
  45 +
  46 +/**
  47 + * \brief Compute the area of the union of b1 and b2
  48 + */
  49 +template <typename Box>
  50 +inline typename default_area_result<Box>::type compute_union_area(Box const& b1, Box const& b2)
  51 +{
  52 + Box enlarged_box = enlarge_box(b1, b2);
  53 + return geometry::area(enlarged_box);
  54 +}
  55 +
  56 +/**
  57 + * \brief Checks if boxes intersects
  58 + */
  59 +// TODO: move to geometry::intersects
  60 +template <typename Box>
  61 +inline bool is_overlapping(Box const& b1, Box const& b2)
  62 +{
  63 + return ! geometry::disjoint(b1, b2);
  64 +}
  65 +
  66 +}}} // namespace boost::geometry::index
  67 +
  68 +#endif // BOOST_GEOMETRY_GGL_INDEX_RTREE_HELPERS_HPP
774 boost/geometry/extensions/index/rtree/rtree.hpp
... ... @@ -0,0 +1,774 @@
  1 +// Boost.Geometry (aka GGL, Generic Geometry Library)
  2 +
  3 +// Boost.SpatialIndex - rtree implementation
  4 +//
  5 +// Copyright 2008 Federico J. Fernandez.
  6 +// Use, modification and distribution is subject to the Boost Software License,
  7 +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  8 +// http://www.boost.org/LICENSE_1_0.txt)
  9 +
  10 +#ifndef BOOST_GEOMETRY_EXTENSIONS_INDEX_RTREE_RTREE_HPP
  11 +#define BOOST_GEOMETRY_EXTENSIONS_INDEX_RTREE_RTREE_HPP
  12 +
  13 +#include <cstddef>
  14 +#include <iostream> // TODO: Remove if print() is removed
  15 +#include <stdexcept>
  16 +#include <utility>
  17 +#include <vector>
  18 +
  19 +#include <boost/concept_check.hpp>
  20 +#include <boost/shared_ptr.hpp>
  21 +
  22 +#include <boost/geometry/algorithms/area.hpp>
  23 +
  24 +#include <boost/geometry/extensions/index/rtree/rtree_node.hpp>
  25 +#include <boost/geometry/extensions/index/rtree/rtree_leaf.hpp>
  26 +
  27 +namespace boost { namespace geometry { namespace index
  28 +{
  29 +
  30 +template <typename Box, typename Value >
  31 +class rtree
  32 +{
  33 +public:
  34 +
  35 + typedef boost::shared_ptr<rtree_node<Box, Value> > node_pointer;
  36 + typedef boost::shared_ptr<rtree_leaf<Box, Value> > leaf_pointer;
  37 +
  38 + /**
  39 + * \brief Creates a rtree with 'maximum' elements per node and 'minimum'.
  40 + */
  41 + rtree(unsigned int const& maximum, unsigned int const& minimum)
  42 + : m_count(0)
  43 + , m_min_elems_per_node(minimum)
  44 + , m_max_elems_per_node(maximum)
  45 + , m_root(new rtree_node<Box, Value>(node_pointer(), 1))
  46 + {
  47 + }
  48 +
  49 + /**
  50 + * \brief Creates a rtree with maximum elements per node
  51 + * and minimum (box is ignored).
  52 + */
  53 + rtree(Box const& box, unsigned int const& maximum, unsigned int const& minimum)
  54 + : m_count(0)
  55 + , m_min_elems_per_node(minimum)
  56 + , m_max_elems_per_node(maximum)
  57 + , m_root(new rtree_node<Box, Value>(node_pointer(), 1))
  58 + {
  59 + boost::ignore_unused_variable_warning(box);
  60 + }
  61 +
  62 + /**
  63 + * \brief destructor (virtual because we have virtual functions)
  64 + */
  65 + virtual ~rtree() {}
  66 +
  67 +
  68 + /**
  69 + * \brief Remove elements inside the 'box'
  70 + */
  71 + inline void remove(Box const& box)
  72 + {
  73 + try
  74 + {
  75 + node_pointer leaf(choose_exact_leaf(box));
  76 + typename rtree_leaf<Box, Value>::leaf_map q_leaves;
  77 +
  78 + leaf->remove(box);
  79 +
  80 + if (leaf->elements() < m_min_elems_per_node && elements() > m_min_elems_per_node)
  81 + {
  82 + q_leaves = leaf->get_leaves();
  83 +
  84 + // we remove the leaf_node in the parent node because now it's empty
  85 + leaf->get_parent()->remove(leaf->get_parent()->get_box(leaf));
  86 + }
  87 +
  88 + typename rtree_node<Box, Value>::node_map q_nodes;
  89 + condense_tree(leaf, q_nodes);
  90 +
  91 + std::vector<std::pair<Box, Value> > s;
  92 + for (typename rtree_node<Box, Value>::node_map::const_iterator it = q_nodes.begin();
  93 + it != q_nodes.end(); ++it)
  94 + {
  95 + typename rtree_leaf<Box, Value>::leaf_map leaves = it->second->get_leaves();
  96 +
  97 + // reinserting leaves from nodes
  98 + for (typename rtree_leaf<Box, Value>::leaf_map::const_iterator itl = leaves.begin();
  99 + itl != leaves.end(); ++itl)
  100 + {
  101 + s.push_back(*itl);
  102 + }
  103 + }
  104 +
  105 + for (typename std::vector<std::pair<Box, Value> >::const_iterator it = s.begin(); it != s.end(); ++it)
  106 + {
  107 + m_count--;
  108 + insert(it->first, it->second);
  109 + }
  110 +
  111 + // if the root has only one child and the child is not a leaf,
  112 + // make it the root
  113 + if (m_root->elements() == 1)
  114 + {
  115 + if (!m_root->first_element()->is_leaf())
  116 + {
  117 + m_root = m_root->first_element();
  118 + }
  119 + }
  120 + // reinserting leaves
  121 + for (typename rtree_leaf<Box, Value>::leaf_map::const_iterator it = q_leaves.begin();
  122 + it != q_leaves.end(); ++it)
  123 + {
  124 + m_count--;
  125 + insert(it->first, it->second);
  126 + }
  127 +
  128 + m_count--;
  129 + }
  130 + catch(std::logic_error & e)
  131 + {
  132 + // TODO: mloskot - replace with Boost.Geometry exception
  133 +
  134 + // not found
  135 + std::cerr << e.what() << std::endl;
  136 + return;
  137 + }
  138 + }
  139 +
  140 + /**
  141 + * \brief Remove element inside the box with value
  142 + */
  143 + void remove(Box const& box, Value const& value)
  144 + {
  145 + try
  146 + {
  147 + node_pointer leaf;
  148 +
  149 + // find possible leaves
  150 + typedef typename std::vector<node_pointer > node_type;
  151 + node_type nodes;
  152 + m_root->find_leaves(box, nodes);
  153 +
  154 + // refine the result
  155 + for (typename node_type::const_iterator it = nodes.begin(); it != nodes.end(); ++it)
  156 + {
  157 + leaf = *it;
  158 + try
  159 + {
  160 + leaf->remove(value);
  161 + break;
  162 + } catch (...)
  163 + {
  164 + leaf = node_pointer();
  165 + }
  166 + }
  167 +
  168 + if (!leaf)
  169 + return;
  170 +
  171 + typename rtree_leaf < Box, Value >::leaf_map q_leaves;
  172 +
  173 + if (leaf->elements() < m_min_elems_per_node && elements() > m_min_elems_per_node)
  174 + {
  175 + q_leaves = leaf->get_leaves();
  176 +
  177 + // we remove the leaf_node in the parent node because now it's empty
  178 + leaf->get_parent()->remove(leaf->get_parent()->get_box(leaf));
  179 + }
  180 +
  181 + typename rtree_node<Box, Value>::node_map q_nodes;
  182 + condense_tree(leaf, q_nodes);
  183 +
  184 + std::vector<std::pair<Box, Value> > s;
  185 + for (typename rtree_node<Box, Value>::node_map::const_iterator it = q_nodes.begin();
  186 + it != q_nodes.end(); ++it)
  187 + {
  188 + typename rtree_leaf<Box, Value>::leaf_map leaves = it->second->get_leaves();
  189 +
  190 + // reinserting leaves from nodes
  191 + for (typename rtree_leaf<Box, Value>::leaf_map::const_iterator itl = leaves.begin();
  192 + itl != leaves.end(); ++itl)
  193 + {
  194 + s.push_back(*itl);
  195 + }
  196 + }
  197 +
  198 + for (typename std::vector<std::pair<Box, Value> >::const_iterator it = s.begin(); it != s.end(); ++it)
  199 + {
  200 + m_count--;
  201 + insert(it->first, it->second);
  202 + }
  203 +
  204 + // if the root has only one child and the child is not a leaf,
  205 + // make it the root
  206 + if (m_root->elements() == 1)
  207 + {
  208 + if (!m_root->first_element()->is_leaf())
  209 + {
  210 + m_root = m_root->first_element();
  211 + }
  212 + }
  213 +
  214 + // reinserting leaves
  215 + for (typename rtree_leaf<Box, Value>::leaf_map::const_iterator it = q_leaves.begin();
  216 + it != q_leaves.end(); ++it)
  217 + {
  218 + m_count--;
  219 + insert(it->first, it->second);
  220 + }
  221 +
  222 + m_count--;
  223 +
  224 + }
  225 + catch(std::logic_error & e)
  226 + {
  227 + // TODO: mloskot - ggl exception
  228 +
  229 + // not found
  230 + std::cerr << e.what() << std::endl;
  231 + return;
  232 + }
  233 + }
  234 +
  235 + /**
  236 + * \brief Returns the number of elements.
  237 + */
  238 + inline unsigned int elements() const
  239 + {
  240 + return m_count;
  241 + }
  242 +
  243 +
  244 + /**
  245 + * \brief Inserts an element with 'box' as key with value.
  246 + */
  247 + inline void insert(Box const& box, Value const& value)
  248 + {
  249 + m_count++;
  250 +
  251 + node_pointer leaf(choose_corresponding_leaf(box));
  252 +
  253 + // check if the selected leaf is full to do the split if necessary
  254 + if (leaf->elements() >= m_max_elems_per_node)
  255 + {
  256 + leaf->insert(box, value);
  257 +
  258 + // split!
  259 + node_pointer n1(new rtree_leaf<Box, Value>(leaf->get_parent()));
  260 + node_pointer n2(new rtree_leaf<Box, Value>(leaf->get_parent()));
  261 +
  262 + split_node(leaf, n1, n2);
  263 + adjust_tree(leaf, n1, n2);
  264 + }
  265 + else
  266 + {
  267 + leaf->insert(box, value);
  268 + adjust_tree(leaf);
  269 + }
  270 + }
  271 +
  272 +
  273 + /**
  274 + * \brief Returns all the values inside 'box'
  275 + */
  276 + inline std::deque<Value> find(Box const& box) const
  277 + {
  278 + std::deque<Value> result;
  279 + m_root->find(box, result, false);
  280 + return result;
  281 + }
  282 +
  283 + /**
  284 + * \brief Print Rtree (mainly for debug)
  285 + */
  286 + inline void print()
  287 + {
  288 + std::cerr << "===================================" << std::endl;
  289 + std::cerr << " Min/Max: " << m_min_elems_per_node << " / " << m_max_elems_per_node << std::endl;
  290 + std::cerr << "Leaves: " << m_root->get_leaves().size() << std::endl;
  291 + m_root->print();
  292 + std::cerr << "===================================" << std::endl;
  293 + }
  294 +
  295 +private:
  296 +
  297 + /// number of elements
  298 + unsigned int m_count;
  299 +
  300 + /// minimum number of elements per node
  301 + unsigned int m_min_elems_per_node;
  302 +
  303 + /// maximum number of elements per node
  304 + unsigned int m_max_elems_per_node;
  305 +
  306 + /// tree root
  307 + node_pointer m_root;
  308 +
  309 + /**
  310 + * \brief Reorganize the tree after a removal. It tries to
  311 + * join nodes with less elements than m.
  312 + */
  313 + void condense_tree(node_pointer const& leaf,
  314 + typename rtree_node<Box, Value>::node_map& q_nodes)
  315 + {
  316 + if (leaf.get() == m_root.get())
  317 + {
  318 + // if it's the root we are done
  319 + return;
  320 + }
  321 +
  322 + node_pointer parent = leaf->get_parent();
  323 + parent->adjust_box(leaf);
  324 +
  325 + if (parent->elements() < m_min_elems_per_node)
  326 + {
  327 + if (parent.get() == m_root.get())
  328 + {
  329 + // if the parent is underfull and it's the root we just exit
  330 + return;
  331 + }
  332 +
  333 + // get the nodes that we should reinsert
  334 + typename rtree_node<Box, Value>::node_map this_nodes = parent->get_nodes();
  335 + for(typename rtree_node<Box, Value>::node_map::const_iterator it = this_nodes.begin();
  336 + it != this_nodes.end(); ++it)
  337 + {
  338 + q_nodes.push_back(*it);
  339 + }
  340 +
  341 + // we remove the node in the parent node because now it should be
  342 + // re inserted
  343 + parent->get_parent()->remove(parent->get_parent()->get_box(parent));
  344 + }
  345 +
  346 + condense_tree(parent, q_nodes);
  347 + }
  348 +
  349 + /**
  350 + * \brief After an insertion splits nodes with more than 'maximum' elements.
  351 + */
  352 + inline void adjust_tree(node_pointer& node)
  353 + {
  354 + if (node.get() == m_root.get())
  355 + {
  356 + // we finished the adjust
  357 + return;
  358 + }
  359 +
  360 + // as there are no splits just adjust the box of the parent and go on
  361 + node_pointer parent = node->get_parent();
  362 + parent->adjust_box(node);
  363 + adjust_tree(parent);
  364 + }
  365 +
  366 + /**
  367 + * \brief After an insertion splits nodes with more than maximum elements
  368 + * (recursive step with subtrees 'n1' and 'n2' to be joined).
  369 + */
  370 + void adjust_tree(node_pointer& leaf, node_pointer& n1, node_pointer& n2)
  371 + {
  372 + // check if we are in the root and do the split
  373 + if (leaf.get() == m_root.get())
  374 + {
  375 + node_pointer new_root(new rtree_node<Box,Value>(node_pointer (), leaf->get_level() + 1));
  376 + new_root->add_node(n1->compute_box(), n1);
  377 + new_root->add_node(n2->compute_box(), n2);
  378 +
  379 + n1->set_parent(new_root);
  380 + n2->set_parent(new_root);
  381 +
  382 + n1->update_parent(n1);
  383 + n2->update_parent(n2);
  384 +
  385 + m_root = new_root;
  386 + return;
  387 + }
  388 +
  389 + node_pointer parent = leaf->get_parent();
  390 +
  391 + parent->replace_node(leaf, n1);
  392 + parent->add_node(n2->compute_box(), n2);
  393 +
  394 + // if parent is full, split and readjust
  395 + if (parent->elements() > m_max_elems_per_node)
  396 + {
  397 + node_pointer p1(new rtree_node<Box, Value>(parent->get_parent(), parent->get_level()));
  398 + node_pointer p2(new rtree_node<Box, Value>(parent->get_parent(), parent->get_level()));
  399 +
  400 + split_node(parent, p1, p2);
  401 + adjust_tree(parent, p1, p2);
  402 + }
  403 + else
  404 + {
  405 + adjust_tree(parent);
  406 + }
  407 + }
  408 +
  409 + /**
  410 + * \brief Splits 'n' in 'n1' and 'n2'
  411 + */
  412 + void split_node(node_pointer const& n, node_pointer& n1, node_pointer& n2) const
  413 + {
  414 + unsigned int seed1 = 0;
  415 + unsigned int seed2 = 0;
  416 + std::vector<Box> boxes = n->get_boxes();
  417 +
  418 + n1->set_parent(n->get_parent());
  419 + n2->set_parent(n->get_parent());
  420 +
  421 + linear_pick_seeds(n, seed1, seed2);
  422 +
  423 + if (n->is_leaf())
  424 + {
  425 + n1->add_value(boxes[seed1], n->get_value(seed1));
  426 + n2->add_value(boxes[seed2], n->get_value(seed2));
  427 + }
  428 + else
  429 + {
  430 + n1->add_node(boxes[seed1], n->get_node(seed1));
  431 + n2->add_node(boxes[seed2], n->get_node(seed2));
  432 + }
  433 +
  434 + unsigned int index = 0;
  435 +
  436 + if (n->is_leaf())
  437 + {
  438 + // TODO: mloskot - add assert(node.size() >= 2); or similar
  439 +
  440 + typename rtree_leaf<Box, Value>::leaf_map nodes = n->get_leaves();
  441 + unsigned int remaining = nodes.size() - 2;
  442 +
  443 + for (typename rtree_leaf<Box, Value>::leaf_map::const_iterator it = nodes.begin();
  444 + it != nodes.end(); ++it, index++)
  445 + {
  446 + if (index != seed1 && index != seed2)
  447 + {
  448 + if (n1->elements() + remaining == m_min_elems_per_node)
  449 + {
  450 + n1->add_value(it->first, it->second);
  451 + continue;
  452 + }
  453 + if (n2->elements() + remaining == m_min_elems_per_node)
  454 + {
  455 + n2->add_value(it->first, it->second);
  456 + continue;
  457 + }
  458 +
  459 + remaining--;
  460 +
  461 + /// current boxes of each group
  462 + Box b1, b2;
  463 +
  464 + /// enlarged boxes of each group
  465 + Box eb1, eb2;
  466 + b1 = n1->compute_box();
  467 + b2 = n2->compute_box();
  468 +
  469 + /// areas
  470 + typedef typename coordinate_type<Box>::type coordinate_type;
  471 + coordinate_type b1_area, b2_area;
  472 + coordinate_type eb1_area, eb2_area;
  473 + b1_area = geometry::area(b1);
  474 + b2_area = geometry::area(b2);
  475 + eb1_area = compute_union_area(b1, it->first);
  476 + eb2_area = compute_union_area(b2, it->first);
  477 +
  478 + if (eb1_area - b1_area > eb2_area - b2_area)
  479 + {
  480 + n2->add_value(it->first, it->second);
  481 + }
  482 + if (eb1_area - b1_area < eb2_area - b2_area)
  483 + {
  484 + n1->add_value(it->first, it->second);
  485 + }
  486 + if (eb1_area - b1_area == eb2_area - b2_area)
  487 + {
  488 + if (b1_area < b2_area)
  489 + {
  490 + n1->add_value(it->first, it->second);
  491 + }
  492 + if (b1_area > b2_area)
  493 + {
  494 + n2->add_value(it->first, it->second);
  495 + }
  496 + if (b1_area == b2_area)
  497 + {
  498 + if (n1->elements() > n2->elements())
  499 + {
  500 + n2->add_value(it->first, it->second);
  501 + }
  502 + else
  503 + {
  504 + n1->add_value(it->first, it->second);
  505 + }
  506 + }
  507 + }
  508 + }
  509 + }
  510 + }
  511 + else
  512 + {
  513 + // TODO: mloskot - add assert(node.size() >= 2); or similar
  514 +
  515 + typename rtree_node<Box, Value>::node_map nodes = n->get_nodes();
  516 + unsigned int remaining = nodes.size() - 2;
  517 + for(typename rtree_node<Box, Value>::node_map::const_iterator it = nodes.begin();
  518 + it != nodes.end(); ++it, index++)
  519 + {
  520 +
  521 + if (index != seed1 && index != seed2)
  522 + {
  523 +
  524 + if (n1->elements() + remaining == m_min_elems_per_node)
  525 + {
  526 + n1->add_node(it->first, it->second);
  527 + continue;
  528 + }
  529 + if (n2->elements() + remaining == m_min_elems_per_node)
  530 + {
  531 + n2->add_node(it->first, it->second);
  532 + continue;
  533 + }
  534 +
  535 + remaining--;
  536 +
  537 + /// current boxes of each group
  538 + Box b1, b2;
  539 +
  540 + /// enlarged boxes of each group
  541 + Box eb1, eb2;
  542 + b1 = n1->compute_box();
  543 + b2 = n2->compute_box();
  544 +
  545 + /// areas
  546 + typedef typename coordinate_type<Box>::type coordinate_type;
  547 + coordinate_type b1_area, b2_area;
  548 + coordinate_type eb1_area, eb2_area;
  549 + b1_area = geometry::area(b1);
  550 + b2_area = geometry::area(b2);
  551 +
  552 + eb1_area = compute_union_area(b1, it->first);
  553 + eb2_area = compute_union_area(b2, it->first);
  554 +
  555 + if (eb1_area - b1_area > eb2_area - b2_area)
  556 + {
  557 + n2->add_node(it->first, it->second);
  558 + }
  559 + if (eb1_area - b1_area < eb2_area - b2_area)
  560 + {
  561 + n1->add_node(it->first, it->second);
  562 + }
  563 + if (eb1_area - b1_area == eb2_area - b2_area)
  564 + {
  565 + if (b1_area < b2_area)
  566 + {
  567 + n1->add_node(it->first, it->second);
  568 + }
  569 + if (b1_area > b2_area)
  570 + {
  571 + n2->add_node(it->first, it->second);
  572 + }
  573 + if (b1_area == b2_area)
  574 + {
  575 + if (n1->elements() > n2->elements())
  576 + {
  577 + n2->add_node(it->first, it->second);
  578 + }
  579 + else
  580 + {
  581 + n1->add_node(it->first, it->second);
  582 + }
  583 + }
  584 + }
  585 +
  586 + }
  587 + }
  588 + }
  589 + }
  590 +
  591 + /**
  592 + * \brief Choose initial values for the split algorithm (linear version)
  593 + */
  594 + void linear_pick_seeds(node_pointer const& n, unsigned int &seed1, unsigned int &seed2) const
  595 + {
  596 + // get boxes from the node
  597 + std::vector<Box>boxes = n->get_boxes();
  598 + if (boxes.size() == 0)
  599 + {
  600 + // TODO: mloskot - throw ggl exception
  601 + throw std::logic_error("Empty Node trying to Pick Seeds");
  602 + }
  603 +
  604 + // only two dim for now
  605 + // unsigned int dimensions =
  606 + // geometry::point_traits<Point>::coordinate_count;
  607 +
  608 + // find the first two elements
  609 + typedef typename coordinate_type<Box>::type coordinate_type;
  610 + coordinate_type separation_x, separation_y;
  611 + unsigned int first_x, second_x;
  612 + unsigned int first_y, second_y;
  613 + find_normalized_separations<0u>(boxes, separation_x, first_x, second_x);
  614 + find_normalized_separations<1u>(boxes, separation_y, first_y, second_y);
  615 +
  616 + if (separation_x > separation_y)
  617 + {
  618 + seed1 = first_x;
  619 + seed2 = second_x;
  620 + }
  621 + else
  622 + {
  623 + seed1 = first_y;
  624 + seed2 = second_y;
  625 + }
  626 + }
  627 +
  628 + /**
  629 + * \brief Find distances between possible initial values for the
  630 + * pick_seeds algorithm.
  631 + */
  632 + template <std::size_t D, typename T>
  633 + void find_normalized_separations(std::vector<Box> const& boxes, T& separation,
  634 + unsigned int& first, unsigned int& second) const
  635 + {
  636 + if (boxes.size() < 2)
  637 + {
  638 + throw std::logic_error("At least two boxes needed to split");
  639 + }
  640 +
  641 + // find the lowest high
  642 + typename std::vector<Box>::const_iterator it = boxes.begin();
  643 + typedef typename coordinate_type<Box>::type coordinate_type;
  644 + coordinate_type lowest_high = geometry::get<max_corner, D>(*it);
  645 + unsigned int lowest_high_index = 0;
  646 + unsigned int index = 1;
  647 + ++it;
  648 + for(; it != boxes.end(); ++it)
  649 + {
  650 + if (geometry::get<max_corner, D>(*it) < lowest_high)
  651 + {
  652 + lowest_high = geometry::get<max_corner, D>(*it);
  653 + lowest_high_index = index;
  654 + }
  655 + index++;
  656 + }
  657 +
  658 + // find the highest low
  659 + coordinate_type highest_low = 0;
  660 + unsigned int highest_low_index = 0;
  661 + if (lowest_high_index == 0)
  662 + {
  663 + highest_low = geometry::get<min_corner, D>(boxes[1]);
  664 + highest_low_index = 1;
  665 + }
  666 + else
  667 + {
  668 + highest_low = geometry::get<min_corner, D>(boxes[0]);
  669 + highest_low_index = 0;
  670 + }
  671 +
  672 + index = 0;
  673 + for (typename std::vector<Box>::const_iterator it = boxes.begin();
  674 + it != boxes.end(); ++it, index++)
  675 + {
  676 + if (geometry::get<min_corner, D>(*it) >= highest_low && index != lowest_high_index)
  677 + {
  678 + highest_low = geometry::get<min_corner, D>(*it);
  679 + highest_low_index = index;
  680 + }
  681 + }
  682 +
  683 + // find the lowest low
  684 + it = boxes.begin();
  685 + coordinate_type lowest_low = geometry::get<min_corner, D>(*it);
  686 + ++it;
  687 + for(; it != boxes.end(); ++it)
  688 + {
  689 + if (geometry::get<min_corner, D>(*it) < lowest_low)
  690 + {
  691 + lowest_low = geometry::get<min_corner, D>(*it);
  692 + }
  693 + }
  694 +
  695 + // find the highest high
  696 + it = boxes.begin();
  697 + coordinate_type highest_high = geometry::get<max_corner, D>(*it);
  698 + ++it;
  699 + for(; it != boxes.end(); ++it)
  700 + {
  701 + if (geometry::get<max_corner, D>(*it) > highest_high)
  702 + {
  703 + highest_high = geometry::get<max_corner, D>(*it);
  704 + }
  705 + }
  706 +
  707 + coordinate_type const width = highest_high - lowest_low;
  708 +
  709 + separation = (highest_low - lowest_high) / width;
  710 + first = highest_low_index;
  711 + second = lowest_high_index;
  712 + }
  713 +
  714 + /**
  715 + * \brief Choose one of the possible leaves to make an insertion
  716 + */
  717 + inline node_pointer choose_corresponding_leaf(Box const& e)
  718 + {
  719 + node_pointer node = m_root;
  720 +
  721 + // if the tree is empty add an initial leaf
  722 + if (m_root->elements() == 0)
  723 + {
  724 + leaf_pointer new_leaf(new rtree_leaf<Box, Value>(m_root));
  725 + m_root->add_leaf_node(Box (), new_leaf);
  726 +
  727 + return new_leaf;
  728 + }
  729 +
  730 + while (!node->is_leaf())
  731 + {
  732 + /// traverse node's map to see which node we should select
  733 + node = node->choose_node(e);
  734 + }
  735 + return node;
  736 + }
  737 +
  738 + /**
  739 + * \brief Choose the exact leaf where an insertion should be done
  740 + */
  741 + node_pointer choose_exact_leaf(Box const&e) const
  742 + {
  743 + // find possible leaves
  744 + typedef typename std::vector<node_pointer> node_type;
  745 + node_type nodes;
  746 + m_root->find_leaves(e, nodes);
  747 +
  748 + // refine the result
  749 + for (typename node_type::const_iterator it = nodes.begin(); it != nodes.end(); ++it)
  750 + {
  751 + typedef std::vector<std::pair<Box, Value> > leaves_type;
  752 + leaves_type leaves = (*it)->get_leaves();
  753 +
  754 + for (typename leaves_type::const_iterator itl = leaves.begin();
  755 + itl != leaves.end(); ++itl)
  756 + {
  757 +
  758 + if (itl->first.max_corner() == e.max_corner()
  759 + && itl->first.min_corner() == e.min_corner())
  760 + {
  761 + return *it;
  762 + }
  763 + }
  764 + }
  765 +
  766 + // TODO: mloskot - ggl exception
  767 + throw std::logic_error("Leaf not found");
  768 + }
  769 +};
  770 +
  771 +}}} // namespace boost::geometry::index
  772 +
  773 +#endif // BOOST_GEOMETRY_EXTENSIONS_INDEX_RTREE_RTREE_HPP
  774 +
253 boost/geometry/extensions/index/rtree/rtree_leaf.hpp
... ... @@ -0,0 +1,253 @@
  1 +// Boost.Geometry (aka GGL, Generic Geometry Library)
  2 +
  3 +// Boost.SpatialIndex - rtree leaf implementation
  4 +//
  5 +// Copyright 2008 Federico J. Fernandez.
  6 +// Use, modification and distribution is subject to the Boost Software License,
  7 +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  8 +// http://www.boost.org/LICENSE_1_0.txt)
  9 +
  10 +#ifndef BOOST_GEOMETRY_EXTENSIONS_INDEX_RTREE_RTREE_LEAF_HPP
  11 +#define BOOST_GEOMETRY_EXTENSIONS_INDEX_RTREE_RTREE_LEAF_HPP
  12 +
  13 +#include <deque>
  14 +#include <iostream> // TODO: Remove if print() is removed
  15 +#include <stdexcept>
  16 +#include <utility>
  17 +#include <vector>
  18 +
  19 +#include <boost/shared_ptr.hpp>
  20 +
  21 +#include <boost/geometry/algorithms/area.hpp>
  22 +#include <boost/geometry/algorithms/assign.hpp>
  23 +#include <boost/geometry/algorithms/expand.hpp>
  24 +
  25 +#include <boost/geometry/extensions/index/rtree/rtree_node.hpp>
  26 +
  27 +namespace boost { namespace geometry { namespace index
  28 +{
  29 +
  30 +template <typename Box, typename Value >
  31 +class rtree_leaf : public rtree_node<Box, Value>
  32 +{
  33 +public:
  34 +
  35 + /// container type for the leaves
  36 + typedef boost::shared_ptr<rtree_node<Box, Value> > node_pointer;
  37 + typedef std::vector<std::pair<Box, Value> > leaf_map;
  38 +
  39 + /**
  40 + * \brief Creates an empty leaf
  41 + */
  42 + inline rtree_leaf()
  43 + {
  44 + }
  45 +
  46 + /**
  47 + * \brief Creates a new leaf, with 'parent' as parent
  48 + */
  49 + inline rtree_leaf(node_pointer const& parent)
  50 + : rtree_node<Box, Value> (parent, 0)
  51 + {
  52 + }
  53 +
  54 + /**
  55 + * \brief Search for elements in 'box' in the Rtree. Add them to 'result'.
  56 + * If exact_match is true only return the elements having as
  57 + * key the 'box'. Otherwise return everything inside 'box'.
  58 + */
  59 + virtual void find(Box const& box, std::deque<Value>& result, bool const exact_match)
  60 + {
  61 + for (typename leaf_map::const_iterator it = m_nodes.begin();
  62 + it != m_nodes.end(); ++it)
  63 + {
  64 + if (exact_match)
  65 + {
  66 + if (geometry::equals(it->first, box))
  67 + {
  68 + result.push_back(it->second);
  69 + }
  70 + }
  71 + else
  72 + {
  73 + if (is_overlapping(it->first, box))
  74 + {
  75 + result.push_back(it->second);
  76 + }
  77 + }
  78 + }
  79 + }
  80 +
  81 + /**
  82 + * \brief Compute bounding box for this leaf
  83 + */
  84 + virtual Box compute_box() const
  85 + {
  86 + if (m_nodes.empty())
  87 + {
  88 + return Box ();
  89 + }
  90 +
  91 + Box r;
  92 + geometry::assign_inverse(r);
  93 + for(typename leaf_map::const_iterator it = m_nodes.begin(); it != m_nodes.end(); ++it)
  94 + {
  95 + geometry::expand(r, it->first);
  96 + }
  97 + return r;
  98 + }
  99 +
  100 + /**
  101 + * \brief True if we are a leaf
  102 + */
  103 + virtual bool is_leaf() const
  104 + {
  105 + return true;
  106 + }
  107 +
  108 + /**
  109 + * \brief Number of elements in the tree
  110 + */
  111 + virtual unsigned int elements() const
  112 + {
  113 + return m_nodes.size();
  114 + }
  115 +
  116 + /**
  117 + * \brief Insert a new element, with key 'box' and value 'v'
  118 + */
  119 + virtual void insert(Box const& box, Value const& v)
  120 + {
  121 + m_nodes.push_back(std::make_pair(box, v));
  122 + }
  123 +
  124 + /**
  125 + * \brief Proyect leaves of this node.
  126 + */
  127 + virtual std::vector< std::pair<Box, Value> > get_leaves() const
  128 + {
  129 + return m_nodes;
  130 + }
  131 +
  132 + /**
  133 + * \brief Add a new child (node) to this node
  134 + */
  135 + virtual void add_node(Box const&, node_pointer const&)
  136 + {
  137 + // TODO: mloskot - define & use GGL exception
  138 + throw std::logic_error("Can't add node to leaf node.");
  139 + }
  140 +
  141 + /**
  142 + * \brief Add a new leaf to this node
  143 + */
  144 + virtual void add_value(Box const& box, Value const& v)
  145 + {
  146 + m_nodes.push_back(std::make_pair(box, v));
  147 + }
  148 +
  149 +
  150 + /**
  151 + * \brief Proyect value in position 'index' in the nodes container
  152 + */
  153 + virtual Value get_value(unsigned int index) const
  154 + {
  155 + return m_nodes[index].second;
  156 + }
  157 +
  158 + /**
  159 + * \brief Box projector for leaf
  160 + */
  161 + virtual Box get_box(unsigned int index) const
  162 + {
  163 + return m_nodes[index].first;
  164 + }
  165 +
  166 + /**
  167 + * \brief Remove value with key 'box' in this leaf
  168 + */
  169 + virtual void remove(Box const& box)
  170 + {
  171 +
  172 + for (typename leaf_map::iterator it = m_nodes.begin();
  173 + it != m_nodes.end(); ++it)
  174 + {
  175 + if (geometry::equals(it->first, box))
  176 + {
  177 + m_nodes.erase(it);
  178 + return;
  179 + }
  180 + }
  181 +
  182 + // TODO: mloskot - use GGL exception
  183 + throw std::logic_error("Node not found.");
  184 + }
  185 +
  186 + /**
  187 + * \brief Remove value in this leaf
  188 + */
  189 + virtual void remove(Value const& v)
  190 + {
  191 + for (typename leaf_map::iterator it = m_nodes.begin();
  192 + it != m_nodes.end(); ++it)
  193 + {
  194 + if (it->second == v)
  195 + {
  196 + m_nodes.erase(it);
  197 + return;
  198 + }
  199 + }
  200 +
  201 + // TODO: mloskot - use GGL exception
  202 + throw std::logic_error("Node not found.");
  203 + }
  204 +
  205 + /**
  206 + * \brief Proyect boxes from this node
  207 + */
  208 + virtual std::vector<Box> get_boxes() const
  209 + {
  210 + std::vector<Box> result;
  211 + for (typename leaf_map::const_iterator it = m_nodes.begin();
  212 + it != m_nodes.end(); ++it)
  213 + {
  214 + result.push_back(it->first);
  215 + }
  216 +
  217 + return result;
  218 + }
  219 +
  220 + /**
  221 + * \brief Print leaf (mainly for debug)
  222 + */
  223 + virtual void print() const
  224 + {
  225 + std::cerr << "\t" << " --> Leaf --------" << std::endl;
  226 + std::cerr << "\t" << " Size: " << m_nodes.size() << std::endl;
  227 + for (typename leaf_map::const_iterator it = m_nodes.begin();
  228 + it != m_nodes.end(); ++it)
  229 + {
  230 + std::cerr << "\t" << " | ";
  231 + std::cerr << "( " << geometry::get<min_corner, 0>
  232 + (it->first) << " , " << geometry::get<min_corner, 1>
  233 + (it->first) << " ) x ";
  234 + std::cerr << "( " << geometry::get<max_corner, 0>
  235 + (it->first) << " , " << geometry::get<max_corner, 1>
  236 + (it->first) << " )";
  237 + std::cerr << " -> ";
  238 + std::cerr << it->second;
  239 + std::cerr << " | " << std::endl;;
  240 + }
  241 + std::cerr << "\t" << " --< Leaf --------" << std::endl;
  242 + }
  243 +
  244 +private:
  245 +
  246 + /// leaves of this node
  247 + leaf_map m_nodes;
  248 +};
  249 +
  250 +}}} // namespace boost::geometry::index
  251 +
  252 +#endif // BOOST_GEOMETRY_EXTENSIONS_INDEX_RTREE_RTREE_LEAF_HPP
  253 +
493 boost/geometry/extensions/index/rtree/rtree_node.hpp
... ... @@ -0,0 +1,493 @@
  1 +// Boost.Geometry (aka GGL, Generic Geometry Library)
  2 +
  3 +// Boost.SpatialIndex - rtree node implementation
  4 +//
  5 +// Copyright 2008 Federico J. Fernandez.
  6 +// Use, modification and distribution is subject to the Boost Software License,
  7 +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  8 +// http://www.boost.org/LICENSE_1_0.txt)
  9 +
  10 +#ifndef BOOST_GEOMETRY_EXTENSIONS_INDEX_RTREE_RTREE_NODE_HPP
  11 +#define BOOST_GEOMETRY_EXTENSIONS_INDEX_RTREE_RTREE_NODE_HPP
  12 +
  13 +#include <deque>
  14 +#include <iostream> // TODO: Remove if print() is removed
  15 +#include <stdexcept>
  16 +#include <utility>