106 using Internal_vertex_handle = std::ptrdiff_t;
110 using Geom_traits = std::conditional_t<Weighted, CGAL::Regular_triangulation_traits_adapter<Kernel>, Kernel>;
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>>;
118 using TDS = CGAL::Triangulation_data_structure<
typename Geom_traits::Dimension,
119 CGAL::Triangulation_vertex<Geom_traits, Internal_vertex_handle>,
120 Triangulation_full_cell >;
123 using Triangulation = std::conditional_t<Weighted, CGAL::Regular_triangulation<Kernel, TDS>,
124 CGAL::Delaunay_triangulation<Kernel, TDS>>;
130 using FT =
typename A_kernel_d::FT;
135 using Sphere =
typename A_kernel_d::Sphere;
138 using Point_d =
typename Geom_traits::Point_d;
142 using CGAL_vertex_iterator =
typename Triangulation::Vertex_iterator;
145 using Vector_vertex_iterator = std::vector< CGAL_vertex_iterator >;
150 Vector_vertex_iterator vertex_handle_to_iterator_;
152 std::unique_ptr<Triangulation> triangulation_;
158 std::vector<Internal_vertex_handle> vertices_;
161 std::vector<Sphere> cache_, old_cache_;
176 std::cerr <<
"Alpha_complex - Unable to read file " << off_file_name <<
"\n";
193 template<
typename InputPo
intRange >
195 init_from_range(points);
210 template <
typename InputPo
intRange,
typename WeightRange>
212 static_assert(Weighted,
"This constructor is not available for non-weighted versions of Alpha_complex");
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);
230 if (triangulation_ ==
nullptr)
233 return triangulation_->number_of_vertices();
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.");
249 template<
typename InputPo
intRange >
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."
257#if CGAL_VERSION_NR < 1050101000
259 static_assert(!Weighted,
"Weighted Alpha_complex is only available for CGAL >= 5.1.0");
262 auto first = std::begin(points);
263 auto last = std::end(points);
267 triangulation_ = std::make_unique<Triangulation>(kernel_.get_dimension(*first));
269 std::vector<Point_d> point_cloud(first, last);
272 std::vector<Internal_vertex_handle> indices(boost::counting_iterator<Internal_vertex_handle>(0),
273 boost::counting_iterator<Internal_vertex_handle>(point_cloud.size()));
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>;
279 CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud)));
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) {
287 hint = pos->full_cell();
293 vertex_handle_to_iterator_.resize(point_cloud.size());
296 vertices_.reserve(triangulation_->number_of_vertices());
298 for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) {
299 if (!triangulation_->is_infinite(*vit)) {
301 std::clog <<
"Vertex insertion - " << vit->data() <<
" -> " << vit->point() << std::endl;
303 vertex_handle_to_iterator_[vit->data()] = vit;
304 vertices_.push_back(vit->data());
307 std::sort(vertices_.begin(), vertices_.end());
318 const Point_d& get_point_(std::size_t vertex)
const {
319 return vertex_handle_to_iterator_[vertex]->point();
323 template<
class SimplicialComplexForAlpha>
325 auto k = cplx.key(s);
326 if(k==cplx.null_key()){
328 cplx.assign_key(s, k);
330 thread_local std::vector<Point_d> v;
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()));
340 template<
class SimplicialComplexForAlpha>
342 auto k = cplx.key(s);
343 if(k!=cplx.null_key())
344 return kernel_.get_squared_radius(old_cache_[k]);
346 thread_local std::vector<Point_d> v;
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());
377 template <
typename SimplicialComplexForAlpha,
380 Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(),
382 bool default_filtration_value =
false) {
384 static_assert(std::numeric_limits<Filtration_value>::has_quiet_NaN);
391 using Vector_vertex = std::vector<Vertex_handle>;
393 if (triangulation_ ==
nullptr) {
394 std::cerr <<
"Alpha_complex cannot create_complex from a NULL triangulation\n";
397 if (triangulation_->maximal_dimension() < 1) {
398 std::cerr <<
"Alpha_complex cannot create_complex from a zero-dimension triangulation\n";
402 std::cerr <<
"Alpha_complex create_complex - complex is not empty\n";
409 std::vector<Vertex_handle> one_vertex(1);
410 for (
auto vertex : vertices_) {
412 std::clog <<
"SimplicialComplex insertion " << vertex << std::endl;
414 one_vertex[0] = vertex;
418 for (
auto cit = triangulation_->finite_full_cells_begin();
419 cit != triangulation_->finite_full_cells_end();
421 Vector_vertex vertexVector;
423 std::clog <<
"SimplicialComplex insertion ";
425 for (
auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) {
426 if (*vit !=
nullptr) {
428 std::clog <<
" " << (*vit)->data();
431 vertexVector.push_back((*vit)->data());
435 std::clog << std::endl;
443 if (!default_filtration_value) {
444 CGAL::NT_converter<FT, Filtration_value> cgal_converter;
447 for (
int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) {
450 int f_simplex_dim = complex.
dimension(f_simplex);
451 if (decr_dim == f_simplex_dim) {
453 if (isnan(complex.filtration(f_simplex))) {
456 if (Weighted || f_simplex_dim > 0) {
457 auto const& sqrad = radius(complex, f_simplex);
458#if CGAL_VERSION_NR >= 1050000000
459 if(exact) CGAL::exact(sqrad);
461 alpha_complex_filtration = cgal_converter(sqrad);
465 std::clog <<
"filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl;
469 if (decr_dim > !Weighted)
470 propagate_alpha_filtration(complex, f_simplex);
473 old_cache_ = std::move(cache_);
491 template <
typename SimplicialComplexForAlpha,
typename Simplex_handle>
499 for (
auto face_opposite_vertex : complex.boundary_opposite_vertex_simplex_range(f_simplex)) {
500 auto f_boundary = face_opposite_vertex.first;
502 std::clog <<
" | --------------------------------------------------\n";
503 std::clog <<
" | Tau ";
504 for (
auto vertex : complex.simplex_vertex_range(f_boundary)) {
505 std::clog << vertex <<
" ";
507 std::clog <<
"is a face of Sigma\n";
508 std::clog <<
" | isnan(complex.filtration(Tau)=" << isnan(complex.filtration(f_boundary)) << std::endl;
511 if (!isnan(complex.filtration(f_boundary))) {
513 Filtration_value alpha_complex_filtration = fmin(complex.filtration(f_boundary),
514 complex.filtration(f_simplex));
517 std::clog <<
" | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << complex.filtration(f_boundary) << std::endl;
521 auto const& cache=get_cache(complex, f_boundary);
522 bool is_gab = kernel_.is_gabriel(cache, get_point_(face_opposite_vertex.second));
524 std::clog <<
" | Tau is_gabriel(Sigma)=" << is_gab <<
" - vertexForGabriel=" << face_opposite_vertex.second << std::endl;
527 if (
false == is_gab) {
532 std::clog <<
" | filt(Tau) = filt(Sigma) = " << complex.filtration(f_boundary) << std::endl;