delaunay off reader/writer fix.

alpha complex algo seems ok. tests are Nok.


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

Former-commit-id: 37cea3b33d783a77a23fa3a54a1bee106508d981

3 files changed, 115 insertions, 68 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 " <<