Modify the algorithm to get the minimal enclosing ball

1 files changed, 102 insertions, 65 deletions

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 <iostream>
 #include <vector>
+#include <set>
+#include <map>
+#include <algorithm>
 #include <cmath>  // 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<Point_d, FT>;
-  
-  
-  // Add an int in TDS to save point index in the structure
-//   using TDS = CGAL::Triangulation_data_structure<typename Kernel::Dimension,
-//                                                  CGAL::Triangulation_vertex<Kernel, std::ptrdiff_t>,
-//                                                  CGAL::Triangulation_full_cell<Kernel> >;
-// 
-//   /** \brief A (Weighted or not) Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/
-//   using Triangulation = CGAL::Delaunay_triangulation<Kernel, TDS>;
-//   // 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<typename Sphere>
+  class CompareSpheresRadii
+  {
+    public:
+      CGAL::NT_converter<FT, double> 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_d>;
-    Point_cloud points;
+    CGAL::NT_converter<FT, double> 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 <Filtration_value> enclosing_ball_radii;
+    std::set <Sphere, CompareSpheresRadii<Sphere>> 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<FT, double> cast_to_double;
-//     CGAL::NT_converter<typename Cech_complex::Point, Point_d> 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<std::string, Sphere> cache_;
 };
 
 }  // namespace cech_complex