Add Cech complex. Do not compile yet.

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/cechcomplex_vincent@3250 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: bef87ed8038444685b964175ea65860300917380

9 files changed, 1020 insertions, 0 deletions

diff --git a/src/Cech_complex/doc/Intro_cech_complex.h b/src/Cech_complex/doc/Intro_cech_complex.h
new file mode 100644
index 00000000..e69de29b
--- /dev/null
+++ b/src/Cech_complex/doc/Intro_cech_complex.h
diff --git a/src/Cech_complex/example/CMakeLists.txt b/src/Cech_complex/example/CMakeLists.txt
new file mode 100644
index 00000000..8097871f
--- /dev/null
+++ b/src/Cech_complex/example/CMakeLists.txt
@@ -0,0 +1,14 @@
+cmake_minimum_required(VERSION 2.6)
+project(Cech_complex_examples)
+
+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_LIBRARY})
+if (TBB_FOUND)
+  target_link_libraries(Cech_complex_example_step_by_step ${TBB_LIBRARIES})
+endif()
+
+add_executable ( Cech_complex_example_one_skeleton_from_points example_one_skeleton_cech_from_points.cpp)
+if (TBB_FOUND)
+  target_link_libraries(Cech_complex_example_one_skeleton_from_points ${TBB_LIBRARIES})
+endif()
+add_test(NAME Cech_complex_example_one_skeleton_from_points COMMAND $<TARGET_FILE:Cech_complex_example_one_skeleton_from_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
new file mode 100644
index 00000000..e71086b6
--- /dev/null
+++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp
@@ -0,0 +1,179 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2018 Inria
+ *
+ *    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 <http://www.gnu.org/licenses/>.
+ */
+
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Points_off_io.h>
+
+#include <Miniball/Miniball.hpp>
+
+#include <boost/program_options.hpp>
+
+#include <string>
+#include <vector>
+#include <limits>  // infinity
+#include <utility>  // for pair
+#include <map>
+
+// ----------------------------------------------------------------------------
+// rips_persistence_step_by_step is an example of each step that is required to
+// build a Rips over a Simplex_tree. Please refer to rips_persistence to see
+// how to do the same thing with the Rips_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 Point = std::vector<double>;
+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 =  Miniball::Miniball <Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>;
+ public:
+  bool operator()(Simplex_handle sh) {
+    std::vector<Point> points;
+    for (auto vertex : simplex_tree_.simplex_vertex_range(sh)) {
+      points.push_back(point_cloud_[vertex]);
+#ifdef DEBUG_TRACES
+      std::cout << "#(" << vertex << ")#";
+#endif  // DEBUG_TRACES
+    }
+    Min_sphere ms(dimension_, points.begin(),points.end());
+    Filtration_value radius = std::sqrt(ms.squared_radius());
+#ifdef DEBUG_TRACES
+    std::cout << "radius = " << radius << " - " << (radius > threshold_) << std::endl;
+#endif  // DEBUG_TRACES
+    simplex_tree_.assign_filtration(sh, radius);
+    return (radius > threshold_);
+  }
+  Cech_blocker(Simplex_tree& simplex_tree, Filtration_value threshold, const std::vector<Point>& point_cloud)
+    : simplex_tree_(simplex_tree),
+      threshold_(threshold),
+      point_cloud_(point_cloud) {
+    dimension_ = point_cloud_[0].size();
+  }
+ private:
+  Simplex_tree simplex_tree_;
+  Filtration_value threshold_;
+  std::vector<Point> point_cloud_;
+  int dimension_;
+};
+
+
+void program_options(int argc, char * argv[]
+                     , std::string & off_file_points
+                     , Filtration_value & threshold
+                     , int & dim_max);
+
+
+int main(int argc, char * argv[]) {
+  std::string off_file_points;
+  Filtration_value threshold;
+  int dim_max;
+
+  program_options(argc, argv, off_file_points, threshold, 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<Simplex_tree>(off_reader.get_point_cloud(),
+                                                                            threshold,
+                                                                            Gudhi::Euclidean_distance());
+
+  // Construct the Rips 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, threshold, off_reader.get_point_cloud()));
+
+  std::cout << "The complex contains " << st.num_simplices() << " simplices \n";
+  std::cout << "   and has dimension " << st.dimension() << " \n";
+
+  // Sort the simplices in the order of the filtration
+  st.initialize_filtration();
+
+#if DEBUG_TRACES
+  std::cout << "********************************************************************\n";
+  // Display the Simplex_tree - Can not be done in the middle of 2 inserts
+  std::cout << "* The complex contains " << st.num_simplices() << " simplices - dimension=" << st.dimension() << "\n";
+  std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
+  for (auto f_simplex : st.filtration_simplex_range()) {
+    std::cout << "   " << "[" << st.filtration(f_simplex) << "] ";
+    for (auto vertex : st.simplex_vertex_range(f_simplex)) {
+      std::cout << static_cast<int>(vertex) << " ";
+    }
+    std::cout << std::endl;
+  }
+#endif  // DEBUG_TRACES
+
+  return 0;
+}
+
+void program_options(int argc, char * argv[]
+                     , std::string & off_file_points
+                     , Filtration_value & threshold
+                     , int & dim_max) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()
+      ("input-file", po::value<std::string>(&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-edge-length,r",
+       po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()),
+       "Maximal length of an edge for the Rips complex construction.")
+      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
+       "Maximal dimension of the Rips 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::cout << std::endl;
+    std::cout << "Construct a Cech complex defined on a set of input points.\n \n";
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
diff --git a/src/Cech_complex/example/example_one_skeleton_cech_from_points.cpp b/src/Cech_complex/example/example_one_skeleton_cech_from_points.cpp
new file mode 100644
index 00000000..9b03616c
--- /dev/null
+++ b/src/Cech_complex/example/example_one_skeleton_cech_from_points.cpp
@@ -0,0 +1,74 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2018 Inria
+ *
+ *    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 <http://www.gnu.org/licenses/>.
+ */
+
+#include <gudhi/Cech_complex.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/distance_functions.h>
+
+#include <iostream>
+#include <string>
+#include <vector>
+#include <limits>  // for std::numeric_limits
+
+int main() {
+  // Type definitions
+  using Point_cloud = std::vector<std::vector<double>>;
+  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>;
+
+  Point_cloud points;
+  points.push_back({1.0, 1.0});
+  points.push_back({7.0, 0.0});
+  points.push_back({4.0, 6.0});
+  points.push_back({9.0, 6.0});
+  points.push_back({0.0, 14.0});
+  points.push_back({2.0, 19.0});
+  points.push_back({9.0, 17.0});
+
+  // ----------------------------------------------------------------------------
+  // Init of a Rips complex from points
+  // ----------------------------------------------------------------------------
+  Filtration_value threshold = 12.0;
+  Cech_complex cech_complex_from_points(points, threshold, Gudhi::Euclidean_distance());
+
+  Simplex_tree stree;
+  cech_complex_from_points.create_complex(stree, 2);
+  // ----------------------------------------------------------------------------
+  // Display information about the one skeleton Rips complex
+  // ----------------------------------------------------------------------------
+  std::cout << "Cech complex is of dimension " << stree.dimension() <<
+               " - " << stree.num_simplices() << " simplices - " <<
+               stree.num_vertices() << " vertices." << std::endl;
+
+  std::cout << "Iterator on Cech complex simplices in the filtration order, with [filtration value]:" <<
+               std::endl;
+  for (auto f_simplex : stree.filtration_simplex_range()) {
+    std::cout << "   ( ";
+    for (auto vertex : stree.simplex_vertex_range(f_simplex)) {
+      std::cout << vertex << " ";
+    }
+    std::cout << ") -> " << "[" << stree.filtration(f_simplex) << "] ";
+    std::cout << std::endl;
+  }
+  return 0;
+}
diff --git a/src/Cech_complex/include/Miniball/Miniball.COPYRIGHT b/src/Cech_complex/include/Miniball/Miniball.COPYRIGHT
new file mode 100644
index 00000000..dbe4c553
--- /dev/null
+++ b/src/Cech_complex/include/Miniball/Miniball.COPYRIGHT
@@ -0,0 +1,4 @@
+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/Miniball/Miniball.README b/src/Cech_complex/include/Miniball/Miniball.README
new file mode 100644
index 00000000..86a96f08
--- /dev/null
+++ b/src/Cech_complex/include/Miniball/Miniball.README
@@ -0,0 +1,23 @@
+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. 
+
+
diff --git a/src/Cech_complex/include/Miniball/Miniball.hpp b/src/Cech_complex/include/Miniball/Miniball.hpp
new file mode 100644
index 00000000..cb76c534
--- /dev/null
+++ b/src/Cech_complex/include/Miniball/Miniball.hpp
@@ -0,0 +1,515 @@
+//    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 <http://www.gnu.org/licenses/>.
+//
+//    Contact:
+//    --------
+//    Bernd Gaertner
+//    Institute of Theoretical Computer Science 
+//    ETH Zuerich
+//    CAB G31.1
+//    CH-8092 Zuerich, Switzerland
+//    http://www.inf.ethz.ch/personal/gaertner
+
+#include <cassert>
+#include <algorithm>
+#include <list>
+#include <ctime>
+#include <limits>
+
+namespace Miniball {
+
+  // Global Functions
+  // ================
+  template <typename NT>
+  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<Pit_, Cit_*> {
+    typedef Pit_                                  Pit;      
+    typedef Cit_*                                 Cit;
+    inline  Cit operator() (Pit it) const { return *it; }
+  };
+
+  // Class Declaration
+  // =================
+
+  template <typename CoordAccessor>
+  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<Cit>::value_type      NT;  
+    // The iterator to go through the support points
+    typedef typename std::list<Pit>::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<Pit> 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<Pit>::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<NT>::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 <typename CoordAccessor>
+  Miniball<CoordAccessor>::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<d; ++j) c[0][j] = nt0;
+    current_c = c[0];
+
+    // compute miniball
+    pivot_mb (points_end);
+
+    // update time
+    time = (clock() - time) / CLOCKS_PER_SEC;
+  }
+  
+  template <typename CoordAccessor>
+  Miniball<CoordAccessor>::~Miniball()
+  {
+    delete_arrays();
+  }
+
+  template <typename CoordAccessor>
+  void Miniball<CoordAccessor>::create_arrays() 
+  {
+    c = new NT*[d+1]; 
+    v = new NT*[d+1]; 
+    a = new NT*[d+1];
+    for (int i=0; i<d+1; ++i) {
+      c[i] = new NT[d];
+      v[i] = new NT[d];
+      a[i] = new NT[d];
+    }
+    sqr_r = new NT[d+1];
+    q0 = new NT[d];
+    z = new NT[d+1];
+    f = new NT[d+1];
+  } 
+
+  template <typename CoordAccessor>
+  void Miniball<CoordAccessor>::delete_arrays() 
+  {
+    delete[] f;
+    delete[] z;
+    delete[] q0;
+    delete[] sqr_r;
+    for (int i=0; i<d+1; ++i) {
+      delete[] a[i];
+      delete[] v[i];
+      delete[] c[i];
+    }
+    delete[] a;
+    delete[] v;
+    delete[] c;
+  }
+
+  template <typename CoordAccessor>
+  const typename Miniball<CoordAccessor>::NT* 
+  Miniball<CoordAccessor>::center () const 
+  {
+    return current_c;
+  }
+ 
+  template <typename CoordAccessor>
+  typename Miniball<CoordAccessor>::NT 
+  Miniball<CoordAccessor>::squared_radius () const 
+  {
+    return current_sqr_r;
+  }
+
+  template <typename CoordAccessor>
+  int Miniball<CoordAccessor>::nr_support_points () const 
+  {
+    assert (ssize < d+2);
+    return ssize;
+  } 
+
+  template <typename CoordAccessor>  
+  typename Miniball<CoordAccessor>::SupportPointIterator 
+  Miniball<CoordAccessor>::support_points_begin () const 
+  {
+    return L.begin();
+  }
+
+  template <typename CoordAccessor>  
+  typename Miniball<CoordAccessor>::SupportPointIterator 
+  Miniball<CoordAccessor>::support_points_end () const  
+  {
+    return support_end;
+  }
+
+  template <typename CoordAccessor>  
+  typename Miniball<CoordAccessor>::NT 
+  Miniball<CoordAccessor>::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 <typename CoordAccessor>  
+  bool Miniball<CoordAccessor>::is_valid (NT tol) const
+  {
+    NT suboptimality;
+    return ( (relative_error (suboptimality) <= tol) && (suboptimality == 0) );
+  }
+
+  template <typename CoordAccessor>  
+  double Miniball<CoordAccessor>::get_time() const 
+  {
+    return time;
+  }
+
+  template <typename CoordAccessor>  
+  void Miniball<CoordAccessor>::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 <typename CoordAccessor>  
+  void Miniball<CoordAccessor>::mtf_move_to_front (Sit j) 
+  {
+    if (support_end == j)
+      support_end++;
+    L.splice (L.begin(), L, j);
+  }
+
+  template <typename CoordAccessor>  
+  void Miniball<CoordAccessor>::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<d; ++j)
+	  e += mb_sqr<NT>(*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 <typename CoordAccessor>  
+  void Miniball<CoordAccessor>::pivot_move_to_front (Pit j) 
+  {
+    L.push_front(j);
+    if (std::distance(L.begin(), support_end) == d+2)
+      support_end--;
+  }
+
+  template <typename CoordAccessor>  
+  inline typename Miniball<CoordAccessor>::NT 
+  Miniball<CoordAccessor>::excess (Pit pit) const 
+  {
+    Cit p = coord_accessor(pit);
+    NT e = -current_sqr_r;
+    NT* c = current_c;
+    for (int k=0; k<d; ++k){
+      e += mb_sqr<NT>(*p++-*c++);
+    }
+    return e;
+  }
+
+  template <typename CoordAccessor>  
+  void Miniball<CoordAccessor>::pop ()  
+  {
+    --fsize;
+  }
+
+  template <typename CoordAccessor>  
+  bool Miniball<CoordAccessor>::push (Pit pit) 
+  {
+    int i, j;
+    NT eps = mb_sqr<NT>(std::numeric_limits<NT>::epsilon());
+   
+    Cit cit = coord_accessor(pit);
+    Cit p = cit;
+ 
+    if (fsize==0) {
+      for (i=0; i<d; ++i)
+	q0[i] = *p++;
+      for (i=0; i<d; ++i)
+	c[0][i] = q0[i];
+      sqr_r[0] = nt0;
+    } 
+    else {
+      // set v_fsize to Q_fsize
+      for (i=0; i<d; ++i)
+	//v[fsize][i] = p[i]-q0[i];
+	v[fsize][i] = *p++-q0[i];
+      
+      // compute the a_{fsize,i}, i< fsize
+      for (i=1; i<fsize; ++i) {
+	a[fsize][i] = nt0;
+	for (j=0; j<d; ++j)
+	  a[fsize][i] += v[i][j] * v[fsize][j];
+	a[fsize][i]*=(2/z[i]);
+      }
+      
+      // update v_fsize to Q_fsize-\bar{Q}_fsize
+      for (i=1; i<fsize; ++i) {
+	for (j=0; j<d; ++j)
+	  v[fsize][j] -= a[fsize][i]*v[i][j];
+      }
+      
+      // compute z_fsize
+      z[fsize]=nt0;
+      for (j=0; j<d; ++j)
+	z[fsize] += mb_sqr<NT>(v[fsize][j]);
+      z[fsize]*=2;
+      
+      // reject push if z_fsize too small
+      if (z[fsize]<eps*current_sqr_r) {
+	return false;
+      }
+      
+      // update c, sqr_r
+      p=cit;
+      NT e = -sqr_r[fsize-1];
+      for (i=0; i<d; ++i)
+	e += mb_sqr<NT>(*p++-c[fsize-1][i]);
+      f[fsize]=e/z[fsize];
+      
+      for (i=0; i<d; ++i)
+	c[fsize][i] = c[fsize-1][i]+f[fsize]*v[fsize][i];
+      sqr_r[fsize] = sqr_r[fsize-1] + e*f[fsize]/2;
+    }
+    current_c = c[fsize];
+    current_sqr_r = sqr_r[fsize];
+    ssize = ++fsize;
+    return true;
+  }
+ 
+  template <typename CoordAccessor>  
+  typename Miniball<CoordAccessor>::NT 
+  Miniball<CoordAccessor>::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;
+  }
+
+} // end Namespace Miniball
diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h
new file mode 100644
index 00000000..3a0d828a
--- /dev/null
+++ b/src/Cech_complex/include/gudhi/Cech_complex.h
@@ -0,0 +1,126 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2018 Inria
+ *
+ *    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 <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef CECH_COMPLEX_H_
+#define CECH_COMPLEX_H_
+
+#include <gudhi/graph_simplicial_complex.h>  // for Gudhi::Proximity_graph
+#include <gudhi/Debug_utils.h>               // for GUDHI_CHECK
+#include <gudhi/Cech_complex_blocker.h>      // for Gudhi::cech_complex::Cech_blocker
+
+#include <iostream>
+#include <vector>
+#include <cstddef>
+#include <iterator>  // for std::size
+
+namespace Gudhi {
+
+namespace cech_complex {
+
+/**
+ * \class Cech_complex
+ * \brief Cech complex data structure.
+ * 
+ * \ingroup Cech_complex
+ * 
+ * \details
+ * The data structure is a one skeleton graph, or Rips graph, containing edges when the edge length is less or equal
+ * to a given threshold. Edge length is computed from a user given point cloud with a given distance function, or a
+ * distance matrix.
+ * 
+ * \tparam Filtration_value is the type used to store the filtration values of the simplicial complex.
+ */
+template<typename SimplicialComplexForCechComplex, typename ForwardPointRange>
+class Cech_complex {
+ private:
+  using Vertex_handle = typename SimplicialComplexForCechComplex::Vertex_handle;
+  using Filtration_value = typename SimplicialComplexForCechComplex::Filtration_value;
+  using Proximity_graph = Gudhi::Proximity_graph<SimplicialComplexForCechComplex>;
+
+ public:
+  /** \brief Cech_complex constructor from a list of points.
+   *
+   * @param[in] points Range of points.
+   * @param[in] threshold Rips value.
+   * @param[in] distance distance function that returns a `Filtration_value` from 2 given points.
+   * 
+   * \tparam ForwardPointRange must be a range for which `std::begin` and `std::end` return input iterators on a
+   * point.
+   *
+   * \tparam Distance furnishes `operator()(const Point& p1, const Point& p2)`, where
+   * `Point` is a point from the `ForwardPointRange`, and that returns a `Filtration_value`.
+   */
+  template<typename Distance >
+  Cech_complex(const ForwardPointRange& points, Filtration_value threshold, Distance distance)
+    : threshold_(threshold),
+      point_cloud_(points) {
+    GUDHI_CHECK(std::size(points) > 0,
+                std::invalid_argument("Cech_complex::create_complex - point cloud is empty"));
+    dimension_ = points[0].size();
+    cech_skeleton_graph_ = Gudhi::compute_proximity_graph<SimplicialComplexForCechComplex>(point_cloud_, threshold_, distance);
+  }
+
+  /** \brief Initializes the simplicial complex from the Rips graph and expands it until a given maximal
+   * dimension.
+   *
+   * @param[in] complex SimplicialComplexForCech to be created.
+   * @param[in] dim_max graph expansion for Rips until this given maximal dimension.
+   * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0.
+   *
+   */
+  void create_complex(SimplicialComplexForCechComplex& complex, int dim_max) {
+    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<SimplicialComplexForCechComplex, ForwardPointRange>(complex, this));
+  }
+
+  Filtration_value threshold() const {
+    return threshold_;
+  }
+
+  std::size_t dimension() const {
+    return dimension_;
+  }
+
+  auto point(std::size_t vertex) const -> decltype(point_cloud_.begin() + vertex) {
+    GUDHI_CHECK((point_cloud_.begin() + vertex) < point_cloud_.end(),
+                std::invalid_argument("Cech_complex::point - simplicial complex is not empty"));
+    return (point_cloud_.begin() + vertex);
+  }
+
+ private:
+  Proximity_graph cech_skeleton_graph_;
+  Filtration_value threshold_;
+  ForwardPointRange point_cloud_;
+  std::size_t dimension_;
+};
+
+}  // namespace cech_complex
+
+}  // namespace Gudhi
+
+#endif  // CECH_COMPLEX_H_
diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h
new file mode 100644
index 00000000..647bf0b7
--- /dev/null
+++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h
@@ -0,0 +1,85 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2018 Inria
+ *
+ *    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 <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef CECH_COMPLEX_BLOCKER_H_
+#define CECH_COMPLEX_BLOCKER_H_
+
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/Cech_complex.h>
+
+#include <Miniball/Miniball.hpp>
+
+#include <iostream>
+#include <vector>
+#include <cmath>  // for std::sqrt
+
+namespace Gudhi {
+
+namespace cech_complex {
+
+// Just declaring Cech_complex class because used and not yet defined.
+template<typename SimplicialComplexForCechComplex, typename ForwardPointRange>
+class Cech_complex;
+
+template <typename SimplicialComplexForCech, typename ForwardPointRange>
+class Cech_blocker {
+ private:
+  using Point = std::vector<double>;
+  using Point_cloud = std::vector<Point>;
+  using Point_iterator = Point_cloud::const_iterator;
+  using Coordinate_iterator = Point::const_iterator;
+  using Min_sphere = Miniball::Miniball<Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>;
+  using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle;
+  using Filtration_value = typename SimplicialComplexForCech::Filtration_value;
+  using Cech_complex = Gudhi::cech_complex::Cech_complex<SimplicialComplexForCech, ForwardPointRange>;
+
+ public:
+  bool operator()(Simplex_handle sh) {
+    Point_cloud points;
+    for (auto vertex : simplicial_complex_.simplex_vertex_range(sh)) {
+      points.push_back(cc_ptr_->point(vertex));
+#ifdef DEBUG_TRACES
+      std::cout << "#(" << vertex << ")#";
+#endif  // DEBUG_TRACES
+    }
+    Min_sphere ms(cc_ptr_->dimension(), points.begin(),points.end());
+    Filtration_value radius = std::sqrt(ms.squared_radius());
+#ifdef DEBUG_TRACES
+    std::cout << "radius = " << radius << " - " << (radius > cc_ptr_->threshold()) << std::endl;
+#endif  // DEBUG_TRACES
+    simplicial_complex_.assign_filtration(sh, radius);
+    return (radius > cc_ptr_->threshold());
+  }
+  Cech_blocker(SimplicialComplexForCech& simplicial_complex, Cech_complex* cc_ptr)
+    : simplicial_complex_(simplicial_complex),
+      cc_ptr_(cc_ptr) {
+  }
+ private:
+  SimplicialComplexForCech simplicial_complex_;
+  Cech_complex* cc_ptr_;
+};
+
+}  // namespace cech_complex
+
+}  // namespace Gudhi
+
+#endif  // CECH_COMPLEX_BLOCKER_H_