Alpha_shapes filtration value

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

Former-commit-id: 303ed34da0718312e4bab3f50b7196acded94cec

7 files changed, 112 insertions, 141 deletions

diff --git a/src/Alpha_shapes/example/CMakeLists.txt b/src/Alpha_shapes/example/CMakeLists.txt
index dab3861c..582a9322 100644
--- a/src/Alpha_shapes/example/CMakeLists.txt
+++ b/src/Alpha_shapes/example/CMakeLists.txt
@@ -1,10 +1,10 @@
 cmake_minimum_required(VERSION 2.6)
 project(GUDHIAlphaShapesExample)
 
-# need CGAL 4.6
-# cmake -DCGAL_DIR=~/workspace/CGAL-4.6-beta1 ../../..
+# need CGAL 4.7
+# cmake -DCGAL_DIR=~/workspace/CGAL-4.7-Ic-41 ../../..
 if(CGAL_FOUND)
-  if (NOT CGAL_VERSION VERSION_LESS 4.6.0)
+  if (NOT CGAL_VERSION VERSION_LESS 4.7.0)
     message(STATUS "CGAL version: ${CGAL_VERSION}.")
 
     include( ${CGAL_USE_FILE} )
@@ -13,12 +13,13 @@ if(CGAL_FOUND)
     if (EIGEN3_FOUND)
       message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.")
       include( ${EIGEN3_USE_FILE} )
+      
+      add_definitions(-DDEBUG_TRACES)
 
       add_executable ( dtoffrw Delaunay_triangulation_off_rw.cpp )
       target_link_libraries(dtoffrw  ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY})
       add_test(dtoffrw_tore3D ${CMAKE_CURRENT_BINARY_DIR}/dtoffrw ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off 3)
 
-      #add_definitions(-DDEBUG_TRACES)
       add_executable ( stfromdt Simplex_tree_from_delaunay_triangulation.cpp )
       target_link_libraries(stfromdt  ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY})
     else()
