GUDHI Tutorial

In this tutorial we will learn how to employ the C++ library GUDHI (Geometric understanding in higher dimensions) in order to compute the Delaunay complex and alpha shapes of given point cloud data. For a more complete picture of GUDHI we refer to the project homepage and recommend reading the tutorials provided there.

This tutorial describes the numerics leading to results as in What is Topological Data Analysis? - A Primer.

The entire script

/* This script computes via GUDHI alpha shapes and persistent homology of points sampled from a circle with noise added
Arguments to this program are:
(1)    Radius of the circle
(2)    Number of points to sample
(3)    Sigma of Gaussian noise

For further information on the alpha shape construction functions employed, cf. e.g.
http://gudhi.gforge.inria.fr/doc/latest/_alpha_complex_2alpha_complex_persistence_8cpp-example.html

To compile use e.g.:
g++ alpha_shapes.cpp -std=c++11 -lgmp -lCGAL -I /usr/local/include -I /usr/local/include/Eigen/ -lboost_system -I /home/daniel/Documents/2018-09-04-14-25-00_GUDHI_2.3.0/include -o alpha_shapes  
*/

#define M_PI 3.14159265358979323846
#include <stdio.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <random>
#include <limits>              // for numeric limits
#include <boost/program_options.hpp>
#define CGAL_EIGEN3_ENABLED    // auxiliary setting
#include <CGAL/Epick_d.h>
#include <gudhi/Alpha_complex.h>
#include <gudhi/Persistent_cohomology.h>
#include <gudhi/Simplex_tree.h>

using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<2> >;
using Point = Kernel::Point_d;
using Vector_of_points = std::vector<Point>;
using Simplex_tree = Gudhi::Simplex_tree<>;
using Filtration_value = Simplex_tree::Filtration_value;

int main(int argc, char *argv[])
{
  // Get spherical point cloud parameters from input
  double r, sigma;
  long n;
  if (argc==4)
  {
    r = std::atof(argv[1]);
    n = std::atol(argv[2]);
    sigma = std::atof(argv[3]);
  } 
  else 
  {
    r = 1.;
    n = 100;
    sigma = 0.1;
    std::cout << "alpha_shapes.cpp main(): Using default input parameters..." << std::endl;
  }  
  // Generate the circular point cloud with noise
  Vector_of_points pts;
  std::random_device rd;  // obtain a seed for the random number engine
  std::mt19937 gen(rd()); // mersenne_twister_engine seeded with rd()
  std::uniform_real_distribution<> dist_uniform(0.,2.*M_PI);
  std::normal_distribution<double> dist_normal(0.,sigma);
  double phase, x, y;
  std::ofstream out_pts("pts.dat", std::ofstream::out);
  for (long i=0; i<n; i++)
  {
    phase = dist_uniform(gen);
    x = r * cos(phase) + dist_normal(gen);
    y = r * sin(phase) + dist_normal(gen);
    pts.push_back(Point(x,y));
    out_pts << x << "\t" << y << std::endl;
  }
  out_pts.close();
    
  // Set maximum alpha radius to infinity, restoring the Delaunay complex in the limiting case
  double alpha_square_max_value {std::numeric_limits<double>::infinity()};
  
  // Initialize the alpha complex of the point cloud
  Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_points(pts);
  
  // Initialize the corresponding simplex tree to store the complex
  Simplex_tree simplex_tree;
  
  // Create the complex
  char name[1000];
  if (alpha_complex_from_points.create_complex(simplex_tree, alpha_square_max_value))
  {
    std::cout << "Alpha complex is of dimension " << simplex_tree.dimension() << " - " << simplex_tree.num_simplices() << " simplices - " << simplex_tree.num_vertices() << " vertices." << std::endl;
    
    // print complex to file
    for (int dim=0; dim<3; dim++)
    {
      sprintf(name, "cplx_dim_%d.dat", dim);
      std::ofstream out_cplx(name, std::ofstream::out);
      for (auto f_simplex : simplex_tree.filtration_simplex_range()) {
        if (simplex_tree.dimension(f_simplex) == dim)    {
          for (auto vertex : simplex_tree.simplex_vertex_range(f_simplex)) {
            out_cplx << vertex << "\t";
          }
          out_cplx << simplex_tree.filtration(f_simplex) << "\n";
        }
      }
      out_cplx.close();
    }
  }
  
  // Sort simplices in filtration order
  simplex_tree.initialize_filtration();
  
  // Compute the persistence diagram of the complex
  int coeff_field_characteristic = 2;
  Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(simplex_tree);
  Filtration_value min_persistence = 0.;
  pcoh.init_coefficients(coeff_field_characteristic);
  pcoh.compute_persistent_cohomology(min_persistence);
  
  // Print persistence diagram to file
  std::ofstream out_dgm("dgm.dat", std::ofstream::out);
  pcoh.output_diagram(out_dgm);
  out_dgm.close();
  
  return 0;
}