Modification for Porquerolles Workshop

Bar_code_persistence and alpha_complex_persistence


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

Former-commit-id: 63d2d21d5ff7c9a94a67eafc326bf8a561bf8166

5 files changed, 207 insertions, 2 deletions

diff --git a/src/GudhUI/CMakeLists.txt b/src/GudhUI/CMakeLists.txt
index 71f4fd1a..5c2afdd0 100644
--- a/src/GudhUI/CMakeLists.txt
+++ b/src/GudhUI/CMakeLists.txt
@@ -78,6 +78,46 @@ if ( CGAL_FOUND AND QT4_FOUND AND OPENGL_FOUND AND QGLVIEWER_FOUND )
   target_link_libraries( GudhUI ${QT_LIBRARIES} ${QGLVIEWER_LIBRARIES} )
   target_link_libraries( GudhUI ${OPENGL_gl_LIBRARY} ${OPENGL_glu_LIBRARY} )	
 
+###############################################################################
+  if (NOT CGAL_VERSION VERSION_LESS 4.7.0)
+    message(STATUS "CGAL version: ${CGAL_VERSION}.")
+
+    include( ${CGAL_USE_FILE} )
+     # In CMakeLists.txt, when include(${CGAL_USE_FILE}), CXX_FLAGS are overwritten.
+    # cf. http://doc.cgal.org/latest/Manual/installation.html#title40
+    # A workaround is to add "-std=c++11" again.
+    # A fix would be to use https://cmake.org/cmake/help/v3.1/prop_gbl/CMAKE_CXX_KNOWN_FEATURES.html 
+    # or even better https://cmake.org/cmake/help/v3.1/variable/CMAKE_CXX_STANDARD.html
+    # but it implies to use cmake version 3.1 at least.
+    if(NOT MSVC)
+      include(CheckCXXCompilerFlag)
+      CHECK_CXX_COMPILER_FLAG(-std=c++11 COMPILER_SUPPORTS_CXX11)
+      if(COMPILER_SUPPORTS_CXX11)
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+      endif()
+    endif()
+    # - End of workaround
+
+    find_package(Eigen3 3.1.0)
+    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_executable (acp alpha_complex_persistence.cpp)
+      target_link_libraries(acp ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY} ${QT_LIBRARIES})
+
+    else()
+      message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.")
+    endif()
+  else()
+    message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile Alpha shapes feature. Version 4.6.0 is required.")
+  endif ()
+###############################################################################
+
 else()
   message(STATUS "NOTICE: GudhUI requires CGAL, the QGLViewer, OpenGL and Qt4, and will not be compiled.")
 endif()
