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/Alpha_complex/Alpha_kernel_d.h>
16 #include <gudhi/Debug_utils.h>
17 // to construct Alpha_complex from a OFF file of points
18 #include <gudhi/Points_off_io.h>
19 
20 #include <stdlib.h>
21 #include <math.h> // isnan, fmax
22 #include <memory> // for std::unique_ptr
23 #include <cstddef> // for std::size_t
24 
25 #include <CGAL/Delaunay_triangulation.h>
26 #include <CGAL/Regular_triangulation.h> // aka. Weighted Delaunay triangulation
27 #include <CGAL/Epeck_d.h> // For EXACT or SAFE version
28 #include <CGAL/Epick_d.h> // For FAST version
29 #include <CGAL/Spatial_sort_traits_adapter_d.h>
30 #include <CGAL/property_map.h> // for CGAL::Identity_property_map
31 #include <CGAL/version.h> // for CGAL_VERSION_NR
32 #include <CGAL/NT_converter.h>
33 
34 #include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST
35 
36 #include <boost/range/size.hpp>
37 #include <boost/range/combine.hpp>
38 #include <boost/range/adaptor/transformed.hpp>
39 
40 #include <iostream>
41 #include <vector>
42 #include <string>
43 #include <limits> // NaN
44 #include <map>
45 #include <utility> // std::pair
46 #include <stdexcept>
47 #include <numeric> // for std::iota
48 
49 // Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
50 #if CGAL_VERSION_NR < 1041101000
51 # error Alpha_complex is only available for CGAL >= 4.11
52 #endif
53 
54 #if !EIGEN_VERSION_AT_LEAST(3,1,0)
55 # error Alpha_complex is only available for Eigen3 >= 3.1.0 installed with CGAL
56 #endif
57 
58 namespace Gudhi {
59 
60 namespace alpha_complex {
61 
62 template<typename D> struct Is_Epeck_D { static const bool value = false; };
63 template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool value = true; };
64 
102 template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>, bool Weighted = false>
104  public:
106  using Geom_traits = std::conditional_t<Weighted, CGAL::Regular_triangulation_traits_adapter<Kernel>, Kernel>;
107 
108  // Add an int in TDS to save point index in the structure
109  using TDS = CGAL::Triangulation_data_structure<typename Geom_traits::Dimension,
110  CGAL::Triangulation_vertex<Geom_traits, std::ptrdiff_t>,
111  CGAL::Triangulation_full_cell<Geom_traits> >;
112 
114  using Triangulation = std::conditional_t<Weighted, CGAL::Regular_triangulation<Kernel, TDS>,
115  CGAL::Delaunay_triangulation<Kernel, TDS>>;
116 
119 
120  // Numeric type of coordinates in the kernel
121  using FT = typename A_kernel_d::FT;
122 
126  using Sphere = typename A_kernel_d::Sphere;
127 
129  using Point_d = typename Geom_traits::Point_d;
130 
131  private:
132  // Vertex_iterator type from CGAL.
133  using CGAL_vertex_iterator = typename Triangulation::Vertex_iterator;
134 
135  // size_type type from CGAL.
136  using size_type = typename Triangulation::size_type;
137 
138  // Structure to switch from simplex tree vertex handle to CGAL vertex iterator.
139  using Vector_vertex_iterator = std::vector< CGAL_vertex_iterator >;
140 
141  private:
144  Vector_vertex_iterator vertex_handle_to_iterator_;
146  std::unique_ptr<Triangulation> triangulation_;
148  A_kernel_d kernel_;
149 
151  std::vector<Sphere> cache_, old_cache_;
152 
153  public:
163  Alpha_complex(const std::string& off_file_name) {
164  Gudhi::Points_off_reader<Point_d> off_reader(off_file_name);
165  if (!off_reader.is_valid()) {
166  std::cerr << "Alpha_complex - Unable to read file " << off_file_name << "\n";
167  exit(-1); // ----- >>
168  }
169 
170  init_from_range(off_reader.get_point_cloud());
171  }
172 
183  template<typename InputPointRange >
184  Alpha_complex(const InputPointRange& points) {
185  init_from_range(points);
186  }
187 
200  template <typename InputPointRange, typename WeightRange>
201  Alpha_complex(const InputPointRange& points, WeightRange weights) {
202  static_assert(Weighted, "This constructor is not available for non-weighted versions of Alpha_complex");
203  // FIXME: this test is only valid if we have a forward range
204  GUDHI_CHECK(boost::size(weights) == boost::size(points),
205  std::invalid_argument("Points number in range different from weights range number"));
206  auto weighted_points = boost::range::combine(points, weights)
207  | boost::adaptors::transformed([](auto const&t){return Point_d(boost::get<0>(t), boost::get<1>(t));});
208  init_from_range(weighted_points);
209  }
210 
211  // Forbid copy/move constructor/assignment operator
212  Alpha_complex(const Alpha_complex& other) = delete;
213  Alpha_complex& operator= (const Alpha_complex& other) = delete;
214  Alpha_complex (Alpha_complex&& other) = delete;
215  Alpha_complex& operator= (Alpha_complex&& other) = delete;
216 
219  std::size_t num_vertices() const {
220  if (triangulation_ == nullptr)
221  return 0;
222  else
223  return triangulation_->number_of_vertices();
224  }
225 
232  const Point_d& get_point(std::size_t vertex) const {
233  return vertex_handle_to_iterator_.at(vertex)->point();
234  }
235 
236  private:
237  template<typename InputPointRange >
238  void init_from_range(const InputPointRange& points) {
239  #if CGAL_VERSION_NR < 1050000000
240  if (Is_Epeck_D<Kernel>::value)
241  std::cerr << "It is strongly advised to use a CGAL version >= 5.0 with Epeck_d Kernel for performance reasons."
242  << std::endl;
243  #endif
244 
245 #if CGAL_VERSION_NR < 1050101000
246  // Make compilation fail if weighted and CGAL < 5.1
247  static_assert(!Weighted, "Weighted Alpha_complex is only available for CGAL >= 5.1");
248 #endif
249 
250  auto first = std::begin(points);
251  auto last = std::end(points);
252 
253  if (first != last) {
254  // Delaunay triangulation init with point dimension.
255  triangulation_ = std::make_unique<Triangulation>(kernel_.get_dimension(*first));
256 
257  std::vector<Point_d> point_cloud(first, last);
258 
259  // Creates a vector {0, 1, ..., N-1}
260  std::vector<std::ptrdiff_t> indices(boost::counting_iterator<std::ptrdiff_t>(0),
261  boost::counting_iterator<std::ptrdiff_t>(point_cloud.size()));
262 
263  using Point_property_map = boost::iterator_property_map<typename std::vector<Point_d>::iterator,
264  CGAL::Identity_property_map<std::ptrdiff_t>>;
265  using Search_traits_d = CGAL::Spatial_sort_traits_adapter_d<Geom_traits, Point_property_map>;
266 
267  CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud)));
268 
269  typename Triangulation::Full_cell_handle hint;
270  for (auto index : indices) {
271  typename Triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint);
272  if (pos != nullptr) {
273  // Save index value as data to retrieve it after insertion
274  pos->data() = index;
275  hint = pos->full_cell();
276  }
277  }
278  // --------------------------------------------------------------------------------------------
279  // structure to retrieve CGAL points from vertex handle - one vertex handle per point.
280  // Needs to be constructed before as vertex handles arrives in no particular order.
281  vertex_handle_to_iterator_.resize(point_cloud.size());
282  // Loop on triangulation vertices list
283  for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) {
284  if (!triangulation_->is_infinite(*vit)) {
285 #ifdef DEBUG_TRACES
286  std::clog << "Vertex insertion - " << vit->data() << " -> " << vit->point() << std::endl;
287 #endif // DEBUG_TRACES
288  vertex_handle_to_iterator_[vit->data()] = vit;
289  }
290  }
291  // --------------------------------------------------------------------------------------------
292  }
293  }
294 
301  const Point_d& get_point_(std::size_t vertex) const {
302  return vertex_handle_to_iterator_[vertex]->point();
303  }
304 
306  template<class SimplicialComplexForAlpha>
307  auto& get_cache(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) {
308  auto k = cplx.key(s);
309  if(k==cplx.null_key()){
310  k = cache_.size();
311  cplx.assign_key(s, k);
312  // Using a transform_range is slower, currently.
313  thread_local std::vector<Point_d> v;
314  v.clear();
315  for (auto vertex : cplx.simplex_vertex_range(s))
316  v.push_back(get_point_(vertex));
317  cache_.emplace_back(kernel_.get_sphere(v.cbegin(), v.cend()));
318  }
319  return cache_[k];
320  }
321 
323  template<class SimplicialComplexForAlpha>
324  auto radius(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) {
325  auto k = cplx.key(s);
326  if(k!=cplx.null_key())
327  return kernel_.get_squared_radius(old_cache_[k]);
328  // Using a transform_range is slower, currently.
329  thread_local std::vector<Point_d> v;
330  v.clear();
331  for (auto vertex : cplx.simplex_vertex_range(s))
332  v.push_back(get_point_(vertex));
333  return kernel_.get_squared_radius(v.cbegin(), v.cend());
334  }
335 
336  public:
360  template <typename SimplicialComplexForAlpha,
363  Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(),
364  bool exact = false,
365  bool default_filtration_value = false) {
366  // From SimplicialComplexForAlpha type required to insert into a simplicial complex (with or without subfaces).
368  using Simplex_handle = typename SimplicialComplexForAlpha::Simplex_handle;
369  using Vector_vertex = std::vector<Vertex_handle>;
370 
371  if (triangulation_ == nullptr) {
372  std::cerr << "Alpha_complex cannot create_complex from a NULL triangulation\n";
373  return false; // ----- >>
374  }
375  if (triangulation_->maximal_dimension() < 1) {
376  std::cerr << "Alpha_complex cannot create_complex from a zero-dimension triangulation\n";
377  return false; // ----- >>
378  }
379  if (complex.num_vertices() > 0) {
380  std::cerr << "Alpha_complex create_complex - complex is not empty\n";
381  return false; // ----- >>
382  }
383 
384  // --------------------------------------------------------------------------------------------
385  // Simplex_tree construction from loop on triangulation finite full cells list
386  if (num_vertices() > 0) {
387  for (auto cit = triangulation_->finite_full_cells_begin();
388  cit != triangulation_->finite_full_cells_end();
389  ++cit) {
390  Vector_vertex vertexVector;
391 #ifdef DEBUG_TRACES
392  std::clog << "Simplex_tree insertion ";
393 #endif // DEBUG_TRACES
394  for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
395  if (*vit != nullptr) {
396 #ifdef DEBUG_TRACES
397  std::clog << " " << (*vit)->data();
398 #endif // DEBUG_TRACES
399  // Vector of vertex construction for simplex_tree structure
400  vertexVector.push_back((*vit)->data());
401  }
402  }
403 #ifdef DEBUG_TRACES
404  std::clog << std::endl;
405 #endif // DEBUG_TRACES
406  // Insert each simplex and its subfaces in the simplex tree - filtration is NaN
407  complex.insert_simplex_and_subfaces(vertexVector, std::numeric_limits<double>::quiet_NaN());
408  }
409  }
410  // --------------------------------------------------------------------------------------------
411 
412  if (!default_filtration_value) {
413  CGAL::NT_converter<FT, Filtration_value> cgal_converter;
414  // --------------------------------------------------------------------------------------------
415  // ### For i : d -> 0
416  for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) {
417  // ### Foreach Sigma of dim i
418  for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) {
419  int f_simplex_dim = complex.dimension(f_simplex);
420  if (decr_dim == f_simplex_dim) {
421  // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma)
422  if (std::isnan(complex.filtration(f_simplex))) {
423  Filtration_value alpha_complex_filtration = 0.0;
424  // No need to compute squared_radius on a non-weighted single point - alpha is 0.0
425  if (Weighted || f_simplex_dim > 0) {
426  auto const& sqrad = radius(complex, f_simplex);
427 #if CGAL_VERSION_NR >= 1050000000
428  if(exact) CGAL::exact(sqrad);
429 #endif
430  alpha_complex_filtration = cgal_converter(sqrad);
431  }
432  complex.assign_filtration(f_simplex, alpha_complex_filtration);
433 #ifdef DEBUG_TRACES
434  std::clog << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl;
435 #endif // DEBUG_TRACES
436  }
437  // No need to propagate further, unweighted points all have value 0
438  if (decr_dim > !Weighted)
439  propagate_alpha_filtration(complex, f_simplex);
440  }
441  }
442  old_cache_ = std::move(cache_);
443  cache_.clear();
444  }
445  // --------------------------------------------------------------------------------------------
446 
447  // --------------------------------------------------------------------------------------------
448  if (!exact)
449  // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension
450  // Only in not exact version, cf. https://github.com/GUDHI/gudhi-devel/issues/57
452  // Remove all simplices that have a filtration value greater than max_alpha_square
453  complex.prune_above_filtration(max_alpha_square);
454  // --------------------------------------------------------------------------------------------
455  }
456  return true;
457  }
458 
459  private:
460  template <typename SimplicialComplexForAlpha, typename Simplex_handle>
461  void propagate_alpha_filtration(SimplicialComplexForAlpha& complex, Simplex_handle f_simplex) {
462  // From SimplicialComplexForAlpha type required to assign filtration values.
464 
465  // ### Foreach Tau face of Sigma
466  for (auto face_opposite_vertex : complex.boundary_opposite_vertex_simplex_range(f_simplex)) {
467  auto f_boundary = face_opposite_vertex.first;
468 #ifdef DEBUG_TRACES
469  std::clog << " | --------------------------------------------------\n";
470  std::clog << " | Tau ";
471  for (auto vertex : complex.simplex_vertex_range(f_boundary)) {
472  std::clog << vertex << " ";
473  }
474  std::clog << "is a face of Sigma\n";
475  std::clog << " | isnan(complex.filtration(Tau)=" << std::isnan(complex.filtration(f_boundary)) << std::endl;
476 #endif // DEBUG_TRACES
477  // ### If filt(Tau) is not NaN
478  if (!std::isnan(complex.filtration(f_boundary))) {
479  // ### filt(Tau) = fmin(filt(Tau), filt(Sigma))
480  Filtration_value alpha_complex_filtration = fmin(complex.filtration(f_boundary),
481  complex.filtration(f_simplex));
482  complex.assign_filtration(f_boundary, alpha_complex_filtration);
483 #ifdef DEBUG_TRACES
484  std::clog << " | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << complex.filtration(f_boundary) << std::endl;
485 #endif // DEBUG_TRACES
486  // ### Else
487  } else {
488  auto const& cache=get_cache(complex, f_boundary);
489  bool is_gab = kernel_.is_gabriel(cache, get_point_(face_opposite_vertex.second));
490 #ifdef DEBUG_TRACES
491  std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << face_opposite_vertex.second << std::endl;
492 #endif // DEBUG_TRACES
493  // ### If Tau is not Gabriel of Sigma
494  if (false == is_gab) {
495  // ### filt(Tau) = filt(Sigma)
496  Filtration_value alpha_complex_filtration = complex.filtration(f_simplex);
497  complex.assign_filtration(f_boundary, alpha_complex_filtration);
498 #ifdef DEBUG_TRACES
499  std::clog << " | filt(Tau) = filt(Sigma) = " << complex.filtration(f_boundary) << std::endl;
500 #endif // DEBUG_TRACES
501  }
502  }
503  }
504  }
505 };
506 
507 } // namespace alpha_complex
508 
509 namespace alphacomplex = alpha_complex;
510 
511 } // namespace Gudhi
512 
513 #endif // ALPHA_COMPLEX_H_
OFF file reader implementation in order to read points from an OFF file.
Definition: Points_off_io.h:122
const std::vector< Point_d > & get_point_cloud() const
Point cloud getter.
Definition: Points_off_io.h:158
bool is_valid() const
Returns if the OFF file read operation was successful or not.
Definition: Points_off_io.h:150
Alpha complex data structure.
Definition: Alpha_complex.h:103
bool create_complex(SimplicialComplexForAlpha &complex, Filtration_value max_alpha_square=std::numeric_limits< Filtration_value >::infinity(), bool exact=false, bool default_filtration_value=false)
Inserts all Delaunay triangulation into the simplicial complex. It also computes the filtration value...
Definition: Alpha_complex.h:362
std::conditional_t< Weighted, CGAL::Regular_triangulation< Kernel, TDS >, CGAL::Delaunay_triangulation< Kernel, TDS > > Triangulation
A (Weighted or not) Delaunay triangulation of a set of points in .
Definition: Alpha_complex.h:115
typename A_kernel_d::Sphere Sphere
Sphere is a std::pair<Kernel::Point_d, Kernel::FT> (aka. circurmcenter and squared radius)....
Definition: Alpha_complex.h:126
Alpha_complex(const InputPointRange &points, WeightRange weights)
Alpha_complex constructor from a list of points and weights.
Definition: Alpha_complex.h:201
std::conditional_t< Weighted, CGAL::Regular_triangulation_traits_adapter< Kernel >, Kernel > Geom_traits
Geometric traits class that provides the geometric types and predicates needed by the triangulations.
Definition: Alpha_complex.h:106
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:232
Alpha_complex(const std::string &off_file_name)
Alpha_complex constructor from an OFF file name.
Definition: Alpha_complex.h:163
typename Geom_traits::Point_d Point_d
A point, or a weighted point in Euclidean space.
Definition: Alpha_complex.h:129
std::size_t num_vertices() const
Returns the number of finite vertices in the triangulation.
Definition: Alpha_complex.h:219
Alpha_complex(const InputPointRange &points)
Alpha_complex constructor from a list of points.
Definition: Alpha_complex.h:184
Alpha complex kernel container.
Definition: Alpha_kernel_d.h:42
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:20
The concept SimplicialComplexForAlpha describes the requirements for a type to implement a simplicial...
Definition: SimplicialComplexForAlpha.h:21
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
int assign_filtration(Simplex_handle simplex, Filtration_value filtration)
void prune_above_filtration(Filtration_value filtration)
Simplex_vertex_range simplex_vertex_range(Simplex_handle const &simplex)
Returns a range over vertices of a given simplex.
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...
unspecified Simplex_handle
Definition: SimplicialComplexForAlpha.h:23
unspecified Vertex_handle
Definition: SimplicialComplexForAlpha.h:25
unspecified Filtration_value
Definition: SimplicialComplexForAlpha.h:27
Handle type for the vertices of a cell complex.
Definition: VertexHandle.h:15