Alpha_complex.h
1 /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2  * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3  * Author(s): Vincent Rouvreau
4  *
5  * Copyright (C) 2015 Inria
6  *
7  * Modification(s):
8  * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL and Eigen3
9  * - YYYY/MM Author: Description of the modification
10  */
11 
12 #ifndef ALPHA_COMPLEX_H_
13 #define ALPHA_COMPLEX_H_
14 
15 #include <gudhi/Debug_utils.h>
16 // to construct Alpha_complex from a OFF file of points
17 #include <gudhi/Points_off_io.h>
18 
19 #include <stdlib.h>
20 #include <math.h> // isnan, fmax
21 
22 #include <CGAL/Delaunay_triangulation.h>
23 #include <CGAL/Epeck_d.h> // For EXACT or SAFE version
24 #include <CGAL/Epick_d.h> // For FAST version
25 #include <CGAL/Spatial_sort_traits_adapter_d.h>
26 #include <CGAL/property_map.h> // for CGAL::Identity_property_map
27 #include <CGAL/version.h> // for CGAL_VERSION_NR
28 #include <CGAL/NT_converter.h>
29 
30 #include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST
31 
32 #include <iostream>
33 #include <vector>
34 #include <string>
35 #include <limits> // NaN
36 #include <map>
37 #include <utility> // std::pair
38 #include <stdexcept>
39 #include <numeric> // for std::iota
40 
41 // Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
42 #if CGAL_VERSION_NR < 1041101000
43 # error Alpha_complex is only available for CGAL >= 4.11
44 #endif
45 
46 #if !EIGEN_VERSION_AT_LEAST(3,1,0)
47 # error Alpha_complex is only available for Eigen3 >= 3.1.0 installed with CGAL
48 #endif
49 
50 namespace Gudhi {
51 
52 namespace alpha_complex {
53 
54 template<typename D> struct Is_Epeck_D { static const bool value = false; };
55 template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool value = true; };
56 
94 template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>>
96  public:
97  // Add an int in TDS to save point index in the structure
98  typedef CGAL::Triangulation_data_structure<typename Kernel::Dimension,
99  CGAL::Triangulation_vertex<Kernel, std::ptrdiff_t>,
100  CGAL::Triangulation_full_cell<Kernel> > TDS;
102  typedef CGAL::Delaunay_triangulation<Kernel, TDS> Delaunay_triangulation;
103 
105  typedef typename Kernel::Point_d Point_d;
108  typedef Kernel Geom_traits;
109 
110  private:
111  typedef typename Kernel::Compute_squared_radius_d Squared_Radius;
112  typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel;
113  typedef typename Kernel::Point_dimension_d Point_Dimension;
114 
115  // Type required to compute squared radius, or side of bounded sphere on a vector of points.
116  typedef typename std::vector<Point_d> Vector_of_CGAL_points;
117 
118  // Vertex_iterator type from CGAL.
119  typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator;
120 
121  // size_type type from CGAL.
122  typedef typename Delaunay_triangulation::size_type size_type;
123 
124  // Structure to switch from simplex tree vertex handle to CGAL vertex iterator.
125  typedef typename std::vector< CGAL_vertex_iterator > Vector_vertex_iterator;
126 
127  private:
130  Vector_vertex_iterator vertex_handle_to_iterator_;
132  Delaunay_triangulation* triangulation_;
134  Kernel kernel_;
135 
136  public:
146  Alpha_complex(const std::string& off_file_name)
147  : triangulation_(nullptr) {
148  Gudhi::Points_off_reader<Point_d> off_reader(off_file_name);
149  if (!off_reader.is_valid()) {
150  std::cerr << "Alpha_complex - Unable to read file " << off_file_name << "\n";
151  exit(-1); // ----- >>
152  }
153 
154  init_from_range(off_reader.get_point_cloud());
155  }
156 
166  template<typename InputPointRange >
167  Alpha_complex(const InputPointRange& points)
168  : triangulation_(nullptr) {
169  init_from_range(points);
170  }
171 
175  delete triangulation_;
176  }
177 
178  // Forbid copy/move constructor/assignment operator
179  Alpha_complex(const Alpha_complex& other) = delete;
180  Alpha_complex& operator= (const Alpha_complex& other) = delete;
181  Alpha_complex (Alpha_complex&& other) = delete;
182  Alpha_complex& operator= (Alpha_complex&& other) = delete;
183 
190  const Point_d& get_point(std::size_t vertex) const {
191  return vertex_handle_to_iterator_.at(vertex)->point();
192  }
193 
194  private:
195  template<typename InputPointRange >
196  void init_from_range(const InputPointRange& points) {
197  #if CGAL_VERSION_NR < 1050000000
198  if (Is_Epeck_D<Kernel>::value)
199  std::cerr << "It is strongly advised to use a CGAL version >= 5.0 with Epeck_d Kernel for performance reasons."
200  << std::endl;
201  #endif
202 
203  auto first = std::begin(points);
204  auto last = std::end(points);
205 
206  if (first != last) {
207  // point_dimension function initialization
208  Point_Dimension point_dimension = kernel_.point_dimension_d_object();
209 
210  // Delaunay triangulation is point dimension.
211  triangulation_ = new Delaunay_triangulation(point_dimension(*first));
212 
213  std::vector<Point_d> point_cloud(first, last);
214 
215  // Creates a vector {0, 1, ..., N-1}
216  std::vector<std::ptrdiff_t> indices(boost::counting_iterator<std::ptrdiff_t>(0),
217  boost::counting_iterator<std::ptrdiff_t>(point_cloud.size()));
218 
219  typedef boost::iterator_property_map<typename std::vector<Point_d>::iterator,
220  CGAL::Identity_property_map<std::ptrdiff_t>> Point_property_map;
221  typedef CGAL::Spatial_sort_traits_adapter_d<Kernel, Point_property_map> Search_traits_d;
222 
223  CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud)));
224 
225  typename Delaunay_triangulation::Full_cell_handle hint;
226  for (auto index : indices) {
227  typename Delaunay_triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint);
228  // Save index value as data to retrieve it after insertion
229  pos->data() = index;
230  hint = pos->full_cell();
231  }
232  // --------------------------------------------------------------------------------------------
233  // structure to retrieve CGAL points from vertex handle - one vertex handle per point.
234  // Needs to be constructed before as vertex handles arrives in no particular order.
235  vertex_handle_to_iterator_.resize(point_cloud.size());
236  // Loop on triangulation vertices list
237  for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) {
238  if (!triangulation_->is_infinite(*vit)) {
239 #ifdef DEBUG_TRACES
240  std::cout << "Vertex insertion - " << vit->data() << " -> " << vit->point() << std::endl;
241 #endif // DEBUG_TRACES
242  vertex_handle_to_iterator_[vit->data()] = vit;
243  }
244  }
245  // --------------------------------------------------------------------------------------------
246  }
247  }
248 
249  public:
268  template <typename SimplicialComplexForAlpha,
271  Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(),
272  bool exact = false) {
273  // From SimplicialComplexForAlpha type required to insert into a simplicial complex (with or without subfaces).
274  typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle;
275  typedef typename SimplicialComplexForAlpha::Simplex_handle Simplex_handle;
276  typedef std::vector<Vertex_handle> Vector_vertex;
277 
278  if (triangulation_ == nullptr) {
279  std::cerr << "Alpha_complex cannot create_complex from a NULL triangulation\n";
280  return false; // ----- >>
281  }
282  if (triangulation_->maximal_dimension() < 1) {
283  std::cerr << "Alpha_complex cannot create_complex from a zero-dimension triangulation\n";
284  return false; // ----- >>
285  }
286  if (complex.num_vertices() > 0) {
287  std::cerr << "Alpha_complex create_complex - complex is not empty\n";
288  return false; // ----- >>
289  }
290 
291  // --------------------------------------------------------------------------------------------
292  // Simplex_tree construction from loop on triangulation finite full cells list
293  if (triangulation_->number_of_vertices() > 0) {
294  for (auto cit = triangulation_->finite_full_cells_begin();
295  cit != triangulation_->finite_full_cells_end();
296  ++cit) {
297  Vector_vertex vertexVector;
298 #ifdef DEBUG_TRACES
299  std::cout << "Simplex_tree insertion ";
300 #endif // DEBUG_TRACES
301  for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
302  if (*vit != nullptr) {
303 #ifdef DEBUG_TRACES
304  std::cout << " " << (*vit)->data();
305 #endif // DEBUG_TRACES
306  // Vector of vertex construction for simplex_tree structure
307  vertexVector.push_back((*vit)->data());
308  }
309  }
310 #ifdef DEBUG_TRACES
311  std::cout << std::endl;
312 #endif // DEBUG_TRACES
313  // Insert each simplex and its subfaces in the simplex tree - filtration is NaN
314  complex.insert_simplex_and_subfaces(vertexVector, std::numeric_limits<double>::quiet_NaN());
315  }
316  }
317  // --------------------------------------------------------------------------------------------
318 
319  // --------------------------------------------------------------------------------------------
320  // Will be re-used many times
321  Vector_of_CGAL_points pointVector;
322  // ### For i : d -> 0
323  for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) {
324  // ### Foreach Sigma of dim i
325  for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) {
326  int f_simplex_dim = complex.dimension(f_simplex);
327  if (decr_dim == f_simplex_dim) {
328  pointVector.clear();
329 #ifdef DEBUG_TRACES
330  std::cout << "Sigma of dim " << decr_dim << " is";
331 #endif // DEBUG_TRACES
332  for (auto vertex : complex.simplex_vertex_range(f_simplex)) {
333  pointVector.push_back(get_point(vertex));
334 #ifdef DEBUG_TRACES
335  std::cout << " " << vertex;
336 #endif // DEBUG_TRACES
337  }
338 #ifdef DEBUG_TRACES
339  std::cout << std::endl;
340 #endif // DEBUG_TRACES
341  // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma)
342  if (std::isnan(complex.filtration(f_simplex))) {
343  Filtration_value alpha_complex_filtration = 0.0;
344  // No need to compute squared_radius on a single point - alpha is 0.0
345  if (f_simplex_dim > 0) {
346  // squared_radius function initialization
347  Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object();
348 
349  CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv;
350  auto sqrad = squared_radius(pointVector.begin(), pointVector.end());
351 #if CGAL_VERSION_NR >= 1050000000
352  if(exact) CGAL::exact(sqrad);
353 #endif
354  alpha_complex_filtration = cv(sqrad);
355  }
356  complex.assign_filtration(f_simplex, alpha_complex_filtration);
357 #ifdef DEBUG_TRACES
358  std::cout << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl;
359 #endif // DEBUG_TRACES
360  }
361  // No need to propagate further, unweighted points all have value 0
362  if (decr_dim > 1)
363  propagate_alpha_filtration(complex, f_simplex);
364  }
365  }
366  }
367  // --------------------------------------------------------------------------------------------
368 
369  // --------------------------------------------------------------------------------------------
370  // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension
372  // Remove all simplices that have a filtration value greater than max_alpha_square
373  complex.prune_above_filtration(max_alpha_square);
374  // --------------------------------------------------------------------------------------------
375  return true;
376  }
377 
378  private:
379  template <typename SimplicialComplexForAlpha, typename Simplex_handle>
380  void propagate_alpha_filtration(SimplicialComplexForAlpha& complex, Simplex_handle f_simplex) {
381  // From SimplicialComplexForAlpha type required to assign filtration values.
382  typedef typename SimplicialComplexForAlpha::Filtration_value Filtration_value;
383 #ifdef DEBUG_TRACES
384  typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle;
385 #endif // DEBUG_TRACES
386 
387  // ### Foreach Tau face of Sigma
388  for (auto f_boundary : complex.boundary_simplex_range(f_simplex)) {
389 #ifdef DEBUG_TRACES
390  std::cout << " | --------------------------------------------------\n";
391  std::cout << " | Tau ";
392  for (auto vertex : complex.simplex_vertex_range(f_boundary)) {
393  std::cout << vertex << " ";
394  }
395  std::cout << "is a face of Sigma\n";
396  std::cout << " | isnan(complex.filtration(Tau)=" << std::isnan(complex.filtration(f_boundary)) << std::endl;
397 #endif // DEBUG_TRACES
398  // ### If filt(Tau) is not NaN
399  if (!std::isnan(complex.filtration(f_boundary))) {
400  // ### filt(Tau) = fmin(filt(Tau), filt(Sigma))
401  Filtration_value alpha_complex_filtration = fmin(complex.filtration(f_boundary),
402  complex.filtration(f_simplex));
403  complex.assign_filtration(f_boundary, alpha_complex_filtration);
404 #ifdef DEBUG_TRACES
405  std::cout << " | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << complex.filtration(f_boundary) << std::endl;
406 #endif // DEBUG_TRACES
407  // ### Else
408  } else {
409  // insert the Tau points in a vector for is_gabriel function
410  Vector_of_CGAL_points pointVector;
411 #ifdef DEBUG_TRACES
412  Vertex_handle vertexForGabriel = Vertex_handle();
413 #endif // DEBUG_TRACES
414  for (auto vertex : complex.simplex_vertex_range(f_boundary)) {
415  pointVector.push_back(get_point(vertex));
416  }
417  // Retrieve the Sigma point that is not part of Tau - parameter for is_gabriel function
418  Point_d point_for_gabriel;
419  for (auto vertex : complex.simplex_vertex_range(f_simplex)) {
420  point_for_gabriel = get_point(vertex);
421  if (std::find(pointVector.begin(), pointVector.end(), point_for_gabriel) == pointVector.end()) {
422 #ifdef DEBUG_TRACES
423  // vertex is not found in Tau
424  vertexForGabriel = vertex;
425 #endif // DEBUG_TRACES
426  // No need to continue loop
427  break;
428  }
429  }
430  // is_gabriel function initialization
431  Is_Gabriel is_gabriel = kernel_.side_of_bounded_sphere_d_object();
432  bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), point_for_gabriel)
433  != CGAL::ON_BOUNDED_SIDE;
434 #ifdef DEBUG_TRACES
435  std::cout << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << vertexForGabriel << std::endl;
436 #endif // DEBUG_TRACES
437  // ### If Tau is not Gabriel of Sigma
438  if (false == is_gab) {
439  // ### filt(Tau) = filt(Sigma)
440  Filtration_value alpha_complex_filtration = complex.filtration(f_simplex);
441  complex.assign_filtration(f_boundary, alpha_complex_filtration);
442 #ifdef DEBUG_TRACES
443  std::cout << " | filt(Tau) = filt(Sigma) = " << complex.filtration(f_boundary) << std::endl;
444 #endif // DEBUG_TRACES
445  }
446  }
447  }
448  }
449 };
450 
451 } // namespace alpha_complex
452 
453 namespace alphacomplex = alpha_complex;
454 
455 } // namespace Gudhi
456 
457 #endif // ALPHA_COMPLEX_H_
OFF file reader implementation in order to read points from an OFF file.
Definition: Points_off_io.h:122
void insert_simplex_and_subfaces(std::vector< Vertex_handle > const &vertex_range, Filtration_value filtration)
Inserts a simplex with vertices from a given simplex (represented by a vector of Vertex_handle) in th...
Simplex_vertex_range simplex_vertex_range(Simplex_handle const &simplex)
Returns a range over vertices of a given simplex.
Definition: SimplicialComplexForAlpha.h:14
const std::vector< Point_d > & get_point_cloud() const
Point cloud getter.
Definition: Points_off_io.h:158
~Alpha_complex()
Alpha_complex destructor deletes the Delaunay triangulation.
Definition: Alpha_complex.h:174
Kernel::Point_d Point_d
A point in Euclidean space.
Definition: Alpha_complex.h:105
Alpha_complex(const InputPointRange &points)
Alpha_complex constructor from a list of points.
Definition: Alpha_complex.h:167
Alpha_complex(const std::string &off_file_name)
Alpha_complex constructor from an OFF file name.
Definition: Alpha_complex.h:146
unspecified Vertex_handle
Definition: SimplicialComplexForAlpha.h:25
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:20
bool is_valid() const
Returns if the OFF file read operation was successful or not.
Definition: Points_off_io.h:150
The concept SimplicialComplexForAlpha describes the requirements for a type to implement a simplicial...
Definition: SimplicialComplexForAlpha.h:21
void prune_above_filtration(Filtration_value filtration)
Kernel Geom_traits
Geometric traits class that provides the geometric types and predicates needed by Delaunay triangulat...
Definition: Alpha_complex.h:108
Alpha complex data structure.
Definition: Alpha_complex.h:95
Skeleton_simplex_range skeleton_simplex_range
Returns a range over the simplices of the skeleton of the simplicial complex, for a given dimension...
Definition: SimplicialComplexForAlpha.h:70
CGAL::Delaunay_triangulation< Kernel, TDS > Delaunay_triangulation
A Delaunay triangulation of a set of points in .
Definition: Alpha_complex.h:102
int assign_filtration(Simplex_handle simplex, Filtration_value filtration)
unspecified Simplex_handle
Definition: SimplicialComplexForAlpha.h:23
bool create_complex(SimplicialComplexForAlpha &complex, Filtration_value max_alpha_square=std::numeric_limits< Filtration_value >::infinity(), bool exact=false)
Inserts all Delaunay triangulation into the simplicial complex. It also computes the filtration value...
Definition: Alpha_complex.h:270
Boundary_simplex_range boundary_simplex_range(Simplex_handle const &simplex)
Returns a range over boundaries of a given simplex.
unspecified Filtration_value
Definition: SimplicialComplexForAlpha.h:27
const Point_d & get_point(std::size_t vertex) const
get_point returns the point corresponding to the vertex given as parameter.
Definition: Alpha_complex.h:190
GUDHI  Version 3.1.1  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : MIT Generated on Fri Feb 7 2020 16:35:36 for GUDHI by Doxygen 1.8.13