12 #ifndef ALPHA_COMPLEX_H_
13 #define ALPHA_COMPLEX_H_
15 #include <gudhi/Debug_utils.h>
17 #include <gudhi/Points_off_io.h>
22 #include <CGAL/Delaunay_triangulation.h>
23 #include <CGAL/Epeck_d.h>
24 #include <CGAL/Epick_d.h>
25 #include <CGAL/Spatial_sort_traits_adapter_d.h>
26 #include <CGAL/property_map.h>
27 #include <CGAL/version.h>
28 #include <CGAL/NT_converter.h>
30 #include <Eigen/src/Core/util/Macros.h>
42 #if CGAL_VERSION_NR < 1041101000
43 # error Alpha_complex is only available for CGAL >= 4.11
46 #if !EIGEN_VERSION_AT_LEAST(3,1,0)
47 # error Alpha_complex is only available for Eigen3 >= 3.1.0 installed with CGAL
52 namespace alpha_complex {
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; };
94 template<
class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>>
98 typedef CGAL::Triangulation_data_structure<
typename Kernel::Dimension,
99 CGAL::Triangulation_vertex<Kernel, std::ptrdiff_t>,
100 CGAL::Triangulation_full_cell<Kernel> > TDS;
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;
116 typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator;
119 typedef typename Delaunay_triangulation::size_type size_type;
122 typedef typename std::vector< CGAL_vertex_iterator > Vector_vertex_iterator;
125 typedef typename Kernel::FT FT;
130 Vector_vertex_iterator vertex_handle_to_iterator_;
136 std::vector<std::pair<Point_d, FT>> cache_, old_cache_;
149 : triangulation_(nullptr) {
152 std::cerr <<
"Alpha_complex - Unable to read file " << off_file_name <<
"\n";
168 template<
typename InputPo
intRange >
170 : triangulation_(nullptr) {
171 init_from_range(points);
177 delete triangulation_;
193 return vertex_handle_to_iterator_.at(vertex)->point();
197 template<
typename InputPo
intRange >
198 void init_from_range(
const InputPointRange& points) {
199 #if CGAL_VERSION_NR < 1050000000
200 if (Is_Epeck_D<Kernel>::value)
201 std::cerr <<
"It is strongly advised to use a CGAL version >= 5.0 with Epeck_d Kernel for performance reasons."
205 auto first = std::begin(points);
206 auto last = std::end(points);
210 Point_Dimension point_dimension = kernel_.point_dimension_d_object();
215 std::vector<Point_d> point_cloud(first, last);
218 std::vector<std::ptrdiff_t> indices(boost::counting_iterator<std::ptrdiff_t>(0),
219 boost::counting_iterator<std::ptrdiff_t>(point_cloud.size()));
221 typedef boost::iterator_property_map<typename std::vector<Point_d>::iterator,
222 CGAL::Identity_property_map<std::ptrdiff_t>> Point_property_map;
223 typedef CGAL::Spatial_sort_traits_adapter_d<Kernel, Point_property_map> Search_traits_d;
225 CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud)));
227 typename Delaunay_triangulation::Full_cell_handle hint;
228 for (
auto index : indices) {
229 typename Delaunay_triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint);
232 hint = pos->full_cell();
237 vertex_handle_to_iterator_.resize(point_cloud.size());
239 for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) {
240 if (!triangulation_->is_infinite(*vit)) {
242 std::clog <<
"Vertex insertion - " << vit->data() <<
" -> " << vit->point() << std::endl;
243 #endif // DEBUG_TRACES
244 vertex_handle_to_iterator_[vit->data()] = vit;
257 const Point_d& get_point_(std::size_t vertex)
const {
258 return vertex_handle_to_iterator_[vertex]->point();
262 template<
class SimplicialComplexForAlpha>
264 auto k = cplx.key(s);
265 if(k==cplx.null_key()){
267 cplx.assign_key(s, k);
269 thread_local std::vector<Point_d> v;
271 for (
auto vertex : cplx.simplex_vertex_range(s))
272 v.push_back(get_point_(vertex));
273 Point_d c = kernel_.construct_circumcenter_d_object()(v.cbegin(), v.cend());
274 FT r = kernel_.squared_distance_d_object()(c, v[0]);
275 cache_.emplace_back(std::move(c), std::move(r));
281 template<
class SimplicialComplexForAlpha>
283 auto k = cplx.key(s);
284 if(k!=cplx.null_key())
285 return old_cache_[k].second;
287 thread_local std::vector<Point_d> v;
289 for (
auto vertex : cplx.simplex_vertex_range(s))
290 v.push_back(get_point_(vertex));
291 return kernel_.compute_squared_radius_d_object()(v.cbegin(), v.cend());
318 template <
typename SimplicialComplexForAlpha,
321 Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(),
323 bool default_filtration_value =
false) {
327 typedef std::vector<Vertex_handle> Vector_vertex;
329 if (triangulation_ ==
nullptr) {
330 std::cerr <<
"Alpha_complex cannot create_complex from a NULL triangulation\n";
333 if (triangulation_->maximal_dimension() < 1) {
334 std::cerr <<
"Alpha_complex cannot create_complex from a zero-dimension triangulation\n";
338 std::cerr <<
"Alpha_complex create_complex - complex is not empty\n";
344 if (triangulation_->number_of_vertices() > 0) {
345 for (
auto cit = triangulation_->finite_full_cells_begin();
346 cit != triangulation_->finite_full_cells_end();
348 Vector_vertex vertexVector;
350 std::clog <<
"Simplex_tree insertion ";
351 #endif // DEBUG_TRACES
352 for (
auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
353 if (*vit !=
nullptr) {
355 std::clog <<
" " << (*vit)->data();
356 #endif // DEBUG_TRACES
358 vertexVector.push_back((*vit)->data());
362 std::clog << std::endl;
363 #endif // DEBUG_TRACES
370 if (!default_filtration_value) {
373 for (
int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) {
376 int f_simplex_dim = complex.
dimension(f_simplex);
377 if (decr_dim == f_simplex_dim) {
379 if (std::isnan(complex.filtration(f_simplex))) {
382 if (f_simplex_dim > 0) {
383 auto const& sqrad = radius(complex, f_simplex);
384 #if CGAL_VERSION_NR >= 1050000000
385 if(exact) CGAL::exact(sqrad);
387 CGAL::NT_converter<FT, Filtration_value> cv;
388 alpha_complex_filtration = cv(sqrad);
392 std::clog <<
"filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl;
393 #endif // DEBUG_TRACES
397 propagate_alpha_filtration(complex, f_simplex);
400 old_cache_ = std::move(cache_);
416 template <
typename SimplicialComplexForAlpha,
typename Simplex_handle>
425 std::clog <<
" | --------------------------------------------------\n";
426 std::clog <<
" | Tau ";
428 std::clog << vertex <<
" ";
430 std::clog <<
"is a face of Sigma\n";
431 std::clog <<
" | isnan(complex.filtration(Tau)=" << std::isnan(complex.filtration(f_boundary)) << std::endl;
432 #endif // DEBUG_TRACES
434 if (!std::isnan(complex.filtration(f_boundary))) {
436 Filtration_value alpha_complex_filtration = fmin(complex.filtration(f_boundary),
437 complex.filtration(f_simplex));
440 std::clog <<
" | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << complex.filtration(f_boundary) << std::endl;
441 #endif // DEBUG_TRACES
447 auto longiter = std::begin(longlist);
448 auto shortiter = std::begin(shortlist);
449 auto enditer = std::end(shortlist);
450 while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; }
451 Vertex_handle extra = *longiter;
452 auto const& cache=get_cache(complex, f_boundary);
453 bool is_gab = kernel_.squared_distance_d_object()(cache.first, get_point_(extra)) >= cache.second;
455 std::clog <<
" | Tau is_gabriel(Sigma)=" << is_gab <<
" - vertexForGabriel=" << extra << std::endl;
458 if (
false == is_gab) {
463 std::clog <<
" | filt(Tau) = filt(Sigma) = " << complex.filtration(f_boundary) << std::endl;
464 #endif // DEBUG_TRACES
473 namespace alphacomplex = alpha_complex;
477 #endif // ALPHA_COMPLEX_H_