diff options
4 files changed, 60 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..f13a98f7 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> @@ -87,6 +88,10 @@ class Kd_tree_search { std::ptrdiff_t, Point_property_map, Traits_base> STraits; + typedef CGAL::Distance_adapter< + std::ptrdiff_t, + Point_property_map, + CGAL::Euclidean_distance<Traits_base> > Orthogonal_distance; typedef CGAL::Orthogonal_k_neighbor_search<STraits> K_neighbor_search; typedef typename K_neighbor_search::Tree Tree; @@ -104,6 +109,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 +170,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 +185,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 +193,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 +206,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 +216,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 +231,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 +239,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 +252,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`. + /// Note: `it` is used this way: `*it++ = each_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); } |