From 91b0ff839b8058d3f5767e6ed80b93c23be2c98a Mon Sep 17 00:00:00 2001 From: Hind-M Date: Mon, 9 Aug 2021 16:21:24 +0200 Subject: First version of cech enhancement --- src/Cech_complex/include/gudhi/Cech_complex.h | 26 +++-- .../include/gudhi/Cech_complex_blocker.h | 110 ++++++++++++++++++++- 2 files changed, 123 insertions(+), 13 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index b0871e10..2c6f202a 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -40,7 +40,7 @@ namespace cech_complex { * \tparam ForwardPointRange must be a range for which `std::begin()` and `std::end()` methods return input * iterators on a point. `std::begin()` and `std::end()` methods are also required for a point. */ -template +template class Cech_complex { private: // Required by compute_proximity_graph @@ -49,14 +49,15 @@ class Cech_complex { using Proximity_graph = Gudhi::Proximity_graph; // Retrieve Coordinate type from ForwardPointRange - using Point_from_range_iterator = typename boost::range_const_iterator::type; - using Point_from_range = typename std::iterator_traits::value_type; - using Coordinate_iterator = typename boost::range_const_iterator::type; - using Coordinate = typename std::iterator_traits::value_type; +// using Point_from_range_iterator = typename boost::range_const_iterator::type; +// using Point_from_range = typename std::iterator_traits::value_type; +// using Coordinate_iterator = typename boost::range_const_iterator::type; +// using Coordinate = typename std::iterator_traits::value_type; public: // Point and Point_cloud type definition - using Point = std::vector; + //using Point = std::vector; + using Point = typename Kernel::Point_d; using Point_cloud = std::vector; public: @@ -70,9 +71,13 @@ class Cech_complex { */ Cech_complex(const ForwardPointRange& points, Filtration_value max_radius) : max_radius_(max_radius) { // Point cloud deep copy - point_cloud_.reserve(boost::size(points)); - for (auto&& point : points) point_cloud_.emplace_back(std::begin(point), std::end(point)); +// point_cloud_.reserve(boost::size(points)); +// for (auto&& point : points) point_cloud_.emplace_back(point.cartesian_begin(), point.cartesian_end()); + + point_cloud_.assign(points.begin(), points.end()); + + std::clog << "Hind: Just before the graph compute" << std::endl; cech_skeleton_graph_ = Gudhi::compute_proximity_graph( point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); } @@ -87,14 +92,17 @@ class Cech_complex { */ template void create_complex(SimplicialComplexForCechComplex& complex, int dim_max) { + std::clog << "Hind: in create complex" << std::endl; GUDHI_CHECK(complex.num_vertices() == 0, std::invalid_argument("Cech_complex::create_complex - simplicial complex is not empty")); // insert the proximity graph in the simplicial complex + std::clog << "Hind: before insert_graph" << std::endl; complex.insert_graph(cech_skeleton_graph_); // expand the graph until dimension dim_max + std::clog << "Hind: before expansion_with_blockers" << std::endl; complex.expansion_with_blockers(dim_max, - Cech_blocker(&complex, this)); + Cech_blocker(&complex, this)); } /** @return max_radius value given at construction. */ diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 31b9aab5..f2cf5ccc 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -11,8 +11,15 @@ #ifndef CECH_COMPLEX_BLOCKER_H_ #define CECH_COMPLEX_BLOCKER_H_ +// TODO to remove #include // for Gudhi::Minimal_enclosing_ball_radius +#include // for CGAL::to_double +#include // +// #include +#include // For EXACT or SAFE version +#include + #include #include #include // for std::sqrt @@ -35,28 +42,120 @@ namespace cech_complex { * * \tparam Chech_complex is required by the blocker. */ -template +template class Cech_blocker { private: - using Point_cloud = typename Cech_complex::Point_cloud; +// using Point_cloud = typename Cech_complex::Point_cloud; using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; using Filtration_value = typename SimplicialComplexForCech::Filtration_value; public: + + using Point_d = typename Kernel::Point_d; + // Numeric type of coordinates in the kernel + using FT = typename Kernel::FT; + // Sphere is a pair of point and squared radius. + using Sphere = typename std::pair; + + + // Add an int in TDS to save point index in the structure +// using TDS = CGAL::Triangulation_data_structure, +// CGAL::Triangulation_full_cell >; +// +// /** \brief A (Weighted or not) Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/ +// using Triangulation = CGAL::Delaunay_triangulation; +// // Vertex_iterator type from CGAL. +// using CGAL_vertex_iterator = typename Triangulation::Vertex_iterator; + + // Structure to switch from simplex tree vertex handle to CGAL vertex iterator. + //using Vector_vertex_iterator = std::vector< CGAL_vertex_iterator >; + + + + /** \brief get_point_ returns the point corresponding to the vertex given as parameter. + * Only for internal use for faster access. + * + * @param[in] vertex Vertex handle of the point to retrieve. + * @return The point found. + */ +/* const Point_d& get_point_(std::size_t vertex) const { + return vertex_handle_to_iterator_[vertex]->point(); + } */ + + /** \internal \brief TODO + * \param[in] + * \return */ + template + FT get_squared_radius(PointIterator begin, PointIterator end) const { + return kernel_.compute_squared_radius_d_object()(begin, end); + } + + /** \internal \brief TODO + * \param[in] + * \return */ + template + Sphere get_sphere(PointIterator begin, PointIterator end) const { + Point_d c = kernel_.construct_circumcenter_d_object()(begin, end); + FT r = kernel_.squared_distance_d_object()(c, *begin); + return std::make_pair(std::move(c), std::move(r)); + } + + /** \internal \brief Čech complex blocker operator() - the oracle - assigns the filtration value from the simplex * radius and returns if the simplex expansion must be blocked. * \param[in] sh The Simplex_handle. * \return true if the simplex radius is greater than the Cech_complex max_radius*/ bool operator()(Simplex_handle sh) { + using Point_cloud = std::vector; Point_cloud points; + + // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices + for (auto face : sc_ptr_->simplex_vertex_range(sh)) { + ///////////////////////////////////////////////////////////////////// + + + + + /////////////////////////////////// + } + for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { - points.push_back(cc_ptr_->get_point(vertex)); + points.push_back(cc_ptr_->get_point(vertex)); +// points.push_back(get_point_(vertex)); #ifdef DEBUG_TRACES std::clog << "#(" << vertex << ")#"; #endif // DEBUG_TRACES } - Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); + // TODO to remove + //Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); + // Hind: Here change the algo of the enclosing Minimal_enclosing_ball_radius + auto point_to_be_inserted = points.back(); + Sphere sph = get_sphere(points.cbegin(), points.cend()-1); + +// Sphere sph = get_sphere(points.cbegin(), points.cend()-1); + CGAL::NT_converter cast_to_double; +// CGAL::NT_converter cast_point_d_to_double; + + std::clog << "circumcenter: " << sph.first << ", radius: " << std::sqrt(cast_to_double(sph.second))<< std::endl; + // TODO to remove + // Filtration_value test = std::sqrt(CGAL::to_double(sph.second)); + + + // Check that the point to be inserted is already included in the sphere of the simplex containing the preceding points + // TODO instead of Euclidean_distance ; use kernel_.squared_distance_d_object()(c, *begin); + // Add a loop on the three faces to check sphere before computing the circumsphere + // Add the computed sphere as cache; a vector of spheres depending on the number of faces ? + // +// if (Gudhi::Euclidean_distance()(cast_point_d_to_double(sph.first), point_to_be_inserted) > std::sqrt(cast_to_double(sph.second))) +// FT r = kernel_.squared_distance_d_object()(sph.first, sph.first); //*(points.cend()-1)); + if (kernel_.squared_distance_d_object()(sph.first, point_to_be_inserted) > sph.second) + sph = get_sphere(points.cbegin(), points.cend()); + + Filtration_value radius = std::sqrt(cast_to_double(sph.second)); + + #ifdef DEBUG_TRACES if (radius > cc_ptr_->max_radius()) std::clog << "radius > max_radius => expansion is blocked\n"; #endif // DEBUG_TRACES @@ -70,6 +169,9 @@ class Cech_blocker { private: SimplicialComplexForCech* sc_ptr_; Cech_complex* cc_ptr_; + Kernel kernel_; + //Vector_vertex_iterator vertex_handle_to_iterator_; + }; } // namespace cech_complex -- cgit v1.2.3 From bc28892cbae3d9a9fcc19a0fbcfcc98bb9195ff7 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Wed, 18 Aug 2021 17:45:51 +0200 Subject: Modify the algorithm to get the minimal enclosing ball --- .../include/gudhi/Cech_complex_blocker.h | 167 +++++++++++++-------- 1 file changed, 102 insertions(+), 65 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index f2cf5ccc..acb53143 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -22,6 +22,9 @@ #include #include +#include +#include +#include #include // for std::sqrt namespace Gudhi { @@ -57,32 +60,7 @@ class Cech_blocker { using FT = typename Kernel::FT; // Sphere is a pair of point and squared radius. using Sphere = typename std::pair; - - - // Add an int in TDS to save point index in the structure -// using TDS = CGAL::Triangulation_data_structure, -// CGAL::Triangulation_full_cell >; -// -// /** \brief A (Weighted or not) Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/ -// using Triangulation = CGAL::Delaunay_triangulation; -// // Vertex_iterator type from CGAL. -// using CGAL_vertex_iterator = typename Triangulation::Vertex_iterator; - - // Structure to switch from simplex tree vertex handle to CGAL vertex iterator. - //using Vector_vertex_iterator = std::vector< CGAL_vertex_iterator >; - - - /** \brief get_point_ returns the point corresponding to the vertex given as parameter. - * Only for internal use for faster access. - * - * @param[in] vertex Vertex handle of the point to retrieve. - * @return The point found. - */ -/* const Point_d& get_point_(std::size_t vertex) const { - return vertex_handle_to_iterator_[vertex]->point(); - } */ /** \internal \brief TODO * \param[in] @@ -102,6 +80,19 @@ class Cech_blocker { return std::make_pair(std::move(c), std::move(r)); } + /** \internal \brief TODO + * \param[in] + * \return */ + template + class CompareSpheresRadii + { + public: + CGAL::NT_converter cast_to_double; + bool operator()(const Sphere& firstSphere, const Sphere& secondSphere) + { + return cast_to_double(firstSphere.second) < cast_to_double(secondSphere.second); + } + }; /** \internal \brief Čech complex blocker operator() - the oracle - assigns the filtration value from the simplex * radius and returns if the simplex expansion must be blocked. @@ -109,51 +100,98 @@ class Cech_blocker { * \return true if the simplex radius is greater than the Cech_complex max_radius*/ bool operator()(Simplex_handle sh) { using Point_cloud = std::vector; - Point_cloud points; + CGAL::NT_converter cast_to_double; + Filtration_value radius = 0.; + std::string key_to_permute; // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices - for (auto face : sc_ptr_->simplex_vertex_range(sh)) { - ///////////////////////////////////////////////////////////////////// - - + // std::set enclosing_ball_radii; + std::set > enclosing_ball_spheres; + for (auto face : sc_ptr_->boundary_simplex_range(sh)) { + // Find which vertex of sh is missing in face. We rely on the fact that simplex_vertex_range is sorted. + auto longlist = sc_ptr_->simplex_vertex_range(sh); + auto shortlist = sc_ptr_->simplex_vertex_range(face); + + std::clog << "Hind debug: within FACE loop "<< std::endl; + // TODO to remove + for (auto i = std::begin(longlist); i != std::end(longlist);++i) + std::clog << "Hind debug: longlist: " << cc_ptr_->get_point(*i) << std::endl; + for (auto i = std::begin(shortlist); i != std::end(shortlist);++i) + std::clog << "Hind debug: shortlist: " << cc_ptr_->get_point(*i) << std::endl; + + auto longiter = std::begin(longlist); + auto shortiter = std::begin(shortlist); + auto enditer = std::end(shortlist); + while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; } + auto extra = *longiter; // Vertex_handle + + std::clog << "Hind debug: extra vertex: " << cc_ptr_->get_point(extra) << std::endl; - /////////////////////////////////// + Point_cloud face_points; + std::string key, key_extra; + for (auto vertex : sc_ptr_->simplex_vertex_range(face)) { + face_points.push_back(cc_ptr_->get_point(vertex)); + key.append(std::to_string(vertex)); + #ifdef DEBUG_TRACES + std::clog << "#(" << vertex << ")#"; + #endif // DEBUG_TRACES + } + key_extra = key; + key_extra.append(std::to_string(extra)); + key_to_permute = key_extra; + std::clog << "END OF VERTICES " << std::endl; + std::clog << "KEY is: " << key << std::endl; + std::clog << "KEY extra is: " << key_extra << std::endl; + Sphere sph; + auto it = cache_.find(key); + if(it != cache_.end()) + sph = it->second; + else { + sph = get_sphere(face_points.cbegin(), face_points.cend()); + } + if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) <= sph.second) { + radius = std::sqrt(cast_to_double(sph.second)); + #ifdef DEBUG_TRACES + std::clog << "circumcenter: " << sph.first << ", radius: " << radius << std::endl; + #endif // DEBUG_TRACES + std::clog << "distance FYI: " << kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " < " << cast_to_double(sph.second) << std::endl; + // enclosing_ball_radii.insert(radius); + enclosing_ball_spheres.insert(sph); + cache_[key_extra] = sph; + } + else {// TODO to remove + std::clog << "vertex not included BECAUSE DISTANCE: "<< kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " AND RAD SPHERE: " << sph.second << std::endl; + } } - - for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { - points.push_back(cc_ptr_->get_point(vertex)); -// points.push_back(get_point_(vertex)); -#ifdef DEBUG_TRACES - std::clog << "#(" << vertex << ")#"; -#endif // DEBUG_TRACES + // Get the minimal radius of all faces enclosing balls if exists + if (!enclosing_ball_spheres.empty()) { + // radius = *enclosing_ball_radii.begin(); + Sphere sph_min = *enclosing_ball_spheres.begin(); + radius = std::sqrt(cast_to_double(sph_min.second)); + // std::clog << "CHECK that radius of min sphere is min radius: " << std::sqrt(cast_to_double(sph_min.second)) << "; and RADIUS min: " << radius << std::endl; + // Set all key_to_permute permutations to min sphere in cache + do + { + cache_[key_to_permute] = sph_min; + } while(std::next_permutation(key_to_permute.begin(), key_to_permute.end())); } - // TODO to remove - //Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); - // Hind: Here change the algo of the enclosing Minimal_enclosing_ball_radius - auto point_to_be_inserted = points.back(); - Sphere sph = get_sphere(points.cbegin(), points.cend()-1); - -// Sphere sph = get_sphere(points.cbegin(), points.cend()-1); - CGAL::NT_converter cast_to_double; -// CGAL::NT_converter cast_point_d_to_double; + std::clog << "END OF FACES ; radius = " << radius << std::endl; - std::clog << "circumcenter: " << sph.first << ", radius: " << std::sqrt(cast_to_double(sph.second))<< std::endl; - // TODO to remove - // Filtration_value test = std::sqrt(CGAL::to_double(sph.second)); - - - // Check that the point to be inserted is already included in the sphere of the simplex containing the preceding points - // TODO instead of Euclidean_distance ; use kernel_.squared_distance_d_object()(c, *begin); - // Add a loop on the three faces to check sphere before computing the circumsphere - // Add the computed sphere as cache; a vector of spheres depending on the number of faces ? - // -// if (Gudhi::Euclidean_distance()(cast_point_d_to_double(sph.first), point_to_be_inserted) > std::sqrt(cast_to_double(sph.second))) -// FT r = kernel_.squared_distance_d_object()(sph.first, sph.first); //*(points.cend()-1)); - if (kernel_.squared_distance_d_object()(sph.first, point_to_be_inserted) > sph.second) - sph = get_sphere(points.cbegin(), points.cend()); - - Filtration_value radius = std::sqrt(cast_to_double(sph.second)); + if (radius == 0.) { // Spheres of each face don't contain the whole simplex + Point_cloud points; + for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { + points.push_back(cc_ptr_->get_point(vertex)); + } + Sphere sph = get_sphere(points.cbegin(), points.cend()); + radius = std::sqrt(cast_to_double(sph.second)); + std::clog << "GLOBAL SPHERE radius = " << radius << std::endl; + // Set all key_to_permute permutations to sphere in cache + do + { + cache_[key_to_permute] = sph; + } while(std::next_permutation(key_to_permute.begin(), key_to_permute.end())); + } #ifdef DEBUG_TRACES @@ -170,8 +208,7 @@ class Cech_blocker { SimplicialComplexForCech* sc_ptr_; Cech_complex* cc_ptr_; Kernel kernel_; - //Vector_vertex_iterator vertex_handle_to_iterator_; - + std::map cache_; }; } // namespace cech_complex -- cgit v1.2.3 From 839093da012bda0ee16744f1e340a8a8eb04f0af Mon Sep 17 00:00:00 2001 From: Hind-M Date: Mon, 6 Sep 2021 17:19:32 +0200 Subject: Remove unnecessary key permutations storage in cache --- .../example/cech_complex_example_from_points.cpp | 16 ++++- .../include/gudhi/Cech_complex_blocker.h | 70 ++++++++++++++-------- 2 files changed, 59 insertions(+), 27 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/example/cech_complex_example_from_points.cpp b/src/Cech_complex/example/cech_complex_example_from_points.cpp index ac17fc73..e78ad51d 100644 --- a/src/Cech_complex/example/cech_complex_example_from_points.cpp +++ b/src/Cech_complex/example/cech_complex_example_from_points.cpp @@ -43,6 +43,18 @@ int main() { std::vector point5({3. + std::sqrt(3.), 3.}); points.emplace_back(point5.begin(), point5.end()); +/* + std::vector point6({1., 4.}); + points.emplace_back(point6.begin(), point6.end()); + std::vector point7({3., 4.}); + points.emplace_back(point7.begin(), point7.end()); + std::vector point8({2., 4. + std::sqrt(3.)}); + points.emplace_back(point8.begin(), point8.end()); + std::vector point9({0., 4.}); + points.emplace_back(point9.begin(), point9.end()); + std::vector point10({-0.5, 2.}); + points.emplace_back(point10.begin(), point10.end()); +*/ // points.emplace_back(Point(std::vector({1., 0.}))); // points.emplace_back(Point(std::vector({0., 1.}))); // points.emplace_back(Point(std::vector({2., 1.}))); @@ -68,13 +80,13 @@ int main() { // ---------------------------------------------------------------------------- // Init of a Cech complex from points // ---------------------------------------------------------------------------- - Filtration_value max_radius = 10.; + Filtration_value max_radius = 10.; //100.; std::clog << "Hind: Just before the Cech constructor" << std::endl; Cech_complex cech_complex_from_points(points, max_radius); std::clog << "Hind: Just after the Cech constructor" << std::endl; Simplex_tree stree; - cech_complex_from_points.create_complex(stree, 3); + cech_complex_from_points.create_complex(stree, 3); //6 // ---------------------------------------------------------------------------- // Display information about the one skeleton Cech complex // ---------------------------------------------------------------------------- diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index acb53143..0fc76c6d 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -102,7 +102,8 @@ class Cech_blocker { using Point_cloud = std::vector; CGAL::NT_converter cast_to_double; Filtration_value radius = 0.; - std::string key_to_permute; +// std::string key_to_permute; + std::vector faces_keys; // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices // std::set enclosing_ball_radii; @@ -113,12 +114,12 @@ class Cech_blocker { auto longlist = sc_ptr_->simplex_vertex_range(sh); auto shortlist = sc_ptr_->simplex_vertex_range(face); - std::clog << "Hind debug: within FACE loop "<< std::endl; +// std::clog << "Hind debug: within FACE loop "<< std::endl; // TODO to remove - for (auto i = std::begin(longlist); i != std::end(longlist);++i) - std::clog << "Hind debug: longlist: " << cc_ptr_->get_point(*i) << std::endl; - for (auto i = std::begin(shortlist); i != std::end(shortlist);++i) - std::clog << "Hind debug: shortlist: " << cc_ptr_->get_point(*i) << std::endl; +// for (auto i = std::begin(longlist); i != std::end(longlist);++i) +// std::clog << "Hind debug: longlist: " << cc_ptr_->get_point(*i) << std::endl; +// for (auto i = std::begin(shortlist); i != std::end(shortlist);++i) +// std::clog << "Hind debug: shortlist: " << cc_ptr_->get_point(*i) << std::endl; auto longiter = std::begin(longlist); auto shortiter = std::begin(shortlist); @@ -126,7 +127,7 @@ class Cech_blocker { while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; } auto extra = *longiter; // Vertex_handle - std::clog << "Hind debug: extra vertex: " << cc_ptr_->get_point(extra) << std::endl; +// std::clog << "Hind debug: extra vertex: " << cc_ptr_->get_point(extra) << std::endl; Point_cloud face_points; std::string key, key_extra; @@ -139,10 +140,11 @@ class Cech_blocker { } key_extra = key; key_extra.append(std::to_string(extra)); - key_to_permute = key_extra; - std::clog << "END OF VERTICES " << std::endl; - std::clog << "KEY is: " << key << std::endl; - std::clog << "KEY extra is: " << key_extra << std::endl; + faces_keys.push_back(key_extra); +// key_to_permute = key_extra; +// std::clog << "END OF VERTICES " << std::endl; +// std::clog << "KEY is: " << key << std::endl; +// std::clog << "KEY extra is: " << key_extra << std::endl; Sphere sph; auto it = cache_.find(key); if(it != cache_.end()) @@ -155,14 +157,14 @@ class Cech_blocker { #ifdef DEBUG_TRACES std::clog << "circumcenter: " << sph.first << ", radius: " << radius << std::endl; #endif // DEBUG_TRACES - std::clog << "distance FYI: " << kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " < " << cast_to_double(sph.second) << std::endl; +// std::clog << "distance FYI: " << kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " < " << cast_to_double(sph.second) << std::endl; // enclosing_ball_radii.insert(radius); enclosing_ball_spheres.insert(sph); cache_[key_extra] = sph; } - else {// TODO to remove - std::clog << "vertex not included BECAUSE DISTANCE: "<< kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " AND RAD SPHERE: " << sph.second << std::endl; - } +// else {// TODO to remove +// std::clog << "vertex not included BECAUSE DISTANCE: "<< kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " AND RAD SPHERE: " << sph.second << std::endl; +// } } // Get the minimal radius of all faces enclosing balls if exists if (!enclosing_ball_spheres.empty()) { @@ -171,12 +173,21 @@ class Cech_blocker { radius = std::sqrt(cast_to_double(sph_min.second)); // std::clog << "CHECK that radius of min sphere is min radius: " << std::sqrt(cast_to_double(sph_min.second)) << "; and RADIUS min: " << radius << std::endl; // Set all key_to_permute permutations to min sphere in cache - do - { - cache_[key_to_permute] = sph_min; - } while(std::next_permutation(key_to_permute.begin(), key_to_permute.end())); +// do +// { +// if (cache_.find(key_to_permute) != cache_.end()) { +// if (cast_to_double(cache_[key_to_permute].second) > cast_to_double(sph_min.second)) +// cache_[key_to_permute] = sph_min; +// } +// else { +// cache_[key_to_permute] = sph_min; +// } +// } while(std::next_permutation(key_to_permute.begin(), key_to_permute.end())); + for (auto k : faces_keys) { + cache_[k] = sph_min; + } } - std::clog << "END OF FACES ; radius = " << radius << std::endl; +// std::clog << "END OF FACES ; radius = " << radius << std::endl; if (radius == 0.) { // Spheres of each face don't contain the whole simplex Point_cloud points; @@ -185,12 +196,21 @@ class Cech_blocker { } Sphere sph = get_sphere(points.cbegin(), points.cend()); radius = std::sqrt(cast_to_double(sph.second)); - std::clog << "GLOBAL SPHERE radius = " << radius << std::endl; +// std::clog << "GLOBAL SPHERE radius = " << radius << std::endl; // Set all key_to_permute permutations to sphere in cache - do - { - cache_[key_to_permute] = sph; - } while(std::next_permutation(key_to_permute.begin(), key_to_permute.end())); +// do +// { +// // if (cache_.find(key_to_permute) != cache_.end()) { +// // if (cast_to_double(cache_[key_to_permute].second) > cast_to_double(sph.second)) +// // cache_[key_to_permute] = sph; +// // } +// // else { +// // cache_[key_to_permute] = sph; +// // } +// } while(std::next_permutation(key_to_permute.begin(), key_to_permute.end())); +// for (auto k : faces_keys) { +// cache_[k] = sph; +// } } -- cgit v1.2.3 From 44f754ee58aeee043891f4494892798b9807374b Mon Sep 17 00:00:00 2001 From: Hind-M Date: Wed, 8 Sep 2021 14:55:48 +0200 Subject: Change cache so that the index of the stored sphere is used as key --- .../include/gudhi/Cech_complex_blocker.h | 44 ++++++++++++++-------- 1 file changed, 29 insertions(+), 15 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 0fc76c6d..3cac9ee2 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -103,7 +103,7 @@ class Cech_blocker { CGAL::NT_converter cast_to_double; Filtration_value radius = 0.; // std::string key_to_permute; - std::vector faces_keys; +// std::vector faces_keys; // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices // std::set enclosing_ball_radii; @@ -130,28 +130,35 @@ class Cech_blocker { // std::clog << "Hind debug: extra vertex: " << cc_ptr_->get_point(extra) << std::endl; Point_cloud face_points; - std::string key, key_extra; +// std::string key, key_extra; for (auto vertex : sc_ptr_->simplex_vertex_range(face)) { face_points.push_back(cc_ptr_->get_point(vertex)); - key.append(std::to_string(vertex)); +// key.append(std::to_string(vertex)); #ifdef DEBUG_TRACES std::clog << "#(" << vertex << ")#"; #endif // DEBUG_TRACES } - key_extra = key; - key_extra.append(std::to_string(extra)); - faces_keys.push_back(key_extra); +// key_extra = key; +// key_extra.append(std::to_string(extra)); +// faces_keys.push_back(key_extra); // key_to_permute = key_extra; // std::clog << "END OF VERTICES " << std::endl; // std::clog << "KEY is: " << key << std::endl; // std::clog << "KEY extra is: " << key_extra << std::endl; Sphere sph; - auto it = cache_.find(key); - if(it != cache_.end()) - sph = it->second; + auto k = sc_ptr_->key(sh); + if(k != sc_ptr_->null_key()) + sph = cache_[k]; else { sph = get_sphere(face_points.cbegin(), face_points.cend()); } +// auto it = cache_.find(key); +// if(it != cache_.end()) +// sph = it->second; +// else { +// sph = get_sphere(face_points.cbegin(), face_points.cend()); +// } + if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) <= sph.second) { radius = std::sqrt(cast_to_double(sph.second)); #ifdef DEBUG_TRACES @@ -160,7 +167,9 @@ class Cech_blocker { // std::clog << "distance FYI: " << kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " < " << cast_to_double(sph.second) << std::endl; // enclosing_ball_radii.insert(radius); enclosing_ball_spheres.insert(sph); - cache_[key_extra] = sph; +// cache_[key_extra] = sph; +// sc_ptr_->assign_key(sh, cache_.size()); +// cache_.push_back(sph); } // else {// TODO to remove // std::clog << "vertex not included BECAUSE DISTANCE: "<< kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " AND RAD SPHERE: " << sph.second << std::endl; @@ -174,7 +183,7 @@ class Cech_blocker { // std::clog << "CHECK that radius of min sphere is min radius: " << std::sqrt(cast_to_double(sph_min.second)) << "; and RADIUS min: " << radius << std::endl; // Set all key_to_permute permutations to min sphere in cache // do -// { +// { // if (cache_.find(key_to_permute) != cache_.end()) { // if (cast_to_double(cache_[key_to_permute].second) > cast_to_double(sph_min.second)) // cache_[key_to_permute] = sph_min; @@ -183,9 +192,11 @@ class Cech_blocker { // cache_[key_to_permute] = sph_min; // } // } while(std::next_permutation(key_to_permute.begin(), key_to_permute.end())); - for (auto k : faces_keys) { - cache_[k] = sph_min; - } +// for (auto k : faces_keys) { +// cache_[k] = sph_min; +// } + sc_ptr_->assign_key(sh, cache_.size()); + cache_.push_back(sph_min); } // std::clog << "END OF FACES ; radius = " << radius << std::endl; @@ -211,6 +222,8 @@ class Cech_blocker { // for (auto k : faces_keys) { // cache_[k] = sph; // } +// sc_ptr_->assign_key(sh, cache_.size()); +// cache_.push_back(sph); } @@ -228,7 +241,8 @@ class Cech_blocker { SimplicialComplexForCech* sc_ptr_; Cech_complex* cc_ptr_; Kernel kernel_; - std::map cache_; +// std::map cache_; + std::vector cache_; }; } // namespace cech_complex -- cgit v1.2.3 From 767d9fca5da2d3dd9698a5c27e9bedc159271f67 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Thu, 30 Sep 2021 11:05:17 +0200 Subject: Move minimal enclosing balls cache to Cech_complex.h instead of the blocker Modify cech test to be more relevant regarding the current algorithm Do some cleaning --- .../benchmark/cech_complex_benchmark.cpp | 2 +- .../example/cech_complex_example_from_points.cpp | 42 +------ .../example/cech_complex_step_by_step.cpp | 8 -- src/Cech_complex/include/gudhi/Cech_complex.h | 38 +++---- .../include/gudhi/Cech_complex_blocker.h | 121 +++------------------ src/Cech_complex/test/test_cech_complex.cpp | 37 ++----- src/Cech_complex/utilities/cech_persistence.cpp | 4 +- src/Simplex_tree/include/gudhi/Simplex_tree.h | 4 - src/common/include/gudhi/distance_functions.h | 20 +--- .../include/gudhi/graph_simplicial_complex.h | 3 - 10 files changed, 54 insertions(+), 225 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index c332c656..4a1aa06e 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -33,7 +33,7 @@ using Points_off_reader = Gudhi::Points_off_reader; using Proximity_graph = Gudhi::Proximity_graph; using Rips_complex = Gudhi::rips_complex::Rips_complex; using Kernel = CGAL::Epeck_d; -using Cech_complex = Gudhi::cech_complex::Cech_complex; +using Cech_complex = Gudhi::cech_complex::Cech_complex; class Minimal_enclosing_ball_radius { public: diff --git a/src/Cech_complex/example/cech_complex_example_from_points.cpp b/src/Cech_complex/example/cech_complex_example_from_points.cpp index e78ad51d..78861951 100644 --- a/src/Cech_complex/example/cech_complex_example_from_points.cpp +++ b/src/Cech_complex/example/cech_complex_example_from_points.cpp @@ -10,25 +10,15 @@ int main() { // Type definitions -// using Point_cloud = std::vector>; using Simplex_tree = Gudhi::Simplex_tree; using Filtration_value = Simplex_tree::Filtration_value; using Kernel = CGAL::Epeck_d; using FT = typename Kernel::FT; using Point = typename Kernel::Point_d; using Point_cloud = std::vector; - using Cech_complex = Gudhi::cech_complex::Cech_complex; + using Cech_complex = Gudhi::cech_complex::Cech_complex; Point_cloud points; -// points.push_back({1., 0.}); // 0 -// points.push_back({0., 1.}); // 1 -// points.push_back({2., 1.}); // 2 -// points.push_back({3., 2.}); // 3 -// points.push_back({0., 3.}); // 4 -// points.push_back({3. + std::sqrt(3.), 3.}); // 5 - -// std::vector point({0.0, 0.0, 0.0, 0.0}); -// points.emplace_back(point.begin(), point.end()); std::vector point0({1., 0.}); points.emplace_back(point0.begin(), point0.end()); @@ -42,8 +32,6 @@ int main() { points.emplace_back(point4.begin(), point4.end()); std::vector point5({3. + std::sqrt(3.), 3.}); points.emplace_back(point5.begin(), point5.end()); - -/* std::vector point6({1., 4.}); points.emplace_back(point6.begin(), point6.end()); std::vector point7({3., 4.}); @@ -54,39 +42,15 @@ int main() { points.emplace_back(point9.begin(), point9.end()); std::vector point10({-0.5, 2.}); points.emplace_back(point10.begin(), point10.end()); -*/ -// points.emplace_back(Point(std::vector({1., 0.}))); -// points.emplace_back(Point(std::vector({0., 1.}))); -// points.emplace_back(Point(std::vector({2., 1.}))); -// points.emplace_back(Point(std::vector({3., 2.}))); -// points.emplace_back(Point(std::vector({0., 3.}))); -// points.emplace_back(Point(std::vector({3. + std::sqrt(3.), 3.}))); - - -// points.push_back(Point(1.0, 0.0)); -// points.push_back(Point(0.0, 1.0)); -// points.push_back(Point(2.0, 1.0)); -// points.push_back(Point(3.0, 2.0)); -// points.push_back(Point(0.0, 3.0)); -// points.push_back(Point(3.0 + std::sqrt(3.0), 3.0)); - - -// points.push_back({1., 4.}); // 6 -// points.push_back({3., 4.}); // 7 -// points.push_back({2., 4. + std::sqrt(3.)}); // 8 -// points.push_back({0., 4.}); // 9 -// points.push_back({-0.5, 2.}); // 10 // ---------------------------------------------------------------------------- // Init of a Cech complex from points // ---------------------------------------------------------------------------- - Filtration_value max_radius = 10.; //100.; - std::clog << "Hind: Just before the Cech constructor" << std::endl; + Filtration_value max_radius = 100.; //100.; Cech_complex cech_complex_from_points(points, max_radius); - std::clog << "Hind: Just after the Cech constructor" << std::endl; Simplex_tree stree; - cech_complex_from_points.create_complex(stree, 3); //6 + cech_complex_from_points.create_complex(stree, 6); //6 // ---------------------------------------------------------------------------- // Display information about the one skeleton Cech complex // ---------------------------------------------------------------------------- diff --git a/src/Cech_complex/example/cech_complex_step_by_step.cpp b/src/Cech_complex/example/cech_complex_step_by_step.cpp index ac08e6cc..44e7f945 100644 --- a/src/Cech_complex/example/cech_complex_step_by_step.cpp +++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp @@ -13,8 +13,6 @@ #include #include -#include // TODO to remove ? - #include #include @@ -36,7 +34,6 @@ using Simplex_tree = Gudhi::Simplex_tree<>; using Simplex_handle = Simplex_tree::Simplex_handle; using Filtration_value = Simplex_tree::Filtration_value; -// using Point = std::vector; using Kernel = CGAL::Epeck_d; using Point = typename Kernel::Point_d; using Points_off_reader = Gudhi::Points_off_reader; @@ -45,9 +42,6 @@ using Proximity_graph = Gudhi::Proximity_graph; class Cech_blocker { private: using Point_cloud = std::vector; -// using Point_iterator = Point_cloud::const_iterator; -// using Coordinate_iterator = Point::const_iterator; -// using Min_sphere = Gudhi::Miniball::Miniball>; public: bool operator()(Simplex_handle sh) { @@ -67,14 +61,12 @@ class Cech_blocker { } Cech_blocker(Simplex_tree& simplex_tree, Filtration_value max_radius, const std::vector& point_cloud) : simplex_tree_(simplex_tree), max_radius_(max_radius), point_cloud_(point_cloud) { -// dimension_ = point_cloud_[0].size(); } private: Simplex_tree simplex_tree_; Filtration_value max_radius_; std::vector point_cloud_; -// int dimension_; }; void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& max_radius, int& dim_max); diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 2c6f202a..32a78aec 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -40,7 +40,7 @@ namespace cech_complex { * \tparam ForwardPointRange must be a range for which `std::begin()` and `std::end()` methods return input * iterators on a point. `std::begin()` and `std::end()` methods are also required for a point. */ -template +template class Cech_complex { private: // Required by compute_proximity_graph @@ -48,17 +48,17 @@ class Cech_complex { using Filtration_value = typename SimplicialComplexForProximityGraph::Filtration_value; using Proximity_graph = Gudhi::Proximity_graph; - // Retrieve Coordinate type from ForwardPointRange -// using Point_from_range_iterator = typename boost::range_const_iterator::type; -// using Point_from_range = typename std::iterator_traits::value_type; -// using Coordinate_iterator = typename boost::range_const_iterator::type; -// using Coordinate = typename std::iterator_traits::value_type; - public: - // Point and Point_cloud type definition - //using Point = std::vector; - using Point = typename Kernel::Point_d; - using Point_cloud = std::vector; + + using cech_blocker = Cech_blocker; + + using Point_d = typename cech_blocker::Point_d; + using Point_cloud = std::vector; + + // Numeric type of coordinates in the kernel + using FT = typename cech_blocker::FT; + // Sphere is a pair of point and squared radius. + using Sphere = typename std::pair; public: /** \brief Cech_complex constructor from a list of points. @@ -77,7 +77,6 @@ class Cech_complex { point_cloud_.assign(points.begin(), points.end()); - std::clog << "Hind: Just before the graph compute" << std::endl; cech_skeleton_graph_ = Gudhi::compute_proximity_graph( point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); } @@ -90,19 +89,14 @@ class Cech_complex { * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0. * */ - template void create_complex(SimplicialComplexForCechComplex& complex, int dim_max) { - std::clog << "Hind: in create complex" << std::endl; GUDHI_CHECK(complex.num_vertices() == 0, std::invalid_argument("Cech_complex::create_complex - simplicial complex is not empty")); // insert the proximity graph in the simplicial complex - std::clog << "Hind: before insert_graph" << std::endl; complex.insert_graph(cech_skeleton_graph_); // expand the graph until dimension dim_max - std::clog << "Hind: before expansion_with_blockers" << std::endl; - complex.expansion_with_blockers(dim_max, - Cech_blocker(&complex, this)); + complex.expansion_with_blockers(dim_max, cech_blocker(&complex, this)); } /** @return max_radius value given at construction. */ @@ -111,12 +105,18 @@ class Cech_complex { /** @param[in] vertex Point position in the range. * @return The point. */ - const Point& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; } + const Point_d& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; } + + /** + * @return Vector of cached spheres. + */ + std::vector & get_cache() { return cache_; } private: Proximity_graph cech_skeleton_graph_; Filtration_value max_radius_; Point_cloud point_cloud_; + std::vector cache_; }; } // namespace cech_complex diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 3cac9ee2..5edd005d 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -11,20 +11,11 @@ #ifndef CECH_COMPLEX_BLOCKER_H_ #define CECH_COMPLEX_BLOCKER_H_ -// TODO to remove -#include // for Gudhi::Minimal_enclosing_ball_radius - -#include // for CGAL::to_double -#include // -// #include -#include // For EXACT or SAFE version -#include +#include // for casting from FT to double #include #include #include -#include -#include #include // for std::sqrt namespace Gudhi { @@ -48,7 +39,6 @@ namespace cech_complex { template class Cech_blocker { private: -// using Point_cloud = typename Cech_complex::Point_cloud; using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; using Filtration_value = typename SimplicialComplexForCech::Filtration_value; @@ -61,28 +51,18 @@ class Cech_blocker { // Sphere is a pair of point and squared radius. using Sphere = typename std::pair; - - /** \internal \brief TODO - * \param[in] - * \return */ template FT get_squared_radius(PointIterator begin, PointIterator end) const { return kernel_.compute_squared_radius_d_object()(begin, end); } - - /** \internal \brief TODO - * \param[in] - * \return */ + template Sphere get_sphere(PointIterator begin, PointIterator end) const { Point_d c = kernel_.construct_circumcenter_d_object()(begin, end); FT r = kernel_.squared_distance_d_object()(c, *begin); return std::make_pair(std::move(c), std::move(r)); } - - /** \internal \brief TODO - * \param[in] - * \return */ + template class CompareSpheresRadii { @@ -93,7 +73,7 @@ class Cech_blocker { return cast_to_double(firstSphere.second) < cast_to_double(secondSphere.second); } }; - + /** \internal \brief Čech complex blocker operator() - the oracle - assigns the filtration value from the simplex * radius and returns if the simplex expansion must be blocked. * \param[in] sh The Simplex_handle. @@ -102,104 +82,54 @@ class Cech_blocker { using Point_cloud = std::vector; CGAL::NT_converter cast_to_double; Filtration_value radius = 0.; -// std::string key_to_permute; -// std::vector faces_keys; - + // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices - // std::set enclosing_ball_radii; std::set > enclosing_ball_spheres; for (auto face : sc_ptr_->boundary_simplex_range(sh)) { - // Find which vertex of sh is missing in face. We rely on the fact that simplex_vertex_range is sorted. auto longlist = sc_ptr_->simplex_vertex_range(sh); auto shortlist = sc_ptr_->simplex_vertex_range(face); -// std::clog << "Hind debug: within FACE loop "<< std::endl; - // TODO to remove -// for (auto i = std::begin(longlist); i != std::end(longlist);++i) -// std::clog << "Hind debug: longlist: " << cc_ptr_->get_point(*i) << std::endl; -// for (auto i = std::begin(shortlist); i != std::end(shortlist);++i) -// std::clog << "Hind debug: shortlist: " << cc_ptr_->get_point(*i) << std::endl; - auto longiter = std::begin(longlist); auto shortiter = std::begin(shortlist); auto enditer = std::end(shortlist); while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; } auto extra = *longiter; // Vertex_handle -// std::clog << "Hind debug: extra vertex: " << cc_ptr_->get_point(extra) << std::endl; - Point_cloud face_points; -// std::string key, key_extra; for (auto vertex : sc_ptr_->simplex_vertex_range(face)) { face_points.push_back(cc_ptr_->get_point(vertex)); -// key.append(std::to_string(vertex)); #ifdef DEBUG_TRACES std::clog << "#(" << vertex << ")#"; #endif // DEBUG_TRACES } -// key_extra = key; -// key_extra.append(std::to_string(extra)); -// faces_keys.push_back(key_extra); -// key_to_permute = key_extra; -// std::clog << "END OF VERTICES " << std::endl; -// std::clog << "KEY is: " << key << std::endl; -// std::clog << "KEY extra is: " << key_extra << std::endl; Sphere sph; - auto k = sc_ptr_->key(sh); - if(k != sc_ptr_->null_key()) - sph = cache_[k]; + auto face_sh = sc_ptr_->find(sc_ptr_->simplex_vertex_range(face)); + auto k = sc_ptr_->key(face_sh); + if(k != sc_ptr_->null_key()) { + sph = cc_ptr_->get_cache().at(k); + } else { sph = get_sphere(face_points.cbegin(), face_points.cend()); } -// auto it = cache_.find(key); -// if(it != cache_.end()) -// sph = it->second; -// else { -// sph = get_sphere(face_points.cbegin(), face_points.cend()); -// } if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) <= sph.second) { radius = std::sqrt(cast_to_double(sph.second)); #ifdef DEBUG_TRACES std::clog << "circumcenter: " << sph.first << ", radius: " << radius << std::endl; #endif // DEBUG_TRACES -// std::clog << "distance FYI: " << kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " < " << cast_to_double(sph.second) << std::endl; - // enclosing_ball_radii.insert(radius); enclosing_ball_spheres.insert(sph); -// cache_[key_extra] = sph; -// sc_ptr_->assign_key(sh, cache_.size()); -// cache_.push_back(sph); } -// else {// TODO to remove -// std::clog << "vertex not included BECAUSE DISTANCE: "<< kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " AND RAD SPHERE: " << sph.second << std::endl; -// } } // Get the minimal radius of all faces enclosing balls if exists if (!enclosing_ball_spheres.empty()) { - // radius = *enclosing_ball_radii.begin(); Sphere sph_min = *enclosing_ball_spheres.begin(); radius = std::sqrt(cast_to_double(sph_min.second)); - // std::clog << "CHECK that radius of min sphere is min radius: " << std::sqrt(cast_to_double(sph_min.second)) << "; and RADIUS min: " << radius << std::endl; - // Set all key_to_permute permutations to min sphere in cache -// do -// { -// if (cache_.find(key_to_permute) != cache_.end()) { -// if (cast_to_double(cache_[key_to_permute].second) > cast_to_double(sph_min.second)) -// cache_[key_to_permute] = sph_min; -// } -// else { -// cache_[key_to_permute] = sph_min; -// } -// } while(std::next_permutation(key_to_permute.begin(), key_to_permute.end())); -// for (auto k : faces_keys) { -// cache_[k] = sph_min; -// } - sc_ptr_->assign_key(sh, cache_.size()); - cache_.push_back(sph_min); + + sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); + cc_ptr_->get_cache().push_back(sph_min); } -// std::clog << "END OF FACES ; radius = " << radius << std::endl; - + if (radius == 0.) { // Spheres of each face don't contain the whole simplex Point_cloud points; for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { @@ -207,25 +137,10 @@ class Cech_blocker { } Sphere sph = get_sphere(points.cbegin(), points.cend()); radius = std::sqrt(cast_to_double(sph.second)); -// std::clog << "GLOBAL SPHERE radius = " << radius << std::endl; - // Set all key_to_permute permutations to sphere in cache -// do -// { -// // if (cache_.find(key_to_permute) != cache_.end()) { -// // if (cast_to_double(cache_[key_to_permute].second) > cast_to_double(sph.second)) -// // cache_[key_to_permute] = sph; -// // } -// // else { -// // cache_[key_to_permute] = sph; -// // } -// } while(std::next_permutation(key_to_permute.begin(), key_to_permute.end())); -// for (auto k : faces_keys) { -// cache_[k] = sph; -// } -// sc_ptr_->assign_key(sh, cache_.size()); -// cache_.push_back(sph); - } + sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); + cc_ptr_->get_cache().push_back(sph); + } #ifdef DEBUG_TRACES if (radius > cc_ptr_->max_radius()) std::clog << "radius > max_radius => expansion is blocked\n"; @@ -241,8 +156,6 @@ class Cech_blocker { SimplicialComplexForCech* sc_ptr_; Cech_complex* cc_ptr_; Kernel kernel_; -// std::map cache_; - std::vector cache_; }; } // namespace cech_complex diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 81efd6ae..51b466da 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -24,26 +24,19 @@ #include #include #include -#include #include // For EXACT or SAFE version - // Type definitions using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; -//using Point = std::vector; using Kernel = CGAL::Epeck_d; using FT = typename Kernel::FT; using Point = typename Kernel::Point_d; using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; -using Cech_complex = Gudhi::cech_complex::Cech_complex; - -// using Point_iterator = Point_cloud::const_iterator; -// using Coordinate_iterator = Point::const_iterator; -// using Min_sphere = Gudhi::Miniball::Miniball>; +using Cech_complex = Gudhi::cech_complex::Cech_complex; BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { // ---------------------------------------------------------------------------- @@ -52,18 +45,6 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { // // ---------------------------------------------------------------------------- Point_cloud points; -// points.push_back({1., 0.}); // 0 -// points.push_back({0., 1.}); // 1 -// points.push_back({2., 1.}); // 2 -// points.push_back({3., 2.}); // 3 -// points.push_back({0., 3.}); // 4 -// points.push_back({3. + std::sqrt(3.), 3.}); // 5 -// points.push_back({1., 4.}); // 6 -// points.push_back({3., 4.}); // 7 -// points.push_back({2., 4. + std::sqrt(3.)}); // 8 -// points.push_back({0., 4.}); // 9 -// points.push_back({-0.5, 2.}); // 10 - std::vector point0({1., 0.}); points.emplace_back(point0.begin(), point0.end()); @@ -156,36 +137,32 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { for (std::size_t vertex = 0; vertex <= 2; vertex++) { points012.push_back(cech_complex_for_doc.get_point(vertex)); } -// std::size_t dimension = points[0].end() - points[0].begin(); -// Min_sphere ms012(dimension, points012.begin(), points012.end()); Kernel kern; Simplex_tree::Filtration_value f012 = st2.filtration(st2.find({0, 1, 2})); - CGAL::NT_converter cast_to_double; - std::clog << "f012= " << f012 << " | points012_radius= " << std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points012.begin(), points012.end()))) << std::endl; - + std::clog << "f012= " << f012 << std::endl; + CGAL::NT_converter cast_to_double; GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points012.begin(), points012.end())))); Point_cloud points1410; points1410.push_back(cech_complex_for_doc.get_point(1)); points1410.push_back(cech_complex_for_doc.get_point(4)); points1410.push_back(cech_complex_for_doc.get_point(10)); -// Min_sphere ms1410(dimension, points1410.begin(), points1410.end()); Simplex_tree::Filtration_value f1410 = st2.filtration(st2.find({1, 4, 10})); - std::clog << "f1410= " << f1410 << " | points1410_radius= " << std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points1410.begin(), points1410.end()))) << std::endl; + std::clog << "f1410= " << f1410 << std::endl; - GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points1410.begin(), points1410.end())))); + // In this case, the computed sphere using CGAL kernel does not match the minimal enclosing ball; the filtration value check is therefore done against a hardcoded value + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, 1.); Point_cloud points469; points469.push_back(cech_complex_for_doc.get_point(4)); points469.push_back(cech_complex_for_doc.get_point(6)); points469.push_back(cech_complex_for_doc.get_point(9)); -// Min_sphere ms469(dimension, points469.begin(), points469.end()); Simplex_tree::Filtration_value f469 = st2.filtration(st2.find({4, 6, 9})); - std::clog << "f469= " << f469 << " | points469_radius= " << std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points469.begin(), points469.end()))) << std::endl; + std::clog << "f469= " << f469 << std::endl; GUDHI_TEST_FLOAT_EQUALITY_CHECK(f469, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points469.begin(), points469.end())))); diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp index ccd7d453..0c945cad 100644 --- a/src/Cech_complex/utilities/cech_persistence.cpp +++ b/src/Cech_complex/utilities/cech_persistence.cpp @@ -25,14 +25,12 @@ // Types definition using Simplex_tree = Gudhi::Simplex_tree; using Filtration_value = Simplex_tree::Filtration_value; -// using Point = std::vector; -// using Point_cloud = std::vector; using Kernel = CGAL::Epeck_d; using Point = typename Kernel::Point_d; using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; -using Cech_complex = Gudhi::cech_complex::Cech_complex; +using Cech_complex = Gudhi::cech_complex::Cech_complex; using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index f69ed6ec..85790baf 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1278,10 +1278,8 @@ class Simplex_tree { intersection.emplace_back(next->first, Node(nullptr, filt)); } } - std::clog << "Hind: after intersection insertion" << std::endl; if (intersection.size() != 0) { // Reverse the order to insert - std::clog << "Hind: declare new siblings" << std::endl; Siblings * new_sib = new Siblings(siblings, // oncles simplex->first, // parent boost::adaptors::reverse(intersection)); // boost::container::ordered_unique_range_t @@ -1290,12 +1288,10 @@ class Simplex_tree { for (auto new_sib_member = new_sib->members().begin(); new_sib_member != new_sib->members().end(); new_sib_member++) { - std::clog << "Hind: check the blocker result" << std::endl; bool blocker_result = block_simplex(new_sib_member); // new_sib member has been blocked by the blocker function // add it to the list to be removed - do not perform it while looping on it if (blocker_result) { - std::clog << "Hind: add to list of blocked sib to be removed" << std::endl; blocked_new_sib_vertex_list.push_back(new_sib_member->first); } } diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index ae5168aa..a8ee4a75 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -65,19 +65,17 @@ class Euclidean_distance { * The points are assumed to have the same dimension. */ class Minimal_enclosing_ball_radius { public: - /** \brief TODO + /** \brief Enclosing ball radius from two points using CGAL. * * @param[in] point_1 * @param[in] point_2 - * @return - * \tparam Point + * @return Enclosing ball radius for the two points. + * \tparam Point must be a Kernel::Point_d from CGAL. * */ - //typename FT = typename Kernel::FT, template< typename Kernel = CGAL::Epeck_d, typename Point= typename Kernel::Point_d> double operator()(const Point& point_1, const Point& point_2) const { - std::clog << "Added template: distance betw points 1 and 2" << std::endl; Kernel kernel_; return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; } @@ -94,25 +92,21 @@ class Minimal_enclosing_ball_radius { template< typename Point > typename std::iterator_traits::type>::value_type operator()(const Point& point_1, const Point& point_2) const { - std::clog << "Hind: Minimal_enclosing_ball_radius point1 et 2; Euclidean" << std::endl; std::clog << "#" << *point_1.begin() << "##" << *point_2.begin() << std::endl; return Euclidean_distance()(point_1, point_2) / 2.; } - - /** \brief TODO + /** \brief Enclosing ball radius from a point cloud using CGAL. * * @param[in] point_cloud The points. - * @return - * \tparam Point_cloud + * @return Enclosing ball radius for the points. + * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. * */ - //typename FT = typename Kernel::FT, template< typename Kernel = CGAL::Epeck_d, typename Point= typename Kernel::Point_d, typename Point_cloud = std::vector> double operator()(const Point_cloud& point_cloud) const { - std::clog << "Added template: distance in point cloud" << std::endl; Kernel kernel_; return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); } @@ -133,8 +127,6 @@ class Minimal_enclosing_ball_radius { typename Coordinate = typename std::iterator_traits::value_type> Coordinate operator()(const Point_cloud& point_cloud) const { - std::clog << "Hind: Minimal_enclosing_ball_radius point cloud; Miniball" << std::endl; - using Min_sphere = Miniball::Miniball>; Min_sphere ms(boost::size(*point_cloud.begin()), point_cloud.begin(), point_cloud.end()); diff --git a/src/common/include/gudhi/graph_simplicial_complex.h b/src/common/include/gudhi/graph_simplicial_complex.h index 9190182c..da9dee7d 100644 --- a/src/common/include/gudhi/graph_simplicial_complex.h +++ b/src/common/include/gudhi/graph_simplicial_complex.h @@ -18,8 +18,6 @@ #include #include // for std::tie -#include - namespace Gudhi { /** @file * @brief Graph simplicial complex methods @@ -78,7 +76,6 @@ Proximity_graph compute_proximity_graph( for (auto it_u = points.begin(); it_u != points.end(); ++it_u) { idx_v = idx_u + 1; for (auto it_v = it_u + 1; it_v != points.end(); ++it_v, ++idx_v) { - std::clog << "#idx_u" << idx_u << "#idx_v " << idx_v << std::endl; fil = distance(*it_u, *it_v); if (fil <= threshold) { edges.emplace_back(idx_u, idx_v); -- cgit v1.2.3 From d3970dbbc16993d348092899eb8fcd1ea1aceb8d Mon Sep 17 00:00:00 2001 From: Hind-M Date: Fri, 1 Oct 2021 15:47:58 +0200 Subject: Separate Minimal_enclosing_ball_radius class which uses CGAL from common distance_functions.h file --- .../benchmark/cech_complex_benchmark.cpp | 2 +- .../example/cech_complex_step_by_step.cpp | 2 +- src/Cech_complex/include/gudhi/Cech_complex.h | 2 +- .../include/gudhi/Cech_complex/Cech_kernel.h | 149 +++++++++++++++++++++ src/Cech_complex/test/test_cech_complex.cpp | 2 +- src/Cech_complex/utilities/cech_persistence.cpp | 2 +- src/common/include/gudhi/distance_functions.h | 83 ------------ 7 files changed, 154 insertions(+), 88 deletions(-) create mode 100644 src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 4a1aa06e..cfeb0725 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -9,7 +9,7 @@ */ #include -#include +#include #include #include #include diff --git a/src/Cech_complex/example/cech_complex_step_by_step.cpp b/src/Cech_complex/example/cech_complex_step_by_step.cpp index 44e7f945..2d8321b1 100644 --- a/src/Cech_complex/example/cech_complex_step_by_step.cpp +++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp @@ -9,7 +9,7 @@ */ #include -#include +#include #include #include diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 32a78aec..7bbf97d1 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -11,7 +11,7 @@ #ifndef CECH_COMPLEX_H_ #define CECH_COMPLEX_H_ -#include // for Gudhi::Minimal_enclosing_ball_radius +#include // for Gudhi::Minimal_enclosing_ball_radius #include // for Gudhi::Proximity_graph #include // for GUDHI_CHECK #include // for Gudhi::cech_complex::Cech_blocker diff --git a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h new file mode 100644 index 00000000..93af90d2 --- /dev/null +++ b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h @@ -0,0 +1,149 @@ +/* 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): Hind Montassif + * + * Copyright (C) 2021 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CECH_KERNEL_H_ +#define CECH_KERNEL_H_ + +#include + +#include // for std::sqrt + +namespace Gudhi { + +// namespace cech_complex { + +/** @brief Compute the radius of the minimal enclosing ball between Points given by a range of coordinates. + * The points are assumed to have the same dimension. */ +class Minimal_enclosing_ball_radius { + public: + /** \brief Enclosing ball radius from two points using CGAL. + * + * @param[in] point_1 + * @param[in] point_2 + * @return Enclosing ball radius for the two points. + * \tparam Point must be a Kernel::Point_d from CGAL. + * + */ + template< typename Kernel = CGAL::Epeck_d, + typename Point= typename Kernel::Point_d> + double operator()(const Point& point_1, const Point& point_2) const { + Kernel kernel_; + return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; + } + + /** \brief Minimal_enclosing_ball_radius from two points. + * + * @param[in] point_1 First point. + * @param[in] point_2 second point. + * @return The minimal enclosing ball radius for the two points (aka. Euclidean distance / 2.). + * + * \tparam Point must be a range of Cartesian coordinates. + * + */ +// template< typename Point > +// typename std::iterator_traits::type>::value_type +// operator()(const Point& point_1, const Point& point_2) const { +// std::clog << "#" << *point_1.begin() << "##" << *point_2.begin() << std::endl; +// return Euclidean_distance()(point_1, point_2) / 2.; +// } + + /** \brief Enclosing ball radius from a point cloud using CGAL. + * + * @param[in] point_cloud The points. + * @return Enclosing ball radius for the points. + * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. + * + */ + template< typename Kernel = CGAL::Epeck_d, + typename Point= typename Kernel::Point_d, + typename Point_cloud = std::vector> + double operator()(const Point_cloud& point_cloud) const { + Kernel kernel_; + return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); + } + + /** \brief Minimal_enclosing_ball_radius from a point cloud. + * + * @param[in] point_cloud The points. + * @return The minimal enclosing ball radius for the points. + * + * \tparam Point_cloud must be a range of points with Cartesian coordinates. + * Point_cloud is a range over a range of Coordinate. + * + */ +// template< typename Point_cloud, +// typename Point_iterator = typename boost::range_const_iterator::type, +// typename Point = typename std::iterator_traits::value_type, +// typename Coordinate_iterator = typename boost::range_const_iterator::type, +// typename Coordinate = typename std::iterator_traits::value_type> +// Coordinate +// operator()(const Point_cloud& point_cloud) const { +// using Min_sphere = Miniball::Miniball>; +// +// Min_sphere ms(boost::size(*point_cloud.begin()), point_cloud.begin(), point_cloud.end()); +// #ifdef DEBUG_TRACES +// std::clog << "Minimal_enclosing_ball_radius = " << std::sqrt(ms.squared_radius()) << " | nb points = " +// << boost::size(point_cloud) << " | dimension = " +// << boost::size(*point_cloud.begin()) << std::endl; +// #endif // DEBUG_TRACES +// +// return std::sqrt(ms.squared_radius()); +// } +}; + +/** + * \class Cech_kernel + * \brief Cech complex kernel container. + * + * \details + * The Cech complex kernel container stores CGAL Kernel and dispatch basic computations. + */ + +// template < typename Kernel > +// class Cech_kernel { +// private: +// // Kernel for functions access. +// Kernel kernel_; +// public: +// using Point_d = typename Kernel::Point_d; +// // Numeric type of coordinates in the kernel +// using FT = typename Kernel::FT; +// // Sphere is a pair of point and squared radius. +// using Sphere = typename std::pair; +// +// int get_dimension(const Point_d& p0) const { +// return kernel_.point_dimension_d_object()(p0); +// } +// +// template +// Sphere get_sphere(PointIterator begin, PointIterator end) const { +// Point_d c = kernel_.construct_circumcenter_d_object()(begin, end); +// FT r = kernel_.squared_distance_d_object()(c, *begin); +// return std::make_pair(std::move(c), std::move(r)); +// } +// +// template +// FT get_squared_radius(PointIterator begin, PointIterator end) const { +// return kernel_.compute_squared_radius_d_object()(begin, end); +// } +// +// FT get_squared_radius(const Sphere& sph) const { +// return sph.second; +// } +// }; + + +//} // namespace cech_complex + +// namespace cechcomplex = cech_complex; + +} // namespace Gudhi + +#endif // CECH_KERNEL_H_ diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 51b466da..7d8c3c22 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -22,7 +22,7 @@ // to construct Cech_complex from a OFF file of points #include #include -#include +#include #include #include // For EXACT or SAFE version diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp index 0c945cad..ccf63e3e 100644 --- a/src/Cech_complex/utilities/cech_persistence.cpp +++ b/src/Cech_complex/utilities/cech_persistence.cpp @@ -9,7 +9,7 @@ */ #include -#include +#include #include #include #include diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index a8ee4a75..5e5a1e31 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -13,13 +13,9 @@ #include -#include - #include #include -#include - #include // for std::sqrt #include // for std::decay #include // for std::begin, std::end @@ -61,85 +57,6 @@ class Euclidean_distance { } }; -/** @brief Compute the radius of the minimal enclosing ball between Points given by a range of coordinates. - * The points are assumed to have the same dimension. */ -class Minimal_enclosing_ball_radius { - public: - /** \brief Enclosing ball radius from two points using CGAL. - * - * @param[in] point_1 - * @param[in] point_2 - * @return Enclosing ball radius for the two points. - * \tparam Point must be a Kernel::Point_d from CGAL. - * - */ - template< typename Kernel = CGAL::Epeck_d, - typename Point= typename Kernel::Point_d> - double operator()(const Point& point_1, const Point& point_2) const { - Kernel kernel_; - return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; - } - - /** \brief Minimal_enclosing_ball_radius from two points. - * - * @param[in] point_1 First point. - * @param[in] point_2 second point. - * @return The minimal enclosing ball radius for the two points (aka. Euclidean distance / 2.). - * - * \tparam Point must be a range of Cartesian coordinates. - * - */ - template< typename Point > - typename std::iterator_traits::type>::value_type - operator()(const Point& point_1, const Point& point_2) const { - std::clog << "#" << *point_1.begin() << "##" << *point_2.begin() << std::endl; - return Euclidean_distance()(point_1, point_2) / 2.; - } - - /** \brief Enclosing ball radius from a point cloud using CGAL. - * - * @param[in] point_cloud The points. - * @return Enclosing ball radius for the points. - * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. - * - */ - template< typename Kernel = CGAL::Epeck_d, - typename Point= typename Kernel::Point_d, - typename Point_cloud = std::vector> - double operator()(const Point_cloud& point_cloud) const { - Kernel kernel_; - return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); - } - - /** \brief Minimal_enclosing_ball_radius from a point cloud. - * - * @param[in] point_cloud The points. - * @return The minimal enclosing ball radius for the points. - * - * \tparam Point_cloud must be a range of points with Cartesian coordinates. - * Point_cloud is a range over a range of Coordinate. - * - */ - template< typename Point_cloud, - typename Point_iterator = typename boost::range_const_iterator::type, - typename Point = typename std::iterator_traits::value_type, - typename Coordinate_iterator = typename boost::range_const_iterator::type, - typename Coordinate = typename std::iterator_traits::value_type> - Coordinate - operator()(const Point_cloud& point_cloud) const { - using Min_sphere = Miniball::Miniball>; - - Min_sphere ms(boost::size(*point_cloud.begin()), point_cloud.begin(), point_cloud.end()); -#ifdef DEBUG_TRACES - std::clog << "Minimal_enclosing_ball_radius = " << std::sqrt(ms.squared_radius()) << " | nb points = " - << boost::size(point_cloud) << " | dimension = " - << boost::size(*point_cloud.begin()) << std::endl; -#endif // DEBUG_TRACES - - return std::sqrt(ms.squared_radius()); - } -}; - } // namespace Gudhi #endif // DISTANCE_FUNCTIONS_H_ -- cgit v1.2.3 From 7845631fc2b4a32d4ca7c3e37bf49d40c22e226e Mon Sep 17 00:00:00 2001 From: Hind-M Date: Tue, 12 Oct 2021 10:08:10 +0200 Subject: Remove old implementation of Minimal_enclosing_ball_radius operator () --- .../include/gudhi/Cech_complex/Cech_kernel.h | 42 ---------------------- 1 file changed, 42 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h index 93af90d2..348bb57d 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h +++ b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h @@ -38,21 +38,6 @@ class Minimal_enclosing_ball_radius { return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; } - /** \brief Minimal_enclosing_ball_radius from two points. - * - * @param[in] point_1 First point. - * @param[in] point_2 second point. - * @return The minimal enclosing ball radius for the two points (aka. Euclidean distance / 2.). - * - * \tparam Point must be a range of Cartesian coordinates. - * - */ -// template< typename Point > -// typename std::iterator_traits::type>::value_type -// operator()(const Point& point_1, const Point& point_2) const { -// std::clog << "#" << *point_1.begin() << "##" << *point_2.begin() << std::endl; -// return Euclidean_distance()(point_1, point_2) / 2.; -// } /** \brief Enclosing ball radius from a point cloud using CGAL. * @@ -69,33 +54,6 @@ class Minimal_enclosing_ball_radius { return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); } - /** \brief Minimal_enclosing_ball_radius from a point cloud. - * - * @param[in] point_cloud The points. - * @return The minimal enclosing ball radius for the points. - * - * \tparam Point_cloud must be a range of points with Cartesian coordinates. - * Point_cloud is a range over a range of Coordinate. - * - */ -// template< typename Point_cloud, -// typename Point_iterator = typename boost::range_const_iterator::type, -// typename Point = typename std::iterator_traits::value_type, -// typename Coordinate_iterator = typename boost::range_const_iterator::type, -// typename Coordinate = typename std::iterator_traits::value_type> -// Coordinate -// operator()(const Point_cloud& point_cloud) const { -// using Min_sphere = Miniball::Miniball>; -// -// Min_sphere ms(boost::size(*point_cloud.begin()), point_cloud.begin(), point_cloud.end()); -// #ifdef DEBUG_TRACES -// std::clog << "Minimal_enclosing_ball_radius = " << std::sqrt(ms.squared_radius()) << " | nb points = " -// << boost::size(point_cloud) << " | dimension = " -// << boost::size(*point_cloud.begin()) << std::endl; -// #endif // DEBUG_TRACES -// -// return std::sqrt(ms.squared_radius()); -// } }; /** -- cgit v1.2.3 From 9db268b5ecf056b87ee2f66c6d3f83de93a8681f Mon Sep 17 00:00:00 2001 From: Hind-M Date: Wed, 29 Dec 2021 15:38:59 +0100 Subject: Get min sphere without using a set Move face_points declaration Remove face_sh redundant variable --- .../include/gudhi/Cech_complex_blocker.h | 30 ++++++++-------------- 1 file changed, 11 insertions(+), 19 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 5edd005d..f7f86534 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -63,16 +63,6 @@ class Cech_blocker { return std::make_pair(std::move(c), std::move(r)); } - template - class CompareSpheresRadii - { - public: - CGAL::NT_converter cast_to_double; - bool operator()(const Sphere& firstSphere, const Sphere& secondSphere) - { - return cast_to_double(firstSphere.second) < cast_to_double(secondSphere.second); - } - }; /** \internal \brief Čech complex blocker operator() - the oracle - assigns the filtration value from the simplex * radius and returns if the simplex expansion must be blocked. @@ -84,7 +74,10 @@ class Cech_blocker { Filtration_value radius = 0.; // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices - std::set > enclosing_ball_spheres; + Sphere min_enclos_ball; + CGAL::NT_converter cast_to_FT; + min_enclos_ball.second = cast_to_FT(std::numeric_limits::max()); + Point_cloud face_points; for (auto face : sc_ptr_->boundary_simplex_range(sh)) { // Find which vertex of sh is missing in face. We rely on the fact that simplex_vertex_range is sorted. auto longlist = sc_ptr_->simplex_vertex_range(sh); @@ -96,7 +89,6 @@ class Cech_blocker { while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; } auto extra = *longiter; // Vertex_handle - Point_cloud face_points; for (auto vertex : sc_ptr_->simplex_vertex_range(face)) { face_points.push_back(cc_ptr_->get_point(vertex)); #ifdef DEBUG_TRACES @@ -104,30 +96,30 @@ class Cech_blocker { #endif // DEBUG_TRACES } Sphere sph; - auto face_sh = sc_ptr_->find(sc_ptr_->simplex_vertex_range(face)); - auto k = sc_ptr_->key(face_sh); + auto k = sc_ptr_->key(face); if(k != sc_ptr_->null_key()) { sph = cc_ptr_->get_cache().at(k); } else { sph = get_sphere(face_points.cbegin(), face_points.cend()); } + face_points.clear(); if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) <= sph.second) { radius = std::sqrt(cast_to_double(sph.second)); #ifdef DEBUG_TRACES std::clog << "circumcenter: " << sph.first << ", radius: " << radius << std::endl; #endif // DEBUG_TRACES - enclosing_ball_spheres.insert(sph); + if (cast_to_double(sph.second) < cast_to_double(min_enclos_ball.second)) + min_enclos_ball = sph; } } // Get the minimal radius of all faces enclosing balls if exists - if (!enclosing_ball_spheres.empty()) { - Sphere sph_min = *enclosing_ball_spheres.begin(); - radius = std::sqrt(cast_to_double(sph_min.second)); + if(cast_to_double(min_enclos_ball.second) != std::numeric_limits::max()) { + radius = std::sqrt(cast_to_double(min_enclos_ball.second)); sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); - cc_ptr_->get_cache().push_back(sph_min); + cc_ptr_->get_cache().push_back(min_enclos_ball); } if (radius == 0.) { // Spheres of each face don't contain the whole simplex -- cgit v1.2.3 From b1f40dd2c4397c1975533c54a54538160c727d55 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Thu, 6 Jan 2022 11:39:05 +0100 Subject: Make kernel a parameter of Minimal_enclosing_ball_radius class Use Epick in cech benchmark instead of Epeck --- src/Cech_complex/benchmark/cech_complex_benchmark.cpp | 8 ++++---- src/Cech_complex/include/gudhi/Cech_complex.h | 3 +-- .../include/gudhi/Cech_complex/Cech_kernel.h | 16 ++++++++-------- src/Cech_complex/test/test_cech_complex.cpp | 4 ++-- 4 files changed, 15 insertions(+), 16 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 06d90757..e715b513 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -17,7 +17,7 @@ #include #include -#include // For EXACT or SAFE version +#include #include "boost/filesystem.hpp" // includes all needed Boost.Filesystem declarations @@ -32,7 +32,7 @@ using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; using Proximity_graph = Gudhi::Proximity_graph; using Rips_complex = Gudhi::rips_complex::Rips_complex; -using Kernel = CGAL::Epeck_d; +using Kernel = CGAL::Epick_d>; using Point_cgal = typename Kernel::Point_d; using Point_cloud_cgal = std::vector; using Points_off_reader_cgal = Gudhi::Points_off_reader; @@ -86,7 +86,7 @@ int main(int argc, char* argv[]) { Gudhi::Clock cgal_miniball_clock("Gudhi::Minimal_enclosing_ball_radius_cgal()"); // Compute the proximity graph of the points Proximity_graph cgal_miniball_prox_graph = Gudhi::compute_proximity_graph( - off_reader_cgal.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius()); + off_reader_cgal.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius()); std::clog << cgal_miniball_clock << std::endl; boost::filesystem::path full_path(boost::filesystem::current_path()); @@ -109,7 +109,7 @@ int main(int argc, char* argv[]) { std::clog << radius << ";"; Gudhi::Clock rips_clock("Rips computation"); Rips_complex rips_complex_from_points(off_reader_cgal.get_point_cloud(), radius, - Gudhi::Minimal_enclosing_ball_radius()); + Gudhi::Minimal_enclosing_ball_radius()); Simplex_tree rips_stree; rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); // ------------------------------------------ diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 7bbf97d1..0031d861 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -18,7 +18,6 @@ #include #include // for exception management -#include namespace Gudhi { @@ -78,7 +77,7 @@ class Cech_complex { point_cloud_.assign(points.begin(), points.end()); cech_skeleton_graph_ = Gudhi::compute_proximity_graph( - point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); + point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); } /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal diff --git a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h index 348bb57d..89012206 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h +++ b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h @@ -11,9 +11,10 @@ #ifndef CECH_KERNEL_H_ #define CECH_KERNEL_H_ -#include +#include // for #include #include // for std::sqrt +#include namespace Gudhi { @@ -21,8 +22,14 @@ namespace Gudhi { /** @brief Compute the radius of the minimal enclosing ball between Points given by a range of coordinates. * The points are assumed to have the same dimension. */ +template class Minimal_enclosing_ball_radius { + private: + Kernel kernel_; public: + using Point = typename Kernel::Point_d; + using Point_cloud = typename std::vector; + /** \brief Enclosing ball radius from two points using CGAL. * * @param[in] point_1 @@ -31,10 +38,7 @@ class Minimal_enclosing_ball_radius { * \tparam Point must be a Kernel::Point_d from CGAL. * */ - template< typename Kernel = CGAL::Epeck_d, - typename Point= typename Kernel::Point_d> double operator()(const Point& point_1, const Point& point_2) const { - Kernel kernel_; return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; } @@ -46,11 +50,7 @@ class Minimal_enclosing_ball_radius { * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. * */ - template< typename Kernel = CGAL::Epeck_d, - typename Point= typename Kernel::Point_d, - typename Point_cloud = std::vector> double operator()(const Point_cloud& point_cloud) const { - Kernel kernel_; return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); } diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 7d8c3c22..ca7a9778 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -108,11 +108,11 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { std::clog << vertex << ","; vp.push_back(points.at(vertex)); } - std::clog << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1)) + std::clog << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1)) << " - filtration =" << st.filtration(f_simplex) << std::endl; BOOST_CHECK(vp.size() == 2); GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), - Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1))); + Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1))); } } -- cgit v1.2.3 From 2d1fb6b63f0ca0c7e027cc298fc16198a6283df1 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Fri, 11 Feb 2022 15:00:27 +0100 Subject: Do some code clean up/renaming --- .../benchmark/cech_complex_benchmark.cpp | 13 ++- .../example/cech_complex_example_from_points.cpp | 2 +- .../example/cech_complex_step_by_step.cpp | 6 +- src/Cech_complex/include/gudhi/Cech_complex.h | 40 +++----- .../include/gudhi/Cech_complex/Cech_kernel.h | 107 --------------------- .../include/gudhi/Cech_complex_blocker.h | 22 ++--- .../include/gudhi/sphere_circumradius.h | 62 ++++++++++++ src/Cech_complex/test/test_cech_complex.cpp | 10 +- src/Cech_complex/utilities/cech_persistence.cpp | 5 +- 9 files changed, 106 insertions(+), 161 deletions(-) delete mode 100644 src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h create mode 100644 src/Cech_complex/include/gudhi/sphere_circumradius.h (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 94c5fa4f..b283e1a8 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -34,9 +34,8 @@ using Proximity_graph = Gudhi::Proximity_graph; using Rips_complex = Gudhi::rips_complex::Rips_complex; using Kernel = CGAL::Epick_d>; using Point_cgal = typename Kernel::Point_d; -using Point_cloud_cgal = std::vector; using Points_off_reader_cgal = Gudhi::Points_off_reader; -using Cech_complex = Gudhi::cech_complex::Cech_complex; +using Cech_complex = Gudhi::cech_complex::Cech_complex; class Minimal_enclosing_ball_radius { public: @@ -83,11 +82,11 @@ int main(int argc, char* argv[]) { off_reader.get_point_cloud(), threshold, Minimal_enclosing_ball_radius()); std::clog << miniball_clock << std::endl; - Gudhi::Clock cgal_miniball_clock("Gudhi::Minimal_enclosing_ball_radius_cgal()"); + Gudhi::Clock cgal_circumsphere_clock("Gudhi::cech_complex::Sphere_circumradius_cgal()"); // Compute the proximity graph of the points - Proximity_graph cgal_miniball_prox_graph = Gudhi::compute_proximity_graph( - off_reader_cgal.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius()); - std::clog << cgal_miniball_clock << std::endl; + Proximity_graph cgal_circumsphere_prox_graph = Gudhi::compute_proximity_graph( + off_reader_cgal.get_point_cloud(), threshold, Gudhi::cech_complex::Sphere_circumradius()); + std::clog << cgal_circumsphere_clock << std::endl; boost::filesystem::path full_path(boost::filesystem::current_path()); std::clog << "Current path is : " << full_path << std::endl; @@ -109,7 +108,7 @@ int main(int argc, char* argv[]) { std::clog << radius << ";"; Gudhi::Clock rips_clock("Rips computation"); Rips_complex rips_complex_from_points(off_reader_cgal.get_point_cloud(), radius, - Gudhi::Minimal_enclosing_ball_radius()); + Gudhi::cech_complex::Sphere_circumradius()); Simplex_tree rips_stree; rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); // ------------------------------------------ diff --git a/src/Cech_complex/example/cech_complex_example_from_points.cpp b/src/Cech_complex/example/cech_complex_example_from_points.cpp index 78861951..38021e4a 100644 --- a/src/Cech_complex/example/cech_complex_example_from_points.cpp +++ b/src/Cech_complex/example/cech_complex_example_from_points.cpp @@ -16,7 +16,7 @@ int main() { using FT = typename Kernel::FT; using Point = typename Kernel::Point_d; using Point_cloud = std::vector; - using Cech_complex = Gudhi::cech_complex::Cech_complex; + using Cech_complex = Gudhi::cech_complex::Cech_complex; Point_cloud points; diff --git a/src/Cech_complex/example/cech_complex_step_by_step.cpp b/src/Cech_complex/example/cech_complex_step_by_step.cpp index c8dd1585..4401f6af 100644 --- a/src/Cech_complex/example/cech_complex_step_by_step.cpp +++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp @@ -9,7 +9,7 @@ */ #include -#include +#include #include #include @@ -52,7 +52,7 @@ class Cech_blocker { std::clog << "#(" << vertex << ")#"; #endif // DEBUG_TRACES } - Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); + Filtration_value radius = Gudhi::cech_complex::Sphere_circumradius()(points); #ifdef DEBUG_TRACES std::clog << "radius = " << radius << " - " << (radius > max_radius_) << std::endl; #endif // DEBUG_TRACES @@ -83,7 +83,7 @@ int main(int argc, char* argv[]) { // Compute the proximity graph of the points Proximity_graph prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), max_radius, - Gudhi::Minimal_enclosing_ball_radius()); + Gudhi::cech_complex::Sphere_circumradius()); // Construct the Cech complex in a Simplex Tree Simplex_tree st; diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 0031d861..375be1d2 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -11,7 +11,7 @@ #ifndef CECH_COMPLEX_H_ #define CECH_COMPLEX_H_ -#include // for Gudhi::Minimal_enclosing_ball_radius +#include // for Gudhi::cech_complex::Sphere_circumradius #include // for Gudhi::Proximity_graph #include // for GUDHI_CHECK #include // for Gudhi::cech_complex::Cech_blocker @@ -31,23 +31,21 @@ namespace cech_complex { * * \details * The data structure is a proximity graph, containing edges when the edge length is less or equal - * to a given max_radius. Edge length is computed from `Gudhi::Minimal_enclosing_ball_radius` distance function. + * to a given max_radius. Edge length is computed from `Gudhi::cech_complex::Sphere_circumradius` distance function. * - * \tparam SimplicialComplexForProximityGraph furnishes `Vertex_handle` and `Filtration_value` type definition required - * by `Gudhi::Proximity_graph`. + * \tparam Kernel CGAL kernel. + * + * \tparam SimplicialComplexForCechComplex furnishes `Vertex_handle` and `Filtration_value` type definition required + * by `Gudhi::Proximity_graph` and Cech blocker. * - * \tparam ForwardPointRange must be a range for which `std::begin()` and `std::end()` methods return input - * iterators on a point. `std::begin()` and `std::end()` methods are also required for a point. */ -template +template class Cech_complex { private: // Required by compute_proximity_graph - using Vertex_handle = typename SimplicialComplexForProximityGraph::Vertex_handle; - using Filtration_value = typename SimplicialComplexForProximityGraph::Filtration_value; - using Proximity_graph = Gudhi::Proximity_graph; - - public: + using Vertex_handle = typename SimplicialComplexForCechComplex::Vertex_handle; + using Filtration_value = typename SimplicialComplexForCechComplex::Filtration_value; + using Proximity_graph = Gudhi::Proximity_graph; using cech_blocker = Cech_blocker; @@ -57,27 +55,21 @@ class Cech_complex { // Numeric type of coordinates in the kernel using FT = typename cech_blocker::FT; // Sphere is a pair of point and squared radius. - using Sphere = typename std::pair; + using Sphere = typename cech_blocker::Sphere; - public: + public: /** \brief Cech_complex constructor from a list of points. * - * @param[in] points Range of points. + * @param[in] points Vector of points where each point is defined as `kernel::Point_d`. * @param[in] max_radius Maximal radius value. * - * \tparam ForwardPointRange must be a range of Point. Point must be a range of copyable Cartesian coordinates. - * */ - Cech_complex(const ForwardPointRange& points, Filtration_value max_radius) : max_radius_(max_radius) { - // Point cloud deep copy - -// point_cloud_.reserve(boost::size(points)); -// for (auto&& point : points) point_cloud_.emplace_back(point.cartesian_begin(), point.cartesian_end()); + Cech_complex(const Point_cloud & points, Filtration_value max_radius) : max_radius_(max_radius) { point_cloud_.assign(points.begin(), points.end()); - cech_skeleton_graph_ = Gudhi::compute_proximity_graph( - point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); + cech_skeleton_graph_ = Gudhi::compute_proximity_graph( + point_cloud_, max_radius_, Sphere_circumradius()); } /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal diff --git a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h deleted file mode 100644 index 89012206..00000000 --- a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h +++ /dev/null @@ -1,107 +0,0 @@ -/* 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): Hind Montassif - * - * Copyright (C) 2021 Inria - * - * Modification(s): - * - YYYY/MM Author: Description of the modification - */ - -#ifndef CECH_KERNEL_H_ -#define CECH_KERNEL_H_ - -#include // for #include - -#include // for std::sqrt -#include - -namespace Gudhi { - -// namespace cech_complex { - -/** @brief Compute the radius of the minimal enclosing ball between Points given by a range of coordinates. - * The points are assumed to have the same dimension. */ -template -class Minimal_enclosing_ball_radius { - private: - Kernel kernel_; - public: - using Point = typename Kernel::Point_d; - using Point_cloud = typename std::vector; - - /** \brief Enclosing ball radius from two points using CGAL. - * - * @param[in] point_1 - * @param[in] point_2 - * @return Enclosing ball radius for the two points. - * \tparam Point must be a Kernel::Point_d from CGAL. - * - */ - double operator()(const Point& point_1, const Point& point_2) const { - return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; - } - - - /** \brief Enclosing ball radius from a point cloud using CGAL. - * - * @param[in] point_cloud The points. - * @return Enclosing ball radius for the points. - * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. - * - */ - double operator()(const Point_cloud& point_cloud) const { - return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); - } - -}; - -/** - * \class Cech_kernel - * \brief Cech complex kernel container. - * - * \details - * The Cech complex kernel container stores CGAL Kernel and dispatch basic computations. - */ - -// template < typename Kernel > -// class Cech_kernel { -// private: -// // Kernel for functions access. -// Kernel kernel_; -// public: -// using Point_d = typename Kernel::Point_d; -// // Numeric type of coordinates in the kernel -// using FT = typename Kernel::FT; -// // Sphere is a pair of point and squared radius. -// using Sphere = typename std::pair; -// -// int get_dimension(const Point_d& p0) const { -// return kernel_.point_dimension_d_object()(p0); -// } -// -// template -// Sphere get_sphere(PointIterator begin, PointIterator end) const { -// Point_d c = kernel_.construct_circumcenter_d_object()(begin, end); -// FT r = kernel_.squared_distance_d_object()(c, *begin); -// return std::make_pair(std::move(c), std::move(r)); -// } -// -// template -// FT get_squared_radius(PointIterator begin, PointIterator end) const { -// return kernel_.compute_squared_radius_d_object()(begin, end); -// } -// -// FT get_squared_radius(const Sphere& sph) const { -// return sph.second; -// } -// }; - - -//} // namespace cech_complex - -// namespace cechcomplex = cech_complex; - -} // namespace Gudhi - -#endif // CECH_KERNEL_H_ diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index f7f86534..1a696422 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -31,17 +31,15 @@ namespace cech_complex { * \details * Čech blocker is an oracle constructed from a Cech_complex and a simplicial complex. * - * \tparam SimplicialComplexForProximityGraph furnishes `Simplex_handle` and `Filtration_value` type definition, + * \tparam SimplicialComplexForCech furnishes `Simplex_handle` and `Filtration_value` type definition, * `simplex_vertex_range(Simplex_handle sh)`and `assign_filtration(Simplex_handle sh, Filtration_value filt)` methods. * - * \tparam Chech_complex is required by the blocker. + * \tparam Cech_complex is required by the blocker. + * + * \tparam Kernel CGAL kernel. */ template class Cech_blocker { - private: - - using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; - using Filtration_value = typename SimplicialComplexForCech::Filtration_value; public: @@ -51,10 +49,10 @@ class Cech_blocker { // Sphere is a pair of point and squared radius. using Sphere = typename std::pair; - template - FT get_squared_radius(PointIterator begin, PointIterator end) const { - return kernel_.compute_squared_radius_d_object()(begin, end); - } + private: + + using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; + using Filtration_value = typename SimplicialComplexForCech::Filtration_value; template Sphere get_sphere(PointIterator begin, PointIterator end) const { @@ -63,6 +61,7 @@ class Cech_blocker { return std::make_pair(std::move(c), std::move(r)); } + public: /** \internal \brief Čech complex blocker operator() - the oracle - assigns the filtration value from the simplex * radius and returns if the simplex expansion must be blocked. @@ -108,10 +107,11 @@ class Cech_blocker { if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) <= sph.second) { radius = std::sqrt(cast_to_double(sph.second)); #ifdef DEBUG_TRACES - std::clog << "circumcenter: " << sph.first << ", radius: " << radius << std::endl; + std::clog << "center: " << sph.first << ", radius: " << radius << std::endl; #endif // DEBUG_TRACES if (cast_to_double(sph.second) < cast_to_double(min_enclos_ball.second)) min_enclos_ball = sph; + break; } } // Get the minimal radius of all faces enclosing balls if exists diff --git a/src/Cech_complex/include/gudhi/sphere_circumradius.h b/src/Cech_complex/include/gudhi/sphere_circumradius.h new file mode 100644 index 00000000..a6dec3dc --- /dev/null +++ b/src/Cech_complex/include/gudhi/sphere_circumradius.h @@ -0,0 +1,62 @@ +/* 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): Hind Montassif + * + * Copyright (C) 2021 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SPHERE_CIRCUMRADIUS_H_ +#define SPHERE_CIRCUMRADIUS_H_ + +#include // for #include + +#include // for std::sqrt +#include + +namespace Gudhi { + +namespace cech_complex { + +/** @brief Compute the circumradius of the sphere passing through points given by a range of coordinates. + * The points are assumed to have the same dimension. */ +template +class Sphere_circumradius { + private: + Kernel kernel_; + public: + using Point = typename Kernel::Point_d; + using Point_cloud = typename std::vector; + + /** \brief Circumradius of sphere passing through two points using CGAL. + * + * @param[in] point_1 + * @param[in] point_2 + * @return Sphere circumradius passing through two points. + * \tparam Point must be a Kernel::Point_d from CGAL. + * + */ + double operator()(const Point& point_1, const Point& point_2) const { + return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; + } + + /** \brief Circumradius of sphere passing through point cloud using CGAL. + * + * @param[in] point_cloud The points. + * @return Sphere circumradius passing through the points. + * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. + * + */ + double operator()(const Point_cloud& point_cloud) const { + return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); + } + +}; + +} // namespace cech_complex + +} // namespace Gudhi + +#endif // SPHERE_CIRCUMRADIUS_H_ diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index ca7a9778..4cf8b68f 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -22,7 +22,7 @@ // to construct Cech_complex from a OFF file of points #include #include -#include +#include #include #include // For EXACT or SAFE version @@ -36,7 +36,7 @@ using Point = typename Kernel::Point_d; using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; -using Cech_complex = Gudhi::cech_complex::Cech_complex; +using Cech_complex = Gudhi::cech_complex::Cech_complex; BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { // ---------------------------------------------------------------------------- @@ -108,11 +108,11 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { std::clog << vertex << ","; vp.push_back(points.at(vertex)); } - std::clog << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1)) + std::clog << ") - distance =" << Gudhi::cech_complex::Sphere_circumradius()(vp.at(0), vp.at(1)) << " - filtration =" << st.filtration(f_simplex) << std::endl; BOOST_CHECK(vp.size() == 2); GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), - Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1))); + Gudhi::cech_complex::Sphere_circumradius()(vp.at(0), vp.at(1))); } } @@ -153,7 +153,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { Simplex_tree::Filtration_value f1410 = st2.filtration(st2.find({1, 4, 10})); std::clog << "f1410= " << f1410 << std::endl; - // In this case, the computed sphere using CGAL kernel does not match the minimal enclosing ball; the filtration value check is therefore done against a hardcoded value + // In this case, the computed circumsphere using CGAL kernel does not match the minimal enclosing ball; the filtration value check is therefore done against a hardcoded value GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, 1.); Point_cloud points469; diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp index ccf63e3e..82992f2d 100644 --- a/src/Cech_complex/utilities/cech_persistence.cpp +++ b/src/Cech_complex/utilities/cech_persistence.cpp @@ -9,7 +9,7 @@ */ #include -#include +#include #include #include #include @@ -28,9 +28,8 @@ using Filtration_value = Simplex_tree::Filtration_value; using Kernel = CGAL::Epeck_d; using Point = typename Kernel::Point_d; -using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; -using Cech_complex = Gudhi::cech_complex::Cech_complex; +using Cech_complex = Gudhi::cech_complex::Cech_complex; using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; -- cgit v1.2.3 From 8730db2e8d1a8663358168ff6a20881c97773002 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Mon, 25 Apr 2022 14:57:26 +0200 Subject: Remove cech_complex_step_by_step example and Miniball --- src/Cech_complex/doc/Intro_cech_complex.h | 4 - src/Cech_complex/example/CMakeLists.txt | 10 - .../example/cech_complex_step_by_step.cpp | 150 ------ src/Cech_complex/include/gudhi/Miniball.COPYRIGHT | 4 - src/Cech_complex/include/gudhi/Miniball.README | 26 - src/Cech_complex/include/gudhi/Miniball.hpp | 523 --------------------- src/common/doc/examples.h | 1 - src/common/doc/main_page.md | 2 +- 8 files changed, 1 insertion(+), 719 deletions(-) delete mode 100644 src/Cech_complex/example/cech_complex_step_by_step.cpp delete mode 100644 src/Cech_complex/include/gudhi/Miniball.COPYRIGHT delete mode 100644 src/Cech_complex/include/gudhi/Miniball.README delete mode 100644 src/Cech_complex/include/gudhi/Miniball.hpp (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/doc/Intro_cech_complex.h b/src/Cech_complex/doc/Intro_cech_complex.h index 644fd6cc..595fb64b 100644 --- a/src/Cech_complex/doc/Intro_cech_complex.h +++ b/src/Cech_complex/doc/Intro_cech_complex.h @@ -62,10 +62,6 @@ namespace cech_complex { * This radius computation is the reason why the Cech_complex is taking much more time to be computed than the * \ref rips_complex but it offers more topological guarantees. * - * If the Cech_complex interfaces are not detailed enough for your need, please refer to - * - * cech_complex_step_by_step.cpp example, where the graph construction over the Simplex_tree is more detailed. - * * \subsection cechpointscloudexample Example from a point cloud * * This example builds the proximity graph from the given points, and maximal radius values. diff --git a/src/Cech_complex/example/CMakeLists.txt b/src/Cech_complex/example/CMakeLists.txt index 4d11ace2..7d52ed5e 100644 --- a/src/Cech_complex/example/CMakeLists.txt +++ b/src/Cech_complex/example/CMakeLists.txt @@ -1,16 +1,6 @@ project(Cech_complex_examples) if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.0.1) - if (TARGET Boost::program_options) - add_executable ( Cech_complex_example_step_by_step cech_complex_step_by_step.cpp ) - target_link_libraries(Cech_complex_example_step_by_step Boost::program_options) - if (TBB_FOUND) - target_link_libraries(Cech_complex_example_step_by_step ${TBB_LIBRARIES}) - endif() - add_test(NAME Cech_complex_utility_from_rips_on_tore_3D COMMAND $ - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-d" "3") - endif() - add_executable ( Cech_complex_example_from_points cech_complex_example_from_points.cpp) if (TBB_FOUND) target_link_libraries(Cech_complex_example_from_points ${TBB_LIBRARIES}) diff --git a/src/Cech_complex/example/cech_complex_step_by_step.cpp b/src/Cech_complex/example/cech_complex_step_by_step.cpp deleted file mode 100644 index 4401f6af..00000000 --- a/src/Cech_complex/example/cech_complex_step_by_step.cpp +++ /dev/null @@ -1,150 +0,0 @@ -/* 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): Vincent Rouvreau - * - * Copyright (C) 2018 Inria - * - * Modification(s): - * - YYYY/MM Author: Description of the modification - */ - -#include -#include -#include -#include - -#include - -#include - -#include -#include -#include // infinity -#include // for pair -#include - -// ---------------------------------------------------------------------------- -// cech_complex_step_by_step is an example of each step that is required to -// build a Cech over a Simplex_tree. Please refer to cech_complex_example_from_points to see -// how to do the same thing with the Cech complex wrapper for less detailed -// steps. -// ---------------------------------------------------------------------------- - -// Types definition -using Simplex_tree = Gudhi::Simplex_tree<>; -using Simplex_handle = Simplex_tree::Simplex_handle; -using Filtration_value = Simplex_tree::Filtration_value; -using Kernel = CGAL::Epeck_d; -using Point = typename Kernel::Point_d; -using Points_off_reader = Gudhi::Points_off_reader; -using Proximity_graph = Gudhi::Proximity_graph; - -class Cech_blocker { - private: - using Point_cloud = std::vector; - - public: - bool operator()(Simplex_handle sh) { - std::vector points; - for (auto vertex : simplex_tree_.simplex_vertex_range(sh)) { - points.push_back(point_cloud_[vertex]); -#ifdef DEBUG_TRACES - std::clog << "#(" << vertex << ")#"; -#endif // DEBUG_TRACES - } - Filtration_value radius = Gudhi::cech_complex::Sphere_circumradius()(points); -#ifdef DEBUG_TRACES - std::clog << "radius = " << radius << " - " << (radius > max_radius_) << std::endl; -#endif // DEBUG_TRACES - simplex_tree_.assign_filtration(sh, radius); - return (radius > max_radius_); - } - Cech_blocker(Simplex_tree& simplex_tree, Filtration_value max_radius, const std::vector& point_cloud) - : simplex_tree_(simplex_tree), max_radius_(max_radius), point_cloud_(point_cloud) { - } - - private: - Simplex_tree simplex_tree_; - Filtration_value max_radius_; - std::vector point_cloud_; -}; - -void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& max_radius, int& dim_max); - -int main(int argc, char* argv[]) { - std::string off_file_points; - Filtration_value max_radius; - int dim_max; - - program_options(argc, argv, off_file_points, max_radius, dim_max); - - // Extract the points from the file filepoints - Points_off_reader off_reader(off_file_points); - - // Compute the proximity graph of the points - Proximity_graph prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), max_radius, - Gudhi::cech_complex::Sphere_circumradius()); - - // Construct the Cech complex in a Simplex Tree - Simplex_tree st; - // insert the proximity graph in the simplex tree - st.insert_graph(prox_graph); - // expand the graph until dimension dim_max - st.expansion_with_blockers(dim_max, Cech_blocker(st, max_radius, off_reader.get_point_cloud())); - - std::clog << "The complex contains " << st.num_simplices() << " simplices \n"; - std::clog << " and has dimension " << st.dimension() << " \n"; - - // Sort the simplices in the order of the filtration - st.initialize_filtration(); - -#if DEBUG_TRACES - std::clog << "********************************************************************\n"; - std::clog << "* The complex contains " << st.num_simplices() << " simplices - dimension=" << st.dimension() << "\n"; - std::clog << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; - for (auto f_simplex : st.filtration_simplex_range()) { - std::clog << " " - << "[" << st.filtration(f_simplex) << "] "; - for (auto vertex : st.simplex_vertex_range(f_simplex)) { - std::clog << static_cast(vertex) << " "; - } - std::clog << std::endl; - } -#endif // DEBUG_TRACES - - return 0; -} - -void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& max_radius, int& dim_max) { - namespace po = boost::program_options; - po::options_description hidden("Hidden options"); - hidden.add_options()("input-file", po::value(&off_file_points), - "Name of an OFF file containing a point set.\n"); - - po::options_description visible("Allowed options", 100); - visible.add_options()("help,h", "produce help message")( - "max-radius,r", - po::value(&max_radius)->default_value(std::numeric_limits::infinity()), - "Maximal length of an edge for the Cech complex construction.")( - "cpx-dimension,d", po::value(&dim_max)->default_value(1), - "Maximal dimension of the Cech complex we want to compute."); - - po::positional_options_description pos; - pos.add("input-file", 1); - - po::options_description all; - all.add(visible).add(hidden); - - po::variables_map vm; - po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm); - po::notify(vm); - - if (vm.count("help") || !vm.count("input-file")) { - std::clog << std::endl; - std::clog << "Construct a Cech complex defined on a set of input points.\n \n"; - - std::clog << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; - std::clog << visible << std::endl; - exit(-1); - } -} diff --git a/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT b/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT deleted file mode 100644 index dbe4c553..00000000 --- a/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT +++ /dev/null @@ -1,4 +0,0 @@ -The miniball software is available under the GNU General Public License (GPLv3 - https://www.gnu.org/copyleft/gpl.html). -If your intended use is not compliant with this license, please buy a commercial license (EUR 500 - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/license.html). -You need a license if the software that you develop using Miniball V3.0 is not open source. - diff --git a/src/Cech_complex/include/gudhi/Miniball.README b/src/Cech_complex/include/gudhi/Miniball.README deleted file mode 100644 index 033d8953..00000000 --- a/src/Cech_complex/include/gudhi/Miniball.README +++ /dev/null @@ -1,26 +0,0 @@ -https://people.inf.ethz.ch/gaertner/subdir/software/miniball.html - -Smallest Enclosing Balls of Points - Fast and Robust in C++. -(high-quality software for smallest enclosing balls of balls is available in the computational geometry algorithms library CGAL) - - -This is the miniball software (V3.0) for computing smallest enclosing balls of points in arbitrary dimensions. It consists of a C++ header file Miniball.hpp (around 500 lines of code) and two example programs miniball_example.cpp and miniball_example_containers.cpp that demonstrate the usage. The first example stores the coordinates of the input points in a two-dimensional array, the second example uses a list of vectors to show how generic containers can be used. - -Credits: Aditya Gupta and Alexandros Konstantinakis-Karmis have significantly contributed to this version of the software. - -Changes - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/changes.txt - from previous versions. - -The theory - https://people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf - behind the miniball software (Proc. 7th Annual European Symposium on Algorithms (ESA), Lecture Notes in Computer Science 1643, Springer-Verlag, pp.325-338, 1999). - -Main Features: - - Very fast in low dimensions. 1 million points in 5-space are processed within 0.05 seconds on any recent machine. - - High numerical stability. Almost all input degeneracies (cospherical points, multiple points, points very close together) are routinely handled. - - Easily integrates into your code. You can freely choose the coordinate type of your points and the container to store the points. If you still need to adapt the code, the header is small and readable and contains documentation for all major methods. - - -Changes done for the GUDHI version of MiniBall: - - Add include guard - - Move Miniball namespace inside a new Gudhi namespace diff --git a/src/Cech_complex/include/gudhi/Miniball.hpp b/src/Cech_complex/include/gudhi/Miniball.hpp deleted file mode 100644 index ce6cbb5b..00000000 --- a/src/Cech_complex/include/gudhi/Miniball.hpp +++ /dev/null @@ -1,523 +0,0 @@ -// Copright (C) 1999-2013, Bernd Gaertner -// $Rev: 3581 $ -// -// This program is free software: 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. - -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. - -// You should have received a copy of the GNU General Public License -// along with this program. If not, see . -// -// Contact: -// -------- -// Bernd Gaertner -// Institute of Theoretical Computer Science -// ETH Zuerich -// CAB G31.1 -// CH-8092 Zuerich, Switzerland -// http://www.inf.ethz.ch/personal/gaertner - -#ifndef MINIBALL_HPP_ -#define MINIBALL_HPP_ - -#include -#include -#include -#include -#include - -namespace Gudhi { - -namespace Miniball { - - // Global Functions - // ================ - template - inline NT mb_sqr (NT r) {return r*r;} - - // Functors - // ======== - - // functor to map a point iterator to the corresponding coordinate iterator; - // generic version for points whose coordinate containers have begin() - template < typename Pit_, typename Cit_ > - struct CoordAccessor { - typedef Pit_ Pit; - typedef Cit_ Cit; - inline Cit operator() (Pit it) const { return (*it).begin(); } - }; - - // partial specialization for points whose coordinate containers are arrays - template < typename Pit_, typename Cit_ > - struct CoordAccessor { - typedef Pit_ Pit; - typedef Cit_* Cit; - inline Cit operator() (Pit it) const { return *it; } - }; - - // Class Declaration - // ================= - - template - class Miniball { - private: - // types - // The iterator type to go through the input points - typedef typename CoordAccessor::Pit Pit; - // The iterator type to go through the coordinates of a single point. - typedef typename CoordAccessor::Cit Cit; - // The coordinate type - typedef typename std::iterator_traits::value_type NT; - // The iterator to go through the support points - typedef typename std::list::iterator Sit; - - // data members... - const int d; // dimension - Pit points_begin; - Pit points_end; - CoordAccessor coord_accessor; - double time; - const NT nt0; // NT(0) - - //...for the algorithms - std::list L; - Sit support_end; - int fsize; // number of forced points - int ssize; // number of support points - - // ...for the ball updates - NT* current_c; - NT current_sqr_r; - NT** c; - NT* sqr_r; - - // helper arrays - NT* q0; - NT* z; - NT* f; - NT** v; - NT** a; - - public: - // The iterator type to go through the support points - typedef typename std::list::const_iterator SupportPointIterator; - - // PRE: [begin, end) is a nonempty range - // POST: computes the smallest enclosing ball of the points in the range - // [begin, end); the functor a maps a point iterator to an iterator - // through the d coordinates of the point - Miniball (int d_, Pit begin, Pit end, CoordAccessor ca = CoordAccessor()); - - // POST: returns a pointer to the first element of an array that holds - // the d coordinates of the center of the computed ball - const NT* center () const; - - // POST: returns the squared radius of the computed ball - NT squared_radius () const; - - // POST: returns the number of support points of the computed ball; - // the support points form a minimal set with the same smallest - // enclosing ball as the input set; in particular, the support - // points are on the boundary of the computed ball, and their - // number is at most d+1 - int nr_support_points () const; - - // POST: returns an iterator to the first support point - SupportPointIterator support_points_begin () const; - - // POST: returns a past-the-end iterator for the range of support points - SupportPointIterator support_points_end () const; - - // POST: returns the maximum excess of any input point w.r.t. the computed - // ball, divided by the squared radius of the computed ball. The - // excess of a point is the difference between its squared distance - // from the center and the squared radius; Ideally, the return value - // is 0. subopt is set to the absolute value of the most negative - // coefficient in the affine combination of the support points that - // yields the center. Ideally, this is a convex combination, and there - // is no negative coefficient in which case subopt is set to 0. - NT relative_error (NT& subopt) const; - - // POST: return true if the relative error is at most tol, and the - // suboptimality is 0; the default tolerance is 10 times the - // coordinate type's machine epsilon - bool is_valid (NT tol = NT(10) * std::numeric_limits::epsilon()) const; - - // POST: returns the time in seconds taken by the constructor call for - // computing the smallest enclosing ball - double get_time() const; - - // POST: deletes dynamically allocated arrays - ~Miniball(); - - private: - void mtf_mb (Sit n); - void mtf_move_to_front (Sit j); - void pivot_mb (Pit n); - void pivot_move_to_front (Pit j); - NT excess (Pit pit) const; - void pop (); - bool push (Pit pit); - NT suboptimality () const; - void create_arrays(); - void delete_arrays(); - }; - - // Class Definition - // ================ - template - Miniball::Miniball (int d_, Pit begin, Pit end, - CoordAccessor ca) - : d (d_), - points_begin (begin), - points_end (end), - coord_accessor (ca), - time (clock()), - nt0 (NT(0)), - L(), - support_end (L.begin()), - fsize(0), - ssize(0), - current_c (NULL), - current_sqr_r (NT(-1)), - c (NULL), - sqr_r (NULL), - q0 (NULL), - z (NULL), - f (NULL), - v (NULL), - a (NULL) - { - assert (points_begin != points_end); - create_arrays(); - - // set initial center - for (int j=0; j - Miniball::~Miniball() - { - delete_arrays(); - } - - template - void Miniball::create_arrays() - { - c = new NT*[d+1]; - v = new NT*[d+1]; - a = new NT*[d+1]; - for (int i=0; i - void Miniball::delete_arrays() - { - delete[] f; - delete[] z; - delete[] q0; - delete[] sqr_r; - for (int i=0; i - const typename Miniball::NT* - Miniball::center () const - { - return current_c; - } - - template - typename Miniball::NT - Miniball::squared_radius () const - { - return current_sqr_r; - } - - template - int Miniball::nr_support_points () const - { - assert (ssize < d+2); - return ssize; - } - - template - typename Miniball::SupportPointIterator - Miniball::support_points_begin () const - { - return L.begin(); - } - - template - typename Miniball::SupportPointIterator - Miniball::support_points_end () const - { - return support_end; - } - - template - typename Miniball::NT - Miniball::relative_error (NT& subopt) const - { - NT e, max_e = nt0; - // compute maximum absolute excess of support points - for (SupportPointIterator it = support_points_begin(); - it != support_points_end(); ++it) { - e = excess (*it); - if (e < nt0) e = -e; - if (e > max_e) { - max_e = e; - } - } - // compute maximum excess of any point - for (Pit i = points_begin; i != points_end; ++i) - if ((e = excess (i)) > max_e) - max_e = e; - - subopt = suboptimality(); - assert (current_sqr_r > nt0 || max_e == nt0); - return (current_sqr_r == nt0 ? nt0 : max_e / current_sqr_r); - } - - template - bool Miniball::is_valid (NT tol) const - { - NT suboptimality; - return ( (relative_error (suboptimality) <= tol) && (suboptimality == 0) ); - } - - template - double Miniball::get_time() const - { - return time; - } - - template - void Miniball::mtf_mb (Sit n) - { - // Algorithm 1: mtf_mb (L_{n-1}, B), where L_{n-1} = [L.begin, n) - // B: the set of forced points, defining the current ball - // S: the superset of support points computed by the algorithm - // -------------------------------------------------------------- - // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, - // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf - - // PRE: B = S - assert (fsize == ssize); - - support_end = L.begin(); - if ((fsize) == d+1) return; - - // incremental construction - for (Sit i = L.begin(); i != n;) - { - // INV: (support_end - L.begin() == |S|-|B|) - assert (std::distance (L.begin(), support_end) == ssize - fsize); - - Sit j = i++; - if (excess(*j) > nt0) - if (push(*j)) { // B := B + p_i - mtf_mb (j); // mtf_mb (L_{i-1}, B + p_i) - pop(); // B := B - p_i - mtf_move_to_front(j); - } - } - // POST: the range [L.begin(), support_end) stores the set S\B - } - - template - void Miniball::mtf_move_to_front (Sit j) - { - if (support_end == j) - support_end++; - L.splice (L.begin(), L, j); - } - - template - void Miniball::pivot_mb (Pit n) - { - // Algorithm 2: pivot_mb (L_{n-1}), where L_{n-1} = [L.begin, n) - // -------------------------------------------------------------- - // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, - // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf - NT old_sqr_r; - const NT* c; - Pit pivot, k; - NT e, max_e, sqr_r; - Cit p; - do { - old_sqr_r = current_sqr_r; - sqr_r = current_sqr_r; - - pivot = points_begin; - max_e = nt0; - for (k = points_begin; k != n; ++k) { - p = coord_accessor(k); - e = -sqr_r; - c = current_c; - for (int j=0; j(*p++-*c++); - if (e > max_e) { - max_e = e; - pivot = k; - } - } - - if (max_e > nt0) { - // check if the pivot is already contained in the support set - if (std::find(L.begin(), support_end, pivot) == support_end) { - assert (fsize == 0); - if (push (pivot)) { - mtf_mb(support_end); - pop(); - pivot_move_to_front(pivot); - } - } - } - } while (old_sqr_r < current_sqr_r); - } - - template - void Miniball::pivot_move_to_front (Pit j) - { - L.push_front(j); - if (std::distance(L.begin(), support_end) == d+2) - support_end--; - } - - template - inline typename Miniball::NT - Miniball::excess (Pit pit) const - { - Cit p = coord_accessor(pit); - NT e = -current_sqr_r; - NT* c = current_c; - for (int k=0; k(*p++-*c++); - } - return e; - } - - template - void Miniball::pop () - { - --fsize; - } - - template - bool Miniball::push (Pit pit) - { - int i, j; - NT eps = mb_sqr(std::numeric_limits::epsilon()); - - Cit cit = coord_accessor(pit); - Cit p = cit; - - if (fsize==0) { - for (i=0; i(v[fsize][j]); - z[fsize]*=2; - - // reject push if z_fsize too small - if (z[fsize](*p++-c[fsize-1][i]); - f[fsize]=e/z[fsize]; - - for (i=0; i - typename Miniball::NT - Miniball::suboptimality () const - { - NT* l = new NT[d+1]; - NT min_l = nt0; - l[0] = NT(1); - for (int i=ssize-1; i>0; --i) { - l[i] = f[i]; - for (int k=ssize-1; k>i; --k) - l[i]-=a[k][i]*l[k]; - if (l[i] < min_l) min_l = l[i]; - l[0] -= l[i]; - } - if (l[0] < min_l) min_l = l[0]; - delete[] l; - if (min_l < nt0) - return -min_l; - return nt0; - } -} // namespace Miniball - -} // namespace Gudhi - -#endif // MINIBALL_HPP_ diff --git a/src/common/doc/examples.h b/src/common/doc/examples.h index 879fb96a..0c4320f6 100644 --- a/src/common/doc/examples.h +++ b/src/common/doc/examples.h @@ -40,7 +40,6 @@ * @example edge_collapse_basic_example.cpp * \section Cech_complex_example_section Cech_complex * @example cech_persistence.cpp - * @example cech_complex_step_by_step.cpp * @example cech_complex_example_from_points.cpp * \section Bitmap_cubical_complex_example_section Bitmap_cubical_complex * @example periodic_cubical_complex_persistence.cpp diff --git a/src/common/doc/main_page.md b/src/common/doc/main_page.md index 17354179..6f995fee 100644 --- a/src/common/doc/main_page.md +++ b/src/common/doc/main_page.md @@ -181,7 +181,7 @@ Author: Vincent Rouvreau
Introduced in: GUDHI 2.2.0
Copyright: MIT [(GPL v3)](../../licensing/)
- Includes: [Miniball](https://people.inf.ethz.ch/gaertner/subdir/software/miniball.html)
+ Requires: \ref cgal -- cgit v1.2.3 From aec7ab1737a5284f4b7c2d1f7fe3eb7977df7537 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Tue, 26 Apr 2022 16:52:26 +0200 Subject: Modify cech doc Use Filtration_value instead of double for casting Use a templated range of points instead of vector in cech constructor Capitalize sphere_circumradius.h file name and make it private in doc --- src/Cech_complex/include/gudhi/Cech_complex.h | 19 +++---- .../include/gudhi/Cech_complex_blocker.h | 18 +++---- .../include/gudhi/Sphere_circumradius.h | 62 ++++++++++++++++++++++ .../include/gudhi/sphere_circumradius.h | 62 ---------------------- src/Cech_complex/test/test_cech_complex.cpp | 15 +++--- src/Cech_complex/utilities/cech_persistence.cpp | 1 - 6 files changed, 88 insertions(+), 89 deletions(-) create mode 100644 src/Cech_complex/include/gudhi/Sphere_circumradius.h delete mode 100644 src/Cech_complex/include/gudhi/sphere_circumradius.h (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 375be1d2..fc39f75b 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -11,7 +11,7 @@ #ifndef CECH_COMPLEX_H_ #define CECH_COMPLEX_H_ -#include // for Gudhi::cech_complex::Sphere_circumradius +#include // for Gudhi::cech_complex::Sphere_circumradius #include // for Gudhi::Proximity_graph #include // for GUDHI_CHECK #include // for Gudhi::cech_complex::Cech_blocker @@ -25,15 +25,15 @@ namespace cech_complex { /** * \class Cech_complex - * \brief Cech complex data structure. + * \brief Cech complex class. * * \ingroup cech_complex * * \details - * The data structure is a proximity graph, containing edges when the edge length is less or equal - * to a given max_radius. Edge length is computed from `Gudhi::cech_complex::Sphere_circumradius` distance function. + * Cech complex is a simplicial complex constructed from a proximity graph, where the set of all simplices is filtered + * by the radius of their minimal enclosing ball and bounded by the given max_radius. * - * \tparam Kernel CGAL kernel. + * \tparam Kernel CGAL kernel: either Epick_d or Epeck_d. * * \tparam SimplicialComplexForCechComplex furnishes `Vertex_handle` and `Filtration_value` type definition required * by `Gudhi::Proximity_graph` and Cech blocker. @@ -58,15 +58,16 @@ class Cech_complex { using Sphere = typename cech_blocker::Sphere; public: - /** \brief Cech_complex constructor from a list of points. + /** \brief Cech_complex constructor from a range of points. * - * @param[in] points Vector of points where each point is defined as `kernel::Point_d`. + * @param[in] points Range of points where each point is defined as `kernel::Point_d`. * @param[in] max_radius Maximal radius value. * */ - Cech_complex(const Point_cloud & points, Filtration_value max_radius) : max_radius_(max_radius) { + template + Cech_complex(const InputPointRange & points, Filtration_value max_radius) : max_radius_(max_radius) { - point_cloud_.assign(points.begin(), points.end()); + point_cloud_.assign(std::begin(points), std::end(points)); cech_skeleton_graph_ = Gudhi::compute_proximity_graph( point_cloud_, max_radius_, Sphere_circumradius()); diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 1a696422..1a09f7e1 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -11,7 +11,7 @@ #ifndef CECH_COMPLEX_BLOCKER_H_ #define CECH_COMPLEX_BLOCKER_H_ -#include // for casting from FT to double +#include // for casting from FT to Filtration_value and double to FT #include #include @@ -36,7 +36,7 @@ namespace cech_complex { * * \tparam Cech_complex is required by the blocker. * - * \tparam Kernel CGAL kernel. + * \tparam Kernel CGAL kernel: either Epick_d or Epeck_d. */ template class Cech_blocker { @@ -69,8 +69,8 @@ class Cech_blocker { * \return true if the simplex radius is greater than the Cech_complex max_radius*/ bool operator()(Simplex_handle sh) { using Point_cloud = std::vector; - CGAL::NT_converter cast_to_double; - Filtration_value radius = 0.; + CGAL::NT_converter cast_to_fv; + Filtration_value radius = 0; // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices Sphere min_enclos_ball; @@ -105,18 +105,18 @@ class Cech_blocker { face_points.clear(); if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) <= sph.second) { - radius = std::sqrt(cast_to_double(sph.second)); + radius = std::sqrt(cast_to_fv(sph.second)); #ifdef DEBUG_TRACES std::clog << "center: " << sph.first << ", radius: " << radius << std::endl; #endif // DEBUG_TRACES - if (cast_to_double(sph.second) < cast_to_double(min_enclos_ball.second)) + if (sph.second < min_enclos_ball.second) min_enclos_ball = sph; break; } } // Get the minimal radius of all faces enclosing balls if exists - if(cast_to_double(min_enclos_ball.second) != std::numeric_limits::max()) { - radius = std::sqrt(cast_to_double(min_enclos_ball.second)); + if(min_enclos_ball.second != std::numeric_limits::max()) { + radius = std::sqrt(cast_to_fv(min_enclos_ball.second)); sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); cc_ptr_->get_cache().push_back(min_enclos_ball); @@ -128,7 +128,7 @@ class Cech_blocker { points.push_back(cc_ptr_->get_point(vertex)); } Sphere sph = get_sphere(points.cbegin(), points.cend()); - radius = std::sqrt(cast_to_double(sph.second)); + radius = std::sqrt(cast_to_fv(sph.second)); sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); cc_ptr_->get_cache().push_back(sph); diff --git a/src/Cech_complex/include/gudhi/Sphere_circumradius.h b/src/Cech_complex/include/gudhi/Sphere_circumradius.h new file mode 100644 index 00000000..b0d9f7cc --- /dev/null +++ b/src/Cech_complex/include/gudhi/Sphere_circumradius.h @@ -0,0 +1,62 @@ +/* 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): Hind Montassif + * + * Copyright (C) 2021 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SPHERE_CIRCUMRADIUS_H_ +#define SPHERE_CIRCUMRADIUS_H_ + +#include // for #include + +#include // for std::sqrt +#include + +namespace Gudhi { + +namespace cech_complex { + +/** \private @brief Compute the circumradius of the sphere passing through points given by a range of coordinates. + * The points are assumed to have the same dimension. */ +template +class Sphere_circumradius { + private: + Kernel kernel_; + public: + using Point = typename Kernel::Point_d; + using Point_cloud = typename std::vector; + + /** \brief Circumradius of sphere passing through two points using CGAL. + * + * @param[in] point_1 + * @param[in] point_2 + * @return Sphere circumradius passing through two points. + * \tparam Point must be a Kernel::Point_d from CGAL. + * + */ + double operator()(const Point& point_1, const Point& point_2) const { + return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; + } + + /** \brief Circumradius of sphere passing through point cloud using CGAL. + * + * @param[in] point_cloud The points. + * @return Sphere circumradius passing through the points. + * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. + * + */ + double operator()(const Point_cloud& point_cloud) const { + return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); + } + +}; + +} // namespace cech_complex + +} // namespace Gudhi + +#endif // SPHERE_CIRCUMRADIUS_H_ diff --git a/src/Cech_complex/include/gudhi/sphere_circumradius.h b/src/Cech_complex/include/gudhi/sphere_circumradius.h deleted file mode 100644 index a6dec3dc..00000000 --- a/src/Cech_complex/include/gudhi/sphere_circumradius.h +++ /dev/null @@ -1,62 +0,0 @@ -/* 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): Hind Montassif - * - * Copyright (C) 2021 Inria - * - * Modification(s): - * - YYYY/MM Author: Description of the modification - */ - -#ifndef SPHERE_CIRCUMRADIUS_H_ -#define SPHERE_CIRCUMRADIUS_H_ - -#include // for #include - -#include // for std::sqrt -#include - -namespace Gudhi { - -namespace cech_complex { - -/** @brief Compute the circumradius of the sphere passing through points given by a range of coordinates. - * The points are assumed to have the same dimension. */ -template -class Sphere_circumradius { - private: - Kernel kernel_; - public: - using Point = typename Kernel::Point_d; - using Point_cloud = typename std::vector; - - /** \brief Circumradius of sphere passing through two points using CGAL. - * - * @param[in] point_1 - * @param[in] point_2 - * @return Sphere circumradius passing through two points. - * \tparam Point must be a Kernel::Point_d from CGAL. - * - */ - double operator()(const Point& point_1, const Point& point_2) const { - return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; - } - - /** \brief Circumradius of sphere passing through point cloud using CGAL. - * - * @param[in] point_cloud The points. - * @return Sphere circumradius passing through the points. - * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. - * - */ - double operator()(const Point_cloud& point_cloud) const { - return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); - } - -}; - -} // namespace cech_complex - -} // namespace Gudhi - -#endif // SPHERE_CIRCUMRADIUS_H_ diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 4cf8b68f..ea32f596 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -22,7 +22,6 @@ // to construct Cech_complex from a OFF file of points #include #include -#include #include #include // For EXACT or SAFE version @@ -139,18 +138,18 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { } Kernel kern; - Simplex_tree::Filtration_value f012 = st2.filtration(st2.find({0, 1, 2})); + Filtration_value f012 = st2.filtration(st2.find({0, 1, 2})); std::clog << "f012= " << f012 << std::endl; - CGAL::NT_converter cast_to_double; - GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points012.begin(), points012.end())))); + CGAL::NT_converter cast_to_fv; + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(cast_to_fv(kern.compute_squared_radius_d_object()(points012.begin(), points012.end())))); Point_cloud points1410; points1410.push_back(cech_complex_for_doc.get_point(1)); points1410.push_back(cech_complex_for_doc.get_point(4)); points1410.push_back(cech_complex_for_doc.get_point(10)); - Simplex_tree::Filtration_value f1410 = st2.filtration(st2.find({1, 4, 10})); + Filtration_value f1410 = st2.filtration(st2.find({1, 4, 10})); std::clog << "f1410= " << f1410 << std::endl; // In this case, the computed circumsphere using CGAL kernel does not match the minimal enclosing ball; the filtration value check is therefore done against a hardcoded value @@ -161,10 +160,10 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { points469.push_back(cech_complex_for_doc.get_point(6)); points469.push_back(cech_complex_for_doc.get_point(9)); - Simplex_tree::Filtration_value f469 = st2.filtration(st2.find({4, 6, 9})); + Filtration_value f469 = st2.filtration(st2.find({4, 6, 9})); std::clog << "f469= " << f469 << std::endl; - GUDHI_TEST_FLOAT_EQUALITY_CHECK(f469, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points469.begin(), points469.end())))); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f469, std::sqrt(cast_to_fv(kern.compute_squared_radius_d_object()(points469.begin(), points469.end())))); BOOST_CHECK((st2.find({6, 7, 8}) == st2.null_simplex())); BOOST_CHECK((st2.find({3, 5, 7}) == st2.null_simplex())); @@ -246,7 +245,7 @@ BOOST_AUTO_TEST_CASE(Cech_create_complex_throw) { // // ---------------------------------------------------------------------------- std::string off_file_name("alphacomplexdoc.off"); - double max_radius = 12.0; + Filtration_value max_radius = 12.0; std::clog << "========== OFF FILE NAME = " << off_file_name << " - Cech max_radius=" << max_radius << "==========" << std::endl; diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp index 82992f2d..75d10c0f 100644 --- a/src/Cech_complex/utilities/cech_persistence.cpp +++ b/src/Cech_complex/utilities/cech_persistence.cpp @@ -9,7 +9,6 @@ */ #include -#include #include #include #include -- cgit v1.2.3 From 70c20c20f89e2037544e7906c5743a30a7e3beb7 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Wed, 27 Apr 2022 16:41:22 +0200 Subject: Remove unnecessary code from cech blocker --- .../include/gudhi/Cech_complex_blocker.h | 41 +++++++++------------- 1 file changed, 17 insertions(+), 24 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 1a09f7e1..72876512 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -12,6 +12,7 @@ #define CECH_COMPLEX_BLOCKER_H_ #include // for casting from FT to Filtration_value and double to FT +#include #include #include @@ -73,10 +74,7 @@ class Cech_blocker { Filtration_value radius = 0; // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices - Sphere min_enclos_ball; - CGAL::NT_converter cast_to_FT; - min_enclos_ball.second = cast_to_FT(std::numeric_limits::max()); - Point_cloud face_points; + boost::optional min_enclos_ball; for (auto face : sc_ptr_->boundary_simplex_range(sh)) { // Find which vertex of sh is missing in face. We rely on the fact that simplex_vertex_range is sorted. auto longlist = sc_ptr_->simplex_vertex_range(sh); @@ -88,41 +86,36 @@ class Cech_blocker { while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; } auto extra = *longiter; // Vertex_handle - for (auto vertex : sc_ptr_->simplex_vertex_range(face)) { - face_points.push_back(cc_ptr_->get_point(vertex)); - #ifdef DEBUG_TRACES - std::clog << "#(" << vertex << ")#"; - #endif // DEBUG_TRACES - } Sphere sph; auto k = sc_ptr_->key(face); if(k != sc_ptr_->null_key()) { sph = cc_ptr_->get_cache().at(k); } else { + Point_cloud face_points; + for (auto vertex : sc_ptr_->simplex_vertex_range(face)) { + face_points.push_back(cc_ptr_->get_point(vertex)); + #ifdef DEBUG_TRACES + std::clog << "#(" << vertex << ")#"; + #endif // DEBUG_TRACES + } sph = get_sphere(face_points.cbegin(), face_points.cend()); + face_points.clear(); } - face_points.clear(); - + // Check if the minimal enclosing ball of current face contains the extra point if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) <= sph.second) { - radius = std::sqrt(cast_to_fv(sph.second)); #ifdef DEBUG_TRACES std::clog << "center: " << sph.first << ", radius: " << radius << std::endl; #endif // DEBUG_TRACES - if (sph.second < min_enclos_ball.second) - min_enclos_ball = sph; + radius = std::sqrt(cast_to_fv(sph.second)); + sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); + cc_ptr_->get_cache().push_back(sph); + min_enclos_ball.emplace(cc_ptr_->get_cache().back()); break; } } - // Get the minimal radius of all faces enclosing balls if exists - if(min_enclos_ball.second != std::numeric_limits::max()) { - radius = std::sqrt(cast_to_fv(min_enclos_ball.second)); - - sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); - cc_ptr_->get_cache().push_back(min_enclos_ball); - } - - if (radius == 0.) { // Spheres of each face don't contain the whole simplex + // Spheres of each face don't contain the whole simplex + if(!min_enclos_ball) { Point_cloud points; for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { points.push_back(cc_ptr_->get_point(vertex)); -- cgit v1.2.3 From a14d9becec7361a2559fa239c1ca8f2c1b5c5768 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Wed, 27 Apr 2022 17:08:14 +0200 Subject: Ultimately, we don't really need to store the min enclosing ball in case it contains the extra point, a bool is enough --- src/Cech_complex/include/gudhi/Cech_complex_blocker.h | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 72876512..fb12946a 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -11,8 +11,7 @@ #ifndef CECH_COMPLEX_BLOCKER_H_ #define CECH_COMPLEX_BLOCKER_H_ -#include // for casting from FT to Filtration_value and double to FT -#include +#include // for casting from FT to Filtration_value #include #include @@ -72,9 +71,9 @@ class Cech_blocker { using Point_cloud = std::vector; CGAL::NT_converter cast_to_fv; Filtration_value radius = 0; + bool is_min_enclos_ball = false; // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices - boost::optional min_enclos_ball; for (auto face : sc_ptr_->boundary_simplex_range(sh)) { // Find which vertex of sh is missing in face. We rely on the fact that simplex_vertex_range is sorted. auto longlist = sc_ptr_->simplex_vertex_range(sh); @@ -107,15 +106,15 @@ class Cech_blocker { #ifdef DEBUG_TRACES std::clog << "center: " << sph.first << ", radius: " << radius << std::endl; #endif // DEBUG_TRACES + is_min_enclos_ball = true; radius = std::sqrt(cast_to_fv(sph.second)); sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); cc_ptr_->get_cache().push_back(sph); - min_enclos_ball.emplace(cc_ptr_->get_cache().back()); break; } } // Spheres of each face don't contain the whole simplex - if(!min_enclos_ball) { + if(!is_min_enclos_ball) { Point_cloud points; for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { points.push_back(cc_ptr_->get_point(vertex)); -- cgit v1.2.3 From a6a68c11455a554619d8a5b5d2f92c1ddbf45e99 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Thu, 28 Apr 2022 16:38:34 +0200 Subject: Put edge sphere in cache --- src/Cech_complex/include/gudhi/Cech_complex_blocker.h | 4 ++++ 1 file changed, 4 insertions(+) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index fb12946a..3141d27a 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -99,6 +99,10 @@ class Cech_blocker { #endif // DEBUG_TRACES } sph = get_sphere(face_points.cbegin(), face_points.cend()); + // Put edge sphere in cache + sc_ptr_->assign_key(face, cc_ptr_->get_cache().size()); + cc_ptr_->get_cache().push_back(sph); + // Clear face_points face_points.clear(); } // Check if the minimal enclosing ball of current face contains the extra point -- cgit v1.2.3 From 868369dd61fb6ef475ffa3af724907927121b6bb Mon Sep 17 00:00:00 2001 From: Hind-M Date: Thu, 16 Jun 2022 15:54:21 +0200 Subject: Add exact option for exact cech variant --- .../benchmark/cech_complex_benchmark.cpp | 22 ++++++++++++++-------- src/Cech_complex/include/gudhi/Cech_complex.h | 6 ++++-- .../include/gudhi/Cech_complex_blocker.h | 21 ++++++++++++++------- 3 files changed, 32 insertions(+), 17 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index d2a71879..19142780 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -31,7 +31,7 @@ using Points_off_reader = Gudhi::Points_off_reader; using Rips_complex = Gudhi::rips_complex::Rips_complex; template -Simplex_tree benchmark_cech(const std::string& off_file_points, const Filtration_value& radius, const int& dim_max) { +Simplex_tree benchmark_cech(const std::string& off_file_points, const Filtration_value& radius, const int& dim_max, const bool exact) { using Point_cgal = typename Kernel::Point_d; using Points_off_reader_cgal = Gudhi::Points_off_reader; using Cech_complex = Gudhi::cech_complex::Cech_complex; @@ -42,7 +42,7 @@ Simplex_tree benchmark_cech(const std::string& off_file_points, const Filtration Gudhi::Clock cech_clock("Cech computation"); Cech_complex cech_complex_from_points(off_reader_cgal.get_point_cloud(), radius); Simplex_tree cech_stree; - cech_complex_from_points.create_complex(cech_stree, dim_max); + cech_complex_from_points.create_complex(cech_stree, dim_max, exact); // ------------------------------------------ // Display information about the Cech complex @@ -56,8 +56,9 @@ int main(int argc, char* argv[]) { boost::filesystem::path full_path(boost::filesystem::current_path()); std::clog << "Current path is : " << full_path << std::endl; - std::clog << "File name ; Radius ; Rips time ; Dim-3 Epick Cech time ; Dynamic_dim Epick Cech time ; " - "Dim-3 Epeck Cech time ; Dynamic_dim Epeck Cech time ; Cech nb simplices ; Rips nb simplices;" + std::clog << "File name ; Radius ; Rips time ; Dim-3 Fast Cech time ; Dynamic_dim Fast Cech time ; " + "Dim-3 Safe Cech time ; Dynamic_dim Safe Cech time ; Dim-3 Exact Cech time ; Dynamic_dim Exact Cech time ; " + "Cech nb simplices ; Rips nb simplices;" << std::endl; boost::filesystem::directory_iterator end_itr; // default construction yields past-the-end for (boost::filesystem::directory_iterator itr(boost::filesystem::current_path()); itr != end_itr; ++itr) { @@ -83,10 +84,15 @@ int main(int argc, char* argv[]) { // -------------- // Cech complex // -------------- - benchmark_cech>>(itr->path().string(), radius, p0.size() - 1); - benchmark_cech>(itr->path().string(), radius, p0.size() - 1); - benchmark_cech>>(itr->path().string(), radius, p0.size() - 1); - auto cech_stree = benchmark_cech>(itr->path().string(), radius, p0.size() - 1); + // Fast + benchmark_cech>>(itr->path().string(), radius, p0.size() - 1, false); + benchmark_cech>(itr->path().string(), radius, p0.size() - 1, false); + // Safe + benchmark_cech>>(itr->path().string(), radius, p0.size() - 1, false); + benchmark_cech>(itr->path().string(), radius, p0.size() - 1, false); + // Exact + benchmark_cech>>(itr->path().string(), radius, p0.size() - 1, true); + auto cech_stree = benchmark_cech>(itr->path().string(), radius, p0.size() - 1, true); std::clog << cech_stree.num_simplices() << " ; "; std::clog << rips_stree.num_simplices() << ";" << std::endl; diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index fc39f75b..2c6d3df5 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -78,17 +78,19 @@ class Cech_complex { * * @param[in] complex SimplicialComplexForCech to be created. * @param[in] dim_max graph expansion until this given maximal dimension. + * @param[in] exact Exact filtration values computation. Not exact if `Kernel` is not CGAL::Epeck_d. * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0. * */ - void create_complex(SimplicialComplexForCechComplex& complex, int dim_max) { + void create_complex(SimplicialComplexForCechComplex& complex, int dim_max, const bool exact = false) { GUDHI_CHECK(complex.num_vertices() == 0, std::invalid_argument("Cech_complex::create_complex - simplicial complex is not empty")); // insert the proximity graph in the simplicial complex complex.insert_graph(cech_skeleton_graph_); // expand the graph until dimension dim_max - complex.expansion_with_blockers(dim_max, cech_blocker(&complex, this)); + complex.expansion_with_blockers(dim_max, cech_blocker(&complex, this, exact)); } /** @return max_radius value given at construction. */ diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 3141d27a..087390b6 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -94,9 +94,9 @@ class Cech_blocker { Point_cloud face_points; for (auto vertex : sc_ptr_->simplex_vertex_range(face)) { face_points.push_back(cc_ptr_->get_point(vertex)); - #ifdef DEBUG_TRACES - std::clog << "#(" << vertex << ")#"; - #endif // DEBUG_TRACES +#ifdef DEBUG_TRACES + std::clog << "#(" << vertex << ")#"; +#endif // DEBUG_TRACES } sph = get_sphere(face_points.cbegin(), face_points.cend()); // Put edge sphere in cache @@ -107,10 +107,13 @@ class Cech_blocker { } // Check if the minimal enclosing ball of current face contains the extra point if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) <= sph.second) { - #ifdef DEBUG_TRACES - std::clog << "center: " << sph.first << ", radius: " << radius << std::endl; - #endif // DEBUG_TRACES +#ifdef DEBUG_TRACES + std::clog << "center: " << sph.first << ", radius: " << radius << std::endl; +#endif // DEBUG_TRACES is_min_enclos_ball = true; +#if CGAL_VERSION_NR >= 1050000000 + if(exact_) CGAL::exact(sph.second); +#endif radius = std::sqrt(cast_to_fv(sph.second)); sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); cc_ptr_->get_cache().push_back(sph); @@ -124,6 +127,9 @@ class Cech_blocker { points.push_back(cc_ptr_->get_point(vertex)); } Sphere sph = get_sphere(points.cbegin(), points.cend()); +#if CGAL_VERSION_NR >= 1050000000 + if(exact_) CGAL::exact(sph.second); +#endif radius = std::sqrt(cast_to_fv(sph.second)); sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); @@ -138,12 +144,13 @@ class Cech_blocker { } /** \internal \brief Čech complex blocker constructor. */ - Cech_blocker(SimplicialComplexForCech* sc_ptr, Cech_complex* cc_ptr) : sc_ptr_(sc_ptr), cc_ptr_(cc_ptr) {} + Cech_blocker(SimplicialComplexForCech* sc_ptr, Cech_complex* cc_ptr, const bool exact) : sc_ptr_(sc_ptr), cc_ptr_(cc_ptr), exact_(exact) {} private: SimplicialComplexForCech* sc_ptr_; Cech_complex* cc_ptr_; Kernel kernel_; + const bool exact_; }; } // namespace cech_complex -- cgit v1.2.3 From 3fa972970514333d4db22ec7628c5c1a4de3c6e8 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Tue, 21 Jun 2022 15:04:27 +0200 Subject: -Add/modify some comments -Some other minor changes -Change license to LGPL --- .../benchmark/cech_complex_benchmark.cpp | 20 +++++++++++--------- src/Cech_complex/include/gudhi/Cech_complex.h | 2 +- .../include/gudhi/Cech_complex_blocker.h | 2 +- src/common/doc/main_page.md | 2 +- 4 files changed, 14 insertions(+), 12 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 19142780..a9dc5d0d 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -61,20 +61,22 @@ int main(int argc, char* argv[]) { "Cech nb simplices ; Rips nb simplices;" << std::endl; boost::filesystem::directory_iterator end_itr; // default construction yields past-the-end + // For every ".off" file in the current directory, and for 3 predefined thresholds, compare Rips and various Cech constructions for (boost::filesystem::directory_iterator itr(boost::filesystem::current_path()); itr != end_itr; ++itr) { if (!boost::filesystem::is_directory(itr->status())) { if (itr->path().extension() == ".off") { Points_off_reader off_reader(itr->path().string()); Point p0 = off_reader.get_point_cloud()[0]; - - for (Filtration_value radius = 0.1; radius < 0.4; radius += 0.1) { + // Loop over the different thresholds + for (Filtration_value radius = 0.1; radius < 0.35; radius += 0.1) { std::clog << itr->path().stem() << " ; "; std::clog << radius << " ; "; Gudhi::Clock rips_clock("Rips computation"); Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), radius, Gudhi::Euclidean_distance()); Simplex_tree rips_stree; - rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); + int dim_max = p0.size() - 1; + rips_complex_from_points.create_complex(rips_stree, dim_max); // ------------------------------------------ // Display information about the Rips complex // ------------------------------------------ @@ -85,14 +87,14 @@ int main(int argc, char* argv[]) { // Cech complex // -------------- // Fast - benchmark_cech>>(itr->path().string(), radius, p0.size() - 1, false); - benchmark_cech>(itr->path().string(), radius, p0.size() - 1, false); + benchmark_cech>>(itr->path().string(), radius, dim_max, false); + benchmark_cech>(itr->path().string(), radius, dim_max, false); // Safe - benchmark_cech>>(itr->path().string(), radius, p0.size() - 1, false); - benchmark_cech>(itr->path().string(), radius, p0.size() - 1, false); + benchmark_cech>>(itr->path().string(), radius, dim_max, false); + benchmark_cech>(itr->path().string(), radius, dim_max, false); // Exact - benchmark_cech>>(itr->path().string(), radius, p0.size() - 1, true); - auto cech_stree = benchmark_cech>(itr->path().string(), radius, p0.size() - 1, true); + benchmark_cech>>(itr->path().string(), radius, dim_max, true); + auto cech_stree = benchmark_cech>(itr->path().string(), radius, dim_max, true); std::clog << cech_stree.num_simplices() << " ; "; std::clog << rips_stree.num_simplices() << ";" << std::endl; diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 2c6d3df5..bae21d28 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -30,7 +30,7 @@ namespace cech_complex { * \ingroup cech_complex * * \details - * Cech complex is a simplicial complex constructed from a proximity graph, where the set of all simplices is filtered + * Cech complex is a simplicial complex where the set of all simplices is filtered * by the radius of their minimal enclosing ball and bounded by the given max_radius. * * \tparam Kernel CGAL kernel: either Epick_d or Epeck_d. diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 087390b6..9cd49a52 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -133,7 +133,7 @@ class Cech_blocker { radius = std::sqrt(cast_to_fv(sph.second)); sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); - cc_ptr_->get_cache().push_back(sph); + cc_ptr_->get_cache().push_back(std::move(sph)); } #ifdef DEBUG_TRACES diff --git a/src/common/doc/main_page.md b/src/common/doc/main_page.md index 2cb02e3f..ce903405 100644 --- a/src/common/doc/main_page.md +++ b/src/common/doc/main_page.md @@ -180,7 +180,7 @@ Author: Vincent Rouvreau
Introduced in: GUDHI 2.2.0
- Copyright: MIT [(GPL v3)](../../licensing/)
+ Copyright: MIT [(LGPL v3)](../../licensing/)
Requires: \ref cgal -- cgit v1.2.3 From b829a198e16fbef4c0cb2698b2c723fa353aac55 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Fri, 24 Jun 2022 11:03:22 +0200 Subject: Use CGAL::NT_converter instead of CGAL::to_double in Sphere_circumradius --- src/Cech_complex/include/gudhi/Cech_complex.h | 2 +- src/Cech_complex/include/gudhi/Cech_complex_blocker.h | 1 + src/Cech_complex/include/gudhi/Sphere_circumradius.h | 15 +++++++++------ src/Cech_complex/test/test_cech_complex.cpp | 4 ++-- 4 files changed, 13 insertions(+), 9 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index bae21d28..08b7a72f 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -70,7 +70,7 @@ class Cech_complex { point_cloud_.assign(std::begin(points), std::end(points)); cech_skeleton_graph_ = Gudhi::compute_proximity_graph( - point_cloud_, max_radius_, Sphere_circumradius()); + point_cloud_, max_radius_, Sphere_circumradius()); } /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 9cd49a52..25d9a71f 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -12,6 +12,7 @@ #define CECH_COMPLEX_BLOCKER_H_ #include // for casting from FT to Filtration_value +#include // for CGAL::exact #include #include diff --git a/src/Cech_complex/include/gudhi/Sphere_circumradius.h b/src/Cech_complex/include/gudhi/Sphere_circumradius.h index b0d9f7cc..790f6950 100644 --- a/src/Cech_complex/include/gudhi/Sphere_circumradius.h +++ b/src/Cech_complex/include/gudhi/Sphere_circumradius.h @@ -11,7 +11,7 @@ #ifndef SPHERE_CIRCUMRADIUS_H_ #define SPHERE_CIRCUMRADIUS_H_ -#include // for #include +#include // for #include which is not working/compiling alone #include // for std::sqrt #include @@ -22,14 +22,17 @@ namespace cech_complex { /** \private @brief Compute the circumradius of the sphere passing through points given by a range of coordinates. * The points are assumed to have the same dimension. */ -template +template class Sphere_circumradius { private: Kernel kernel_; public: + using FT = typename Kernel::FT; using Point = typename Kernel::Point_d; using Point_cloud = typename std::vector; + CGAL::NT_converter cast_to_fv; + /** \brief Circumradius of sphere passing through two points using CGAL. * * @param[in] point_1 @@ -38,8 +41,8 @@ class Sphere_circumradius { * \tparam Point must be a Kernel::Point_d from CGAL. * */ - double operator()(const Point& point_1, const Point& point_2) const { - return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; + Filtration_value operator()(const Point& point_1, const Point& point_2) const { + return std::sqrt(cast_to_fv(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; } /** \brief Circumradius of sphere passing through point cloud using CGAL. @@ -49,8 +52,8 @@ class Sphere_circumradius { * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. * */ - double operator()(const Point_cloud& point_cloud) const { - return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); + Filtration_value operator()(const Point_cloud& point_cloud) const { + return std::sqrt(cast_to_fv(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); } }; diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index ea32f596..f5980e6d 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -107,11 +107,11 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { std::clog << vertex << ","; vp.push_back(points.at(vertex)); } - std::clog << ") - distance =" << Gudhi::cech_complex::Sphere_circumradius()(vp.at(0), vp.at(1)) + std::clog << ") - distance =" << Gudhi::cech_complex::Sphere_circumradius()(vp.at(0), vp.at(1)) << " - filtration =" << st.filtration(f_simplex) << std::endl; BOOST_CHECK(vp.size() == 2); GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), - Gudhi::cech_complex::Sphere_circumradius()(vp.at(0), vp.at(1))); + Gudhi::cech_complex::Sphere_circumradius()(vp.at(0), vp.at(1))); } } -- cgit v1.2.3