12 #ifndef ALPHA_COMPLEX_H_ 13 #define ALPHA_COMPLEX_H_ 15 #include <gudhi/Alpha_complex/Alpha_kernel_d.h> 16 #include <gudhi/Debug_utils.h> 18 #include <gudhi/Points_off_io.h> 24 #include <CGAL/Delaunay_triangulation.h> 25 #include <CGAL/Regular_triangulation.h> 26 #include <CGAL/Epeck_d.h> 27 #include <CGAL/Epick_d.h> 28 #include <CGAL/Spatial_sort_traits_adapter_d.h> 29 #include <CGAL/property_map.h> 30 #include <CGAL/version.h> 31 #include <CGAL/NT_converter.h> 33 #include <Eigen/src/Core/util/Macros.h> 35 #include <boost/range/size.hpp> 36 #include <boost/range/combine.hpp> 37 #include <boost/range/adaptor/transformed.hpp> 49 #if CGAL_VERSION_NR < 1041101000 50 # error Alpha_complex is only available for CGAL >= 4.11 53 #if !EIGEN_VERSION_AT_LEAST(3,1,0) 54 # error Alpha_complex is only available for Eigen3 >= 3.1.0 installed with CGAL 59 namespace alpha_complex {
61 template<
typename D>
struct Is_Epeck_D {
static const bool value =
false; };
62 template<
typename D>
struct Is_Epeck_D<CGAL::Epeck_d<D>> {
static const bool value =
true; };
101 template<
class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>,
bool Weighted = false>
105 using Geom_traits = std::conditional_t<Weighted, CGAL::Regular_triangulation_traits_adapter<Kernel>, Kernel>;
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> >;
113 using Triangulation = std::conditional_t<Weighted, CGAL::Regular_triangulation<Kernel, TDS>,
114 CGAL::Delaunay_triangulation<Kernel, TDS>>;
120 using FT =
typename A_kernel_d::FT;
125 using Sphere =
typename A_kernel_d::Sphere;
128 using Point_d =
typename Geom_traits::Point_d;
132 using CGAL_vertex_iterator =
typename Triangulation::Vertex_iterator;
135 using size_type =
typename Triangulation::size_type;
138 using Vector_vertex_iterator = std::vector< CGAL_vertex_iterator >;
143 Vector_vertex_iterator vertex_handle_to_iterator_;
145 std::unique_ptr<Triangulation> triangulation_;
150 std::vector<Sphere> cache_, old_cache_;
165 std::cerr <<
"Alpha_complex - Unable to read file " << off_file_name <<
"\n";
182 template<
typename InputPo
intRange >
184 init_from_range(points);
199 template <
typename InputPo
intRange,
typename WeightRange>
201 static_assert(Weighted,
"This constructor is not available for non-weighted versions of Alpha_complex");
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);
223 return vertex_handle_to_iterator_.at(vertex)->point();
227 template<
typename InputPo
intRange >
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." 235 #if CGAL_VERSION_NR < 1050101000 237 static_assert(!Weighted,
"Weighted Alpha_complex is only available for CGAL >= 5.1");
240 auto first = std::begin(points);
241 auto last = std::end(points);
245 triangulation_ = std::make_unique<Triangulation>(kernel_.get_dimension(*first));
247 std::vector<Point_d> point_cloud(first, last);
250 std::vector<std::ptrdiff_t> indices(boost::counting_iterator<std::ptrdiff_t>(0),
251 boost::counting_iterator<std::ptrdiff_t>(point_cloud.size()));
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>;
257 CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud)));
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) {
265 hint = pos->full_cell();
271 vertex_handle_to_iterator_.resize(point_cloud.size());
273 for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) {
274 if (!triangulation_->is_infinite(*vit)) {
276 std::clog <<
"Vertex insertion - " << vit->data() <<
" -> " << vit->point() << std::endl;
277 #endif // DEBUG_TRACES 278 vertex_handle_to_iterator_[vit->data()] = vit;
291 const Point_d& get_point_(std::size_t vertex)
const {
292 return vertex_handle_to_iterator_[vertex]->point();
296 template<
class SimplicialComplexForAlpha>
298 auto k = cplx.
key(s);
303 thread_local std::vector<Point_d> v;
306 v.push_back(get_point_(vertex));
307 cache_.emplace_back(kernel_.get_sphere(v.cbegin(), v.cend()));
313 template<
class SimplicialComplexForAlpha>
315 auto k = cplx.
key(s);
317 return kernel_.get_squared_radius(old_cache_[k]);
319 thread_local std::vector<Point_d> v;
322 v.push_back(get_point_(vertex));
323 return kernel_.get_squared_radius(v.cbegin(), v.cend());
353 Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(),
355 bool default_filtration_value =
false) {
359 using Vector_vertex = std::vector<Vertex_handle>;
361 if (triangulation_ ==
nullptr) {
362 std::cerr <<
"Alpha_complex cannot create_complex from a NULL triangulation\n";
365 if (triangulation_->maximal_dimension() < 1) {
366 std::cerr <<
"Alpha_complex cannot create_complex from a zero-dimension triangulation\n";
370 std::cerr <<
"Alpha_complex create_complex - complex is not empty\n";
376 if (triangulation_->number_of_vertices() > 0) {
377 for (
auto cit = triangulation_->finite_full_cells_begin();
378 cit != triangulation_->finite_full_cells_end();
380 Vector_vertex vertexVector;
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) {
387 std::clog <<
" " << (*vit)->data();
388 #endif // DEBUG_TRACES 390 vertexVector.push_back((*vit)->data());
394 std::clog << std::endl;
395 #endif // DEBUG_TRACES 402 if (!default_filtration_value) {
403 CGAL::NT_converter<FT, Filtration_value> cgal_converter;
406 for (
int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) {
409 int f_simplex_dim = complex.
dimension(f_simplex);
410 if (decr_dim == f_simplex_dim) {
412 if (std::isnan(complex.filtration(f_simplex))) {
413 Filtration_value alpha_complex_filtration = 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);
420 alpha_complex_filtration = cgal_converter(sqrad);
424 std::clog <<
"filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl;
425 #endif // DEBUG_TRACES 428 if (decr_dim > !Weighted)
429 propagate_alpha_filtration(complex, f_simplex);
432 old_cache_ = std::move(cache_);
448 template <
typename SimplicialComplexForAlpha,
typename Simplex_handle>
457 std::clog <<
" | --------------------------------------------------\n";
458 std::clog <<
" | Tau ";
460 std::clog << vertex <<
" ";
462 std::clog <<
"is a face of Sigma\n";
463 std::clog <<
" | isnan(complex.filtration(Tau)=" << std::isnan(complex.filtration(f_boundary)) << std::endl;
464 #endif // DEBUG_TRACES 466 if (!std::isnan(complex.filtration(f_boundary))) {
468 Filtration_value alpha_complex_filtration = fmin(complex.filtration(f_boundary),
469 complex.filtration(f_simplex));
472 std::clog <<
" | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << complex.filtration(f_boundary) << std::endl;
473 #endif // DEBUG_TRACES 479 auto longiter = std::begin(longlist);
480 auto shortiter = std::begin(shortlist);
481 auto enditer = std::end(shortlist);
482 while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; }
484 auto const& cache=get_cache(complex, f_boundary);
485 bool is_gab = kernel_.is_gabriel(cache, get_point_(extra));
487 std::clog <<
" | Tau is_gabriel(Sigma)=" << is_gab <<
" - vertexForGabriel=" << extra << std::endl;
488 #endif // DEBUG_TRACES 490 if (
false == is_gab) {
492 Filtration_value alpha_complex_filtration = complex.filtration(f_simplex);
495 std::clog <<
" | filt(Tau) = filt(Sigma) = " << complex.filtration(f_boundary) << std::endl;
496 #endif // DEBUG_TRACES 505 namespace alphacomplex = alpha_complex;
509 #endif // ALPHA_COMPLEX_H_ 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
Alpha_complex(const InputPointRange &points)
Alpha_complex constructor from a list of points.
Definition: Alpha_complex.h:183
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
OFF file reader implementation in order to read points from an OFF file.
Definition: Points_off_io.h:122
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...
Simplex_vertex_range simplex_vertex_range(Simplex_handle const &simplex)
Returns a range over vertices of a given simplex.
Simplex_key null_key()
Returns a constant dummy number that is either negative, or at least as large as the number of simpli...
typename Geom_traits::Point_d Point_d
A point, or a weighted point in Euclidean space.
Definition: Alpha_complex.h:128
Definition: SimplicialComplexForAlpha.h:14
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
const std::vector< Point_d > & get_point_cloud() const
Point cloud getter.
Definition: Points_off_io.h:158
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
unspecified Vertex_handle
Definition: SimplicialComplexForAlpha.h:25
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:20
bool is_valid() const
Returns if the OFF file read operation was successful or not.
Definition: Points_off_io.h:150
The concept SimplicialComplexForAlpha describes the requirements for a type to implement a simplicial...
Definition: SimplicialComplexForAlpha.h:21
void make_filtration_non_decreasing()
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
int dimension(Simplex_handle simplex)
Alpha_complex(const InputPointRange &points, WeightRange weights)
Alpha_complex constructor from a list of points and weights.
Definition: Alpha_complex.h:200
void prune_above_filtration(Filtration_value filtration)
Alpha complex data structure.
Definition: Alpha_complex.h:102
Alpha complex kernel container.
Definition: Alpha_kernel_d.h:42
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)
unspecified Simplex_handle
Definition: SimplicialComplexForAlpha.h:23
std::size_t num_vertices()
Simplex_key key(Simplex_handle sh)
Returns the number stored for a simplex by assign_key().
Boundary_simplex_range boundary_simplex_range(Simplex_handle const &simplex)
Returns a range over boundaries of a given simplex.
Handle type for the vertices of a cell complex.
Definition: VertexHandle.h:15
unspecified Filtration_value
Definition: SimplicialComplexForAlpha.h:27
Alpha_complex(const std::string &off_file_name)
Alpha_complex constructor from an OFF file name.
Definition: Alpha_complex.h:162
void assign_key(Simplex_handle sh, Simplex_key n)
Store a number for a simplex, which can later be retrieved with key().