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/Cech_complex_blocker.h') 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