First version of cech enhancement

10 files changed, 311 insertions, 64 deletions

diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
index e489e8a4..c332c656 100644
--- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
+++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp
@@ -17,6 +17,8 @@
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Miniball.hpp>
 
+#include <CGAL/Epeck_d.h>  // For EXACT or SAFE version
+
 #include "boost/filesystem.hpp"  // includes all needed Boost.Filesystem declarations
 
 #include <string>
@@ -30,7 +32,8 @@ using Point_cloud = std::vector<Point>;
 using Points_off_reader = Gudhi::Points_off_reader<Point>;
 using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>;
 using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
-using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>;
+using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel>;
 
 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 1a1f708c..ac17fc73 100644
--- a/src/Cech_complex/example/cech_complex_example_from_points.cpp
+++ b/src/Cech_complex/example/cech_complex_example_from_points.cpp
@@ -1,6 +1,8 @@
 #include <gudhi/Cech_complex.h>
 #include <gudhi/Simplex_tree.h>
 
+#include <CGAL/Epeck_d.h>  // For EXACT or SAFE version
+
 #include <iostream>
 #include <string>
 #include <vector>
@@ -8,32 +10,71 @@
 
 int main() {
   // Type definitions
-  using Point_cloud = std::vector<std::array<double, 2>>;
+//   using Point_cloud = std::vector<std::array<double, 2>>;
   using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
   using Filtration_value = Simplex_tree::Filtration_value;
-  using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>;
+  using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+  using FT = typename Kernel::FT;
+  using Point = typename Kernel::Point_d;
+  using Point_cloud = std::vector<Point>;
+  using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel>;
 
   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
+//   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<FT> point({0.0, 0.0, 0.0, 0.0});
+//   points.emplace_back(point.begin(), point.end());
+
+  std::vector<FT> point0({1., 0.});
+  points.emplace_back(point0.begin(), point0.end());
+  std::vector<FT> point1({0., 1.});
+  points.emplace_back(point1.begin(), point1.end());
+  std::vector<FT> point2({2., 1.});
+  points.emplace_back(point2.begin(), point2.end());
+  std::vector<FT> point3({3., 2.});
+  points.emplace_back(point3.begin(), point3.end());
+  std::vector<FT> point4({0., 3.});
+  points.emplace_back(point4.begin(), point4.end());
+  std::vector<FT> point5({3. + std::sqrt(3.), 3.});
+  points.emplace_back(point5.begin(), point5.end());
+
+//   points.emplace_back(Point(std::vector<FT>({1., 0.})));
+//   points.emplace_back(Point(std::vector<FT>({0., 1.})));
+//   points.emplace_back(Point(std::vector<FT>({2., 1.})));
+//   points.emplace_back(Point(std::vector<FT>({3., 2.})));
+//   points.emplace_back(Point(std::vector<FT>({0., 3.})));
+//   points.emplace_back(Point(std::vector<FT>({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 = 1.;
+  Filtration_value max_radius = 10.;
+  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, 2);
+  cech_complex_from_points.create_complex(stree, 3);
   // ----------------------------------------------------------------------------
   // 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 60ae9712..ac08e6cc 100644
--- a/src/Cech_complex/example/cech_complex_step_by_step.cpp
+++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp
@@ -13,7 +13,9 @@
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Points_off_io.h>
 
-#include <gudhi/Miniball.hpp>
+#include <gudhi/Miniball.hpp> // TODO to remove ?
+
+#include <CGAL/Epeck_d.h>
 
 #include <boost/program_options.hpp>
 
@@ -34,16 +36,18 @@
 using Simplex_tree = Gudhi::Simplex_tree<>;
 using Simplex_handle = Simplex_tree::Simplex_handle;
 using Filtration_value = Simplex_tree::Filtration_value;
