Add radius search

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/Spatial_searching-Add_radius_search@2341 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: c30c32b27b8c3c1c8dd65e45d2bdd49738794062

4 files changed, 56 insertions, 19 deletions

diff --git a/src/Spatial_searching/doc/Intro_spatial_searching.h b/src/Spatial_searching/doc/Intro_spatial_searching.h
index 23705378..9a3c1b65 100644
--- a/src/Spatial_searching/doc/Intro_spatial_searching.h
+++ b/src/Spatial_searching/doc/Intro_spatial_searching.h
@@ -46,7 +46,7 @@ namespace spatial_searching {
  * 
  * \section spatial_searching_examples Example
  * 
- * This example generates 500 random points, then performs queries for nearest and farthest points using different methods.
+ * This example generates 500 random points, then performs radius search, and queries for nearest and farthest points using different methods.
  * 
  * \include Spatial_searching/example_spatial_searching.cpp
  * 
diff --git a/src/Spatial_searching/example/example_spatial_searching.cpp b/src/Spatial_searching/example/example_spatial_searching.cpp
index 14b324ae..9e6a8f32 100644
--- a/src/Spatial_searching/example/example_spatial_searching.cpp
+++ b/src/Spatial_searching/example/example_spatial_searching.cpp
@@ -48,5 +48,13 @@ int main(void) {
   for (auto ifs_iterator = ifn_range.begin(); ifs_iterator->first != 0; ++ifs_iterator)
     std::cout << ifs_iterator->first << " (sq. dist. = " << ifs_iterator->second << ")\n";
 
+  // Radius search
+  std::cout << "Radius search:\n";
+  std::vector<std::size_t> rs_result;
+  points_ds.radius_search(points[45], 0.5, std::back_inserter(rs_result));
+  K k;
+  for (auto const& p_idx : rs_result)
+    std::cout << p_idx << " (sq. dist. = " << k.squared_distance_d_object()(points[p_idx], points[45]) << ")\n";
+
   return 0;
 }
diff --git a/src/Spatial_searching/include/gudhi/Kd_tree_search.h b/src/Spatial_searching/include/gudhi/Kd_tree_search.h
index 6728d56e..374ebed6 100644
--- a/src/Spatial_searching/include/gudhi/Kd_tree_search.h
+++ b/src/Spatial_searching/include/gudhi/Kd_tree_search.h
@@ -27,6 +27,7 @@
 #include <CGAL/Orthogonal_incremental_neighbor_search.h>
 #include <CGAL/Search_traits.h>
 #include <CGAL/Search_traits_adapter.h>
+#include <CGAL/Fuzzy_sphere.h>
 #include <CGAL/property_map.h>
 
 #include <boost/property_map/property_map.hpp>
@@ -104,6 +105,7 @@ class Kd_tree_search {
   /// of a point P and `second` is the squared distance between P and the query point.
   typedef Incremental_neighbor_search                       INS_range;
 
+  typedef CGAL::Fuzzy_sphere<STraits>                       Fuzzy_sphere;
   /// \brief Constructor
   /// @param[in] points Const reference to the point range. This range
   /// is not copied, so it should not be destroyed or modified afterwards.
@@ -164,9 +166,9 @@ class Kd_tree_search {
   /// @param[in] k Number of nearest points to search.
   /// @param[in] sorted Indicates if the computed sequence of k-nearest neighbors needs to be sorted.
   /// @param[in] eps Approximation factor.
-  /// @return A range containing the k-nearest neighbors.
-  KNS_range query_k_nearest_neighbors(const
-    Point &p,
+  /// @return A range (whose `value_type` is `std::size_t`) containing the k-nearest neighbors.
+  KNS_range query_k_nearest_neighbors(
+    Point const& p,
     unsigned int k,
     bool sorted = true,
     FT eps = FT(0)) const {
@@ -179,8 +181,7 @@ class Kd_tree_search {
       k,
       eps,
       true,
-      CGAL::Distance_adapter<std::ptrdiff_t, Point_property_map, CGAL::Euclidean_distance<Traits_base> >(
-        std::begin(m_points)), sorted);
+      Orthogonal_distance(std::begin(m_points)), sorted);
 
     return search;
   }
@@ -188,10 +189,11 @@ class Kd_tree_search {
   /// \brief Search incrementally for the nearest neighbors from a query point.
   /// @param[in] p The query point.
   /// @param[in] eps Approximation factor.
-  /// @return A range containing the neighbors sorted by their distance to p.
+  /// @return A range (whose `value_type` is `std::size_t`) containing the 
+  /// neighbors sorted by their distance to p.
   /// All the neighbors are not computed by this function, but they will be
   /// computed incrementally when the iterator on the range is incremented.
-  INS_range query_incremental_nearest_neighbors(const Point &p, FT eps = FT(0)) const {
+  INS_range query_incremental_nearest_neighbors(Point const& p, FT eps = FT(0)) const {
     // Initialize the search structure, and search all N points
     // Note that we need to pass the Distance explicitly since it needs to
     // know the property map
@@ -200,8 +202,7 @@ class Kd_tree_search {
       p,
       eps,
       true,
-      CGAL::Distance_adapter<std::ptrdiff_t, Point_property_map, CGAL::Euclidean_distance<Traits_base> >(
-        std::begin(m_points)) );
+      Orthogonal_distance(std::begin(m_points)) );
 
     return search;
   }
@@ -211,9 +212,9 @@ class Kd_tree_search {
   /// @param[in] k Number of farthest points to search.
   /// @param[in] sorted Indicates if the computed sequence of k-farthest neighbors needs to be sorted.
   /// @param[in] eps Approximation factor.
-  /// @return A range containing the k-farthest neighbors.
-  KNS_range query_k_farthest_neighbors(const
-    Point &p,
+  /// @return A range (whose `value_type` is `std::size_t`) containing the k-farthest neighbors.
+  KNS_range query_k_farthest_neighbors(
+    Point const& p,
     unsigned int k,
     bool sorted = true,
     FT eps = FT(0)) const {
@@ -226,8 +227,7 @@ class Kd_tree_search {
       k,
       eps,
       false,
-      CGAL::Distance_adapter<std::ptrdiff_t, Point_property_map, CGAL::Euclidean_distance<Traits_base> >(
-        std::begin(m_points)), sorted);
+      Orthogonal_distance(std::begin(m_points)), sorted);
 
     return search;
   }
@@ -235,10 +235,11 @@ class Kd_tree_search {
   /// \brief Search incrementally for the farthest neighbors from a query point.
   /// @param[in] p The query point.
   /// @param[in] eps Approximation factor.
-  /// @return A range containing the neighbors sorted by their distance to p.
+  /// @return A range (whose `value_type` is `std::size_t`) 
+  /// containing the neighbors sorted by their distance to p.
   /// All the neighbors are not computed by this function, but they will be
   /// computed incrementally when the iterator on the range is incremented.
-  INS_range query_incremental_farthest_neighbors(const Point &p, FT eps = FT(0)) const {
+  INS_range query_incremental_farthest_neighbors(Point const& p, FT eps = FT(0)) const {
     // Initialize the search structure, and search all N points
     // Note that we need to pass the Distance explicitly since it needs to
     // know the property map
@@ -247,12 +248,32 @@ class Kd_tree_search {
       p,
       eps,
       false,
-      CGAL::Distance_adapter<std::ptrdiff_t, Point_property_map, CGAL::Euclidean_distance<Traits_base> >(
-        std::begin(m_points)) );
+      Orthogonal_distance(std::begin(m_points)) );
 
     return search;
   }
 
+  /// \brief Search for all the neighbors in a ball.
+  /// @param[in] p The query point.
+  /// @param[in] radius The search radius
+  /// @param[out] it The points that lie inside the sphere of center `p` and radius `radius`.
+  ///                The `value_type` of the iterator must be `Point`.
+  /// @param[in] eps Approximation factor.
+  template <typename OutputIterator>
+  void radius_search(
+    Point const& p,
+    FT radius,
+    OutputIterator it,
+    FT eps = FT(0)) const {
+    
+    m_tree.search(it, Fuzzy_sphere(p, radius, eps, m_tree.traits()));
+  }
+
+  int tree_depth() const
+  {
+    return m_tree.root()->depth();
+  }
+
  private:
   Point_range const& m_points;
   Tree m_tree;
diff --git a/src/Spatial_searching/test/test_Kd_tree_search.cpp b/src/Spatial_searching/test/test_Kd_tree_search.cpp
index 0ef22023..f79114bc 100644
--- a/src/Spatial_searching/test/test_Kd_tree_search.cpp
+++ b/src/Spatial_searching/test/test_Kd_tree_search.cpp
@@ -109,4 +109,12 @@ BOOST_AUTO_TEST_CASE(test_Kd_tree_search) {
 
   // Same result for KFN and IFN?
   BOOST_CHECK(kfn_result == ifn_result);
+
+  // Test radius search
+  Point rs_q(rd.get_double(-1., 1), rd.get_double(-1., 1), rd.get_double(-1., 1), rd.get_double(-1., 1));
+  std::vector<std::size_t> rs_result;
+  points_ds.radius_search(rs_q, 0.5, std::back_inserter(rs_result));
+  K k;
+  for (auto const& p_idx : rs_result)
+    BOOST_CHECK(k.squared_distance_d_object()(points[p_idx], rs_q) <= 0.5);
 }