61 #ifndef OPENMS_DATASTRUCTURES_KDTREE_H
62 #define OPENMS_DATASTRUCTURES_KDTREE_H
64 #ifdef KDTREE_DEFINE_OSTREAM_OPERATORS
102 template <
typename _Val>
116 #ifdef KDTREE_DEFINE_OSTREAM_OPERATORS
118 template <
typename Char,
typename Traits>
120 std::basic_ostream<Char, Traits>&
121 operator<<(
typename std::basic_ostream<Char, Traits>& out,
125 out <<
" parent: " << node._M_parent;
126 out <<
"; left: " << node._M_left;
127 out <<
"; right: " << node._M_right;
131 template <
typename Char,
typename Traits>
133 std::basic_ostream<Char, Traits>&
134 operator<<(
typename std::basic_ostream<Char, Traits>& out,
135 _Node<_Val>
const& node)
138 out <<
' ' << node._M_value;
139 out <<
"; parent: " << node._M_parent;
140 out <<
"; left: " << node._M_left;
141 out <<
"; right: " << node._M_right;
148 template <
typename _Val,
typename _Acc,
typename _Cmp>
173 template <
typename _ValA,
typename _ValB,
typename _Cmp,
178 const _Cmp& __cmp,
const _Acc& __acc,
179 const _ValA& __a,
const _ValB& __b)
181 return __cmp(__acc(__a, __dim), __acc(__b, __dim));
189 template <
typename _ValA,
typename _ValB,
typename _Dist,
192 typename _Dist::distance_type
194 const _Dist& __dist,
const _Acc& __acc,
195 const _ValA& __a,
const _ValB& __b)
197 return __dist(__acc(__a, __dim), __acc(__b, __dim));
206 template <
typename _ValA,
typename _ValB,
typename _Dist,
209 typename _Dist::distance_type
211 const _Dist& __dist,
const _Acc& __acc,
212 const _ValA& __a,
const _ValB& __b)
214 typename _Dist::distance_type d = 0;
215 for (
size_t i=0; i<__dim; ++i)
216 d += __dist(__acc(__a, i), __acc(__b, i));
225 template <
typename _Val,
typename _Cmp,
typename _Acc,
typename NodeType>
229 const _Cmp& __cmp,
const _Acc& __acc,
230 const _Val& __val,
const NodeType* __node)
233 return static_cast<const NodeType *>(__node->_M_left);
234 return static_cast<const NodeType *>(__node->_M_right);
245 template <
class SearchVal,
246 typename NodeType,
typename _Cmp,
247 typename _Acc,
typename _Dist,
250 std::pair<
const NodeType*,
251 std::pair<size_t, typename _Dist::distance_type> >
253 const NodeType* __node,
const _Node_base* __end,
254 const NodeType* __best,
typename _Dist::distance_type __max,
255 const _Cmp& __cmp,
const _Acc& __acc,
const _Dist& __dist,
258 typedef const NodeType* NodePtr;
259 NodePtr pcur = __node;
260 NodePtr cur =
_S_node_descend(__dim % __k, __cmp, __acc, __val, __node);
261 size_t cur_dim = __dim+1;
265 if (__p(cur->_M_value))
267 typename _Dist::distance_type d = 0;
268 for (
size_t i=0; i != __k; ++i)
292 NodePtr pprobe = probe;
295 size_t probe_dim = cur_dim;
296 if (
_S_node_compare(probe_dim % __k, __cmp, __acc, __val, probe->_M_value))
297 near_node = static_cast<NodePtr>(probe->_M_right);
299 near_node = static_cast<NodePtr>(probe->_M_left);
302 && (std::sqrt(
_S_node_distance(probe_dim % __k, __dist, __acc, __val, probe->_M_value)) <= __max))
311 if (
_S_node_compare(probe_dim % __k, __cmp, __acc, __val, probe->_M_value))
313 near_node = static_cast<NodePtr>(probe->_M_left);
314 far_node = static_cast<NodePtr>(probe->_M_right);
318 near_node = static_cast<NodePtr>(probe->_M_right);
319 far_node = static_cast<NodePtr>(probe->_M_left);
321 if (pprobe == probe->_M_parent)
323 if (__p(probe->_M_value))
325 typename _Dist::distance_type d = 0;
326 for (
size_t i=0; i < __k; ++i)
344 std::sqrt(
_S_node_distance(probe_dim % __k, __dist, __acc, __val, probe->_M_value)) <= __max)
351 probe = static_cast<NodePtr>(probe->_M_parent);
357 if (pprobe == near_node && far_node
359 && std::sqrt(
_S_node_distance(probe_dim % __k, __dist, __acc, __val, probe->_M_value)) <= __max)
368 probe = static_cast<NodePtr>(probe->_M_parent);
374 cur = static_cast<NodePtr>(cur->_M_parent);
381 if (pcur == cur->_M_left)
382 near_node = static_cast<NodePtr>(cur->_M_right);
384 near_node = static_cast<NodePtr>(cur->_M_left);
387 && (std::sqrt(
_S_node_distance(cur_dim % __k, __dist, __acc, __val, cur->_M_value)) <= __max))
394 return std::pair<NodePtr,
395 std::pair<size_t, typename _Dist::distance_type> >
396 (__best, std::pair<size_t, typename _Dist::distance_type>
403 #endif // include guard
424 #ifndef INCLUDE_KDTREE_ACCESSOR_HPP
425 #define INCLUDE_KDTREE_ACCESSOR_HPP
431 template <
typename _Val>
443 template <
typename _Tp>
449 template <
typename _Tp,
typename _Dist>
462 template <
typename _Tp,
typename _Dist>
492 #endif // include guard
512 #ifndef INCLUDE_KDTREE_ALLOCATOR_HPP
513 #define INCLUDE_KDTREE_ALLOCATOR_HPP
520 template <
typename _Tp,
typename _Alloc>
574 new (__p)
_Node_(__V, __PARENT, __LEFT, __RIGHT);
586 #endif // include guard
606 #ifndef INCLUDE_KDTREE_ITERATOR_HPP
607 #define INCLUDE_KDTREE_ITERATOR_HPP
613 template <
typename _Val,
typename _Ref,
typename _Ptr>
616 template<
typename _Val,
typename _Ref,
typename _Ptr>
621 template<
typename _Val>
626 template<
typename _Val>
631 template<
typename _Val,
typename _Ref,
typename _Ptr>
636 template<
typename _Val>
641 template<
typename _Val>
668 while (__p &&
_M_node == __p->_M_right)
689 while (x->_M_right) x = x->_M_right;
695 while (__p &&
_M_node == __p->_M_left)
706 template <
size_t const __K,
typename _Val,
typename _Acc,
707 typename _Dist,
typename _Cmp,
typename _Alloc>
711 template <
typename _Val,
typename _Ref,
typename _Ptr>
804 template<
typename _Val,
typename _Ref,
typename _Ptr>
810 template<
typename _Val>
816 template<
typename _Val>
822 template<
typename _Val,
typename _Ref,
typename _Ptr>
828 template<
typename _Val>
834 template<
typename _Val>
842 #endif // include guard
862 #ifndef INCLUDE_KDTREE_REGION_HPP
863 #define INCLUDE_KDTREE_REGION_HPP
870 template <
size_t const __K,
typename _Val,
typename _SubVal,
871 typename _Acc,
typename _Cmp>
882 _Region(_Acc
const& __acc=_Acc(),
const _Cmp& __cmp=_Cmp())
885 template <
typename Val>
887 _Acc
const& __acc=_Acc(),
const _Cmp& __cmp=_Cmp())
890 for (
size_t __i = 0; __i != __K; ++__i)
896 template <
typename Val>
898 _Acc
const& __acc=_Acc(),
const _Cmp& __cmp=_Cmp())
901 for (
size_t __i = 0; __i != __K; ++__i)
911 for (
size_t __i = 0; __i != __K; ++__i)
929 for (
size_t __i = 0; __i != __K; ++__i)
941 for (
size_t __i = 0; __i != __K; ++__i)
971 #endif // include guard
1030 #ifndef INCLUDE_KDTREE_KDTREE_HPP
1031 #define INCLUDE_KDTREE_KDTREE_HPP
1040 #define KDTREE_VERSION 701
1045 #define KDTREE_LIB_VERSION "0_7_1"
1050 #ifdef KDTREE_CHECK_PERFORMANCE_COUNTERS
1053 #include <algorithm>
1054 #include <functional>
1056 #ifdef KDTREE_DEFINE_OSTREAM_OPERATORS
1068 #ifdef KDTREE_CHECK_PERFORMANCE
1069 unsigned long long num_dist_calcs = 0;
1072 template <
size_t const __K,
typename _Val,
1073 typename _Acc = _Bracket_accessor<_Val>,
1074 typename _Dist = squared_difference<
typename _Acc::result_type,
1075 typename _Acc::result_type>,
1076 typename _Cmp = std::less<typename _Acc::result_type>,
1077 typename _Alloc = std::allocator<_Node<_Val> > >
1103 KDTree(_Acc
const& __acc = _Acc(), _Dist
const& __dist = _Dist(),
1105 :
_Base(__a), _M_header(),
1106 _M_count(0), _M_acc(__acc), _M_cmp(__cmp), _M_dist(__dist)
1108 _M_empty_initialise();
1112 :
_Base(__x.get_allocator()), _M_header(), _M_count(0),
1113 _M_acc(__x._M_acc), _M_cmp(__x._M_cmp), _M_dist(__x._M_dist)
1115 _M_empty_initialise();
1124 std::vector<value_type> temp;
1125 temp.reserve(__x.size());
1126 std::copy(__x.begin(),__x.end(),std::back_inserter(temp));
1127 _M_optimise(temp.begin(), temp.end(), 0);
1130 template<
typename _InputIterator>
1131 KDTree(_InputIterator __first, _InputIterator __last,
1132 _Acc
const& acc = _Acc(), _Dist
const& __dist = _Dist(),
1134 :
_Base(__a), _M_header(), _M_count(0),
1135 _M_acc(acc), _M_cmp(__cmp), _M_dist(__dist)
1137 _M_empty_initialise();
1146 std::vector<value_type> temp;
1147 temp.reserve(std::distance(__first,__last));
1148 std::copy(__first,__last,std::back_inserter(temp));
1149 _M_optimise(temp.begin(), temp.end(), 0);
1172 _M_optimise(writable_vector.begin(), writable_vector.end(), 0);
1182 _M_acc = __x._M_acc;
1183 _M_dist = __x._M_dist;
1184 _M_cmp = __x._M_cmp;
1193 std::vector<value_type> temp;
1194 temp.reserve(__x.size());
1195 std::copy(__x.begin(),__x.end(),std::back_inserter(temp));
1196 efficient_replace_and_optimise(temp);
1209 return _Base::get_allocator();
1227 return this->size() == 0;
1233 _M_erase_subtree(_M_get_root());
1234 _M_set_leftmost(&_M_header);
1235 _M_set_rightmost(&_M_header);
1236 _M_set_root(
nullptr);
1290 return this->insert(__V);
1298 _Link_type __n = _M_new_node(__V, &_M_header);
1301 _M_set_leftmost(__n);
1302 _M_set_rightmost(__n);
1305 return _M_insert(_M_get_root(), __V, 0);
1308 template <
class _InputIterator>
1309 void insert(_InputIterator __first, _InputIterator __last) {
1310 for (; __first != __last; ++__first)
1311 this->insert(*__first);
1317 for (; __n > 0; --__n)
1318 this->insert(__pos, __x);
1321 template<
typename _InputIterator>
1324 for (; __first != __last; ++__first)
1325 this->insert(__pos, *__first);
1346 this->erase(this->find_exact(__V));
1353 assert(__IT != this->end());
1357 while ((n = _S_parent(n)) != &_M_header)
1359 _M_erase( const_cast<_Link_type>(target), level );
1360 _M_delete_node( const_cast<_Link_type>(target) );
1384 template <
class SearchVal>
1388 if (!_M_get_root())
return this->end();
1389 return _M_find(_M_get_root(), __V, 0);
1406 template <
class SearchVal>
1410 if (!_M_get_root())
return this->end();
1411 return _M_find_exact(_M_get_root(), __V, 0);
1418 if (!_M_get_root())
return 0;
1419 _Region_ __region(__V, __R, _M_acc, _M_cmp);
1420 return this->count_within_range(__region);
1426 if (!_M_get_root())
return 0;
1429 return _M_count_within_range(_M_get_root(),
1430 __REGION, __bounds, 0);
1434 template <
typename SearchVal,
class Visitor>
1438 if (!_M_get_root())
return visitor;
1440 return this->visit_within_range(region, visitor);
1443 template <
class Visitor>
1450 return _M_visit_within_range(visitor, _M_get_root(), REGION, bounds, 0);
1465 template <
typename SearchVal,
typename _OutputIterator>
1468 _OutputIterator out)
const
1470 if (!_M_get_root())
return out;
1471 _Region_ region(val, range, _M_acc, _M_cmp);
1472 return this->find_within_range(region, out);
1475 template <
typename _OutputIterator>
1478 _OutputIterator out)
const
1483 out = _M_find_within_range(out, _M_get_root(),
1489 template <
class SearchVal>
1490 std::pair<const_iterator, distance_type>
1495 std::pair<const _Node<_Val>*,
1496 std::pair<size_type, typename _Acc::result_type> >
1498 _M_get_root(), &_M_header, _M_get_root(),
1500 (__K, _M_dist, _M_acc, _M_get_root()->_M_value, __val)),
1501 _M_cmp, _M_acc, _M_dist,
1503 return std::pair<const_iterator, distance_type>
1504 (best.first, best.second.second);
1506 return std::pair<const_iterator, distance_type>(end(), 0);
1509 template <
class SearchVal>
1510 std::pair<const_iterator, distance_type>
1515 bool root_is_candidate =
false;
1519 (__K, _M_dist, _M_acc, _M_get_root()->_M_value, __val));
1520 if (root_dist <= __max)
1522 root_is_candidate =
true;
1526 std::pair<const _Node<_Val>*,
1527 std::pair<size_type, typename _Acc::result_type> >
1529 node, __max, _M_cmp, _M_acc, _M_dist,
1532 if (root_is_candidate || best.first != _M_get_root())
1533 return std::pair<const_iterator, distance_type>
1534 (best.first, best.second.second);
1536 return std::pair<const_iterator, distance_type>(end(), __max);
1539 template <
class SearchVal,
class _Predicate>
1540 std::pair<const_iterator, distance_type>
1542 _Predicate __p)
const
1546 bool root_is_candidate =
false;
1548 if (__p(_M_get_root()->_M_value))
1552 (__K, _M_dist, _M_acc, _M_get_root()->_M_value, __val));
1553 if (root_dist <= __max)
1555 root_is_candidate =
true;
1560 std::pair<const _Node<_Val>*,
1561 std::pair<size_type, typename _Acc::result_type> >
1563 node, __max, _M_cmp, _M_acc, _M_dist, __p);
1565 if (root_is_candidate || best.first != _M_get_root())
1566 return std::pair<const_iterator, distance_type>
1567 (best.first, best.second.second);
1569 return std::pair<const_iterator, distance_type>(end(), __max);
1575 std::vector<value_type> __v(this->begin(),this->end());
1577 _M_optimise(__v.begin(), __v.end(), 0);
1588 _M_check_node(_M_get_root(),0);
1605 _M_check_children(_S_left(child),parent,level,to_the_left);
1606 _M_check_children(_S_right(child),parent,level,to_the_left);
1616 _M_check_children( _S_left(node), node, level,
true );
1618 _M_check_children( _S_right(node), node, level,
false );
1620 _M_check_node( _S_left(node), level+1 );
1621 _M_check_node( _S_right(node), level+1 );
1627 _M_set_leftmost(&_M_header);
1628 _M_set_rightmost(&_M_header);
1629 _M_header._M_parent =
nullptr;
1630 _M_set_root(
nullptr);
1636 _S_set_left(__N, _M_new_node(__V)); ++_M_count;
1637 _S_set_parent( _S_left(__N), __N );
1638 if (__N == _M_get_leftmost())
1639 _M_set_leftmost( _S_left(__N) );
1646 _S_set_right(__N, _M_new_node(__V)); ++_M_count;
1647 _S_set_parent( _S_right(__N), __N );
1648 if (__N == _M_get_rightmost())
1649 _M_set_rightmost( _S_right(__N) );
1660 return _M_insert_left(__N, __V);
1661 return _M_insert(_S_left(__N), __V, __L+1);
1665 if (!_S_right(__N) || __N == _M_get_rightmost())
1666 return _M_insert_right(__N, __V);
1667 return _M_insert(_S_right(__N), __V, __L+1);
1675 _Link_type step_dad = _M_get_erase_replacement(dead_dad, level);
1678 if (dead_dad == _M_get_root())
1679 _M_set_root(step_dad);
1680 else if (_S_left(_S_parent(dead_dad)) == dead_dad)
1681 _S_set_left(_S_parent(dead_dad), step_dad);
1683 _S_set_right(_S_parent(dead_dad), step_dad);
1688 if (dead_dad == _M_get_leftmost())
1689 _M_set_leftmost( (step_dad ? step_dad : _S_parent(dead_dad)) );
1690 if (dead_dad == _M_get_rightmost())
1691 _M_set_rightmost( (step_dad ? step_dad : _S_parent(dead_dad)) );
1696 _S_set_parent(step_dad, _S_parent(dead_dad));
1699 if (_S_left(dead_dad))
1700 _S_set_parent(_S_left(dead_dad), step_dad);
1701 if (_S_right(dead_dad))
1702 _S_set_parent(_S_right(dead_dad), step_dad);
1705 _S_set_left(step_dad, _S_left(dead_dad));
1706 _S_set_right(step_dad, _S_right(dead_dad));
1718 if (_S_is_leaf(node))
1721 std::pair<_Link_type,size_type> candidate;
1724 candidate = _M_get_j_min( std::pair<_Link_type,size_type>(_S_right(node),level), level+1);
1726 else if ((!_S_right(node)))
1727 candidate = _M_get_j_max( std::pair<_Link_type,size_type>(_S_left(node),level), level+1);
1740 if (compare(_S_right(node)->_M_value, _S_left(node)->_M_value))
1742 candidate = _M_get_j_min(std::pair<_Link_type,size_type>(_S_right(node),level), level+1);
1744 candidate = _M_get_j_max( std::pair<_Link_type,size_type>(_S_left(node),level), level+1);
1750 _Link_type parent = _S_parent(candidate.first);
1751 if (_S_left(parent) == candidate.first)
1752 _S_set_left(parent, _M_erase(candidate.first, candidate.second));
1754 _S_set_right(parent, _M_erase(candidate.first, candidate.second));
1756 return candidate.first;
1762 std::pair<_Link_type,size_type>
1765 typedef std::pair<_Link_type,size_type> Result;
1766 if (_S_is_leaf(node.first))
1767 return Result(node.first,level);
1770 Result candidate = node;
1771 if (_S_left(node.first))
1773 Result left = _M_get_j_min(Result(_S_left(node.first), node.second), level+1);
1774 if (compare(left.first->_M_value, candidate.first->_M_value))
1777 if (_S_right(node.first))
1779 Result right = _M_get_j_min( Result(_S_right(node.first),node.second), level+1);
1780 if (compare(right.first->_M_value, candidate.first->_M_value))
1783 if (candidate.first == node.first)
1784 return Result(candidate.first,level);
1792 std::pair<_Link_type,size_type>
1795 typedef std::pair<_Link_type,size_type> Result;
1797 if (_S_is_leaf(node.first))
1798 return Result(node.first,level);
1801 Result candidate = node;
1802 if (_S_left(node.first))
1804 Result left = _M_get_j_max( Result(_S_left(node.first),node.second), level+1);
1805 if (compare(candidate.first->_M_value, left.first->_M_value))
1808 if (_S_right(node.first))
1810 Result right = _M_get_j_max(Result(_S_right(node.first),node.second), level+1);
1811 if (compare(candidate.first->_M_value, right.first->_M_value))
1815 if (candidate.first == node.first)
1816 return Result(candidate.first,level);
1827 _M_erase_subtree(_S_right(__n));
1829 _M_delete_node(__n);
1844 if (!compare(node->
_M_value,value))
1847 if (_M_matches_node(node, value, level))
1850 found = _M_find(_S_left(node), value, level+1);
1852 if ( _S_right(node) && found == this->end() && !compare(value,node->
_M_value))
1853 found = _M_find(_S_right(node), value, level+1);
1867 if (!compare(node->
_M_value,value))
1873 found = _M_find_exact(_S_left(node), value, level+1);
1877 if ( _S_right(node) && found == this->end() && !compare(value,node->
_M_value))
1878 found = _M_find_exact(_S_right(node), value, level+1);
1895 while ((__i = (__i + 1) % __K) != __L % __K)
1896 if (!_M_matches_node_in_d(__N, __V, __i))
return false;
1904 return _M_matches_node_in_d(__N, __V, __L)
1905 && _M_matches_node_in_other_ds(__N, __V, __L);
1914 if (__REGION.
encloses(_S_value(__N)))
1923 count += _M_count_within_range(_S_left(__N),
1924 __REGION, __bounds, __L+1);
1931 count += _M_count_within_range(_S_right(__N),
1932 __REGION, __bounds, __L+1);
1939 template <
class Visitor>
1948 visitor(_S_value(N));
1955 visitor = _M_visit_within_range(visitor, _S_left(N),
1956 REGION, bounds, L+1);
1963 visitor = _M_visit_within_range(visitor, _S_right(N),
1964 REGION, bounds, L+1);
1972 template <
typename _OutputIterator>
1979 if (__REGION.
encloses(_S_value(__N)))
1981 *out++ = _S_value(__N);
1988 out = _M_find_within_range(out, _S_left(__N),
1989 __REGION, __bounds, __L+1);
1996 out = _M_find_within_range(out, _S_right(__N),
1997 __REGION, __bounds, __L+1);
2004 template <
typename _Iter>
2009 if (__A == __B)
return;
2011 _Iter __m = __A + (__B - __A) / 2;
2012 std::nth_element(__A, __m, __B, compare);
2014 if (__m != __A) _M_optimise(__A, __m, __L+1);
2015 if (++__m != __B) _M_optimise(__m, __B, __L+1);
2021 return const_cast<_Link_const_type>(_M_root);
2038 return static_cast<_Link_type>(_M_header._M_left);
2044 _M_header._M_left = a;
2050 return static_cast<_Link_type>( _M_header._M_right );
2056 _M_header._M_right = a;
2062 return static_cast<_Link_type>( N->
_M_parent );
2065 static _Link_const_type
2068 return static_cast<_Link_const_type>( N->_M_parent );
2086 return static_cast<_Link_type>( N->
_M_left );
2089 static _Link_const_type
2092 return static_cast<_Link_const_type>( N->_M_left );
2104 return static_cast<_Link_type>( N->
_M_right );
2107 static _Link_const_type
2110 return static_cast<_Link_const_type>( N->_M_right );
2116 return !_S_left(N) && !_S_right(N);
2119 static const_reference
2125 static const_reference
2128 return static_cast<_Link_const_type>(N)->_M_value;
2131 static _Link_const_type
2137 static _Link_const_type
2150 typename _Base::NoLeakAlloc noleak(
this);
2152 _Base::_M_construct_node(new_node, __V, __PARENT, __LEFT, __RIGHT);
2153 noleak.disconnect();
2170 _Base::_M_destroy_node(__p);
2171 _Base::_M_deallocate_node(__p);
2181 #ifdef KDTREE_DEFINE_OSTREAM_OPERATORS
2182 friend std::ostream&
2186 o <<
"meta node: " << tree._M_header << std::endl;
2187 o <<
"root node: " << tree._M_root << std::endl;
2190 return o <<
"[empty " << __K <<
"d-tree " << &tree <<
"]";
2192 o <<
"nodes total: " << tree.size() << std::endl;
2193 o <<
"dimensions: " << __K << std::endl;
2196 typedef typename _Tree::_Link_type
_Link_type;
2198 std::stack<_Link_const_type> s;
2199 s.push(tree._M_get_root());
2205 o << *n << std::endl;
2206 if (_Tree::_S_left(n)) s.push(_Tree::_S_left(n));
2207 if (_Tree::_S_right(n)) s.push(_Tree::_S_right(n));
2219 #endif // include guard