diff --git a/src/GudhUI/alpha_complex_persistence.cpp b/src/GudhUI/alpha_complex_persistence.cpp
new file mode 100644
index 00000000..4f85459a
--- /dev/null
+++ b/src/GudhUI/alpha_complex_persistence.cpp
@@ -0,0 +1,78 @@
+#include <iostream>
+#include <string>
+
+
+#include <QtGui/QApplication>
+
+// to construct a Delaunay_triangulation from a OFF file
+#include <gudhi/Delaunay_triangulation_off_io.h>
+#include <gudhi/Alpha_complex.h>
+#include <gudhi/Persistent_cohomology.h>
+
+#include "utils/Bar_code_persistence.h"
+
+void usage(char * const progName) {
+  std::cerr << "Usage: " << progName << " filename.off " <<  // alpha_square_max_value[double] " <<
+      "coeff_field_characteristic[integer > 0] min_persistence[double >= -1.0]" << std::endl;
+  std::cerr << "       i.e.: " << progName << " ../../data/points/alphacomplexdoc.off 60.0 2 0.02" << std::endl;
+  exit(-1); // ----- >>
+}
+
+int main(int argc, char **argv) {
+  if (argc != 4) {
+    std::cerr << "Error: Number of arguments (" << argc << ") is not correct" << std::endl;
+    usage(argv[0]);
+  }
+  
+  QApplication qtapp(argc, argv);
+
+  std::string off_file_name(argv[1]);
+  // double alpha_square_max_value = atof(argv[2]);
+  double alpha_square_max_value = 1e20;
+  int coeff_field_characteristic = atoi(argv[2]);  // argv[3]
+  double min_persistence = atof(argv[3]);  // argv[4]
+  
+  // ----------------------------------------------------------------------------
+  // Init of an alpha complex from an OFF file
+  // ----------------------------------------------------------------------------
+  typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel;
+  Gudhi::alphacomplex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_name, alpha_square_max_value);
+
+  // ----------------------------------------------------------------------------
+  // Display information about the alpha complex
+  // ----------------------------------------------------------------------------
+  std::cout << "Alpha complex is of dimension " << alpha_complex_from_file.dimension() <<
+      " - " << alpha_complex_from_file.num_simplices() << " simplices - " <<
+      alpha_complex_from_file.num_vertices() << " vertices." << std::endl;
+
+  // Sort the simplices in the order of the filtration
+  alpha_complex_from_file.initialize_filtration();
+
+  std::cout << "Simplex_tree dim: " << alpha_complex_from_file.dimension() << std::endl;
+  // Compute the persistence diagram of the complex
+  Gudhi::persistent_cohomology::Persistent_cohomology< Gudhi::alphacomplex::Alpha_complex<Kernel>,
+      Gudhi::persistent_cohomology::Field_Zp > pcoh(alpha_complex_from_file);
+  
+  std::cout << "coeff_field_characteristic " << coeff_field_characteristic <<
+      " - min_persistence " << min_persistence << std::endl;
+
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(coeff_field_characteristic);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  pcoh.output_diagram();
+  
+  std::vector<std::pair<double, double>> persistence_vector;
+  pcoh.get_persistence(persistence_vector);
+
+  Bar_code_persistence bc_persistence;
+  
+  for (auto persistence : persistence_vector) {
+    bc_persistence.insert(persistence.first, persistence.second);
+  }
+
+  bc_persistence.show();
+
+  return qtapp.exec();
+}
diff --git a/src/GudhUI/model/Model.h b/src/GudhUI/model/Model.h
index d78cbad9..4f8d48ee 100644
--- a/src/GudhUI/model/Model.h
+++ b/src/GudhUI/model/Model.h
@@ -71,7 +71,7 @@ class CGAL_geometric_flag_complex_wrapper {
 
   void maximal_face(std::vector<int> vertices) {
     if (!load_only_points_) {
-      std::cout << "size:" << vertices.size() << std::endl;
+      //std::cout << "size:" << vertices.size() << std::endl;
       for (int i = 0; i < vertices.size(); ++i)
         for (int j = i + 1; j < vertices.size(); ++j)
           complex_.add_edge(Vertex_handle(vertices[i]), Vertex_handle(vertices[j]));
diff --git a/src/GudhUI/utils/Bar_code_persistence.h b/src/GudhUI/utils/Bar_code_persistence.h
new file mode 100644
index 00000000..a1a46ea8
--- /dev/null
+++ b/src/GudhUI/utils/Bar_code_persistence.h
@@ -0,0 +1,84 @@
+#include <math.h>  // isfinite
+
+#include <QtGui/QApplication>
+
+#include <QGraphicsView>
+#include <QGraphicsScene>
+#include <QPointF>
+#include <QVector>
+#include <QGraphicsTextItem>
+
+#include <iostream>
+#include <vector>
+#include <limits>  // NaN, infinity
+#include <utility>  // for pair
+
+class Bar_code_persistence {
+ private:
+  typedef std::vector<std::pair<double, double>> Persistence;
+  Persistence persistence_vector;
+  double min_birth;
+  double max_death;
+
+ public:
+
+  Bar_code_persistence()
+      : min_birth(std::numeric_limits<double>::quiet_NaN()),
+      max_death(std::numeric_limits<double>::quiet_NaN()) { }
+
+  void insert(double birth, double death) {
+    persistence_vector.push_back(std::make_pair(birth, death));
+    if (std::isfinite(birth)) {
+      if ((birth < min_birth) || (std::isnan(min_birth)))
+        min_birth = birth;
+      if ((birth > max_death) || (std::isnan(max_death)))
+        max_death = birth;
+    }
+    if (std::isfinite(death))
+      if ((death > max_death) || (std::isnan(max_death)))
+        max_death = death;
+  }
+
+  void show() {
+    // Create a view, put a scene in it
+    QGraphicsView * view = new QGraphicsView();
+    QGraphicsScene * scene = new QGraphicsScene();
+    view->setScene(scene);
+    double ratio = 600.0 / (max_death - min_birth);
+    //std::cout << "min_birth=" << min_birth << " - max_death=" << max_death << " - ratio=" << ratio << std::endl;
+
+    double height = 0.0, birth = 0.0, death = 0.0;
+    int pers_num = 1;
+    for (auto& persistence : persistence_vector) {
+      height = 5.0 * pers_num;
+      //std::cout << "[" << pers_num << "] birth=" << persistence.first << " - death=" << persistence.second << std::endl;
+      if (std::isfinite(persistence.first))
+        birth = ((persistence.first - min_birth) * ratio) + 50.0;
+      else
+        birth = 0.0;
+
+      if (std::isfinite(persistence.second))
+        death = ((persistence.second - min_birth) * ratio) + 50.0;
+      else
+        death = 700.0;
+
+      scene->addLine(birth, height, death, height, QPen(Qt::blue, 2));
+      pers_num++;
+    }
+    height += 10.0;
+    // scale line
+    scene->addLine(0, height, 700.0, height, QPen(Qt::black, 1));
+    int modulo = 0;
+    for (double scale = 50.0; scale < 700.0; scale += 50.0) {
+      modulo++;
+      // scale small dash
+      scene->addLine(scale, height - 3.0, scale, height + 3.0, QPen(Qt::black, 1));
+      // scale text
+      QString scale_value = QString::number(((scale - 50.0) / ratio) + min_birth);
+      QGraphicsTextItem* dimText = scene->addText(scale_value, QFont("Helvetica", 8));
+      dimText->setPos(scale - (3.0 * scale_value.size()), height + 9.0 * (modulo % 2));
+    }
+    // Show the view
+    view->show();
+  }
+};
diff --git a/src/GudhUI/view/FirstCoordProjector.h b/src/GudhUI/view/FirstCoordProjector.h
index 529d2d42..3ceda3f5 100644
--- a/src/GudhUI/view/FirstCoordProjector.h
+++ b/src/GudhUI/view/FirstCoordProjector.h
@@ -32,8 +32,11 @@ class FirstCoordProjector3D : public Projector3D {
   typedef Projector3D::Point_3 Point_3;
 
   Point_3 operator()(const Point& p) const {
-    assert(p.dimension() >= 3);
+    if (p.dimension() >= 3)
     return Point_3(p.x(), p.y(), p.z());
+    else if  (p.dimension() >= 2)
+    return Point_3(p.x(), p.y(), 0.0);
+
   }
 };