diff --git a/STL_Extension/include/CGAL/Compact_container.h b/STL_Extension/include/CGAL/Compact_container.h index a5004a8674f..86539276015 100644 --- a/STL_Extension/include/CGAL/Compact_container.h +++ b/STL_Extension/include/CGAL/Compact_container.h @@ -35,6 +35,7 @@ #include #include #include +#include #include @@ -87,24 +88,6 @@ namespace CGAL { -#define CGAL_GENERATE_MEMBER_DETECTOR(X) \ -template class has_##X { \ - struct Fallback { int X; }; \ - struct Derived : T, Fallback { }; \ - \ - template struct Check; \ - \ - typedef char ArrayOfOne[1]; \ - typedef char ArrayOfTwo[2]; \ - \ - template static ArrayOfOne & func( \ - Check *); \ - template static ArrayOfTwo & func(...); \ - public: \ - typedef has_##X type; \ - enum { value = sizeof(func(0)) == 2 }; \ -} // semicolon is after the macro call - #define CGAL_INIT_COMPACT_CONTAINER_BLOCK_SIZE 14 #define CGAL_INCREMENT_COMPACT_CONTAINER_BLOCK_SIZE 16 diff --git a/STL_Extension/include/CGAL/Concurrent_compact_container.h b/STL_Extension/include/CGAL/Concurrent_compact_container.h index 6f85eb09c00..d7cfaa6d222 100644 --- a/STL_Extension/include/CGAL/Concurrent_compact_container.h +++ b/STL_Extension/include/CGAL/Concurrent_compact_container.h @@ -33,6 +33,7 @@ #include #include #include +#include #include @@ -40,24 +41,6 @@ namespace CGAL { -#define CGAL_GENERATE_MEMBER_DETECTOR(X) \ -template class has_##X { \ - struct Fallback { int X; }; \ - struct Derived : T, Fallback { }; \ - \ - template struct Check; \ - \ - typedef char ArrayOfOne[1]; \ - typedef char ArrayOfTwo[2]; \ - \ - template static ArrayOfOne & func( \ - Check *); \ - template static ArrayOfTwo & func(...); \ - public: \ - typedef has_##X type; \ - enum { value = sizeof(func(0)) == 2 }; \ -} // semicolon is after the macro call - #define CGAL_INIT_CONCURRENT_COMPACT_CONTAINER_BLOCK_SIZE 14 #define CGAL_INCREMENT_CONCURRENT_COMPACT_CONTAINER_BLOCK_SIZE 16 diff --git a/STL_Extension/include/CGAL/Has_member.h b/STL_Extension/include/CGAL/Has_member.h new file mode 100644 index 00000000000..79cf5a31b35 --- /dev/null +++ b/STL_Extension/include/CGAL/Has_member.h @@ -0,0 +1,46 @@ +// Copyright (c) 2017 Inria (France) +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public License as +// published by the Free Software Foundation; either version 3 of the License, +// or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// +// +// Author(s) : Clement Jamin + +#ifndef CGAL_HAS_MEMBER_H +#define CGAL_HAS_MEMBER_H + +// Macro used to check if a type T has a member named `X` +// It generates a class has_X where has_X::value is a boolean +// See example in Concurrent_compact_container.h +// Note: this macro is similar to BOOST_TTI_HAS_MEMBER_FUNCTION +// which was introduced in Boost 1.54 (CGAL supports version 1.48 and +// later at this time - June 2017) +#define CGAL_GENERATE_MEMBER_DETECTOR(X) \ +template class has_##X { \ + struct Fallback { int X; }; \ + struct Derived : T, Fallback { }; \ + \ + template struct Check; \ + \ + typedef char ArrayOfOne[1]; \ + typedef char ArrayOfTwo[2]; \ + \ + template static ArrayOfOne & func( \ + Check *); \ + template static ArrayOfTwo & func(...); \ + public: \ + typedef has_##X type; \ + enum { value = sizeof(func(0)) == 2 }; \ +} // semicolon is after the macro call + +#endif // CGAL_HAS_MEMBER_H diff --git a/Spatial_searching/include/CGAL/Euclidean_distance.h b/Spatial_searching/include/CGAL/Euclidean_distance.h index 3178170962d..308f4a69680 100644 --- a/Spatial_searching/include/CGAL/Euclidean_distance.h +++ b/Spatial_searching/include/CGAL/Euclidean_distance.h @@ -17,6 +17,7 @@ // // // Author(s) : Hans Tangelder () +// Clement Jamin (clement.jamin.pro@gmail.com) #ifndef CGAL_EUCLIDEAN_DISTANCE_H @@ -60,61 +61,114 @@ namespace CGAL { // default constructor Euclidean_distance(const SearchTraits& traits_=SearchTraits()):traits(traits_) {} - - inline FT transformed_distance(const Query_item& q, const Point_d& p) const { - return transformed_distance(q,p, D()); + inline FT transformed_distance(const Query_item& q, const Point_d& p) const + { + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d p_begin = construct_it(p), p_end = construct_it(p, 0); + return transformed_distance_from_coordinates(q, p_begin, p_end); } - //Dynamic version for runtime dimension - inline FT transformed_distance(const Query_item& q, const Point_d& p, Dynamic_dimension_tag) const { - FT distance = FT(0); - typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); - typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q), - qe = construct_it(q,1), pit = construct_it(p); - for(; qit != qe; qit++, pit++){ - distance += ((*qit)-(*pit))*((*qit)-(*pit)); - } - return distance; + template + inline FT transformed_distance_from_coordinates(const Query_item& q, + Coord_iterator it_coord_begin, Coord_iterator it_coord_end) const + { + return transformed_distance_from_coordinates(q, it_coord_begin, it_coord_end, D()); } - //Generic version for DIM > 3 - template < int DIM > - inline FT transformed_distance(const Query_item& q, const Point_d& p, Dimension_tag) const { - FT distance = FT(0); - typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); - typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q), - qe = construct_it(q,1), pit = construct_it(p); - for(; qit != qe; qit++, pit++){ - distance += ((*qit)-(*pit))*((*qit)-(*pit)); + // Static dim = 2 loop unrolled + template + inline FT transformed_distance_from_coordinates(const Query_item& q, + Coord_iterator it_coord_begin, Coord_iterator /*unused*/, + Dimension_tag<2>) const + { + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q); + FT distance = square(*qit - *it_coord_begin); + qit++; it_coord_begin++; + distance += square(*qit - *it_coord_begin); + return distance; + } + + // Static dim = 3 loop unrolled + template + inline FT transformed_distance_from_coordinates(const Query_item& q, + Coord_iterator it_coord_begin, Coord_iterator /*unused*/, + Dimension_tag<3>) const + { + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q); + FT distance = square(*qit - *it_coord_begin); + qit++; it_coord_begin++; + distance += square(*qit - *it_coord_begin); + qit++; it_coord_begin++; + distance += square(*qit - *it_coord_begin); + return distance; + } + + // Other cases: static dim > 3 or dynamic dim + template + inline FT transformed_distance_from_coordinates(const Query_item& q, + Coord_iterator it_coord_begin, Coord_iterator /*unused*/, + Dim) const + { + FT distance = FT(0); + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q), qe = construct_it(q, 1); + for (; qit != qe; ++qit, ++it_coord_begin) + { + FT diff = (*qit) - (*it_coord_begin); + distance += diff*diff; + } + return distance; + } + + // During the computation, if the partially-computed distance `pcd` gets greater or equal + // to `stop_if_geq_to_this`, the computation is stopped and `pcd` is returned + template + inline FT interruptible_transformed_distance(const Query_item& q, + Coord_iterator it_coord_begin, Coord_iterator /*unused*/, + FT stop_if_geq_to_this) const + { + FT distance = FT(0); + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object(); + typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q), qe = construct_it(q, 1); + if (qe - qit >= 6) + { + // Every 4 coordinates, the current partially-computed distance + // is compared to stop_if_geq_to_this + // Note: the concept SearchTraits specifies that Cartesian_const_iterator_d + // must be a random-access iterator + typename SearchTraits::Cartesian_const_iterator_d qe_minus_5 = qe - 5; + for (;;) + { + FT diff = (*qit) - (*it_coord_begin); + distance += diff*diff; + ++qit; ++it_coord_begin; + diff = (*qit) - (*it_coord_begin); + distance += diff*diff; + ++qit; ++it_coord_begin; + diff = (*qit) - (*it_coord_begin); + distance += diff*diff; + ++qit; ++it_coord_begin; + diff = (*qit) - (*it_coord_begin); + distance += diff*diff; + ++qit, ++it_coord_begin; + + if (distance >= stop_if_geq_to_this) + return distance; + + if (qit >= qe_minus_5) + break; } - return distance; + } + for (; qit != qe; ++qit, ++it_coord_begin) + { + FT diff = (*qit) - (*it_coord_begin); + distance += diff*diff; + } + return distance; } - //DIM = 2 loop unrolled - inline FT transformed_distance(const Query_item& q, const Point_d& p, Dimension_tag<2> ) const { - typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); - typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),pit = construct_it(p); - FT distance = square(*qit - *pit); - qit++;pit++; - distance += square(*qit - *pit); - return distance; - } - - //DIM = 3 loop unrolled - inline FT transformed_distance(const Query_item& q, const Point_d& p, Dimension_tag<3> ) const { - typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object(); - typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),pit = construct_it(p); - FT distance = square(*qit - *pit); - qit++;pit++; - distance += square(*qit - *pit); - qit++;pit++; - distance += square(*qit - *pit); - return distance; - } - - - - inline FT min_distance_to_rectangle(const Query_item& q, const Kd_tree_rectangle& r) const { FT distance = FT(0); diff --git a/Spatial_searching/include/CGAL/Kd_tree.h b/Spatial_searching/include/CGAL/Kd_tree.h index a5f38aed069..3896706dcfa 100644 --- a/Spatial_searching/include/CGAL/Kd_tree.h +++ b/Spatial_searching/include/CGAL/Kd_tree.h @@ -16,7 +16,8 @@ // $Id$ // // Author(s) : Hans Tangelder (), -// : Waqar Khan +// : Waqar Khan , +// Clement Jamin (clement.jamin.pro@gmail.com) #ifndef CGAL_KD_TREE_H #define CGAL_KD_TREE_H @@ -44,7 +45,11 @@ namespace CGAL { //template , class UseExtendedNode = Tag_true > -template , class UseExtendedNode = Tag_true > +template < + class SearchTraits, + class Splitter_=Sliding_midpoint, + class UseExtendedNode = Tag_true, + class EnablePointsCache = Tag_false> class Kd_tree { public: @@ -54,11 +59,11 @@ public: typedef typename Splitter::Container Point_container; typedef typename SearchTraits::FT FT; - typedef Kd_tree_node Node; - typedef Kd_tree_leaf_node Leaf_node; - typedef Kd_tree_internal_node Internal_node; - typedef Kd_tree Tree; - typedef Kd_tree Self; + typedef Kd_tree_node Node; + typedef Kd_tree_leaf_node Leaf_node; + typedef Kd_tree_internal_node Internal_node; + typedef Kd_tree Tree; + typedef Kd_tree Self; typedef Node* Node_handle; typedef const Node* Node_const_handle; @@ -76,6 +81,8 @@ public: typedef typename internal::Get_dimension_tag::Dimension D; + typedef EnablePointsCache Enable_points_cache; + private: SearchTraits traits_; Splitter split; @@ -95,6 +102,10 @@ private: Kd_tree_rectangle* bbox; std::vector pts; + // Store a contiguous copy of the point coordinates + // for faster queries (reduce the number of cache misses) + std::vector points_cache; + // Instead of storing the points in arrays in the Kd_tree_node // we put all the data in a vector in the Kd_tree. // and we only store an iterator range in the Kd_tree_node. @@ -282,9 +293,18 @@ public: //Reorder vector for spatial locality std::vector ptstmp; ptstmp.resize(pts.size()); - for (std::size_t i = 0; i < pts.size(); ++i){ + for (std::size_t i = 0; i < pts.size(); ++i) ptstmp[i] = *data[i]; + + // Cache? + if (Enable_points_cache::value) + { + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits_.construct_cartesian_const_iterator_d_object(); + points_cache.reserve(dim * pts.size()); + for (std::size_t i = 0; i < pts.size(); ++i) + points_cache.insert(points_cache.end(), construct_it(ptstmp[i]), construct_it(ptstmp[i], 0)); } + for(std::size_t i = 0; i < leaf_nodes.size(); ++i){ std::ptrdiff_t tmp = leaf_nodes[i].begin() - pts.begin(); leaf_nodes[i].data = ptstmp.begin() + tmp; @@ -546,6 +566,13 @@ public: return *bbox; } + + typename std::vector::const_iterator + cache_begin() const + { + return points_cache.begin(); + } + const_iterator begin() const { diff --git a/Spatial_searching/include/CGAL/Kd_tree_node.h b/Spatial_searching/include/CGAL/Kd_tree_node.h index 453aee65b59..71fd3d7a8e6 100644 --- a/Spatial_searching/include/CGAL/Kd_tree_node.h +++ b/Spatial_searching/include/CGAL/Kd_tree_node.h @@ -31,27 +31,29 @@ namespace CGAL { - template + template class Kd_tree; - template < class TreeTraits, class Splitter, class UseExtendedNode > + template < class TreeTraits, class Splitter, class UseExtendedNode, class EnablePointsCache > class Kd_tree_node { - friend class Kd_tree; + friend class Kd_tree; - typedef typename Kd_tree::Node_handle Node_handle; - typedef typename Kd_tree::Node_const_handle Node_const_handle; - typedef typename Kd_tree::Internal_node_handle Internal_node_handle; - typedef typename Kd_tree::Internal_node_const_handle Internal_node_const_handle; - typedef typename Kd_tree::Leaf_node_handle Leaf_node_handle; - typedef typename Kd_tree::Leaf_node_const_handle Leaf_node_const_handle; + typedef Kd_tree Kdt; + + typedef typename Kdt::Node_handle Node_handle; + typedef typename Kdt::Node_const_handle Node_const_handle; + typedef typename Kdt::Internal_node_handle Internal_node_handle; + typedef typename Kdt::Internal_node_const_handle Internal_node_const_handle; + typedef typename Kdt::Leaf_node_handle Leaf_node_handle; + typedef typename Kdt::Leaf_node_const_handle Leaf_node_const_handle; typedef typename TreeTraits::Point_d Point_d; typedef typename TreeTraits::FT FT; - typedef typename Kd_tree::Separator Separator; - typedef typename Kd_tree::Point_d_iterator Point_d_iterator; - typedef typename Kd_tree::iterator iterator; - typedef typename Kd_tree::D D; + typedef typename Kdt::Separator Separator; + typedef typename Kdt::Point_d_iterator Point_d_iterator; + typedef typename Kdt::iterator iterator; + typedef typename Kdt::D D; bool leaf; @@ -267,13 +269,13 @@ namespace CGAL { }; - template < class TreeTraits, class Splitter, class UseExtendedNode > - class Kd_tree_leaf_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ + template < class TreeTraits, class Splitter, class UseExtendedNode, class EnablePointsCache > + class Kd_tree_leaf_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode, EnablePointsCache >{ - friend class Kd_tree; + friend class Kd_tree; - typedef typename Kd_tree::iterator iterator; - typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; + typedef typename Kd_tree::iterator iterator; + typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode, EnablePointsCache> Base; typedef typename TreeTraits::Point_d Point_d; private: @@ -332,18 +334,20 @@ namespace CGAL { - template < class TreeTraits, class Splitter, class UseExtendedNode> - class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{ + template < class TreeTraits, class Splitter, class UseExtendedNode, class EnablePointsCache> + class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode, EnablePointsCache >{ - friend class Kd_tree; + friend class Kd_tree; - typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base; - typedef typename Kd_tree::Node_handle Node_handle; - typedef typename Kd_tree::Node_const_handle Node_const_handle; + typedef Kd_tree Kdt; + + typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode, EnablePointsCache> Base; + typedef typename Kdt::Node_handle Node_handle; + typedef typename Kdt::Node_const_handle Node_const_handle; typedef typename TreeTraits::FT FT; - typedef typename Kd_tree::Separator Separator; - typedef typename Kd_tree::D D; + typedef typename Kdt::Separator Separator; + typedef typename Kdt::D D; private: @@ -480,18 +484,21 @@ namespace CGAL { } };//internal node - template < class TreeTraits, class Splitter> - class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, Tag_false >{ + template < class TreeTraits, class Splitter, class EnablePointsCache> + class Kd_tree_internal_node + : public Kd_tree_node< TreeTraits, Splitter, Tag_false, EnablePointsCache > + { + friend class Kd_tree; + + typedef Kd_tree Kdt; - friend class Kd_tree; - - typedef Kd_tree_node< TreeTraits, Splitter, Tag_false> Base; - typedef typename Kd_tree::Node_handle Node_handle; - typedef typename Kd_tree::Node_const_handle Node_const_handle; + typedef Kd_tree_node< TreeTraits, Splitter, Tag_false, EnablePointsCache> Base; + typedef typename Kdt::Node_handle Node_handle; + typedef typename Kdt::Node_const_handle Node_const_handle; typedef typename TreeTraits::FT FT; - typedef typename Kd_tree::Separator Separator; - typedef typename Kd_tree::D D; + typedef typename Kdt::Separator Separator; + typedef typename Kdt::D D; private: diff --git a/Spatial_searching/include/CGAL/Orthogonal_incremental_neighbor_search.h b/Spatial_searching/include/CGAL/Orthogonal_incremental_neighbor_search.h index 27f5dd78bdd..05df9d48739 100644 --- a/Spatial_searching/include/CGAL/Orthogonal_incremental_neighbor_search.h +++ b/Spatial_searching/include/CGAL/Orthogonal_incremental_neighbor_search.h @@ -17,6 +17,7 @@ // // // Author(s) : Hans Tangelder () +// Clement Jamin (clement.jamin.pro@gmail.com) #ifndef CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH #define CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH @@ -32,6 +33,8 @@ #include #include +#include + namespace CGAL { template m_distance_helper; FT multiplication_factor; @@ -91,6 +95,8 @@ namespace CGAL { FT rd; + int m_dim; + Tree const& m_tree; class Priority_higher { public: @@ -146,27 +152,30 @@ namespace CGAL { FT Eps=FT(0.0), bool search_nearest=true) : traits(tree.traits()),number_of_neighbours_computed(0), number_of_internal_nodes_visited(0), number_of_leaf_nodes_visited(0), number_of_items_visited(0), - Orthogonal_distance_instance(tr), multiplication_factor(Orthogonal_distance_instance.transformed_distance(FT(1.0)+Eps)), + orthogonal_distance_instance(tr), + m_distance_helper(orthogonal_distance_instance), + multiplication_factor(orthogonal_distance_instance.transformed_distance(FT(1.0)+Eps)), query_point(q), search_nearest_neighbour(search_nearest), PriorityQueue(Priority_higher(search_nearest)), Item_PriorityQueue(Distance_smaller(search_nearest)), - reference_count(1) + reference_count(1), + m_tree(tree) { - if (tree.empty()) return; + if (m_tree.empty()) return; typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits.construct_cartesian_const_iterator_d_object(); - int dim = static_cast(std::distance(ccci(q), ccci(q,0))); + m_dim = static_cast(std::distance(ccci(q), ccci(q,0))); - dists.resize(dim); - for(int i=0 ; i + bool search_in_leaf(typename Tree::Leaf_node_const_handle node, bool search_furthest); + + // With cache + template<> + bool search_in_leaf(typename Tree::Leaf_node_const_handle node, bool search_furthest) + { + typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end(); + typename std::vector::const_iterator cache_point_begin = m_tree.cache_begin() + m_dim*(it_node_point - m_tree.begin()); + + for (; it_node_point != node->end(); ++it_node_point) + { + number_of_items_visited++; + FT distance_to_query_point = + m_distance_helper.transformed_distance_from_coordinates( + query_point, *it_node_point, cache_point_begin, cache_point_begin + m_dim); + + Point_with_transformed_distance *NN_Candidate = + new Point_with_transformed_distance(*it_node_point, distance_to_query_point); + Item_PriorityQueue.push(NN_Candidate); + + cache_point_begin += m_dim; + } + // old top of PriorityQueue has been processed, + // hence update rd + + bool next_neighbour_found; + if (!(PriorityQueue.empty())) + { + rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); + next_neighbour_found = (search_furthest ? + multiplication_factor*rd < Item_PriorityQueue.top()->second + : multiplication_factor*rd > Item_PriorityQueue.top()->second); + } + else // priority queue empty => last neighbour found + { + next_neighbour_found = true; + } + + number_of_neighbours_computed++; + return next_neighbour_found; + } + + // Without cache + template<> + bool search_in_leaf(typename Tree::Leaf_node_const_handle node, bool search_furthest) + { + typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end(); + + for (; it_node_point != node->end(); ++it_node_point) + { + number_of_items_visited++; + FT distance_to_query_point= + orthogonal_distance_instance.transformed_distance(query_point, *it_node_point); + + Point_with_transformed_distance *NN_Candidate = + new Point_with_transformed_distance(*it_node_point, distance_to_query_point); + Item_PriorityQueue.push(NN_Candidate); + } + // old top of PriorityQueue has been processed, + // hence update rd + + bool next_neighbour_found; + if (!(PriorityQueue.empty())) + { + rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); + next_neighbour_found = (search_furthest ? + multiplication_factor*rd < Item_PriorityQueue.top()->second + : multiplication_factor*rd > Item_PriorityQueue.top()->second); + } + else // priority queue empty => last neighbour found + { + next_neighbour_found=true; + } + + number_of_neighbours_computed++; + return next_neighbour_found; + } + + void Compute_the_next_nearest_neighbour() { @@ -290,7 +380,7 @@ namespace CGAL { FT diff2 = val - node->lower_high_value(); if (diff1 + diff2 < FT(0.0)) { new_rd= - Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); + orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); CGAL_assertion(new_rd >= copy_rd); dists[new_cut_dim] = diff1; @@ -302,7 +392,7 @@ namespace CGAL { } else { // compute new distance - new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); + new_rd=orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); CGAL_assertion(new_rd >= copy_rd); dists[new_cut_dim] = diff2; Node_with_distance *Lower_Child = @@ -316,31 +406,9 @@ namespace CGAL { typename Tree::Leaf_node_const_handle node = static_cast(N); number_of_leaf_nodes_visited++; - if (node->size() > 0) { - for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { - number_of_items_visited++; - FT distance_to_query_point= - Orthogonal_distance_instance.transformed_distance(query_point,*it); - Point_with_transformed_distance *NN_Candidate= - new Point_with_transformed_distance(*it,distance_to_query_point); - Item_PriorityQueue.push(NN_Candidate); - } - // old top of PriorityQueue has been processed, - // hence update rd - - if (!(PriorityQueue.empty())) { - rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); - next_neighbour_found = - (multiplication_factor*rd > - Item_PriorityQueue.top()->second); - } - else // priority queue empty => last neighbour found - { - next_neighbour_found=true; - } - - number_of_neighbours_computed++; - } + if (node->size() > 0) + next_neighbour_found = + search_in_leaf::type::value>::value>(node, false); } // next_neighbour_found or priority queue is empty // in the latter case also the item priority quee is empty } @@ -377,7 +445,7 @@ namespace CGAL { if (diff1 + diff2 < FT(0.0)) { diff1 = val - node->upper_high_value(); new_rd= - Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); + orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim); Node_with_distance *Lower_Child = new Node_with_distance(node->lower(), copy_rd, dists); PriorityQueue.push(Lower_Child); @@ -388,7 +456,7 @@ namespace CGAL { } else { // compute new distance diff2 = val - node->lower_low_value(); - new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); + new_rd=orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim); Node_with_distance *Upper_Child = new Node_with_distance(node->upper(), copy_rd, dists); PriorityQueue.push(Upper_Child); @@ -401,31 +469,9 @@ namespace CGAL { typename Tree::Leaf_node_const_handle node = static_cast(N); number_of_leaf_nodes_visited++; - if (node->size() > 0) { - for (typename Tree::iterator it=node->begin(); it != node->end(); it++) { - number_of_items_visited++; - FT distance_to_query_point= - Orthogonal_distance_instance.transformed_distance(query_point,*it); - Point_with_transformed_distance *NN_Candidate= - new Point_with_transformed_distance(*it,distance_to_query_point); - Item_PriorityQueue.push(NN_Candidate); - } - // old top of PriorityQueue has been processed, - // hence update rd - - if (!(PriorityQueue.empty())) { - rd = CGAL::cpp11::get<1>(*PriorityQueue.top()); - next_neighbour_found = - (multiplication_factor*rd < - Item_PriorityQueue.top()->second); - } - else // priority queue empty => last neighbour found - { - next_neighbour_found=true; - } - - number_of_neighbours_computed++; - } + if (node->size() > 0) + next_neighbour_found = + search_in_leaf::type::value>::value>(node, true); } // next_neighbour_found or priority queue is empty // in the latter case also the item priority quee is empty } diff --git a/Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h b/Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h index f98e11df40c..a5431c7b5ec 100644 --- a/Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h +++ b/Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h @@ -16,7 +16,9 @@ // $Id$ // // -// Author(s) : Gael Guennebaud (gael.guennebaud@inria.fr), Hans Tangelder () +// Author(s) : Gael Guennebaud (gael.guennebaud@inria.fr), +// Hans Tangelder (), +// Clement Jamin (clement.jamin.pro@gmail.com) #ifndef CGAL_ORTHOGONAL_K_NEIGHBOR_SEARCH_H #define CGAL_ORTHOGONAL_K_NEIGHBOR_SEARCH_H @@ -25,6 +27,7 @@ #include +#include namespace CGAL { @@ -32,29 +35,39 @@ template ::type, class Splitter= Sliding_midpoint , class Tree= Kd_tree > -class Orthogonal_k_neighbor_search: public internal::K_neighbor_search { - typedef internal::K_neighbor_search Base; +class Orthogonal_k_neighbor_search: public internal::K_neighbor_search +{ + typedef internal::K_neighbor_search Base; + typedef typename Tree::Point_d Point; - typename SearchTraits::Cartesian_const_iterator_d query_object_it; - - std::vector dists; public: typedef typename Base::FT FT; +private: + typename SearchTraits::Cartesian_const_iterator_d query_object_it; + + internal::Distance_helper m_distance_helper; + std::vector dists; + int m_dim; + Tree const& m_tree; + +public: Orthogonal_k_neighbor_search(const Tree& tree, const typename Base::Query_item& q, unsigned int k=1, FT Eps=FT(0.0), bool Search_nearest=true, const Distance& d=Distance(),bool sorted=true) - : Base(q,k,Eps,Search_nearest,d) + : Base(q,k,Eps,Search_nearest,d), + m_distance_helper(this->distance_instance), + m_tree(tree) { if (tree.empty()) return; typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=tree.traits().construct_cartesian_const_iterator_d_object(); query_object_it = construct_it(this->query_object); - int dim = static_cast(std::distance(query_object_it, construct_it(this->query_object,0))); + m_dim = static_cast(std::distance(query_object_it, construct_it(this->query_object,0))); - dists.resize(dim); - for(int i=0;iqueue.sort(); } private: + template + void search_nearest_in_leaf(typename Tree::Leaf_node_const_handle node); + + // With cache + template<> + void search_nearest_in_leaf(typename Tree::Leaf_node_const_handle node) + { + typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end(); + typename std::vector::const_iterator cache_point_begin = m_tree.cache_begin() + m_dim*(it_node_point - m_tree.begin()); + // As long as the queue is not full, the node is just inserted + for (; !this->queue.full() && it_node_point != it_node_point_end; ++it_node_point) + { + this->number_of_items_visited++; + + FT distance_to_query_object = + m_distance_helper.transformed_distance_from_coordinates( + this->query_object, *it_node_point, cache_point_begin, cache_point_begin + m_dim); + this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object)); + + cache_point_begin += m_dim; + } + // Now that the queue is full, we can gain time by keeping track + // of the current worst distance to interrupt the distance computation earlier + FT worst_dist = this->queue.top().second; + for (; it_node_point != it_node_point_end; ++it_node_point) + { + this->number_of_items_visited++; + + FT distance_to_query_object = + m_distance_helper.interruptible_transformed_distance( + this->query_object, *it_node_point, cache_point_begin, cache_point_begin + m_dim, worst_dist); + + if (distance_to_query_object < worst_dist) + { + this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object)); + worst_dist = this->queue.top().second; + } + + cache_point_begin += m_dim; + } + } + + // Without cache + template<> + void search_nearest_in_leaf(typename Tree::Leaf_node_const_handle node) + { + typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end(); + // As long as the queue is not full, the node is just inserted + for (; !this->queue.full() && it_node_point != it_node_point_end; ++it_node_point) + { + this->number_of_items_visited++; + + FT distance_to_query_object = + this->distance_instance.transformed_distance(this->query_object, *it_node_point); + this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object)); + } + // Now that the queue is full, we can gain time by keeping track + // of the current worst distance to interrupt the distance computation earlier + FT worst_dist = this->queue.top().second; + for (; it_node_point != it_node_point_end; ++it_node_point) + { + this->number_of_items_visited++; + + FT distance_to_query_object = + m_distance_helper.interruptible_transformed_distance( + this->query_object, *it_node_point, worst_dist); + + if (distance_to_query_object < worst_dist) + { + this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object)); + worst_dist = this->queue.top().second; + } + } + } + void compute_nearest_neighbors_orthogonally(typename Base::Node_const_handle N, FT rd) { - if (!(N->is_leaf())) - { - typename Tree::Internal_node_const_handle node = - static_cast(N); - this->number_of_internal_nodes_visited++; - int new_cut_dim=node->cutting_dimension(); - typename Base::Node_const_handle bestChild, otherChild; - FT new_off; - FT val = *(query_object_it + new_cut_dim); - FT diff1 = val - node->upper_low_value(); - FT diff2 = val - node->lower_high_value(); - if ( (diff1 + diff2 < FT(0.0)) ) - { - new_off = diff1; - bestChild = node->lower(); - otherChild = node->upper(); - } - else // compute new distance - { - new_off= diff2; - bestChild = node->upper(); - otherChild = node->lower(); - } - compute_nearest_neighbors_orthogonally(bestChild, rd); - FT dst=dists[new_cut_dim]; - FT new_rd = this->distance_instance.new_distance(rd,dst,new_off,new_cut_dim); - dists[new_cut_dim]=new_off; - if (this->branch_nearest(new_rd)) - { - compute_nearest_neighbors_orthogonally(otherChild, new_rd); - } - dists[new_cut_dim]=dst; - } - else + if (N->is_leaf()) { // n is a leaf typename Tree::Leaf_node_const_handle node = static_cast(N); this->number_of_leaf_nodes_visited++; - bool full = this->queue.full(); - FT worst_dist = this->queue.top().second; if (node->size() > 0) + search_nearest_in_leaf::type::value>::value>(node); + } + else + { + typename Tree::Internal_node_const_handle node = + static_cast(N); + this->number_of_internal_nodes_visited++; + int new_cut_dim = node->cutting_dimension(); + typename Base::Node_const_handle bestChild, otherChild; + FT new_off; + FT val = *(query_object_it + new_cut_dim); + FT diff1 = val - node->upper_low_value(); + FT diff2 = val - node->lower_high_value(); + if ((diff1 + diff2 < FT(0.0))) { - for (typename Tree::iterator it=node->begin(); it != node->end(); it++) - { - this->number_of_items_visited++; - FT distance_to_query_object= - this->distance_instance.transformed_distance(this->query_object,*it); - - if(!full || distance_to_query_object < worst_dist) - this->queue.insert(std::make_pair(&(*it),distance_to_query_object)); - } + new_off = diff1; + bestChild = node->lower(); + otherChild = node->upper(); } + else // compute new distance + { + new_off = diff2; + bestChild = node->upper(); + otherChild = node->lower(); + } + compute_nearest_neighbors_orthogonally(bestChild, rd); + FT dst = dists[new_cut_dim]; + FT new_rd = this->distance_instance.new_distance(rd, dst, new_off, new_cut_dim); + dists[new_cut_dim] = new_off; + if (this->branch_nearest(new_rd)) + { + compute_nearest_neighbors_orthogonally(otherChild, new_rd); + } + dists[new_cut_dim] = dst; } } - void compute_furthest_neighbors_orthogonally(typename Base::Node_const_handle N, FT rd) + template + void search_furthest_in_leaf(typename Tree::Leaf_node_const_handle node); + + // With cache + template<> + void search_furthest_in_leaf(typename Tree::Leaf_node_const_handle node) { - if (!(N->is_leaf())) + typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end(); + typename std::vector::const_iterator cache_point_begin = m_tree.cache_begin() + m_dim*(it_node_point - m_tree.begin()); + // In furthest search mode, the interruptible distance cannot be used to optimize + for (; it_node_point != it_node_point_end; ++it_node_point) + { + this->number_of_items_visited++; + + FT distance_to_query_object = + m_distance_helper.transformed_distance_from_coordinates( + this->query_object, *it_node_point, cache_point_begin, cache_point_begin + m_dim); + this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object)); + + cache_point_begin += m_dim; + } + } + + // Without cache + template<> + void search_furthest_in_leaf(typename Tree::Leaf_node_const_handle node) + { + typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end(); + // In furthest search mode, the interruptible distance cannot be used to optimize + for (; it_node_point != it_node_point_end; ++it_node_point) + { + this->number_of_items_visited++; + + FT distance_to_query_object = + this->distance_instance.transformed_distance(this->query_object, *it_node_point); + this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object)); + } + } + + void compute_furthest_neighbors_orthogonally(typename Base::Node_const_handle N, FT rd) + { + if (N->is_leaf()) + { + // n is a leaf + typename Tree::Leaf_node_const_handle node = + static_cast(N); + this->number_of_leaf_nodes_visited++; + if (node->size() > 0) + search_furthest_in_leaf::type::value>::value>(node); + } + else { typename Tree::Internal_node_const_handle node = static_cast(N); @@ -170,24 +290,7 @@ private: compute_furthest_neighbors_orthogonally(otherChild, new_rd); dists[new_cut_dim]=dst; } - else - { - // n is a leaf - typename Tree::Leaf_node_const_handle node = - static_cast(N); - this->number_of_leaf_nodes_visited++; - if (node->size() > 0) - { - for (typename Tree::iterator it=node->begin(); it != node->end(); it++) - { - this->number_of_items_visited++; - FT distance_to_query_object= - this->distance_instance.transformed_distance(this->query_object,*it); - this->queue.insert(std::make_pair(&(*it),distance_to_query_object)); - } - } - } - } + } }; // class diff --git a/Spatial_searching/include/CGAL/Search_traits_adapter.h b/Spatial_searching/include/CGAL/Search_traits_adapter.h index 60f7703255a..27b568376bb 100644 --- a/Spatial_searching/include/CGAL/Search_traits_adapter.h +++ b/Spatial_searching/include/CGAL/Search_traits_adapter.h @@ -143,7 +143,6 @@ public: template class Distance_adapter : public Base_distance { PointPropertyMap ppmap; - typedef typename Base_distance::FT FT; CGAL_static_assertion( ( boost::is_same< boost::lvalue_property_map_tag, typename boost::property_traits::category @@ -155,7 +154,8 @@ public: ):Base_distance(distance),ppmap(ppmap_){} using Base_distance::transformed_distance; - + + typedef typename Base_distance::FT FT; typedef Point_with_info Point_d; typedef typename Base_distance::Query_item Query_item; diff --git a/Spatial_searching/include/CGAL/internal/Search_helpers.h b/Spatial_searching/include/CGAL/internal/Search_helpers.h new file mode 100644 index 00000000000..782031e2c62 --- /dev/null +++ b/Spatial_searching/include/CGAL/internal/Search_helpers.h @@ -0,0 +1,148 @@ +// Copyright (c) 2002,2011 Utrecht University (The Netherlands). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Clement Jamin (clement.jamin.pro@gmail.com) + +#ifndef CGAL_INTERNAL_SEARCH_HELPERS_H +#define CGAL_INTERNAL_SEARCH_HELPERS_H + +#include + +#include + +namespace CGAL { +namespace internal { + +// Helper struct to know at compile-time if there is a cache of the points +// stored in the tree +template +struct Has_points_cache; + +template +struct Has_points_cache +{ + static const bool value = Tree::Enable_points_cache::value; +}; + +template +struct Has_points_cache +{ + static const bool value = false; +}; + +CGAL_GENERATE_MEMBER_DETECTOR(transformed_distance_from_coordinates); +CGAL_GENERATE_MEMBER_DETECTOR(interruptible_transformed_distance); +BOOST_MPL_HAS_XXX_TRAIT_NAMED_DEF(has_Enable_points_cache, Enable_points_cache, false) + + +template +class Distance_helper +{ + typedef typename Distance::FT FT; + typedef typename Distance::Point_d Point; + typedef typename Distance::Query_item Query_item; + +public: + + Distance_helper(Distance const& distance) + : m_distance(distance) + {} + + // If transformed_distance_from_coordinates does not exist in `Distance` + template ::value> + FT + transformed_distance_from_coordinates( + const typename Query_item& q, + Point const& p, + typename std::vector::const_iterator it_coord_begin, + typename std::vector::const_iterator it_coord_end) + { + return m_distance.transformed_distance(q, p); + } + // ... or if it exists + template <> + FT + transformed_distance_from_coordinates( + const typename Query_item& q, + Point const& p, + typename std::vector::const_iterator it_coord_begin, + typename std::vector::const_iterator it_coord_end) + { + return m_distance.transformed_distance_from_coordinates(q, it_coord_begin, it_coord_end); + } + + // *** Version with cache *** + // If interruptible_transformed_distance does not exist in `Distance` + template ::value> + FT + interruptible_transformed_distance( + const typename Query_item& q, + Point const& p, + typename std::vector::const_iterator it_coord_begin, + typename std::vector::const_iterator it_coord_end, + FT) + { + return transformed_distance_from_coordinates(q, p, it_coord_begin, it_coord_end); + } + // ... or if it exists + template <> + FT + interruptible_transformed_distance( + const typename Query_item& q, + Point const& p, + typename std::vector::const_iterator it_coord_begin, + typename std::vector::const_iterator it_coord_end, + FT stop_if_geq_to_this) + { + return m_distance.interruptible_transformed_distance( + q, it_coord_begin, it_coord_end, stop_if_geq_to_this); + } + + // *** Version without cache *** + // If interruptible_transformed_distance does not exist in `Distance` + template ::value> + FT + interruptible_transformed_distance( + const typename Query_item& q, + Point const& p, + FT) + { + return m_distance.transformed_distance(q, p); + } + // ... or if it exists + template <> + FT + interruptible_transformed_distance( + const typename Query_item& q, + Point const& p, + FT stop_if_geq_to_this) + { + typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = m_tree.traits().construct_cartesian_const_iterator_d_object(); + return m_distance.interruptible_transformed_distance( + q, construct_it(p), construct_it(p, 0), stop_if_geq_to_this); + } + +private: + Distance const & m_distance; + +}; // Distance_helper + + +}} // namespace CGAL::internal + +#endif // CGAL_INTERNAL_SEARCH_HELPERSS_H