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