Debug traces in DEBUG mode

Alpha complex from delaunay triangulation in static dimension


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

Former-commit-id: 9dda5a16874ae009c24be39c5a9c9d4124a29063

7 files changed, 81 insertions, 53 deletions

diff --git a/src/Alpha_complex/example/Alpha_complex_from_off.cpp b/src/Alpha_complex/example/Alpha_complex_from_off.cpp
index ce278419..a17155de 100644
--- a/src/Alpha_complex/example/Alpha_complex_from_off.cpp
+++ b/src/Alpha_complex/example/Alpha_complex_from_off.cpp
@@ -22,7 +22,8 @@ int main(int argc, char **argv) {
   // ----------------------------------------------------------------------------
   // Init of an alpha complex from an OFF file
   // ----------------------------------------------------------------------------
-  Gudhi::alphacomplex::Alpha_complex alpha_complex_from_file(off_file_name);
+  typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel;
+  Gudhi::alphacomplex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_name);
   
   // ----------------------------------------------------------------------------
   // Display information about the alpha complex
diff --git a/src/Alpha_complex/example/Alpha_complex_from_points.cpp b/src/Alpha_complex/example/Alpha_complex_from_points.cpp
index fc0e2460..33680a40 100644
--- a/src/Alpha_complex/example/Alpha_complex_from_points.cpp
+++ b/src/Alpha_complex/example/Alpha_complex_from_points.cpp
@@ -49,7 +49,7 @@ int main(int argc, char **argv) {
   // ----------------------------------------------------------------------------
   // Init of an alpha complex from the list of points
   // ----------------------------------------------------------------------------
-  Gudhi::alphacomplex::Alpha_complex alpha_complex_from_points(3, points.size(), points.begin(), points.end());
+  Gudhi::alphacomplex::Alpha_complex<Kernel> alpha_complex_from_points(3, points.size(), points.begin(), points.end());
 
   // ----------------------------------------------------------------------------
   // Display information about the alpha complex
diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt
index 2e64e4db..04fc34af 100644
--- a/src/Alpha_complex/example/CMakeLists.txt
+++ b/src/Alpha_complex/example/CMakeLists.txt
@@ -13,8 +13,11 @@ if(CGAL_FOUND)
     if (EIGEN3_FOUND)
       message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.")
       include( ${EIGEN3_USE_FILE} )
+      if (CMAKE_BUILD_TYPE MATCHES Debug)
+        # For programs to be more verbose
+        add_definitions(-DDEBUG_TRACES)
+      endif()
       
-      #add_definitions(-DDEBUG_TRACES)
       add_executable ( alphaoffreader Alpha_complex_from_off.cpp )
       target_link_libraries(alphaoffreader  ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY})
 
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h
index 16781563..5834e3df 100644
--- a/src/Alpha_complex/include/gudhi/Alpha_complex.h
+++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h
@@ -43,8 +43,8 @@
 #include <iterator>
 #include <vector>
 #include <string>
