5 files changed, 242 insertions, 95 deletions

diff --git a/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp b/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp
index 6c9fcdab..fe889ec0 100644
--- a/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp
+++ b/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp
@@ -44,12 +44,12 @@ typedef CGAL::Delaunay_triangulation<K> T;
 // TriangulationDataStructure template parameter
 
 void usage(char * const progName) {
-  std::cerr << "Usage: " << progName << " filename.off dimension" << std::endl;
+  std::cerr << "Usage: " << progName << " inputFile.off dimension outputFile.off" << std::endl;
   exit(-1); // ----- >>
 }
 
 int main(int argc, char **argv) {
-  if (argc != 3) {
+  if (argc != 4) {
     std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl;
     usage(argv[0]);
   }
@@ -68,37 +68,11 @@ int main(int argc, char **argv) {
     std::cerr << "Unable to read file " << offFileName << std::endl;
     exit(-1); // ----- >>
   }
+  
+  std::cout << "number_of_finite_full_cells= " << dt.number_of_finite_full_cells() << std::endl;
 
-  std::cout << "number of vertices=" << dt.number_of_vertices() << std::endl;
-  std::cout << "number of full cells=" << dt.number_of_full_cells() << std::endl;
-  std::cout << "number of finite full cells=" << dt.number_of_finite_full_cells() << std::endl;
-
-  // Points list
-  /*for (T::Vertex_iterator vit = dt.vertices_begin(); vit != dt.vertices_end(); ++vit) {
-        for (auto Coord = vit->point().cartesian_begin(); Coord != vit->point().cartesian_end(); ++Coord) {
-          std::cout << *Coord << " ";
-        }
-        std::cout << std::endl;
-  }
-  std::cout << std::endl;*/
- 
-  /*int i = 0, j = 0;
-  typedef T::Full_cell_iterator Full_cell_iterator;
-  typedef T::Facet Facet;
-
-  for (Full_cell_iterator cit = dt.full_cells_begin(); cit != dt.full_cells_end(); ++cit) {
-    if (!dt.is_infinite(cit)) {
-      j++;
-      continue;
-    }
-    Facet fct(cit, cit->index(dt.infinite_vertex()));
-    i++;
-  }
-  std::cout << "There are " << i << " facets on the convex hull." << std::endl;
-  std::cout << "There are " << j << " facets not on the convex hull." << std::endl;
-*/
-
-  std::string offOutputFile("out.off");
+  std::string outFileName(argv[3]);
+  std::string offOutputFile(outFileName);
   Gudhi::alphashapes::Delaunay_triangulation_off_writer<T> off_writer(offOutputFile, dt);
 
   return 0;
diff --git a/src/Alpha_shapes/include/gudhi/Alpha_shapes.h b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h
index 841c883a..2bc8b221 100644
--- a/src/Alpha_shapes/include/gudhi/Alpha_shapes.h
+++ b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h
@@ -32,6 +32,9 @@
 
 #include <stdio.h>
 #include <stdlib.h>
+#include <math.h>  // isnan, fmax
+
+#include <boost/bimap.hpp>
 
 #include <CGAL/Delaunay_triangulation.h>
 #include <CGAL/Epick_d.h>
@@ -208,12 +211,14 @@ class Alpha_shapes {
     Squared_Radius squared_radius Kinit(compute_squared_radius_d_object);
     Is_Gabriel is_gabriel Kinit(side_of_bounded_sphere_d_object);
 
-    std::map<Kernel::Point_d, Vertex_handle> points_to_vh;
+    // bimap to retrieve vertex handles from points and vice versa
+    typedef boost::bimap< Kernel::Point_d, Vertex_handle > bimap_points_vh;
+    bimap_points_vh points_to_vh;
     // Start to insert at handle = 0 - default integer value
     Vertex_handle vertex_handle = Vertex_handle();
     // Loop on triangulation vertices list
     for (auto vit = triangulation.vertices_begin(); vit != triangulation.vertices_end(); ++vit) {
-      points_to_vh[vit->point()] = vertex_handle;
+      points_to_vh.insert(bimap_points_vh::value_type(vit->point(), vertex_handle));
       vertex_handle++;
     }
 
@@ -223,10 +228,10 @@ class Alpha_shapes {
       typeVectorPoint pointVector;
       for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
 #ifdef DEBUG_TRACES
-        std::cout << "points_to_vh=" << points_to_vh[(*vit)->point()] << std::endl;
+        std::cout << "points_to_vh=" << points_to_vh.left.at((*vit)->point()) << std::endl;
 #endif  // DEBUG_TRACES
         // Vector of vertex construction for simplex_tree structure
-        vertexVector.push_back(points_to_vh[(*vit)->point()]);
+        vertexVector.push_back(points_to_vh.left.at((*vit)->point()));
         // Vector of points for alpha_shapes filtration value computation
         pointVector.push_back((*vit)->point());
       }
@@ -244,17 +249,99 @@ class Alpha_shapes {
         filtration_max = fmax(filtration_max, alpha_shapes_filtration);
       }
     }
-    
-    // Loop on triangulation finite full cells list
-    for (auto f_simplex : st_.skeleton_simplex_range(st_.dimension() - 1)) {
-      std::cout << "vertex = [" << st_.filtration(f_simplex) << "] ";
+
+    // ### For i : d -> 0
+    for (int decr_dim = st_.dimension(); decr_dim >= 0; decr_dim--) {
+      // ### Foreach Sigma of dim i
+      for (auto f_simplex : st_.skeleton_simplex_range(decr_dim)) {
+        int f_simplex_dim = st_.dimension(f_simplex);
+#ifdef DEBUG_TRACES
+        std::cout << "f_simplex_dim= " << f_simplex_dim << " - decr_dim= " << decr_dim << std::endl;
+#endif  // DEBUG_TRACES
+        if (decr_dim == f_simplex_dim) {
+          typeVectorPoint pointVector;
+#ifdef DEBUG_TRACES
+          std::cout << "vertex ";
+#endif  // DEBUG_TRACES
+          for (auto vertex : st_.simplex_vertex_range(f_simplex)) {
+            pointVector.push_back(points_to_vh.right.at(vertex));
+#ifdef DEBUG_TRACES
+            std::cout << (int) vertex << " ";
+#endif  // DEBUG_TRACES
+          }
+#ifdef DEBUG_TRACES
+          std::cout << std::endl;
+#endif  // DEBUG_TRACES
+          // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma)
+          if (isnan(st_.filtration(f_simplex))) {
+            Filtration_value alpha_shapes_filtration = 0.0;
+            // No need to compute squared_radius on a single point - alpha is 0.0
+            if (f_simplex_dim > 0) {
+              alpha_shapes_filtration = squared_radius(pointVector.begin(), pointVector.end());
+            }
+            st_.assign_filtration(f_simplex, alpha_shapes_filtration);
+#ifdef DEBUG_TRACES
+            std::cout << "From NaN to alpha_shapes_filtration=" << st_.filtration(f_simplex) << std::endl;
+#endif  // DEBUG_TRACES
+          }
+
+          // ### Foreach Tau face of Sigma
+          for (auto f_boundary : st_.boundary_simplex_range(f_simplex)) {
+#ifdef DEBUG_TRACES
+            std::cout << "Sigma ";
+            for (auto vertex : st_.simplex_vertex_range(f_simplex)) {
+              std::cout << (int) vertex << " ";
+            }
+            std::cout << " - Tau ";
+            for (auto vertex : st_.simplex_vertex_range(f_boundary)) {
+              std::cout << (int) vertex << " ";
+            }
+            std::cout << std::endl;
+#endif  // DEBUG_TRACES
+            // insert the Tau points in a vector for is_gabriel function
+            typeVectorPoint pointVector;
+            Vertex_handle vertexForGabriel = Vertex_handle();
+            for (auto vertex : st_.simplex_vertex_range(f_boundary)) {
+              pointVector.push_back(points_to_vh.right.at(vertex));
+            }
+            // Retrieve the Sigma point that is not part of Tau - parameter for is_gabriel function
+            for (auto vertex : st_.simplex_vertex_range(f_simplex)) {
+              if (std::find(pointVector.begin(), pointVector.end(), points_to_vh.right.at(vertex)) == pointVector.end()) {
+                // vertex is not found in Tau
+                vertexForGabriel = vertex;
+                // No need to continue loop
+                break;
+              }
+            }
+            // ### If filt(Tau) is not NaN
+            // ### or Tau is not Gabriel of Sigma
+            if (!isnan(st_.filtration(f_boundary)) ||
+                !is_gabriel(pointVector.begin(), pointVector.end(), points_to_vh.right.at(vertexForGabriel))
+                ) {
+              // ### filt(Tau) = fmin(filt(Tau), filt(Sigma))
+              Filtration_value alpha_shapes_filtration = fmin(st_.filtration(f_boundary), st_.filtration(f_simplex));
+              st_.assign_filtration(f_boundary, alpha_shapes_filtration);
+#ifdef DEBUG_TRACES
+              std::cout << "From Boundary to alpha_shapes_filtration=" << st_.filtration(f_boundary) << std::endl;
+#endif  // DEBUG_TRACES
+            }
+          }
+        }
+      }
+    }
+
+#ifdef DEBUG_TRACES
+    std::cout << "The complex contains " << st_.num_simplices() << " simplices" << std::endl;
+    std::cout << "   - dimension " << st_.dimension() << "   - filtration " << st_.filtration() << std::endl;
+    std::cout << std::endl << std::endl << "Iterator on Simplices in the filtration, with [filtration value]:" << std::endl;
+    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 << (int) vertex << " ";
       }
       std::cout << std::endl;
-      if (st_.filtration(f_simplex) == filtration_unknown)
-        st_.assign_filtration(f_simplex, filtration_max);  // TODO(VR) Compute filtration value from simplex
     }
+#endif  // DEBUG_TRACES
 
 #ifdef DEBUG_TRACES
     std::cout << "filtration_max=" << filtration_max << std::endl;
@@ -286,7 +373,7 @@ class Alpha_shapes {
     return st_.filtration();
   }
 
-  friend std::ostream& operator<<(std::ostream& os, const Alpha_shapes& alpha_shape) {
+  friend std::ostream& operator<<(std::ostream& os, const Alpha_shapes & alpha_shape) {
     // TODO: Program terminated with signal SIGABRT, Aborted - Maybe because of copy constructor
     Gudhi::Simplex_tree<> st = alpha_shape.st_;
     os << st << std::endl;
diff --git a/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h b/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h
index 1413ad89..a4e5e2fe 100644
--- a/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h
+++ b/src/Alpha_shapes/include/gudhi/Alpha_shapes/Delaunay_triangulation_off_io.h
@@ -25,7 +25,7 @@
 #include <string>
 #include <vector>
 #include <fstream>
-#include <iterator>     // std::distance
+#include <map>
 
 #include "gudhi/Off_reader.h"
 
@@ -130,6 +130,7 @@ class Delaunay_triangulation_off_reader {
 template<typename Complex>
 class Delaunay_triangulation_off_writer {
  public:
+  typedef typename Complex::Point Point;
 
   /**
    * name_file : file where the off will be written
@@ -153,44 +154,29 @@ class Delaunay_triangulation_off_writer {
 
       }
 
+      // bimap to retrieve vertex handles from points and vice versa
+      std::map< Point, int > points_to_vh;
+      // Start to insert at default handle value
+      int vertex_handle = int();
+
       // Points list
       for (auto vit = save_complex.vertices_begin(); vit != save_complex.vertices_end(); ++vit) {
         for (auto Coord = vit->point().cartesian_begin(); Coord != vit->point().cartesian_end(); ++Coord) {
           stream << *Coord << " ";
         }
         stream << std::endl;
+        points_to_vh[vit->point()] = vertex_handle;
+        vertex_handle++;
       }
 
-
-      for (auto cit = save_complex.full_cells_begin(); cit != save_complex.full_cells_end(); ++cit) {
+      for (auto cit = save_complex.finite_full_cells_begin(); cit != save_complex.finite_full_cells_end(); ++cit) {
         std::vector<int> vertexVector;
         stream << std::distance(cit->vertices_begin(), cit->vertices_end()) << " ";
         for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
-          // Vector of vertex construction for simplex_tree structure
-          // Vertex handle is distance - 1
-          //int vertexHdl = std::distance(save_complex.vertices_begin(), *vit);
-          // infinite cell is -1 for infinite
-          //vertexVector.push_back(vertexHdl);
-          // Vector of points for alpha_shapes filtration value computation
+          stream << points_to_vh[(*vit)->point()] << " ";
         }
         stream << std::endl;
       }
-
-
-
-
-
-      /*
-              // Finite cells list
-              for (auto cit = save_complex.full_cells_begin(); cit != save_complex.full_cells_end(); ++cit) {
-                stream << std::distance(cit->vertices_begin(), cit->vertices_end()) << " "; // Dimension
-                for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
-                  //auto vertexHdl = *vit;
-                  auto vertexHdl = std::distance(save_complex.vertices_begin(), vit) - 1;
-                  // stream << std::distance(save_complex.vertices_begin(), *(vit)) - 1 << " ";
-                }
-                stream << std::endl;
-              }*/
       stream.close();
     } else {
       std::cerr << "Delaunay_triangulation_off_writer::Delaunay_triangulation_off_writer could not open file " <<
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index b79e3c8f..32fb2f43 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -37,7 +37,6 @@
 #include <vector>
 
 namespace Gudhi {
-  
 /** \defgroup simplex_tree Filtered Complexes
  *
  * A simplicial complex \f$\mathbf{K}\f$
@@ -84,8 +83,8 @@ namespace Gudhi {
  *
  */
 template<typename IndexingTag = linear_indexing_tag,
-    typename FiltrationValue = double, typename SimplexKey = int  // must be a signed integer type
-    , typename VertexHandle = int  // must be a signed integer type, int convertible to it
+typename FiltrationValue = double, typename SimplexKey = int // must be a signed integer type
+, typename VertexHandle = int // must be a signed integer type, int convertible to it
 //         , bool ContiguousVertexHandles = true   //true is Vertex_handles are exactly the set [0;n)
 >
 class Simplex_tree {
@@ -245,6 +244,7 @@ class Simplex_tree {
   Filtration_simplex_range filtration_simplex_range() {
     return filtration_simplex_range(Indexing_tag());
   }
+
   /** \brief Returns a range over the vertices of a simplex.
    *
    * The order in which the vertices are visited is the decreasing order for < on Vertex_handles,
@@ -316,12 +316,14 @@ class Simplex_tree {
   Simplex_key key(Simplex_handle sh) {
     return sh->second.key();
   }
+
   /** \brief Returns the simplex associated to a key.
    *
    * The filtration must be initialized. */
   Simplex_handle simplex(Simplex_key key) {
     return filtration_vect_[key];
   }
+
   /** \brief Returns the filtration value of a simplex.
    *
    * Called on the null_simplex, returns INFINITY. */
@@ -330,12 +332,23 @@ class Simplex_tree {
       return sh->second.filtration();
     } else {
       return INFINITY;
-    }  // filtration(); }
+    }
+  }
+
+  /** \brief Sets the filtration value of a simplex.
+   *
+   * No action if called on the null_simplex*/
+  void assign_filtration(Simplex_handle sh, Filtration_value fv) {
+    if (sh != null_simplex()) {
+      sh->second.assign_filtration(fv);
+    }
   }
+
   /** \brief Returns an upper bound of the filtration values of the simplices. */
   Filtration_value filtration() {
     return threshold_;
   }
+
   /** \brief Returns a Simplex_handle different from all Simplex_handles
    * associated to the simplices in the simplicial complex.
    *
@@ -343,20 +356,24 @@ class Simplex_tree {
   Simplex_handle null_simplex() {
     return Dictionary_it(NULL);
   }
+
   /** \brief Returns a key different for all keys associated to the
    * simplices of the simplicial complex. */
   Simplex_key null_key() {
     return -1;
   }
+
   /** \brief Returns a Vertex_handle different from all Vertex_handles associated
    * to the vertices of the simplicial complex. */
   Vertex_handle null_vertex() {
     return null_vertex_;
   }
+
   /** \brief Returns the number of vertices in the complex. */
   size_t num_vertices() {
     return root_.members_.size();
   }
+
   /** \brief Returns the number of simplices in the complex.
    *
    * Does not count the empty simplex. */
@@ -376,6 +393,7 @@ class Simplex_tree {
     }
     return dim - 1;
   }
+
   /** \brief Returns an upper bound on the dimension of the simplicial complex. */
   int dimension() {
     return dimension_;
@@ -423,7 +441,6 @@ class Simplex_tree {
   Simplex_handle find_vertex(Vertex_handle v) {
     return root_.members_.begin() + v;
   }
-//{ return root_.members_.find(v); }
 
   /** \brief Insert a simplex, represented by a range of Vertex_handles, in the simplicial complex.
    *
@@ -479,7 +496,6 @@ class Simplex_tree {
     return res_insert;
   }
 
-
   /** \brief Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of
    * Vertex_handles, in the simplicial complex.
    *
@@ -487,33 +503,35 @@ class Simplex_tree {
    * @param[in]  filtration the filtration value assigned to the new N-simplex.
   */
   template<class RandomAccessVertexRange>
-  void insert_simplex_and_subfaces(RandomAccessVertexRange& Nsimplex,
-                                         Filtration_value filtration = 0.0) {
+  std::pair<Simplex_handle, bool> insert_simplex_and_subfaces(RandomAccessVertexRange& Nsimplex,
+                                                              Filtration_value filtration = 0.0) {
+    std::pair<Simplex_handle, bool> returned;
     if (Nsimplex.size() > 1) {
       for (unsigned int NIndex = 0; NIndex < Nsimplex.size(); NIndex++) {
         // insert N (N-1)-Simplex
         RandomAccessVertexRange NsimplexMinusOne;
         for (unsigned int NListIter = 0; NListIter < Nsimplex.size() - 1; NListIter++) {
           // (N-1)-Simplex creation
-          NsimplexMinusOne.push_back( Nsimplex[(NIndex + NListIter) % Nsimplex.size()]);
+          NsimplexMinusOne.push_back(Nsimplex[(NIndex + NListIter) % Nsimplex.size()]);
         }
         // (N-1)-Simplex recursive call
-        insert_simplex_and_subfaces(NsimplexMinusOne, filtration);
+        returned = insert_simplex_and_subfaces(NsimplexMinusOne, filtration);
       }
       // N-Simplex insert
-      std::pair<Simplex_handle, bool> returned = insert_simplex(Nsimplex, filtration);
+      returned = insert_simplex(Nsimplex, filtration);
       if (returned.second == true) {
         num_simplices_++;
       }
     } else if (Nsimplex.size() == 1) {
       // 1-Simplex insert - End of recursivity
-      std::pair<Simplex_handle, bool> returned = insert_simplex(Nsimplex, filtration);
+      returned = insert_simplex(Nsimplex, filtration);
       if (returned.second == true) {
         num_simplices_++;
       }
     } else {
       // Nothing to insert - empty vector
     }
+    return returned;
   }
 
   /** \brief Assign a value 'key' to the key of the simplex
@@ -540,17 +558,6 @@ class Simplex_tree {
       return sh->second.children();
   }
 
-// void display_simplex(Simplex_handle sh)
-// {
-//   std::cout << "   " << "[" << filtration(sh) << "] ";
-//   for( auto vertex : simplex_vertex_range(sh) )
-//     { std::cout << vertex << " "; }
-// }
-
-  // void print(Simplex_handle sh, std::ostream& os = std::cout)
-  // {  for(auto v : simplex_vertex_range(sh)) {os << v << " ";}
-  //    os << std::endl; }
-
  public:
   /** Returns a pointer to the root nodes of the simplex tree. */
   Siblings * root() {
@@ -562,10 +569,12 @@ class Simplex_tree {
   void set_filtration(Filtration_value fil) {
     threshold_ = fil;
   }
+
   /** Set a number of simplices for the simplicial complex. */
   void set_num_simplices(const unsigned int& num_simplices) {
     num_simplices_ = num_simplices;
   }
+
   /** Set a dimension for the simplicial complex. */
   void set_dimension(int dimension) {
     dimension_ = dimension;
@@ -623,6 +632,7 @@ class Simplex_tree {
     }
     return ((it1 == rg1.end()) && (it2 != rg2.end()));
   }
+
   /** \brief StrictWeakOrdering, for the simplices, defined by the filtration.
    *
    * It corresponds to the partial order
@@ -631,8 +641,7 @@ class Simplex_tree {
    * to be smaller. The filtration function must be monotonic. */
   struct is_before_in_filtration {
     explicit is_before_in_filtration(Simplex_tree * st)
-        : st_(st) {
-    }
+        : st_(st) { }
 
     bool operator()(const Simplex_handle sh1, const Simplex_handle sh2) const {
       if (st_->filtration(sh1) != st_->filtration(sh2)) {
@@ -708,6 +717,7 @@ class Simplex_tree {
       }
     }
   }
+
   /** \brief Expands the Simplex_tree containing only its one skeleton
    * until dimension max_dim.
    *
@@ -731,6 +741,7 @@ class Simplex_tree {
   }
 
  private:
+
   /** \brief Recursive expansion of the simplex tree.*/
   void siblings_expansion(Siblings * siblings,  // must contain elements
       int k) {
@@ -769,6 +780,7 @@ class Simplex_tree {
       }
     }
   }
+
   /** \brief Intersects Dictionary 1 [begin1;end1) with Dictionary 2 [begin2,end2)
    * and assigns the maximal possible Filtration_value to the Nodes. */
   void intersection(std::vector<std::pair<Vertex_handle, Node> >& intersection,
@@ -800,6 +812,7 @@ class Simplex_tree {
       }
     }
   }
+
   /** Maximum over 3 values.*/
   Filtration_value maximum(Filtration_value a, Filtration_value b,
                            Filtration_value c) {
@@ -824,6 +837,92 @@ class Simplex_tree {
       os << filtration(sh) << " \n";
     }
   }
+  //----------------------------------------------------------------------------------------------
+  //----------------------------------------------------------------------------------------------
+ private:
+
+  /** Recursive search of cofaces
+   */
+  template <class RandomAccessVertexRange>
+  void rec_coface(RandomAccessVertexRange &vertices, Siblings *curr_sib, Dictionary *curr_res, std::vector<Dictionary>& cofaces, unsigned int length, unsigned long codimension) {
+    for (auto sib = curr_sib->members().begin(); sib != curr_sib->members().end() && (vertices.empty() || sib->first <= vertices[vertices.size() - 1]); ++sib) {
+      bool continueRecursion = (codimension == length || curr_res->size() <= codimension); // dimension of actual simplex <= codimension
+      if (vertices.empty()) {
+        if (curr_res->size() >= length && continueRecursion)
+          // If we reached the end of the vertices, and the simplex has more vertices than the given simplex, we found a coface
+        {
+          curr_res->emplace(sib->first, sib->second);
+          bool egalDim = (codimension == length || curr_res->size() == codimension); // dimension of actual simplex == codimension
+          if (egalDim)
+            cofaces.push_back(*curr_res);
+          if (has_children(sib))
+            rec_coface(vertices, sib->second.children(), curr_res, cofaces, length, codimension);
+          curr_res->erase(curr_res->end() - 1);
+        }
+      } else if (continueRecursion) {
+        if (sib->first == vertices[vertices.size() - 1]) // If curr_sib matches with the top vertex
+        {
+          curr_res->emplace(sib->first, sib->second);
+          bool egalDim = (codimension == length || curr_res->size() == codimension); // dimension of actual simplex == codimension
+          if (vertices.size() == 1 && curr_res->size() > length && egalDim)
+            cofaces.push_back(*curr_res);
+          if (has_children(sib)) { // Rec call
+            Vertex_handle tmp = vertices[vertices.size() - 1];
+            vertices.pop_back();
+            rec_coface(vertices, sib->second.children(), curr_res, cofaces, length, codimension);
+            vertices.push_back(tmp);
+          }
+          curr_res->erase(curr_res->end() - 1);
+        } else // (sib->first < vertices[vertices.size()-1])
+        {
+          if (has_children(sib)) {
+            curr_res->emplace(sib->first, sib->second);
+            rec_coface(vertices, sib->second.children(), curr_res, cofaces, length, codimension);
+            curr_res->erase(curr_res->end() - 1);
+          }
+        }
+      }
+    }
+  }
+
+ public:
+
+  /** \brief Compute the cofaces of a n simplex
+   * \param vertices List of vertices which represent the n simplex.
+   * \param codimension The function returns the n+codimension-simplices. If codimension = 0, return all cofaces
+   * \return Vector of Dictionary, empty vector if no cofaces found. 
+   * \warning n+codimension must be lower than Simplex_tree dimension, otherwise an an empty vector is returned.
+   */
+
+  template<class RandomAccessVertexRange>
+  std::vector<Dictionary> coface(const RandomAccessVertexRange &vertices, int codimension) {
+    RandomAccessVertexRange copy = vertices;
+    std::vector<Dictionary> cofaces;
+    std::sort(copy.begin(), copy.end(), std::greater<Vertex_handle>()); // must be sorted in decreasing order
+    if (root_.members().empty()) {
+      std::cerr << "Simplex_tree::coface - empty Simplex_tree" << std::endl;
+      return cofaces; // ----->>
+    }
+    if (vertices.empty()) {
+      std::cerr << "Simplex_tree::coface - empty vertices list" << std::endl;
+      return cofaces; // ----->>
+    }
+    if (codimension < 0) {
+      std::cerr << "Simplex_tree::coface - codimension is empty" << std::endl;
+      return cofaces; // ----->>
+    }
+    if (codimension + vertices.size() >= (unsigned long) dimension_) {
+      std::cerr << "Simplex_tree::coface - codimension + vertices list size cannot be greater than Simplex_tree dimension" << std::endl;
+      return cofaces; // ----->>
+    }
+    std::sort(copy.begin(), copy.end(), std::greater<Vertex_handle>()); // must be sorted in decreasing order
+    Dictionary res;
+    rec_coface(copy, &root_, &res, cofaces, vertices.size(), codimension + vertices.size());
+    return cofaces;
+  }
+
+  //----------------------------------------------------------------------------------------------
+  //----------------------------------------------------------------------------------------------
 
  private:
   Vertex_handle null_vertex_;
@@ -851,6 +950,7 @@ std::ostream& operator<<(std::ostream & os, Simplex_tree<T1, T2, T3> & st) {
   }
   return os;
 }
+
 template<typename T1, typename T2, typename T3>
 std::istream& operator>>(std::istream & is, Simplex_tree<T1, T2, T3> & st) {
   // assert(st.num_simplices() == 0);
diff --git a/src/common/include/gudhi/Off_reader.h b/src/common/include/gudhi/Off_reader.h
index e29218d8..618d1b4d 100644
--- a/src/common/include/gudhi/Off_reader.h
+++ b/src/common/include/gudhi/Off_reader.h
@@ -7,7 +7,7 @@
  *
  *    Author(s):       David Salinas
  *
- *    Copyright (C) 2014  INRIA Sophia Antipolis-Méditerranée (France)
+ *    Copyright (C) 2014  INRIA Sophia Antipolis-Méditerranée (France)
  *
  *    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
@@ -121,7 +121,7 @@ private:
 
     	if(!goto_next_uncomment_line(line)) return false;
     	std::istringstream iss(line);
-   		if(is_off_file){
+   		if((is_off_file) && (!is_noff_file)) {
 	  		off_info_.dim = 3;
 	  		if(!(iss >> off_info_.num_vertices >> off_info_.num_faces >> off_info_.num_edges)){
 	  			std::cerr << "incorrect number of vertices/faces/edges\n";