All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
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 * - 2024/10 Vincent Rouvreau: Add square root filtration values interface
10 * - YYYY/MM Author: Description of the modification
11 */
12
13#ifndef ALPHA_COMPLEX_H_
14#define ALPHA_COMPLEX_H_
15
16#include <gudhi/Alpha_complex/Alpha_kernel_d.h>
17#include <gudhi/Debug_utils.h>
18// to construct Alpha_complex from a OFF file of points
19#include <gudhi/Points_off_io.h>
20
21#include <CGAL/Delaunay_triangulation.h>
22#include <CGAL/Regular_triangulation.h> // aka. Weighted Delaunay triangulation
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 <boost/range/size.hpp>
33#include <boost/range/combine.hpp>
34#include <boost/range/adaptor/transformed.hpp>
35
36#include <iostream>
37#include <vector>
38#include <string>
39#include <limits> // NaN
40#include <map>
41#include <utility> // std::pair
42#include <stdexcept>
43#include <numeric> // for std::iota
44#include <algorithm> // for std::sort
45#include <type_traits> // for std::is_same_v
46#include <cmath> // isnan, fmax
47#include <memory> // for std::unique_ptr
48#include <cstddef> // for std::size_t
49
50// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
51#if CGAL_VERSION_NR < 1041101000
52# error Alpha_complex is only available for CGAL >= 4.11
53#endif
54
55#if !EIGEN_VERSION_AT_LEAST(3,1,0)
56# error Alpha_complex is only available for Eigen3 >= 3.1.0 installed with CGAL
57#endif
58
59namespace Gudhi {
60
61namespace alpha_complex {
62
63template<typename D> struct Is_Epeck_D { static const bool value = false; };
64template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool value = true; };
65
102template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>, bool Weighted = false>
104 private:
105 // Vertex_handle internal type (required by triangulation_ and vertices_).
106 using Internal_vertex_handle = std::ptrdiff_t;
107
108 public:
110 using Geom_traits = std::conditional_t<Weighted, CGAL::Regular_triangulation_traits_adapter<Kernel>, Kernel>;
111
112 // CGAL::Triangulation_ds_full_cell<void, CGAL::TDS_full_cell_mirror_storage_policy> has been enhanced for CGAL >= 6.0
113 // But faster only with static dimensions
114 using Triangulation_full_cell = std::conditional_t<std::is_same_v<typename Kernel::Dimension, CGAL::Dynamic_dimension_tag>,
115 CGAL::Triangulation_ds_full_cell<>,
116 CGAL::Triangulation_ds_full_cell<void, CGAL::TDS_full_cell_mirror_storage_policy>>;
117 // Add an int in TDS to save point index in the structure
118 using TDS = CGAL::Triangulation_data_structure<typename Geom_traits::Dimension,
119 CGAL::Triangulation_vertex<Geom_traits, Internal_vertex_handle>,
120 Triangulation_full_cell >;
121
123 using Triangulation = std::conditional_t<Weighted, CGAL::Regular_triangulation<Kernel, TDS>,
124 CGAL::Delaunay_triangulation<Kernel, TDS>>;
125
128
129 // Numeric type of coordinates in the kernel
130 using FT = typename A_kernel_d::FT;
131
135 using Sphere = typename A_kernel_d::Sphere;
136
138 using Point_d = typename Geom_traits::Point_d;
139
140 private:
141 // Vertex_iterator type from CGAL.
142 using CGAL_vertex_iterator = typename Triangulation::Vertex_iterator;
143
144 // Structure to switch from simplex tree vertex handle to CGAL vertex iterator.
145 using Vector_vertex_iterator = std::vector< CGAL_vertex_iterator >;
146
147 private:
150 Vector_vertex_iterator vertex_handle_to_iterator_;
152 std::unique_ptr<Triangulation> triangulation_;
154 A_kernel_d kernel_;
158 std::vector<Internal_vertex_handle> vertices_;
159
161 std::vector<Sphere> cache_, old_cache_;
162
163 public:
173 Alpha_complex(const std::string& off_file_name) {
174 Gudhi::Points_off_reader<Point_d> off_reader(off_file_name);
175 if (!off_reader.is_valid()) {
176 std::cerr << "Alpha_complex - Unable to read file " << off_file_name << "\n";
177 exit(-1); // ----- >>
178 }
179
180 init_from_range(off_reader.get_point_cloud());
181 }
182
193 template<typename InputPointRange >
194 Alpha_complex(const InputPointRange& points) {
195 init_from_range(points);
196 }
197
210 template <typename InputPointRange, typename WeightRange>
211 Alpha_complex(const InputPointRange& points, WeightRange weights) {
212 static_assert(Weighted, "This constructor is not available for non-weighted versions of Alpha_complex");
213 // FIXME: this test is only valid if we have a forward range
214 GUDHI_CHECK(boost::size(weights) == boost::size(points),
215 std::invalid_argument("Points number in range different from weights range number"));
216 auto weighted_points = boost::range::combine(points, weights)
217 | boost::adaptors::transformed([](auto const&t){return Point_d(boost::get<0>(t), boost::get<1>(t));});
218 init_from_range(weighted_points);
219 }
220
221 // Forbid copy/move constructor/assignment operator
222 Alpha_complex(const Alpha_complex& other) = delete;
223 Alpha_complex& operator= (const Alpha_complex& other) = delete;
224 Alpha_complex (Alpha_complex&& other) = delete;
225 Alpha_complex& operator= (Alpha_complex&& other) = delete;
226
229 std::size_t num_vertices() const {
230 if (triangulation_ == nullptr)
231 return 0;
232 else
233 return triangulation_->number_of_vertices();
234 }
235
242 const Point_d& get_point(std::size_t vertex) const {
243 auto it = vertex_handle_to_iterator_.at(vertex);
244 if (it == nullptr) throw std::out_of_range("This vertex is missing, maybe hidden by a duplicate or another heavier point.");
245 return it->point();
246 }
247
248 private:
249 template<typename InputPointRange >
250 void init_from_range(const InputPointRange& points) {
251 #if CGAL_VERSION_NR < 1050000000
252 if (Is_Epeck_D<Kernel>::value)
253 std::cerr << "It is strongly advised to use a CGAL version >= 5.0 with Epeck_d Kernel for performance reasons."
254 << std::endl;
255 #endif
256
257#if CGAL_VERSION_NR < 1050101000
258 // Make compilation fail if weighted and CGAL < 5.1.0
259 static_assert(!Weighted, "Weighted Alpha_complex is only available for CGAL >= 5.1.0");
260#endif
261
262 auto first = std::begin(points);
263 auto last = std::end(points);
264
265 if (first != last) {
266 // Delaunay triangulation init with point dimension.
267 triangulation_ = std::make_unique<Triangulation>(kernel_.get_dimension(*first));
268
269 std::vector<Point_d> point_cloud(first, last);
270
271 // Creates a vector {0, 1, ..., N-1}
272 std::vector<Internal_vertex_handle> indices(boost::counting_iterator<Internal_vertex_handle>(0),
273 boost::counting_iterator<Internal_vertex_handle>(point_cloud.size()));
274
275 using Point_property_map = boost::iterator_property_map<typename std::vector<Point_d>::iterator,
276 CGAL::Identity_property_map<Internal_vertex_handle>>;
277 using Search_traits_d = CGAL::Spatial_sort_traits_adapter_d<Geom_traits, Point_property_map>;
278
279 CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud)));
280
281 typename Triangulation::Full_cell_handle hint;
282 for (auto index : indices) {
283 typename Triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint);
284 if (pos != nullptr) {
285 // Save index value as data to retrieve it after insertion
286 pos->data() = index;
287 hint = pos->full_cell();
288 }
289 }
290 // --------------------------------------------------------------------------------------------
291 // structure to retrieve CGAL points from vertex handle - one vertex handle per point.
292 // Needs to be constructed before as vertex handles arrives in no particular order.
293 vertex_handle_to_iterator_.resize(point_cloud.size());
294 // List of sorted unique vertices in the triangulation. We take advantage of the existing loop to construct it
295 // Vertices list avoids quadratic complexity with the Simplex_tree. We should not fill it up with Toplex_map e.g.
296 vertices_.reserve(triangulation_->number_of_vertices());
297 // Loop on triangulation vertices list
298 for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) {
299 if (!triangulation_->is_infinite(*vit)) {
300#ifdef DEBUG_TRACES
301 std::clog << "Vertex insertion - " << vit->data() << " -> " << vit->point() << std::endl;
302#endif // DEBUG_TRACES
303 vertex_handle_to_iterator_[vit->data()] = vit;
304 vertices_.push_back(vit->data());
305 }
306 }
307 std::sort(vertices_.begin(), vertices_.end());
308 // --------------------------------------------------------------------------------------------
309 }
310 }
311
318 const Point_d& get_point_(std::size_t vertex) const {
319 return vertex_handle_to_iterator_[vertex]->point();
320 }
321
323 template<class SimplicialComplexForAlpha>
324 auto& get_cache(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) {
325 auto k = cplx.key(s);
326 if(k==cplx.null_key()){
327 k = cache_.size();
328 cplx.assign_key(s, k);
329 // Using a transform_range is slower, currently.
330 thread_local std::vector<Point_d> v;
331 v.clear();
332 for (auto vertex : cplx.simplex_vertex_range(s))
333 v.push_back(get_point_(vertex));
334 cache_.emplace_back(kernel_.get_sphere(v.cbegin(), v.cend()));
335 }
336 return cache_[k];
337 }
338
340 template<class SimplicialComplexForAlpha>
341 auto radius(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) {
342 auto k = cplx.key(s);
343 if(k!=cplx.null_key())
344 return kernel_.get_squared_radius(old_cache_[k]);
345 // Using a transform_range is slower, currently.
346 thread_local std::vector<Point_d> v;
347 v.clear();
348 for (auto vertex : cplx.simplex_vertex_range(s))
349 v.push_back(get_point_(vertex));
350 return kernel_.get_squared_radius(v.cbegin(), v.cend());
351 }
352
353 public:
383 template <bool Output_squared_values = true,
384 typename SimplicialComplexForAlpha,
387 Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(),
388 bool exact = false,
389 bool default_filtration_value = false) {
390 // Filtration_value must be capable to represent the special value "Not-A-Number"
391 static_assert(std::numeric_limits<Filtration_value>::has_quiet_NaN);
392 // Forbids weighted and output square root filtration values - cannot static_assert for python purpose
393 if constexpr(Weighted && !Output_squared_values)
394 throw std::invalid_argument("Weighted Alpha complex cannot output square root filtration values");
395 // To support more general types for Filtration_value
396 using std::isnan;
397
398 // From SimplicialComplexForAlpha type required to insert into a simplicial complex (with or without subfaces).
400 using Simplex_handle = typename SimplicialComplexForAlpha::Simplex_handle;
401 using Vector_vertex = std::vector<Vertex_handle>;
402
403 // For users to be able to define their own sqrt function on their desired Filtration_value type
404 using std::sqrt;
405
406 if (triangulation_ == nullptr) {
407 std::cerr << "Alpha_complex cannot create_complex from a NULL triangulation\n";
408 return false; // ----- >>
409 }
410 if (triangulation_->maximal_dimension() < 1) {
411 std::cerr << "Alpha_complex cannot create_complex from a zero-dimension triangulation\n";
412 return false; // ----- >>
413 }
414 if (complex.num_vertices() > 0) {
415 std::cerr << "Alpha_complex create_complex - complex is not empty\n";
416 return false; // ----- >>
417 }
418
419 // --------------------------------------------------------------------------------------------
420 // Simplex_tree construction from loop on triangulation finite full cells list
421 if (num_vertices() > 0) {
422 std::vector<Vertex_handle> one_vertex(1);
423 for (auto vertex : vertices_) {
424#ifdef DEBUG_TRACES
425 std::clog << "SimplicialComplex insertion " << vertex << std::endl;
426#endif // DEBUG_TRACES
427 one_vertex[0] = vertex;
428 complex.insert_simplex_and_subfaces(one_vertex, std::numeric_limits<Filtration_value>::quiet_NaN());
429 }
430
431 for (auto cit = triangulation_->finite_full_cells_begin();
432 cit != triangulation_->finite_full_cells_end();
433 ++cit) {
434 Vector_vertex vertexVector;
435#ifdef DEBUG_TRACES
436 std::clog << "SimplicialComplex insertion ";
437#endif // DEBUG_TRACES
438 for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
439 if (*vit != nullptr) {
440#ifdef DEBUG_TRACES
441 std::clog << " " << (*vit)->data();
442#endif // DEBUG_TRACES
443 // Vector of vertex construction for simplex_tree structure
444 vertexVector.push_back((*vit)->data());
445 }
446 }
447#ifdef DEBUG_TRACES
448 std::clog << std::endl;
449#endif // DEBUG_TRACES
450 // Insert each simplex and its subfaces in the simplex tree - filtration is NaN
451 complex.insert_simplex_and_subfaces(vertexVector, std::numeric_limits<Filtration_value>::quiet_NaN());
452 }
453 }
454 // --------------------------------------------------------------------------------------------
455 if (!default_filtration_value) {
456 CGAL::NT_converter<FT, Filtration_value> cgal_converter;
457 // --------------------------------------------------------------------------------------------
458 // ### For i : d -> 0
459 for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) {
460 // ### Foreach Sigma of dim i
461 for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) {
462 int f_simplex_dim = complex.dimension(f_simplex);
463 if (decr_dim == f_simplex_dim) {
464 // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma)
465 if (isnan(complex.filtration(f_simplex))) {
466 Filtration_value alpha_complex_filtration = 0.0;
467 // No need to compute squared_radius on a non-weighted single point - alpha is 0.0
468 if (Weighted || f_simplex_dim > 0) {
469 auto const& sqrad = radius(complex, f_simplex);
470#if CGAL_VERSION_NR >= 1050000000
471 if(exact) CGAL::exact(sqrad);
472#endif
473 alpha_complex_filtration = cgal_converter(sqrad);
474 if constexpr (!Output_squared_values) {
475 alpha_complex_filtration = sqrt(alpha_complex_filtration);
476 }
477 }
478 complex.assign_filtration(f_simplex, alpha_complex_filtration);
479#ifdef DEBUG_TRACES
480 std::clog << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl;
481#endif // DEBUG_TRACES
482 }
483 // No need to propagate further, unweighted points all have value 0
484 if (decr_dim > !Weighted)
485 propagate_alpha_filtration(complex, f_simplex);
486 }
487 }
488 old_cache_ = std::move(cache_);
489 cache_.clear();
490 }
491 // --------------------------------------------------------------------------------------------
492
493 // --------------------------------------------------------------------------------------------
494 if (!exact)
495 // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension
496 // Only in not exact version, cf. https://github.com/GUDHI/gudhi-devel/issues/57
498 // Remove all simplices that have a filtration value greater than max_alpha_square
499 if constexpr (!Output_squared_values) {
500 complex.prune_above_filtration(sqrt(max_alpha_square));
501 } else {
502 complex.prune_above_filtration(max_alpha_square);
503 }
504 // --------------------------------------------------------------------------------------------
505 }
506 return true;
507 }
508
509 private:
510 template <typename SimplicialComplexForAlpha, typename Simplex_handle>
511 void propagate_alpha_filtration(SimplicialComplexForAlpha& complex, Simplex_handle f_simplex) {
512 // From SimplicialComplexForAlpha type required to assign filtration values.
514 // To support more general types for Filtration_value
515 using std::isnan;
516
517 // ### Foreach Tau face of Sigma
518 for (auto face_opposite_vertex : complex.boundary_opposite_vertex_simplex_range(f_simplex)) {
519 auto f_boundary = face_opposite_vertex.first;
520#ifdef DEBUG_TRACES
521 std::clog << " | --------------------------------------------------\n";
522 std::clog << " | Tau ";
523 for (auto vertex : complex.simplex_vertex_range(f_boundary)) {
524 std::clog << vertex << " ";
525 }
526 std::clog << "is a face of Sigma\n";
527 std::clog << " | isnan(complex.filtration(Tau)=" << isnan(complex.filtration(f_boundary)) << std::endl;
528#endif // DEBUG_TRACES
529 // ### If filt(Tau) is not NaN
530 if (!isnan(complex.filtration(f_boundary))) {
531 // ### filt(Tau) = fmin(filt(Tau), filt(Sigma))
532 Filtration_value alpha_complex_filtration = fmin(complex.filtration(f_boundary),
533 complex.filtration(f_simplex));
534 complex.assign_filtration(f_boundary, alpha_complex_filtration);
535#ifdef DEBUG_TRACES
536 std::clog << " | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << complex.filtration(f_boundary) << std::endl;
537#endif // DEBUG_TRACES
538 // ### Else
539 } else {
540 auto const& cache=get_cache(complex, f_boundary);
541 bool is_gab = kernel_.is_gabriel(cache, get_point_(face_opposite_vertex.second));
542#ifdef DEBUG_TRACES
543 std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << face_opposite_vertex.second << std::endl;
544#endif // DEBUG_TRACES
545 // ### If Tau is not Gabriel of Sigma
546 if (false == is_gab) {
547 // ### filt(Tau) = filt(Sigma)
548 Filtration_value alpha_complex_filtration = complex.filtration(f_simplex);
549 complex.assign_filtration(f_boundary, alpha_complex_filtration);
550#ifdef DEBUG_TRACES
551 std::clog << " | filt(Tau) = filt(Sigma) = " << complex.filtration(f_boundary) << std::endl;
552#endif // DEBUG_TRACES
553 }
554 }
555 }
556 }
557};
558
559} // namespace alpha_complex
560
561namespace alphacomplex = alpha_complex;
562
563} // namespace Gudhi
564
565#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
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:124
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:135
Alpha_complex(const InputPointRange &points, WeightRange weights)
Alpha_complex constructor from a list of points and weights.
Definition: Alpha_complex.h:211
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:110
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:386
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:242
Alpha_complex(const std::string &off_file_name)
Alpha_complex constructor from an OFF file name.
Definition: Alpha_complex.h:173
typename Geom_traits::Point_d Point_d
A point, or a weighted point in Euclidean space.
Definition: Alpha_complex.h:138
std::size_t num_vertices() const
Returns the number of finite vertices in the triangulation.
Definition: Alpha_complex.h:229
Alpha_complex(const InputPointRange &points)
Alpha_complex constructor from a list of points.
Definition: Alpha_complex.h:194
Alpha complex kernel container.
Definition: Alpha_kernel_d.h:42
Gudhi namespace.
Definition: SimplicialComplexForAlpha.h:14
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:72
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:29
Handle type for the vertices of a cell complex.
Definition: VertexHandle.h:15