-#include <limits> // NaN
-#include <iterator>     // std::iterator
+#include <limits>  // NaN
+#include <iterator>  // std::iterator
 
 namespace Gudhi {
 
@@ -61,6 +61,7 @@ namespace alphacomplex {
  * Please refer to \ref alpha_complex for examples.
  *
  */
+template<class Kernel>
 class Alpha_complex : public Simplex_tree<> {
  private:
   // From Simplex_tree
@@ -72,35 +73,43 @@ class Alpha_complex : public Simplex_tree<> {
 
   // From CGAL
   // Kernel for the Delaunay_triangulation. Dimension can be set dynamically.
-  typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel;
+  //typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel;
 
   // Delaunay_triangulation type required to create an alpha-complex.
-  typedef CGAL::Delaunay_triangulation<Kernel> Delaunay_triangulation;
+  typedef typename 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;
 
+  typedef typename Kernel::Point_d Point_d;
+
   // Type required to compute squared radius, or side of bounded sphere on a vector of points.
-  typedef std::vector<Kernel::Point_d> Vector_of_CGAL_points;
+  typedef typename std::vector<Point_d> Vector_of_CGAL_points;
 
   // Vertex_iterator type from CGAL.
-  typedef Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator;
+  typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator;
 
-  // Boost bimap type to switch from CGAL vertex iterator to simplex tree vertex handle and vice versa.
-  typedef boost::bimap< CGAL_vertex_iterator, Vertex_handle > Bimap_vertex;
-  
   // size_type type from CGAL.
-  typedef Delaunay_triangulation::size_type size_type;
+  typedef typename Delaunay_triangulation::size_type size_type;
+
+  // Boost bimap type to switch from CGAL vertex iterator to simplex tree vertex handle and vice versa.
+  //typedef typename boost::bimap< CGAL_vertex_iterator, Vertex_handle > Bimap_vertex;
+  //typedef typename Bimap_vertex::value_type value_type;
+  typedef typename std::map< CGAL_vertex_iterator, Vertex_handle > Map_vertex_iterator_to_handle;
+  typedef typename std::map< Vertex_handle, CGAL_vertex_iterator > Map_vertex_handle_to_iterator;
 
  private:
-  /** \brief Boost bimap to switch from CGAL vertex iterator to simplex tree vertex handle and vice versa.*/
-  Bimap_vertex cgal_simplextree;
+  /** \brief Map to switch from CGAL vertex iterator to simplex tree vertex handle.*/
+  Map_vertex_iterator_to_handle vertex_iterator_to_handle;
+  /** \brief Map to switch from simplex tree vertex handle to CGAL vertex iterator.*/
+  Map_vertex_handle_to_iterator vertex_handle_to_iterator;
   /** \brief Pointer on the CGAL Delaunay triangulation.*/
   Delaunay_triangulation* triangulation;
   /** \brief Kernel for triangulation functions access.*/
   Kernel kernel;
 
  public:
+
   /** \brief Alpha_complex constructor from an OFF file name.
    * Uses the Delaunay_triangulation_off_reader to construct the Delaunay triangulation required to initialize 
    * the Alpha_complex.
@@ -140,9 +149,9 @@ class Alpha_complex : public Simplex_tree<> {
   Alpha_complex(int dimension, size_type size, ForwardIterator firstPoint, ForwardIterator lastPoint)
       : triangulation(nullptr) {
     triangulation = new Delaunay_triangulation(dimension);
-    Delaunay_triangulation::size_type inserted = triangulation->insert<ForwardIterator>(firstPoint, lastPoint);
+    size_type inserted = triangulation->insert(firstPoint, lastPoint);
     if (inserted != size) {
-      std::cerr << "Alpha_complex - insertion failed " << inserted << " != " << size<< std::endl;
+      std::cerr << "Alpha_complex - insertion failed " << inserted << " != " << size << std::endl;
       exit(-1); // ----- >>
     }
     init();
@@ -161,18 +170,20 @@ class Alpha_complex : public Simplex_tree<> {
    * @param[in] vertex Vertex handle of the point to retrieve.
    * @return The founded point.
    */
-  Kernel::Point_d get_point(Vertex_handle vertex) {
-    Kernel::Point_d point;
+  Point_d get_point(Vertex_handle vertex) {
+    Point_d point;
     try {
-      point = cgal_simplextree.right.at(vertex)->point();
-    }
-    catch(...) {
+      if (vertex_handle_to_iterator[vertex] != nullptr) {
+        point = vertex_handle_to_iterator[vertex]->point();
+      }
+    }    catch (...) {
       std::cerr << "Alpha_complex - getPoint not found on vertex " << vertex << std::endl;
     }
     return point;
   }
-  
+
  private:
+
   /** \brief Initialize the Alpha_complex from the Delaunay triangulation.
    *
    * @warning Delaunay triangulation must be already constructed with at least one vertex and dimension must be more 
@@ -197,7 +208,7 @@ class Alpha_complex : public Simplex_tree<> {
       std::cerr << "Alpha_complex init - Cannot init twice" << std::endl;
       return; // ----- >>
     }
-    
+
     set_dimension(triangulation->maximal_dimension());
 
     // --------------------------------------------------------------------------------------------
@@ -206,7 +217,11 @@ class Alpha_complex : public Simplex_tree<> {
     Vertex_handle vertex_handle = Vertex_handle();
     // Loop on triangulation vertices list
     for (CGAL_vertex_iterator vit = triangulation->vertices_begin(); vit != triangulation->vertices_end(); ++vit) {
-      cgal_simplextree.insert(Bimap_vertex::value_type(vit, vertex_handle));
+#ifdef DEBUG_TRACES
+      std::cout << "Vertex insertion - " << vertex_handle << " -> " << vit->point() << std::endl;
+#endif  // DEBUG_TRACES
+      vertex_iterator_to_handle[vit] = vertex_handle;
+      vertex_handle_to_iterator[vertex_handle] = vit;
       vertex_handle++;
     }
     // --------------------------------------------------------------------------------------------
@@ -219,18 +234,20 @@ class Alpha_complex : public Simplex_tree<> {
       std::cout << "Simplex_tree insertion ";
 #endif  // DEBUG_TRACES
       for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
+        if (*vit != nullptr) {
 #ifdef DEBUG_TRACES
-        std::cout << " " << cgal_simplextree.left.at(*vit);
+          std::cout << " " << vertex_iterator_to_handle[*vit];
 #endif  // DEBUG_TRACES
-        // Vector of vertex construction for simplex_tree structure
-        vertexVector.push_back(cgal_simplextree.left.at(*vit));
+          // Vector of vertex construction for simplex_tree structure
+          vertexVector.push_back(vertex_iterator_to_handle[*vit]);
+        }
       }
 #ifdef DEBUG_TRACES
       std::cout << std::endl;
 #endif  // DEBUG_TRACES
       // Insert each simplex and its subfaces in the simplex tree - filtration is NaN
       Simplex_result insert_result = insert_simplex_and_subfaces(vertexVector,
-                                                                     std::numeric_limits<double>::quiet_NaN());
+                                                                 std::numeric_limits<double>::quiet_NaN());
       if (!insert_result.second) {
         std::cerr << "Alpha_complex::init insert_simplex_and_subfaces failed" << std::endl;
       }
@@ -250,7 +267,7 @@ class Alpha_complex : public Simplex_tree<> {
           std::cout << "Sigma of dim " << decr_dim << " is";
 #endif  // DEBUG_TRACES
           for (auto vertex : simplex_vertex_range(f_simplex)) {
-            pointVector.push_back((cgal_simplextree.right.at(vertex))->point());
+            pointVector.push_back(get_point(vertex));
 #ifdef DEBUG_TRACES
             std::cout << " " << vertex;
 #endif  // DEBUG_TRACES
@@ -265,7 +282,7 @@ class Alpha_complex : public Simplex_tree<> {
             if (f_simplex_dim > 0) {
               // squared_radius function initialization
               Squared_Radius squared_radius = kernel.compute_squared_radius_d_object();
-              
+
               alpha_complex_filtration = squared_radius(pointVector.begin(), pointVector.end());
             }
             assign_filtration(f_simplex, alpha_complex_filtration);
@@ -317,12 +334,11 @@ class Alpha_complex : public Simplex_tree<> {
           Vector_of_CGAL_points pointVector;
           Vertex_handle vertexForGabriel = Vertex_handle();
           for (auto vertex : simplex_vertex_range(f_boundary)) {
-            pointVector.push_back((cgal_simplextree.right.at(vertex))->point());
+            pointVector.push_back(get_point(vertex));
           }
           // Retrieve the Sigma point that is not part of Tau - parameter for is_gabriel function
           for (auto vertex : simplex_vertex_range(f_simplex)) {
-            if (std::find(pointVector.begin(), pointVector.end(), (cgal_simplextree.right.at(vertex))->point())
-                == pointVector.end()) {
+            if (std::find(pointVector.begin(), pointVector.end(), get_point(vertex)) == pointVector.end()) {
               // vertex is not found in Tau
               vertexForGabriel = vertex;
               // No need to continue loop
@@ -331,14 +347,13 @@ class Alpha_complex : public Simplex_tree<> {
           }
           // is_gabriel function initialization
           Is_Gabriel is_gabriel = kernel.side_of_bounded_sphere_d_object();
-#ifdef DEBUG_TRACES
-          bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), (cgal_simplextree.right.at(vertexForGabriel))->point())
+          bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), get_point(vertexForGabriel))
               != CGAL::ON_BOUNDED_SIDE;
+#ifdef DEBUG_TRACES
           std::cout << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << vertexForGabriel << std::endl;
 #endif  // DEBUG_TRACES
           // ### If Tau is not Gabriel of Sigma
-          if ((is_gabriel(pointVector.begin(), pointVector.end(), (cgal_simplextree.right.at(vertexForGabriel))->point())
-               == CGAL::ON_BOUNDED_SIDE)) {
+          if (false == is_gab) {
             // ### filt(Tau) = filt(Sigma)
             Filtration_value alpha_complex_filtration = filtration(f_simplex);
             assign_filtration(f_boundary, alpha_complex_filtration);
diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
index b55b5e2e..4202c9e9 100644
--- a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
+++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp
@@ -35,11 +35,9 @@
 #include <limits>
 
 // Use dynamic_dimension_tag for the user to be able to set dimension
-typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel;
-typedef Kernel::Point_d Point;
-typedef std::vector<Point> Vector_of_points;
+typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel_d;
 // The triangulation uses the default instantiation of the TriangulationDataStructure template parameter
-
+/*
 BOOST_AUTO_TEST_CASE(S4_100_OFF_file) {
   // ----------------------------------------------------------------------------
   //
@@ -49,7 +47,7 @@ BOOST_AUTO_TEST_CASE(S4_100_OFF_file) {
   std::string off_file_name("S4_100.off");
   std::cout << "========== OFF FILE NAME = " << off_file_name << " ==========" << std::endl;
 
-  Gudhi::alphacomplex::Alpha_complex alpha_complex_from_file(off_file_name);
+  Gudhi::alphacomplex::Alpha_complex<Kernel_d> alpha_complex_from_file(off_file_name);
 
   const int DIMENSION = 4;
   std::cout << "alpha_complex_from_file.dimension()=" << alpha_complex_from_file.dimension() << std::endl;
@@ -74,7 +72,7 @@ BOOST_AUTO_TEST_CASE(S8_10_OFF_file) {
   std::string off_file_name("S8_10.off");
   std::cout << "========== OFF FILE NAME = " << off_file_name << " ==========" << std::endl;
 
-  Gudhi::alphacomplex::Alpha_complex alpha_complex_from_file(off_file_name);
+  Gudhi::alphacomplex::Alpha_complex<Kernel_d> alpha_complex_from_file(off_file_name);
 
   const int DIMENSION = 8;
   std::cout << "alpha_complex_from_file.dimension()=" << alpha_complex_from_file.dimension() << std::endl;
@@ -88,11 +86,17 @@ BOOST_AUTO_TEST_CASE(S8_10_OFF_file) {
   std::cout << "alpha_complex_from_file.num_simplices()=" << alpha_complex_from_file.num_simplices() << std::endl;
   BOOST_CHECK(alpha_complex_from_file.num_simplices() == NUMBER_OF_SIMPLICES);
 }
-
+*/
 bool are_almost_the_same(float a, float b) {
   return std::fabs(a - b) < std::numeric_limits<float>::epsilon();
 }
 
+// Use dynamic_dimension_tag for the user to be able to set dimension
+typedef CGAL::Epick_d< CGAL::Dimension_tag<4> > Kernel_s;
+typedef Kernel_s::Point_d Point;
+typedef std::vector<Point> Vector_of_points;
+
+
 bool is_point_in_list(Vector_of_points points_list, Point point) {
   for (auto& point_in_list : points_list) {
     if (point_in_list == point) {
@@ -101,6 +105,7 @@ bool is_point_in_list(Vector_of_points points_list, Point point) {
   }
   return false; // point not found
 }
+
 BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) {
 
   // ----------------------------------------------------------------------------
@@ -109,6 +114,7 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) {
   Vector_of_points points;
   std::vector<double> coords;
 
+  points.clear();
   coords.clear();
   coords.push_back(0.0);
   coords.push_back(0.0);
@@ -137,12 +143,12 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) {
   // ----------------------------------------------------------------------------
   // Init of an alpha complex from the list of points
   // ----------------------------------------------------------------------------
-  Gudhi::alphacomplex::Alpha_complex alpha_complex_from_points(3, points.size(), points.begin(), points.end());
+  Gudhi::alphacomplex::Alpha_complex<Kernel_s> alpha_complex_from_points(3, points.size(), points.begin(), points.end());
 
   std::cout << "========== Alpha_complex_from_points ==========" << std::endl;
 
   std::cout << "alpha_complex_from_points.dimension()=" << alpha_complex_from_points.dimension() << std::endl;
-  BOOST_CHECK(alpha_complex_from_points.dimension() == 3);
+  BOOST_CHECK(alpha_complex_from_points.dimension() == 4);
   std::cout << "alpha_complex_from_points.num_simplices()=" << alpha_complex_from_points.num_simplices() << std::endl;
   BOOST_CHECK(alpha_complex_from_points.num_simplices() == 15);
   std::cout << "alpha_complex_from_points.num_vertices()=" << alpha_complex_from_points.num_vertices() << std::endl;
@@ -190,17 +196,14 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) {
 
   Point p5 = alpha_complex_from_points.get_point(5);
   std::cout << "alpha_complex_from_points.get_point(5)=" << p5 << std::endl;
-  BOOST_CHECK(0 == p5.dimension());
   BOOST_CHECK(!is_point_in_list(points, p5));
 
   Point p0 = alpha_complex_from_points.get_point(0);
   std::cout << "alpha_complex_from_points.get_point(0)=" << p0 << std::endl;
-  BOOST_CHECK(0 == p0.dimension());
   BOOST_CHECK(!is_point_in_list(points, p0));
 
   Point p1234 = alpha_complex_from_points.get_point(1234);
   std::cout << "alpha_complex_from_points.get_point(1234)=" << p1234.dimension() << std::endl;
-  BOOST_CHECK(0 == p1234.dimension());
   BOOST_CHECK(!is_point_in_list(points, p1234));
 
 }
diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt
index 79300790..87b05a3a 100644
--- a/src/Alpha_complex/test/CMakeLists.txt
+++ b/src/Alpha_complex/test/CMakeLists.txt
@@ -14,8 +14,11 @@ if(CGAL_FOUND)
       message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.")
       include( ${EIGEN3_USE_FILE} )
       include_directories (BEFORE "../../include")
+      if (CMAKE_BUILD_TYPE MATCHES Debug)
+        # For programs to be more verbose
+        add_definitions(-DDEBUG_TRACES)
+      endif()
 
-      #add_definitions(-DDEBUG_TRACES)
       add_executable ( AlphaComplexUT Alpha_complex_unit_test.cpp )
       target_link_libraries(AlphaComplexUT  ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
       add_test(AlphaComplexUT ${CMAKE_CURRENT_BINARY_DIR}/AlphaComplexUT
diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt
index 9487cce6..3276989d 100644
--- a/src/Persistent_cohomology/example/CMakeLists.txt
+++ b/src/Persistent_cohomology/example/CMakeLists.txt
@@ -28,8 +28,11 @@ if (NOT MSVC)
       target_link_libraries(performance_rips_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES})
    
       if(CGAL_FOUND)
-         # uncomment to display debug traces
-         # add_definitions(-DDEBUG_TRACES)
+        if (CMAKE_BUILD_TYPE MATCHES Debug)
+          # For programs to be more verbose
+          add_definitions(-DDEBUG_TRACES)
+        endif()
+
          add_executable(alpha_shapes_persistence alpha_shapes_persistence.cpp)
          target_link_libraries(alpha_shapes_persistence ${Boost_SYSTEM_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES} ${CGAL_LIBRARY})
          add_test(alpha_shapes_persistence_2_0_5 ${CMAKE_CURRENT_BINARY_DIR}/alpha_shapes_persistence ${CMAKE_SOURCE_DIR}/data/points/bunny_5000 2 0.5)