-using Point = std::vector<double>;
+// using Point = std::vector<double>;
+using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+using Point = typename Kernel::Point_d;
 using Points_off_reader = Gudhi::Points_off_reader<Point>;
 using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>;
 
 class Cech_blocker {
  private:
   using Point_cloud = std::vector<Point>;
-  using Point_iterator = Point_cloud::const_iterator;
-  using Coordinate_iterator = Point::const_iterator;
-  using Min_sphere = Gudhi::Miniball::Miniball<Gudhi::Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>;
+//   using Point_iterator = Point_cloud::const_iterator;
+//   using Coordinate_iterator = Point::const_iterator;
+//   using Min_sphere = Gudhi::Miniball::Miniball<Gudhi::Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>;
 
  public:
   bool operator()(Simplex_handle sh) {
@@ -63,14 +67,14 @@ class Cech_blocker {
   }
   Cech_blocker(Simplex_tree& simplex_tree, Filtration_value max_radius, const std::vector<Point>& point_cloud)
       : simplex_tree_(simplex_tree), max_radius_(max_radius), point_cloud_(point_cloud) {
-    dimension_ = point_cloud_[0].size();
+//     dimension_ = point_cloud_[0].size();
   }
 
  private:
   Simplex_tree simplex_tree_;
   Filtration_value max_radius_;
   std::vector<Point> point_cloud_;
-  int dimension_;
+//   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 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 <typename SimplicialComplexForProximityGraph, typename ForwardPointRange>
+template <typename SimplicialComplexForProximityGraph, typename ForwardPointRange, typename Kernel>
 class Cech_complex {
  private:
   // Required by compute_proximity_graph
@@ -49,14 +49,15 @@ class Cech_complex {
   using Proximity_graph = Gudhi::Proximity_graph<SimplicialComplexForProximityGraph>;
 
   // Retrieve Coordinate type from ForwardPointRange
-  using Point_from_range_iterator = typename boost::range_const_iterator<ForwardPointRange>::type;
-  using Point_from_range = typename std::iterator_traits<Point_from_range_iterator>::value_type;
-  using Coordinate_iterator = typename boost::range_const_iterator<Point_from_range>::type;
-  using Coordinate = typename std::iterator_traits<Coordinate_iterator>::value_type;
+//   using Point_from_range_iterator = typename boost::range_const_iterator<ForwardPointRange>::type;
+//   using Point_from_range = typename std::iterator_traits<Point_from_range_iterator>::value_type;
+//   using Coordinate_iterator = typename boost::range_const_iterator<Point_from_range>::type;
+//   using Coordinate = typename std::iterator_traits<Coordinate_iterator>::value_type;
 
  public:
   // Point and Point_cloud type definition
-  using Point = std::vector<Coordinate>;
+  //using Point = std::vector<Coordinate>;
+  using Point = typename Kernel::Point_d;
   using Point_cloud = std::vector<Point>;
 
  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<SimplicialComplexForProximityGraph>(
         point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius());
   }
@@ -87,14 +92,17 @@ class Cech_complex {
    */
   template <typename SimplicialComplexForCechComplex>
   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<SimplicialComplexForCechComplex, Cech_complex>(&complex, this));
+                                    Cech_blocker<SimplicialComplexForCechComplex, Cech_complex, Kernel>(&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 <gudhi/distance_functions.h>  // for Gudhi::Minimal_enclosing_ball_radius
 
+#include <CGAL/number_utils.h> // for CGAL::to_double
+#include <CGAL/NT_converter.h> //
+// #include <CGAL/Delaunay_triangulation.h>
+#include <CGAL/Epeck_d.h>  // For EXACT or SAFE version
+#include <CGAL/Spatial_sort_traits_adapter_d.h>
+
 #include <iostream>
 #include <vector>
 #include <cmath>  // for std::sqrt
@@ -35,28 +42,120 @@ namespace cech_complex {
  *
  * \tparam Chech_complex is required by the blocker.
  */
-template <typename SimplicialComplexForCech, typename Cech_complex>
+template <typename SimplicialComplexForCech, typename Cech_complex, typename Kernel>
 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<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] 
+   *  \return  */
+  template<class PointIterator>
+  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<class PointIterator>
+  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_d>;
     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<FT, double> cast_to_double;
+//     CGAL::NT_converter<typename Cech_complex::Point, Point_d> 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
diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp
index 6e00d7b5..81efd6ae 100644
--- a/src/Cech_complex/test/test_cech_complex.cpp
+++ b/src/Cech_complex/test/test_cech_complex.cpp
@@ -26,17 +26,24 @@
 #include <gudhi/Unitary_tests_utils.h>
 #include <gudhi/Miniball.hpp>
 
+#include <CGAL/Epeck_d.h>  // For EXACT or SAFE version
+
+
 // Type definitions
 using Simplex_tree = Gudhi::Simplex_tree<>;
 using Filtration_value = Simplex_tree::Filtration_value;
-using Point = std::vector<Filtration_value>;
+//using Point = std::vector<Filtration_value>;
+using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+using FT = typename Kernel::FT;
+using Point = typename Kernel::Point_d;
+
 using Point_cloud = std::vector<Point>;
 using Points_off_reader = Gudhi::Points_off_reader<Point>;
-using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>;
+using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel>;
 
-using Point_iterator = Point_cloud::const_iterator;
-using Coordinate_iterator = Point::const_iterator;
-using Min_sphere = Gudhi::Miniball::Miniball<Gudhi::Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>;
+// using Point_iterator = Point_cloud::const_iterator;
+// using Coordinate_iterator = Point::const_iterator;
+// using Min_sphere = Gudhi::Miniball::Miniball<Gudhi::Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>;
 
 BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) {
   // ----------------------------------------------------------------------------
@@ -45,17 +52,41 @@ 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
+//   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<FT> point0({1., 0.});
+  points.emplace_back(point0.begin(), point0.end());
+  std::vector<FT> point1({0., 1.});
+  points.emplace_back(point1.begin(), point1.end());
+  std::vector<FT> point2({2., 1.});
+  points.emplace_back(point2.begin(), point2.end());
+  std::vector<FT> point3({3., 2.});
+  points.emplace_back(point3.begin(), point3.end());
+  std::vector<FT> point4({0., 3.});
+  points.emplace_back(point4.begin(), point4.end());
+  std::vector<FT> point5({3. + std::sqrt(3.), 3.});
+  points.emplace_back(point5.begin(), point5.end());
+  std::vector<FT> point6({1., 4.});
+  points.emplace_back(point6.begin(), point6.end());
+  std::vector<FT> point7({3., 4.});
+  points.emplace_back(point7.begin(), point7.end());
+  std::vector<FT> point8({2., 4. + std::sqrt(3.)});
+  points.emplace_back(point8.begin(), point8.end());
+  std::vector<FT> point9({0., 4.});
+  points.emplace_back(point9.begin(), point9.end());
+  std::vector<FT> point10({-0.5, 2.});
+  points.emplace_back(point10.begin(), point10.end());
 
   Filtration_value max_radius = 1.0;
   std::clog << "========== NUMBER OF POINTS = " << points.size() << " - Cech max_radius = " << max_radius
@@ -125,35 +156,38 @@ 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());
+//   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}));
-  std::clog << "f012= " << f012 << " | ms012_radius= " << std::sqrt(ms012.squared_radius()) << std::endl;
+  CGAL::NT_converter<FT, double> 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;
+
 
-  GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(ms012.squared_radius()));
+  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());
+//   Min_sphere ms1410(dimension, points1410.begin(), points1410.end());
 
   Simplex_tree::Filtration_value f1410 = st2.filtration(st2.find({1, 4, 10}));
-  std::clog << "f1410= " << f1410 << " | ms1410_radius= " << std::sqrt(ms1410.squared_radius()) << std::endl;
+  std::clog << "f1410= " << f1410 << " | points1410_radius= " << std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points1410.begin(), points1410.end()))) << std::endl;
 
