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

GUDHI uses the Computational Geometry Algorithms Library, CGAL, to compute e.g. the Delaunay complex and alpha shapes of a given point cloud. As such, point clouds need to be provided as CGAL data type objects. In the two dimensional case of interest in this tutorial the data type of pts reads std::vector<CGAL::Epick_d< CGAL::Dimension_tag<2> >::Point_d>.

Generating alpha shapes

Using the GUDHI library, the following lines of code compute all alpha shapes of the point cloud pts at once. For details on the data type of pts we refer to the end of the previous subsection.

// 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.

Computing persistent homology

Given the alpha shapes of an input point cloud, stored in the object simplex_tree, we can finally compute the persistent homology of the filtration of alpha shapes. The following lines of code first initialize the alpha shape filtration. Then the coefficient field of interest, [math]\mathbb{Z}/2\mathbb{Z}[/math], is set and the persistent homology structure, called pcoh, is initialized. The last line then actually does the persistent homology computation, using an efficient persistent cohomology algorithm as mentioned in Maria et al 2014.

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

As the nomenclature in this example code already suggests, GUDHI is capable of computing persistent homology with coefficients in different fields of interest. Actually, the GUDHI persistent cohomology algorithm is capable of computing persistent homology classes of different fields simultaneously, providing effective means to unveil torsional features in point clouds. Again, for details we refer to .

The entire C++ script

The following code is the entire script described in the previous subsections. Output generation lines of code have been added here.

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