diff options
Diffstat (limited to 'src/Spatial_searching')
| -rw-r--r-- | src/Spatial_searching/doc/Intro_spatial_searching.h | 48 | ||||
| -rw-r--r-- | src/Spatial_searching/example/CMakeLists.txt | 9 | ||||
| -rw-r--r-- | src/Spatial_searching/example/example_spatial_searching.cpp | 60 | ||||
| -rw-r--r-- | src/Spatial_searching/include/gudhi/Kd_tree_search.h | 287 | ||||
| -rw-r--r-- | src/Spatial_searching/test/CMakeLists.txt | 11 | ||||
| -rw-r--r-- | src/Spatial_searching/test/test_Kd_tree_search.cpp | 108 |
6 files changed, 523 insertions, 0 deletions
diff --git a/src/Spatial_searching/doc/Intro_spatial_searching.h b/src/Spatial_searching/doc/Intro_spatial_searching.h new file mode 100644 index 00000000..30805570 --- /dev/null +++ b/src/Spatial_searching/doc/Intro_spatial_searching.h @@ -0,0 +1,48 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clement Jamin + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef DOC_SPATIAL_SEARCHING_INTRO_SPATIAL_SEARCHING_H_ +#define DOC_SPATIAL_SEARCHING_INTRO_SPATIAL_SEARCHING_H_ + +// needs namespaces for Doxygen to link on classes +namespace Gudhi { +namespace spatial_searching { + +/** \defgroup spatial_searching Spatial_searching + * + * \author Clément Jamin + * + * @{ + * + * \section introduction Introduction + * + * This Gudhi component is a wrapper around + * <a target="_blank" href="http://doc.cgal.org/latest/Spatial_searching/index.html">CGAL dD spatial searching algorithms</a>. + * It provides a simplified API to perform (approximate) neighbor searches. Contrary to CGAL default behavior, the tree + * does not store the points themselves, but stores indices. + * + * For more details about the data structure or the algorithms, or for more advanced usages, reading + * <a target="_blank" href="http://doc.cgal.org/latest/Spatial_searching/index.html">CGAL documentation</a> + * is highly recommended. + * + * \section spatial_searching_examples Example + * + * This example generates 500 random points, then performs all-near-neighbors searches, and queries for nearest and furthest neighbors using different methods. + * + * \include Spatial_searching/example_spatial_searching.cpp + * + */ +/** @} */ // end defgroup spatial_searching + +} // namespace spatial_searching + +} // namespace Gudhi + +#endif // DOC_SPATIAL_SEARCHING_INTRO_SPATIAL_SEARCHING_H_ diff --git a/src/Spatial_searching/example/CMakeLists.txt b/src/Spatial_searching/example/CMakeLists.txt new file mode 100644 index 00000000..eeb3e85f --- /dev/null +++ b/src/Spatial_searching/example/CMakeLists.txt @@ -0,0 +1,9 @@ +project(Spatial_searching_examples) + +if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + add_executable( Spatial_searching_example_spatial_searching example_spatial_searching.cpp ) + target_link_libraries(Spatial_searching_example_spatial_searching ${CGAL_LIBRARY}) + add_test(NAME Spatial_searching_example_spatial_searching + COMMAND $<TARGET_FILE:Spatial_searching_example_spatial_searching>) + install(TARGETS Spatial_searching_example_spatial_searching DESTINATION bin) +endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) diff --git a/src/Spatial_searching/example/example_spatial_searching.cpp b/src/Spatial_searching/example/example_spatial_searching.cpp new file mode 100644 index 00000000..034ad24a --- /dev/null +++ b/src/Spatial_searching/example/example_spatial_searching.cpp @@ -0,0 +1,60 @@ +#include <gudhi/Kd_tree_search.h> + +#include <CGAL/Epick_d.h> +#include <CGAL/Random.h> + +#include <vector> + +namespace gss = Gudhi::spatial_searching; + +int main(void) { + typedef CGAL::Epick_d<CGAL::Dimension_tag<4> > K; + typedef typename K::Point_d Point; + typedef std::vector<Point> Points; + + typedef gss::Kd_tree_search<K, Points> Points_ds; + + CGAL::Random rd; + + Points points; + for (int i = 0; i < 500; ++i) + points.push_back(Point(rd.get_double(-1., 1), rd.get_double(-1., 1), rd.get_double(-1., 1), rd.get_double(-1., 1))); + + Points_ds points_ds(points); + + // 10-nearest neighbor query + std::cout << "10 nearest neighbors from points[20]:\n"; + auto knn_range = points_ds.k_nearest_neighbors(points[20], 10, true); + for (auto const& nghb : knn_range) + std::cout << nghb.first << " (sq. dist. = " << nghb.second << ")\n"; + + // Incremental nearest neighbor query + std::cout << "Incremental nearest neighbors:\n"; + auto inn_range = points_ds.incremental_nearest_neighbors(points[45]); + // Get the neighbors in distance order until we hit the first point + for (auto ins_iterator = inn_range.begin(); ins_iterator->first != 0; ++ins_iterator) + std::cout << ins_iterator->first << " (sq. dist. = " << ins_iterator->second << ")\n"; + + // 10-furthest neighbor query + std::cout << "10 furthest neighbors from points[20]:\n"; + auto kfn_range = points_ds.k_furthest_neighbors(points[20], 10, true); + for (auto const& nghb : kfn_range) + std::cout << nghb.first << " (sq. dist. = " << nghb.second << ")\n"; + + // Incremental furthest neighbor query + std::cout << "Incremental furthest neighbors:\n"; + auto ifn_range = points_ds.incremental_furthest_neighbors(points[45]); + // Get the neighbors in distance reverse order until we hit the first point + for (auto ifs_iterator = ifn_range.begin(); ifs_iterator->first != 0; ++ifs_iterator) + std::cout << ifs_iterator->first << " (sq. dist. = " << ifs_iterator->second << ")\n"; + + // All-near-neighbors search + std::cout << "All-near-neighbors search:\n"; + std::vector<std::size_t> rs_result; + points_ds.all_near_neighbors(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 new file mode 100644 index 00000000..fedbb32e --- /dev/null +++ b/src/Spatial_searching/include/gudhi/Kd_tree_search.h @@ -0,0 +1,287 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clement Jamin + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL and Eigen3 + * - YYYY/MM Author: Description of the modification + */ + +#ifndef KD_TREE_SEARCH_H_ +#define KD_TREE_SEARCH_H_ + +#include <CGAL/Orthogonal_k_neighbor_search.h> +#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 <CGAL/version.h> // for CGAL_VERSION_NR + +#include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST + +#include <boost/property_map/property_map.hpp> +#include <boost/iterator/counting_iterator.hpp> + +#include <cstddef> +#include <vector> + +// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 +#if CGAL_VERSION_NR < 1041101000 +# error Alpha_complex_3d is only available for CGAL >= 4.11 +#endif + +#if !EIGEN_VERSION_AT_LEAST(3,1,0) +# error Alpha_complex_3d is only available for Eigen3 >= 3.1.0 installed with CGAL +#endif + +namespace Gudhi { +namespace spatial_searching { + + + /** + * \class Kd_tree_search Kd_tree_search.h gudhi/Kd_tree_search.h + * \brief Spatial tree data structure to perform (approximate) nearest and furthest neighbor search. + * + * \ingroup spatial_searching + * + * \details + * The class Kd_tree_search is a tree-based data structure, based on + * <a target="_blank" href="http://doc.cgal.org/latest/Spatial_searching/index.html">CGAL dD spatial searching data structures</a>. + * It provides a simplified API to perform (approximate) nearest and furthest neighbor searches. Contrary to CGAL default behavior, the tree + * does not store the points themselves, but stores indices. + * + * There are two types of queries: the <i>k-nearest or k-furthest neighbor query</i>, where <i>k</i> is fixed and the <i>k</i> nearest + * or furthest points are computed right away, + * and the <i>incremental nearest or furthest neighbor query</i>, where no number of neighbors is provided during the call, as the + * neighbors will be computed incrementally when the iterator on the range is incremented. + * + * \tparam Search_traits must be a model of the <a target="_blank" + * href="http://doc.cgal.org/latest/Spatial_searching/classSearchTraits.html">SearchTraits</a> + * concept, such as the <a target="_blank" + * href="http://doc.cgal.org/latest/Kernel_d/classCGAL_1_1Epick__d.html">CGAL::Epick_d</a> class, which + * can be static if you know the ambiant dimension at compile-time, or dynamic if you don't. + * \tparam Point_range is the type of the range that provides the points. + * It must be a range whose iterator type is a `RandomAccessIterator`. + */ +template <typename Search_traits, typename Point_range> +class Kd_tree_search { + typedef boost::iterator_property_map< + typename Point_range::const_iterator, + CGAL::Identity_property_map<std::ptrdiff_t> > Point_property_map; + + public: + /// The Traits. + typedef Search_traits Traits; + /// Number type used for distances. + typedef typename Traits::FT FT; + /// The point type. + typedef typename Point_range::value_type Point; + + typedef CGAL::Search_traits< + FT, Point, + typename Traits::Cartesian_const_iterator_d, + typename Traits::Construct_cartesian_const_iterator_d> Traits_base; + + typedef CGAL::Search_traits_adapter< + 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; + typedef typename K_neighbor_search::Distance Distance; + /// \brief The range returned by a k-nearest or k-furthest neighbor search. + /// Its value type is `std::pair<std::size_t, FT>` where `first` is the index + /// of a point P and `second` is the squared distance between P and the query point. + typedef K_neighbor_search KNS_range; + + typedef CGAL::Orthogonal_incremental_neighbor_search< + STraits, Distance, CGAL::Sliding_midpoint<STraits>, Tree> + Incremental_neighbor_search; + /// \brief The range returned by an incremental nearest or furthest neighbor search. + /// Its value type is `std::pair<std::size_t, FT>` where `first` is the index + /// 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. + Kd_tree_search(Point_range const& points) + : m_points(points), + m_tree(boost::counting_iterator<std::ptrdiff_t>(0), + boost::counting_iterator<std::ptrdiff_t>(points.size()), + typename Tree::Splitter(), + STraits(std::begin(points))) { + // Build the tree now (we don't want to wait for the first query) + m_tree.build(); + } + + /// \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. + /// @param[in] only_these_points Specifies the indices of the points that + /// should be actually inserted into the tree. The other points are ignored. + template <typename Point_indices_range> + Kd_tree_search( + Point_range const& points, + Point_indices_range const& only_these_points) + : m_points(points), + m_tree( + only_these_points.begin(), only_these_points.end(), + typename Tree::Splitter(), + STraits(std::begin(points))) { + // Build the tree now (we don't want to wait for the first query) + m_tree.build(); + } + + /// \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. + /// @param[in] begin_idx, past_the_end_idx Define the subset of the points that + /// should be actually inserted into the tree. The other points are ignored. + Kd_tree_search( + Point_range const& points, + std::size_t begin_idx, std::size_t past_the_end_idx) + : m_points(points), + m_tree( + boost::counting_iterator<std::ptrdiff_t>(begin_idx), + boost::counting_iterator<std::ptrdiff_t>(past_the_end_idx), + typename Tree::Splitter(), + STraits(std::begin(points))) { + // Build the tree now (we don't want to wait for the first query) + m_tree.build(); + } + + // Be careful, this function invalidates the tree, + // which will be recomputed at the next query + void insert(std::ptrdiff_t point_idx) { + m_tree.insert(point_idx); + } + + /// \brief Search for the k-nearest neighbors from a query point. + /// @param[in] p The query point. + /// @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 (whose `value_type` is `std::size_t`) containing the k-nearest neighbors. + KNS_range k_nearest_neighbors( + Point const& p, + unsigned int k, + bool sorted = true, + 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 + K_neighbor_search search( + m_tree, + p, + k, + eps, + true, + Orthogonal_distance(std::begin(m_points)), sorted); + + return 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 (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 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 + Incremental_neighbor_search search( + m_tree, + p, + eps, + true, + Orthogonal_distance(std::begin(m_points)) ); + + return search; + } + + /// \brief Search for the k-furthest points from a query point. + /// @param[in] p The query point. + /// @param[in] k Number of furthest points to search. + /// @param[in] sorted Indicates if the computed sequence of k-furthest neighbors needs to be sorted. + /// @param[in] eps Approximation factor. + /// @return A range (whose `value_type` is `std::size_t`) containing the k-furthest neighbors. + KNS_range k_furthest_neighbors( + Point const& p, + unsigned int k, + bool sorted = true, + 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 + K_neighbor_search search( + m_tree, + p, + k, + eps, + false, + Orthogonal_distance(std::begin(m_points)), sorted); + + return search; + } + + /// \brief Search incrementally for the furthest neighbors from a query point. + /// @param[in] p The query point. + /// @param[in] eps Approximation factor. + /// @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 incremental_furthest_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 + Incremental_neighbor_search search( + m_tree, + p, + eps, + false, + 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 all_near_neighbors(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; +}; + +} // namespace spatial_searching +} // namespace Gudhi + +#endif // KD_TREE_SEARCH_H_ diff --git a/src/Spatial_searching/test/CMakeLists.txt b/src/Spatial_searching/test/CMakeLists.txt new file mode 100644 index 00000000..18f7c6b8 --- /dev/null +++ b/src/Spatial_searching/test/CMakeLists.txt @@ -0,0 +1,11 @@ +project(Spatial_searching_tests) + +if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + include(GUDHI_test_coverage) + + add_executable( Spatial_searching_test_Kd_tree_search test_Kd_tree_search.cpp ) + target_link_libraries(Spatial_searching_test_Kd_tree_search + ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + + gudhi_add_coverage_test(Spatial_searching_test_Kd_tree_search) +endif () diff --git a/src/Spatial_searching/test/test_Kd_tree_search.cpp b/src/Spatial_searching/test/test_Kd_tree_search.cpp new file mode 100644 index 00000000..d6c6fba3 --- /dev/null +++ b/src/Spatial_searching/test/test_Kd_tree_search.cpp @@ -0,0 +1,108 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Clement Jamin + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE Spatial_searching - test Kd_tree_search +#include <boost/test/unit_test.hpp> + +#include <gudhi/Kd_tree_search.h> + +#include <CGAL/Epick_d.h> +#include <CGAL/Random.h> + +#include <vector> + +BOOST_AUTO_TEST_CASE(test_Kd_tree_search) { + typedef CGAL::Epick_d<CGAL::Dimension_tag<4> > K; + typedef K::FT FT; + typedef K::Point_d Point; + typedef std::vector<Point> Points; + + typedef Gudhi::spatial_searching::Kd_tree_search< + K, Points> Points_ds; + + CGAL::Random rd; + + Points points; + for (int i = 0; i < 500; ++i) + points.push_back(Point(rd.get_double(-1., 1), rd.get_double(-1., 1), rd.get_double(-1., 1), rd.get_double(-1., 1))); + + Points_ds points_ds(points); + + // Test k_nearest_neighbors + std::size_t closest_pt_index = + points_ds.k_nearest_neighbors(points[10], 1, false).begin()->first; + BOOST_CHECK(closest_pt_index == 10); + + auto kns_range = points_ds.k_nearest_neighbors(points[20], 10, true); + + std::vector<std::size_t> knn_result; + FT last_dist = -1.; + for (auto const& nghb : kns_range) { + BOOST_CHECK(nghb.second > last_dist); + knn_result.push_back(nghb.second); + last_dist = nghb.second; + } + + // Test incremental_nearest_neighbors + closest_pt_index = + points_ds.incremental_nearest_neighbors(points[10]).begin()->first; + BOOST_CHECK(closest_pt_index == 10); + + auto inn_range = points_ds.incremental_nearest_neighbors(points[20]); + + std::vector<std::size_t> inn_result; + last_dist = -1.; + auto inn_it = inn_range.begin(); + for (int i = 0; i < 10; ++inn_it, ++i) { + auto const& nghb = *inn_it; + BOOST_CHECK(nghb.second > last_dist); + inn_result.push_back(nghb.second); + last_dist = nghb.second; + } + + // Same result for KNN and INN? + BOOST_CHECK(knn_result == inn_result); + + // Test k_furthest_neighbors + auto kfn_range = points_ds.k_furthest_neighbors(points[20], 10, true); + + std::vector<std::size_t> kfn_result; + last_dist = kfn_range.begin()->second; + for (auto const& nghb : kfn_range) { + BOOST_CHECK(nghb.second <= last_dist); + kfn_result.push_back(nghb.second); + last_dist = nghb.second; + } + + // Test k_furthest_neighbors + auto ifn_range = points_ds.incremental_furthest_neighbors(points[20]); + + std::vector<std::size_t> ifn_result; + last_dist = ifn_range.begin()->second; + auto ifn_it = ifn_range.begin(); + for (int i = 0; i < 10; ++ifn_it, ++i) { + auto const& nghb = *ifn_it; + BOOST_CHECK(nghb.second <= last_dist); + ifn_result.push_back(nghb.second); + last_dist = nghb.second; + } + + // Same result for KFN and IFN? + BOOST_CHECK(kfn_result == ifn_result); + + // Test all_near_neighbors + 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.all_near_neighbors(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); +} |