-  GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, std::sqrt(ms1410.squared_radius()));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points1410.begin(), points1410.end()))));
 
   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());
+//   Min_sphere ms469(dimension, points469.begin(), points469.end());
 
   Simplex_tree::Filtration_value f469 = st2.filtration(st2.find({4, 6, 9}));
-  std::clog << "f469= " << f469 << " | ms469_radius= " << std::sqrt(ms469.squared_radius()) << std::endl;
+  std::clog << "f469= " << f469 << " | points469_radius= " << std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points469.begin(), points469.end()))) << std::endl;
 
-  GUDHI_TEST_FLOAT_EQUALITY_CHECK(f469, std::sqrt(ms469.squared_radius()));
+  GUDHI_TEST_FLOAT_EQUALITY_CHECK(f469, std::sqrt(cast_to_double(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()));
diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp
index daea08e2..ccd7d453 100644
--- a/src/Cech_complex/utilities/cech_persistence.cpp
+++ b/src/Cech_complex/utilities/cech_persistence.cpp
@@ -16,6 +16,8 @@
 
 #include <boost/program_options.hpp>
 
+#include <CGAL/Epeck_d.h>  // For EXACT or SAFE version
+
 #include <string>
 #include <vector>
 #include <limits>  // infinity
@@ -23,10 +25,14 @@
 // Types definition
 using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
 using Filtration_value = Simplex_tree::Filtration_value;
-using Point = std::vector<double>;
+// using Point = std::vector<double>;
+// using Point_cloud = std::vector<Point>;
+
+using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+using Point = typename Kernel::Point_d;
 using Point_cloud = std::vector<Point>;
 using Points_off_reader = Gudhi::Points_off_reader<Point>;
-using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>;
+using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel>;
 using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
 using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>;
 
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index 85790baf..f69ed6ec 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -1278,8 +1278,10 @@ 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
@@ -1288,10 +1290,12 @@ 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 9bbc62b7..ae5168aa 100644
--- a/src/common/include/gudhi/distance_functions.h
+++ b/src/common/include/gudhi/distance_functions.h
@@ -18,6 +18,8 @@
 #include <boost/range/metafunctions.hpp>
 #include <boost/range/size.hpp>
 
+#include <CGAL/Epeck_d.h>
+
 #include <cmath>  // for std::sqrt
 #include <type_traits>  // for std::decay
 #include <iterator>  // for std::begin, std::end
@@ -63,6 +65,23 @@ class Euclidean_distance {
  * The points are assumed to have the same dimension. */
 class Minimal_enclosing_ball_radius {
  public:
+   /** \brief TODO
+   *
+   * @param[in] point_1
+   * @param[in] point_2
+   * @return
+   * \tparam Point
+   *
+   */
+  //typename FT = typename Kernel::FT,
+  template< typename Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>,
+            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.;
+  }
+
   /** \brief Minimal_enclosing_ball_radius from two points.
    *
    * @param[in] point_1 First point.
@@ -75,8 +94,29 @@ class Minimal_enclosing_ball_radius {
   template< typename Point >
   typename std::iterator_traits<typename boost::range_iterator<Point>::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
+   *
+   * @param[in] point_cloud The points.
+   * @return
+   * \tparam Point_cloud
+   *
+   */
+  //typename FT = typename Kernel::FT,
+  template< typename Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>,
+            typename Point= typename Kernel::Point_d,
+            typename Point_cloud = std::vector<Point>>
+  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())));
+  }
+
   /** \brief Minimal_enclosing_ball_radius from a point cloud.
    *
    * @param[in] point_cloud The points.
@@ -93,6 +133,8 @@ class Minimal_enclosing_ball_radius {
             typename Coordinate = typename std::iterator_traits<Coordinate_iterator>::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<Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>;
 
     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 da9dee7d..9190182c 100644
--- a/src/common/include/gudhi/graph_simplicial_complex.h
+++ b/src/common/include/gudhi/graph_simplicial_complex.h
@@ -18,6 +18,8 @@
 #include <map>
 #include <tuple>  // for std::tie
 
+#include <iostream>
+
 namespace Gudhi {
 /** @file
  * @brief Graph simplicial complex methods
@@ -76,6 +78,7 @@ Proximity_graph<SimplicialComplexForProximityGraph> 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);