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/Epick_d.h>
24 #include <CGAL/Spatial_sort_traits_adapter_d.h>
25 #include <CGAL/property_map.h>
26 #include <CGAL/NT_converter.h>
27 #include <CGAL/version.h>
29 #include <Eigen/src/Core/util/Macros.h>
41 #if CGAL_VERSION_NR < 1041101000
42 # error Alpha_complex_3d is only available for CGAL >= 4.11
45 #if !EIGEN_VERSION_AT_LEAST(3,1,0)
46 # error Alpha_complex_3d is only available for Eigen3 >= 3.1.0 installed with CGAL
51 namespace alpha_complex {
76 template<
class Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>>
80 typedef CGAL::Triangulation_data_structure<
typename Kernel::Dimension,
81 CGAL::Triangulation_vertex<Kernel, std::ptrdiff_t>,
82 CGAL::Triangulation_full_cell<Kernel> > TDS;
93 typedef typename Kernel::Compute_squared_radius_d Squared_Radius;
94 typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel;
95 typedef typename Kernel::Point_dimension_d Point_Dimension;
98 typedef typename std::vector<Point_d> Vector_of_CGAL_points;
101 typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator;
104 typedef typename Delaunay_triangulation::size_type size_type;
107 typedef typename std::map< std::size_t, CGAL_vertex_iterator > Vector_vertex_iterator;
112 Vector_vertex_iterator vertex_handle_to_iterator_;
129 : triangulation_(nullptr) {
132 std::cerr <<
"Alpha_complex - Unable to read file " << off_file_name <<
"\n";
148 template<
typename InputPo
intRange >
150 : triangulation_(nullptr) {
151 init_from_range(points);
157 delete triangulation_;
173 return vertex_handle_to_iterator_.at(vertex)->point();
181 return vertex_handle_to_iterator_.size();
185 template<
typename InputPo
intRange >
186 void init_from_range(
const InputPointRange& points) {
187 auto first = std::begin(points);
188 auto last = std::end(points);
192 Point_Dimension point_dimension = kernel_.point_dimension_d_object();
197 std::vector<Point_d> point_cloud(first, last);
200 std::vector<std::ptrdiff_t> indices(boost::counting_iterator<std::ptrdiff_t>(0),
201 boost::counting_iterator<std::ptrdiff_t>(point_cloud.size()));
203 typedef boost::iterator_property_map<typename std::vector<Point_d>::iterator,
204 CGAL::Identity_property_map<std::ptrdiff_t>> Point_property_map;
205 typedef CGAL::Spatial_sort_traits_adapter_d<Kernel, Point_property_map> Search_traits_d;
207 CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud)));
209 typename Delaunay_triangulation::Full_cell_handle hint;
210 for (
auto index : indices) {
211 typename Delaunay_triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint);
214 hint = pos->full_cell();
219 for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) {
220 if (!triangulation_->is_infinite(*vit)) {
222 std::cout <<
"Vertex insertion - " << vit->data() <<
" -> " << vit->point() << std::endl;
223 #endif // DEBUG_TRACES
224 vertex_handle_to_iterator_.emplace(vit->data(), vit);
248 template <
typename SimplicialComplexForAlpha,
251 Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity()) {
255 typedef std::vector<Vertex_handle> Vector_vertex;
257 if (triangulation_ ==
nullptr) {
258 std::cerr <<
"Alpha_complex cannot create_complex from a NULL triangulation\n";
261 if (triangulation_->maximal_dimension() < 1) {
262 std::cerr <<
"Alpha_complex cannot create_complex from a zero-dimension triangulation\n";
266 std::cerr <<
"Alpha_complex create_complex - complex is not empty\n";
272 if (triangulation_->number_of_vertices() > 0) {
273 for (
auto cit = triangulation_->finite_full_cells_begin();
274 cit != triangulation_->finite_full_cells_end();
276 Vector_vertex vertexVector;
278 std::cout <<
"Simplex_tree insertion ";
279 #endif // DEBUG_TRACES
280 for (
auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
281 if (*vit !=
nullptr) {
283 std::cout <<
" " << (*vit)->data();
284 #endif // DEBUG_TRACES
286 vertexVector.push_back((*vit)->data());
290 std::cout << std::endl;
291 #endif // DEBUG_TRACES
300 Vector_of_CGAL_points pointVector;
302 for (
int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) {
305 int f_simplex_dim = complex.
dimension(f_simplex);
306 if (decr_dim == f_simplex_dim) {
309 std::cout <<
"Sigma of dim " << decr_dim <<
" is";
310 #endif // DEBUG_TRACES
312 pointVector.push_back(
get_point(vertex));
314 std::cout <<
" " << vertex;
315 #endif // DEBUG_TRACES
318 std::cout << std::endl;
319 #endif // DEBUG_TRACES
321 if (std::isnan(complex.filtration(f_simplex))) {
324 if (f_simplex_dim > 0) {
326 Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object();
327 CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv;
329 alpha_complex_filtration = cv(squared_radius(pointVector.begin(), pointVector.end()));
333 std::cout <<
"filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl;
334 #endif // DEBUG_TRACES
338 propagate_alpha_filtration(complex, f_simplex);
354 template <
typename SimplicialComplexForAlpha,
typename Simplex_handle>
360 #endif // DEBUG_TRACES
365 std::cout <<
" | --------------------------------------------------\n";
366 std::cout <<
" | Tau ";
368 std::cout << vertex <<
" ";
370 std::cout <<
"is a face of Sigma\n";
371 std::cout <<
" | isnan(complex.filtration(Tau)=" << std::isnan(complex.filtration(f_boundary)) << std::endl;
372 #endif // DEBUG_TRACES
374 if (!std::isnan(complex.filtration(f_boundary))) {
376 Filtration_value alpha_complex_filtration = fmin(complex.filtration(f_boundary),
377 complex.filtration(f_simplex));
380 std::cout <<
" | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << complex.filtration(f_boundary) << std::endl;
381 #endif // DEBUG_TRACES
385 Vector_of_CGAL_points pointVector;
387 Vertex_handle vertexForGabriel = Vertex_handle();
388 #endif // DEBUG_TRACES
390 pointVector.push_back(
get_point(vertex));
396 if (std::find(pointVector.begin(), pointVector.end(), point_for_gabriel) == pointVector.end()) {
399 vertexForGabriel = vertex;
400 #endif // DEBUG_TRACES
406 Is_Gabriel is_gabriel = kernel_.side_of_bounded_sphere_d_object();
407 bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), point_for_gabriel)
408 != CGAL::ON_BOUNDED_SIDE;
410 std::cout <<
" | Tau is_gabriel(Sigma)=" << is_gab <<
" - vertexForGabriel=" << vertexForGabriel << std::endl;
411 #endif // DEBUG_TRACES
413 if (
false == is_gab) {
418 std::cout <<
" | filt(Tau) = filt(Sigma) = " << complex.filtration(f_boundary) << std::endl;
419 #endif // DEBUG_TRACES
428 namespace alphacomplex = alpha_complex;
432 #endif // ALPHA_COMPLEX_H_