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.

Creating the example point cloud

The basis of any topological data analysis routine is a point cloud of data. For instance, the following lines of code generate [math]n[/math] points lying on a circle of radius [math]r[/math] with Gaussian noise of width [math]\sigma[/math] added on top:

  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;
  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));
  }

Generating alpha shapes

Using the GUDHI library, the following lines of code compute all alpha shapes of the point cloud pts at once.

  // 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
  alpha_complex_from_points.create_complex(simplex_tree, alpha_square_max_value)

A simplex tree is a particularly efficient and memory-saving data format to store abstract simplicial complexes, providing cheap algorithms to compute for example faces and cofaces of a given simplex. For more details we refer to the original publication by Boissonat and Maria 2012.

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;
}