diff --git a/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp b/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp
index 03f6440d..a31c44ab 100644
--- a/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp
+++ b/src/Alpha_shapes/example/Delaunay_triangulation_off_rw.cpp
@@ -63,7 +63,7 @@ int main(int argc, char **argv) {
 
   T dt(dimension);
   std::string offFileName(argv[1]);
-  Gudhi::alphashapes::Delaunay_triangulation_off_reader<T> off_reader(offFileName, dt, true, true);
+  Gudhi::alphashapes::Delaunay_triangulation_off_reader<T> off_reader(offFileName, dt);
   if (!off_reader.is_valid()) {
     std::cerr << "Unable to read file " << offFileName << std::endl;
     exit(-1); // ----- >>
diff --git a/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp b/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp
index 3a17b519..11079a03 100644
--- a/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp
+++ b/src/Alpha_shapes/example/Simplex_tree_from_delaunay_triangulation.cpp
@@ -48,50 +48,24 @@ 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 << " filename.off" << std::endl;
   exit(-1); // ----- >>
 }
 
 int main(int argc, char **argv) {
-  if (argc != 3) {
+  if (argc != 2) {
     std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl;
     usage(argv[0]);
   }
 
-  int dimension = 0;
-  int returnedScanValue = sscanf(argv[2], "%d", &dimension);
-  if ((returnedScanValue == EOF) || (dimension <= 0)) {
-    std::cerr << "Error: " << argv[2] << " is not correct" << std::endl;
-    usage(argv[0]);
-  }
-
-  // ----------------------------------------------------------------------------
-  //
-  // Init of an alpha-shape from a Delaunay triangulation
-  //
-  // ----------------------------------------------------------------------------
-  T dt(dimension);
   std::string off_file_name(argv[1]);
 
-  Gudhi::alphashapes::Delaunay_triangulation_off_reader<T> off_reader(off_file_name, dt, false, false);
-  if (!off_reader.is_valid()) {
-    std::cerr << "Unable to read file " << off_file_name << std::endl;
-    exit(-1); // ----- >>
-  }
-
-  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;
-
-  Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_dt(dt);
-  //std::cout << alpha_shapes_from_dt << std::endl;
-
   // ----------------------------------------------------------------------------
   //
   // Init of an alpha-shape from a OFF file
   //
   // ----------------------------------------------------------------------------
-  Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_file(off_file_name, dimension);
+  Gudhi::alphashapes::Alpha_shapes alpha_shapes_from_file(off_file_name);
   //std::cout << alpha_shapes_from_file << std::endl;
   
   std::cout << "alpha_shapes_from_file.dimension()=" << alpha_shapes_from_file.dimension() << std::endl;
diff --git a/src/Alpha_shapes/include/gudhi/Alpha_shapes.h b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h
index b8efdb4d..6f290938 100644
--- a/src/Alpha_shapes/include/gudhi/Alpha_shapes.h
+++ b/src/Alpha_shapes/include/gudhi/Alpha_shapes.h
@@ -47,6 +47,8 @@ namespace Gudhi {
 
 namespace alphashapes {
 
+#define Kinit(f) =k.f()
+
 /** \defgroup alpha_shapes Alpha shapes in dimension N
  *
  <DT>Implementations:</DT>
@@ -86,16 +88,22 @@ class Alpha_shapes {
    */
   typedef CGAL::Delaunay_triangulation<Kernel> Delaunay_triangulation;
 
+  typedef typename Kernel::Compute_squared_radius_d Squared_Radius;
+  typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel;
+
+  /** \brief Type required to insert into a simplex_tree (with or without subfaces).*/
+  typedef std::vector<Kernel::Point_d> typeVectorPoint;
+
  private:
   /** \brief Upper bound on the simplex tree of the simplicial complex.*/
   Gudhi::Simplex_tree<> _st;
 
  public:
 
-  Alpha_shapes(std::string off_file_name, int dimension) {
-    Delaunay_triangulation dt(dimension);
-    Gudhi::alphashapes::Delaunay_triangulation_off_reader<Delaunay_triangulation>
-        off_reader(off_file_name, dt, true, true);
+  Alpha_shapes(std::string off_file_name) {
+    // Construct a default Delaunay_triangulation (dim=0) - dim will be set in visitor reader init function
+    Delaunay_triangulation dt(3);
+    Gudhi::alphashapes::Delaunay_triangulation_off_reader<Delaunay_triangulation> off_reader(off_file_name, dt);
     if (!off_reader.is_valid()) {
       std::cerr << "Unable to read file " << off_file_name << std::endl;
       exit(-1); // ----- >>
@@ -120,45 +128,70 @@ class Alpha_shapes {
   template<typename T>
   void init(T triangulation) {
     _st.set_dimension(triangulation.maximal_dimension());
-    _st.set_filtration(0.0);
-    // triangulation points list
-    for (auto vit = triangulation.finite_vertices_begin();
-         vit != triangulation.finite_vertices_end(); ++vit) {
-      typeVectorVertex vertexVector;
-      Vertex_handle vertexHdl = std::distance(triangulation.finite_vertices_begin(), vit);
-      vertexVector.push_back(vertexHdl);
+    Filtration_value filtration_max = 0.0;
 
-      // Insert each point in the simplex tree
-      _st.insert_simplex(vertexVector, 0.0);
+    Kernel k;
+    Squared_Radius squared_radius Kinit(compute_squared_radius_d_object);
+    Is_Gabriel     is_gabriel     Kinit(side_of_bounded_sphere_d_object);
 
+    // triangulation full cells list
+    for (auto cit = triangulation.full_cells_begin(); cit != triangulation.full_cells_end(); ++cit) {
+      typeVectorVertex vertexVector;
+      typeVectorPoint pointVector;
+      for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
+        if (!triangulation.is_infinite(*vit)) {
+          // Vector of vertex construction for simplex_tree structure
+          // Vertex handle is distance - 1
+          Vertex_handle vertexHdl = std::distance(triangulation.vertices_begin(), *vit) - 1;
+          // infinite cell is -1 for infinite
+          vertexVector.push_back(vertexHdl);
+          // Vector of points for alpha_shapes filtration value computation
+          pointVector.push_back((*vit)->point());
 #ifdef DEBUG_TRACES
-      std::cout << "P" << vertexHdl << ":";
-      for (auto Coord = vit->point().cartesian_begin(); Coord != vit->point().cartesian_end(); ++Coord) {
-        std::cout << *Coord << " ";
+          std::cout << "Point ";
+          for (auto Coord = (*vit)->point().cartesian_begin(); Coord != (*vit)->point().cartesian_end(); ++Coord) {
+            std::cout << *Coord << " | ";
+          }
+          std::cout << std::endl;
+#endif  // DEBUG_TRACES
+        }
       }
-      std::cout << std::endl;
+      Filtration_value alpha_shapes_filtration = 0.0;
+
+      if (!triangulation.is_infinite(cit)) {
+        alpha_shapes_filtration = squared_radius(pointVector.begin(), pointVector.end());
+#ifdef DEBUG_TRACES
+        std::cout << "Alpha_shape filtration value = " << alpha_shapes_filtration << std::endl;
 #endif  // DEBUG_TRACES
-    }
-    // triangulation finite full cells list
-    for (auto cit = triangulation.finite_full_cells_begin();
-         cit != triangulation.finite_full_cells_end(); ++cit) {
-      typeVectorVertex vertexVector;
-      for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
-        // Vertex handle is distance - 1
-        Vertex_handle vertexHdl = std::distance(triangulation.vertices_begin(), *vit) - 1;
-        vertexVector.push_back(vertexHdl);
+      } else {
+        Filtration_value tmp_filtration = 0.0;
+        bool is_gab = true;
+        for (auto vit = triangulation.finite_vertices_begin(); vit != triangulation.finite_vertices_end(); ++vit) {
+          if (CGAL::ON_UNBOUNDED_SIDE != is_gabriel(pointVector.begin(), pointVector.end(), vit->point())) {
+            is_gab = false;
+            // TODO(VR) : Compute minimum 
+            
+          }
+        }
+        if (true == is_gab) {
+          alpha_shapes_filtration = squared_radius(pointVector.begin(), pointVector.end());
+#ifdef DEBUG_TRACES
+          std::cout << "Alpha_shape filtration value = " << alpha_shapes_filtration << std::endl;
+#endif  // DEBUG_TRACES
+        }
       }
       // Insert each point in the simplex tree
-      _st.insert_simplex_and_subfaces(vertexVector, 0.0);
+      _st.insert_simplex_and_subfaces(vertexVector, alpha_shapes_filtration);
 
 #ifdef DEBUG_TRACES
-      std::cout << "C" << std::distance(triangulation.finite_full_cells_begin(), cit) << ":";
+      std::cout << "C" << std::distance(triangulation.full_cells_begin(), cit) << ":";
       for (auto value : vertexVector) {
         std::cout << value << ' ';
       }
       std::cout << std::endl;
 #endif  // DEBUG_TRACES
     }
+    _st.set_filtration(filtration_max);
   }
 
  public:
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 693b393e..3215c8f6 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
@@ -36,27 +36,35 @@ namespace alphashapes {
  *@brief Off reader visitor with flag that can be passed to Off_reader to read a Delaunay_triangulation_complex.
  */
 template<typename Complex>
-class Delaunay_triangulation_off_flag_visitor_reader {
+class Delaunay_triangulation_off_visitor_reader {
   Complex& complex_;
   typedef typename Complex::Point Point;
 
-  const bool load_only_points_;
-
  public:
-  explicit Delaunay_triangulation_off_flag_visitor_reader(Complex& complex, bool load_only_points = false) :
-      complex_(complex),
-      load_only_points_(load_only_points) { }
+
+  explicit Delaunay_triangulation_off_visitor_reader(Complex& complex) :
+      complex_(complex) { }
 
   void init(int dim, int num_vertices, int num_faces, int num_edges) {
 #ifdef DEBUG_TRACES
-    std::cout << "init" << std::endl;
+    std::cout << "Delaunay_triangulation_off_visitor_reader::init - dim=" << dim << " - num_vertices=" << 
+        num_vertices << " - num_faces=" << num_faces << " - num_edges=" << num_edges << std::endl;
 #endif  // DEBUG_TRACES
+    if (num_faces > 0) {
+      std::cerr << "Delaunay_triangulation_off_visitor_reader::init faces are not taken into account from OFF " <<
+          "file for Delaunay triangulation - faces are computed." << std::endl;
+    }
+    if (num_edges > 0) {
+      std::cerr << "Delaunay_triangulation_off_visitor_reader::init edges are not taken into account from OFF " <<
+          "file for Delaunay triangulation - edges are computed." << std::endl;
+    }
+    //complex_.set_current_dimension(dim);
   }
 
   void point(const std::vector<double>& point) {
 #ifdef DEBUG_TRACES
-    std::cout << "p ";
-    for (auto coordinate: point) {
+    std::cout << "Delaunay_triangulation_off_visitor_reader::point ";
+    for (auto coordinate : point) {
       std::cout << coordinate << " | ";
     }
     std::cout << std::endl;
@@ -65,62 +73,19 @@ class Delaunay_triangulation_off_flag_visitor_reader {
   }
 
   void maximal_face(const std::vector<int>& face) {
-    // For alpha shapes, only points are read
-  }
-
-  void done() {
+    // For Delaunay Triangulation, only points are read
 #ifdef DEBUG_TRACES
-    std::cout << "done" << std::endl;
-#endif  // DEBUG_TRACES
-  }
-};
-
-/**
- *@brief Off reader visitor that can be passed to Off_reader to read a Delaunay_triangulation_complex.
- */
-template<typename Complex>
-class Delaunay_triangulation_off_visitor_reader {
-  Complex& complex_;
-  // typedef typename Complex::Vertex_handle Vertex_handle;
-  // typedef typename Complex::Simplex_handle Simplex_handle;
-  typedef typename Complex::Point Point;
-
-  const bool load_only_points_;
-  std::vector<Point> points_;
-  // std::vector<Simplex_handle> maximal_faces_;
-
- public:
-  explicit Delaunay_triangulation_off_visitor_reader(Complex& complex, bool load_only_points = false) :
-      complex_(complex),
-      load_only_points_(load_only_points) { }
-
-  void init(int dim, int num_vertices, int num_faces, int num_edges) {
-#ifdef DEBUG_TRACES
-    std::cout << "init - " << num_vertices << std::endl;
-#endif  // DEBUG_TRACES
-    // maximal_faces_.reserve(num_faces);
-    points_.reserve(num_vertices);
-  }
-
-  void point(const std::vector<double>& point) {
-#ifdef DEBUG_TRACES
-    std::cout << "p ";
-    for (auto coordinate: point) {
-      std::cout << coordinate << " | ";
+    std::cout << "Delaunay_triangulation_off_visitor_reader::face ";
+    for (auto vertex : face) {
+      std::cout << vertex << " | ";
     }
     std::cout << std::endl;
 #endif  // DEBUG_TRACES
-    points_.emplace_back(Point(point.size(), point.begin(), point.end()));
-  }
-
-  void maximal_face(const std::vector<int>& face) {
-    // For alpha shapes, only points are read
   }
 
   void done() {
-    complex_.insert(points_.begin(), points_.end());
 #ifdef DEBUG_TRACES
-    std::cout << "done" << std::endl;
+    std::cout << "Delaunay_triangulation_off_visitor_reader::done" << std::endl;
 #endif  // DEBUG_TRACES
   }
 };
@@ -131,27 +96,23 @@ class Delaunay_triangulation_off_visitor_reader {
 template<typename Complex>
 class Delaunay_triangulation_off_reader {
  public:
+
   /**
    * name_file : file to read
    * read_complex : complex that will receive the file content
    * read_only_points : specify true if only the points must be read
    */
-  Delaunay_triangulation_off_reader(const std::string & name_file, Complex& read_complex, bool read_only_points = false,
-                                    bool is_flag = false) : valid_(false) {
+  Delaunay_triangulation_off_reader(const std::string & name_file, Complex& read_complex) : valid_(false) {
     std::ifstream stream(name_file);
     if (stream.is_open()) {
-      if (is_flag) {
-        // For alpha shapes, only points are read
-        Delaunay_triangulation_off_flag_visitor_reader<Complex> off_visitor(read_complex, true);
-        Off_reader off_reader(stream);
-        valid_ = off_reader.read(off_visitor);
-      } else {
-        // For alpha shapes, only points are read
-        Delaunay_triangulation_off_visitor_reader<Complex> off_visitor(read_complex, true);
-        Off_reader off_reader(stream);
-        valid_ = off_reader.read(off_visitor);
-      }
+      Delaunay_triangulation_off_visitor_reader<Complex> off_visitor(read_complex);
+      Off_reader off_reader(stream);
+      valid_ = off_reader.read(off_visitor);
+    } else {
+      std::cerr << "Delaunay_triangulation_off_reader::Delaunay_triangulation_off_reader could not open file " <<
+          name_file << std::endl;
     }
+
   }
 
   /**
@@ -160,7 +121,7 @@ class Delaunay_triangulation_off_reader {
   bool is_valid() const {
     return valid_;
   }
-
+  
  private:
   bool valid_;
 };
@@ -168,6 +129,7 @@ class Delaunay_triangulation_off_reader {
 template<typename Complex>
 class Delaunay_triangulation_off_writer {
  public:
+
   /**
    * name_file : file where the off will be written
    * save_complex : complex that be outputted in the file
@@ -191,7 +153,7 @@ class Delaunay_triangulation_off_writer {
 
       // Finite cells list
       for (auto cit = save_complex.finite_full_cells_begin(); cit != save_complex.finite_full_cells_end(); ++cit) {
-        stream << std::distance(cit->vertices_begin(), cit->vertices_end()) << " ";  // Dimension
+        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;
@@ -201,13 +163,14 @@ class Delaunay_triangulation_off_writer {
       }
       stream.close();
     } else {
-      std::cerr << "could not open file " << name_file << std::endl;
+      std::cerr << "Delaunay_triangulation_off_writer::Delaunay_triangulation_off_writer could not open file " <<
+          name_file << std::endl;
     }
   }
 };
 
-}  // namespace alphashapes
+} // namespace alphashapes
 
-}  // namespace Gudhi
+} // namespace Gudhi
 
 #endif  // SRC_ALPHA_SHAPES_INCLUDE_GUDHI_ALPHA_SHAPES_DELAUNAY_TRIANGULATION_OFF_IO_H_
diff --git a/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp b/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp
index a90704b6..d5db3bfa 100644
--- a/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp
+++ b/src/Alpha_shapes/test/Alpha_shapes_unit_test.cpp
@@ -90,7 +90,7 @@ BOOST_AUTO_TEST_CASE( Delaunay_triangulation ) {
   std::string off_file_name("S8_10.off");
   std::cout << "========== OFF FILE NAME = " << off_file_name << " ==========" << std::endl;
   
-  Gudhi::alphashapes::Delaunay_triangulation_off_reader<T> off_reader(off_file_name, dt, true, true);
+  Gudhi::alphashapes::Delaunay_triangulation_off_reader<T> off_reader(off_file_name, dt);
   std::cout << "off_reader.is_valid()=" << off_reader.is_valid() << std::endl;
   BOOST_CHECK(off_reader.is_valid());
 
diff --git a/src/Alpha_shapes/test/CMakeLists.txt b/src/Alpha_shapes/test/CMakeLists.txt
index a48c1a8f..3cf97b71 100644
--- a/src/Alpha_shapes/test/CMakeLists.txt
+++ b/src/Alpha_shapes/test/CMakeLists.txt
@@ -1,10 +1,10 @@
 cmake_minimum_required(VERSION 2.6)
 project(GUDHIAlphaShapesTest)
 
-# need CGAL 4.6
-# cmake -DCGAL_DIR=~/workspace/CGAL-4.6-beta1 ../../..
+# need CGAL 4.7
+# cmake -DCGAL_DIR=~/workspace/CGAL-4.7-Ic-41 ../../..
 if(CGAL_FOUND)
-  if (NOT CGAL_VERSION VERSION_LESS 4.6.0)
+  if (NOT CGAL_VERSION VERSION_LESS 4.7.0)
     message(STATUS "CGAL version: ${CGAL_VERSION}.")
 
     include( ${CGAL_USE_FILE} )