12 #ifndef ALPHA_COMPLEX_3D_H_ 13 #define ALPHA_COMPLEX_3D_H_ 15 #include <boost/version.hpp> 16 #include <boost/variant.hpp> 17 #include <boost/range/size.hpp> 18 #include <boost/range/combine.hpp> 19 #include <boost/range/adaptor/transformed.hpp> 21 #include <gudhi/Debug_utils.h> 22 #include <gudhi/Alpha_complex_options.h> 24 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h> 25 #include <CGAL/Exact_predicates_exact_constructions_kernel.h> 26 #include <CGAL/Delaunay_triangulation_3.h> 27 #include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h> 28 #include <CGAL/Periodic_3_Delaunay_triangulation_3.h> 29 #include <CGAL/Periodic_3_regular_triangulation_traits_3.h> 30 #include <CGAL/Periodic_3_regular_triangulation_3.h> 31 #include <CGAL/Regular_triangulation_3.h> 32 #include <CGAL/Alpha_shape_3.h> 33 #include <CGAL/Alpha_shape_cell_base_3.h> 34 #include <CGAL/Alpha_shape_vertex_base_3.h> 36 #include <CGAL/Object.h> 37 #include <CGAL/tuple.h> 38 #include <CGAL/iterator.h> 39 #include <CGAL/version.h> 41 #include <boost/container/static_vector.hpp> 45 #include <unordered_map> 49 #include <type_traits> 53 #if CGAL_VERSION_NR < 1041101000 54 # error Alpha_complex_3d is only available for CGAL >= 4.11 59 namespace alpha_complex {
61 thread_local
double RELATIVE_PRECISION_OF_TO_DOUBLE = 0.00001;
68 template <complexity Complexity>
69 struct Value_from_iterator {
70 template <
typename Iterator>
71 static double perform(Iterator it) {
73 return CGAL::to_double(*it);
79 template <
typename Iterator>
80 static double perform(Iterator it) {
81 return CGAL::to_double(it->exact());
120 template <complexity Complexity = complexity::SAFE,
bool Weighted = false,
bool Periodic = false>
135 using Predicates =
typename std::conditional<(Complexity ==
complexity::FAST),
136 CGAL::Exact_predicates_inexact_constructions_kernel,
137 CGAL::Exact_predicates_exact_constructions_kernel>::type;
140 template <
typename Predicates,
bool Weighted_version,
bool Periodic_version>
143 template <
typename Predicates,
bool Is_periodic>
144 struct Kernel_3<Predicates, Is_periodic, false> {
145 using Kernel = Predicates;
148 template <
typename Predicates>
149 struct Kernel_3<Predicates, false, true> {
150 using Kernel = CGAL::Periodic_3_Delaunay_triangulation_traits_3<Predicates>;
152 template <
typename Predicates>
153 struct Kernel_3<Predicates, true, true> {
154 using Kernel = CGAL::Periodic_3_regular_triangulation_traits_3<Predicates>;
158 using Kernel =
typename Kernel_3<Predicates, Weighted, Periodic>::Kernel;
161 using TdsVb =
typename std::conditional<Periodic, CGAL::Periodic_3_triangulation_ds_vertex_base_3<>,
162 CGAL::Triangulation_ds_vertex_base_3<>>::type;
164 using Tvb =
typename std::conditional<Weighted, CGAL::Regular_triangulation_vertex_base_3<Kernel, TdsVb>,
165 CGAL::Triangulation_vertex_base_3<Kernel, TdsVb>>::type;
167 using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel, Tvb>;
169 using TdsCb =
typename std::conditional<Periodic, CGAL::Periodic_3_triangulation_ds_cell_base_3<>,
170 CGAL::Triangulation_ds_cell_base_3<>>::type;
172 using Tcb =
typename std::conditional<Weighted, CGAL::Regular_triangulation_cell_base_3<Kernel, TdsCb>,
173 CGAL::Triangulation_cell_base_3<Kernel, TdsCb>>::type;
175 using Cb = CGAL::Alpha_shape_cell_base_3<Kernel, Tcb>;
176 using Tds = CGAL::Triangulation_data_structure_3<Vb, Cb>;
179 template <
typename Kernel,
typename Tds,
bool Weighted_version,
bool Periodic_version>
180 struct Triangulation_3 {};
182 template <
typename Kernel,
typename Tds>
183 struct Triangulation_3<Kernel, Tds, false, false> {
184 using Dt = CGAL::Delaunay_triangulation_3<Kernel, Tds>;
187 template <
typename Kernel,
typename Tds>
188 struct Triangulation_3<Kernel, Tds, true, false> {
189 using Dt = CGAL::Regular_triangulation_3<Kernel, Tds>;
192 template <
typename Kernel,
typename Tds>
193 struct Triangulation_3<Kernel, Tds, false, true> {
194 using Dt = CGAL::Periodic_3_Delaunay_triangulation_3<Kernel, Tds>;
197 template <
typename Kernel,
typename Tds>
198 struct Triangulation_3<Kernel, Tds, true, true> {
199 using Dt = CGAL::Periodic_3_regular_triangulation_3<Kernel, Tds>;
210 using Dt =
typename Triangulation_3<Kernel, Tds, Weighted, Periodic>::Dt;
222 using FT =
typename Alpha_shape_3::FT;
246 using Weighted_point_3 =
typename Triangulation_3<Kernel, Tds, Weighted, Periodic>::Weighted_point_3;
251 using Point_3 =
typename Alpha_shape_3::Point;
255 CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<CGAL::Object, FT>,
256 CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<CGAL::Object>>,
257 std::back_insert_iterator<std::vector<FT>>>>;
259 using Cell_handle =
typename Alpha_shape_3::Cell_handle;
260 using Facet =
typename Alpha_shape_3::Facet;
261 using Edge =
typename Alpha_shape_3::Edge;
262 using Alpha_vertex_handle =
typename Alpha_shape_3::Vertex_handle;
263 using Vertex_list = boost::container::static_vector<Alpha_vertex_handle, 4>;
275 template <
typename InputPo
intRange>
277 static_assert(!Periodic,
"This constructor is not available for periodic versions of Alpha_complex_3d");
279 alpha_shape_3_ptr_ = std::make_unique<Alpha_shape_3>(
280 std::begin(points), std::end(points), 0, Alpha_shape_3::GENERAL);
297 template <
typename InputPo
intRange,
typename WeightRange>
299 static_assert(Weighted,
"This constructor is not available for non-weighted versions of Alpha_complex_3d");
300 static_assert(!Periodic,
"This constructor is not available for periodic versions of Alpha_complex_3d");
302 GUDHI_CHECK(boost::size(weights) == boost::size(points),
303 std::invalid_argument(
"Points number in range different from weights range number"));
305 auto weighted_points_3 = boost::range::combine(points, weights)
306 | boost::adaptors::transformed([](
auto const&t){
return Weighted_point_3(boost::get<0>(t), boost::get<1>(t));});
308 alpha_shape_3_ptr_ = std::make_unique<Alpha_shape_3>(
309 std::begin(weighted_points_3), std::end(weighted_points_3), 0, Alpha_shape_3::GENERAL);
332 template <
typename InputPo
intRange>
334 FT z_min,
FT x_max,
FT y_max,
FT z_max) {
335 static_assert(Periodic,
"This constructor is not available for non-periodic versions of Alpha_complex_3d");
338 (x_max - x_min == y_max - y_min) && (x_max - x_min == z_max - z_min) && (z_max - z_min == y_max - y_min),
339 std::invalid_argument(
"The size of the cuboid in every directions is not the same."));
342 Dt pdt(
typename Kernel::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max));
344 pdt.insert(std::begin(points), std::end(points),
true);
346 if (!pdt.is_triangulation_in_1_sheet()) {
347 throw std::invalid_argument(
"Unable to construct a triangulation within a single periodic domain.");
349 pdt.convert_to_1_sheeted_covering();
353 alpha_shape_3_ptr_ = std::make_unique<Alpha_shape_3>(pdt, 0, Alpha_shape_3::GENERAL);
380 template <
typename InputPo
intRange,
typename WeightRange>
382 FT z_min,
FT x_max,
FT y_max,
FT z_max) {
383 static_assert(Weighted,
"This constructor is not available for non-weighted versions of Alpha_complex_3d");
384 static_assert(Periodic,
"This constructor is not available for non-periodic versions of Alpha_complex_3d");
386 GUDHI_CHECK(boost::size(weights) == boost::size(points),
387 std::invalid_argument(
"Points number in range different from weights range number"));
390 (x_max - x_min == y_max - y_min) && (x_max - x_min == z_max - z_min) && (z_max - z_min == y_max - y_min),
391 std::invalid_argument(
"The size of the cuboid in every directions is not the same."));
395 FT maximal_possible_weight = 0.015625 * (x_max - x_min) * (x_max - x_min);
398 auto weighted_points_3 = boost::range::combine(points, weights)
399 | boost::adaptors::transformed([=](
auto const&t){
400 auto w = boost::get<1>(t);
401 GUDHI_CHECK((w < maximal_possible_weight) && (w >= 0),
402 std::invalid_argument(
"Invalid weight " + std::to_string(w) +
403 ". Must be non-negative and less than maximal possible weight = 1/64*cuboid length squared."));
408 Dt pdt(
typename Kernel::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max));
410 pdt.insert(std::begin(weighted_points_3), std::end(weighted_points_3),
true);
412 if (!pdt.is_triangulation_in_1_sheet()) {
413 throw std::invalid_argument(
"Unable to construct a triangulation within a single periodic domain.");
415 pdt.convert_to_1_sheeted_covering();
419 alpha_shape_3_ptr_ = std::make_unique<Alpha_shape_3>(pdt, 0, Alpha_shape_3::GENERAL);
439 Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity()) {
441 std::cerr <<
"Alpha_complex_3d create_complex - complex is not empty\n";
446 using Simplex_tree_vector_vertex = std::vector<Complex_vertex_handle>;
449 std::size_t count_vertices = 0;
450 std::size_t count_edges = 0;
451 std::size_t count_facets = 0;
452 std::size_t count_cells = 0;
453 #endif // DEBUG_TRACES 454 std::vector<CGAL::Object> objects;
455 std::vector<FT> alpha_values;
457 Dispatch dispatcher = CGAL::dispatch_output<CGAL::Object, FT>(std::back_inserter(objects),
458 std::back_inserter(alpha_values));
460 alpha_shape_3_ptr_->filtration_with_alpha_values(dispatcher);
462 std::clog <<
"filtration_with_alpha_values returns : " << objects.size() <<
" objects" << std::endl;
463 #endif // DEBUG_TRACES 464 if (objects.size() == 0) {
465 std::cerr <<
"Alpha_complex_3d create_complex - no triangulation as points are on a 2d plane\n";
469 using Alpha_value_iterator =
typename std::vector<FT>::const_iterator;
470 Alpha_value_iterator alpha_value_iterator = alpha_values.begin();
471 for (
auto object_iterator : objects) {
472 Vertex_list vertex_list;
475 if (
const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
476 for (
auto i = 0; i < 4; i++) {
478 std::clog <<
"from cell[" << i <<
"] - Point coordinates (" << (*cell)->vertex(i)->point() <<
")" 480 #endif // DEBUG_TRACES 481 vertex_list.push_back((*cell)->vertex(i));
485 #endif // DEBUG_TRACES 486 }
else if (
const Facet* facet = CGAL::object_cast<Facet>(&object_iterator)) {
487 for (
auto i = 0; i < 4; i++) {
488 if ((*facet).second != i) {
490 std::clog <<
"from facet=[" << i <<
"] - Point coordinates (" << (*facet).first->vertex(i)->point() <<
")" 492 #endif // DEBUG_TRACES 493 vertex_list.push_back((*facet).first->vertex(i));
498 #endif // DEBUG_TRACES 499 }
else if (
const Edge* edge = CGAL::object_cast<Edge>(&object_iterator)) {
500 for (
auto i : {(*edge).second, (*edge).third}) {
502 std::clog <<
"from edge[" << i <<
"] - Point coordinates (" << (*edge).first->vertex(i)->point() <<
")" 504 #endif // DEBUG_TRACES 505 vertex_list.push_back((*edge).first->vertex(i));
509 #endif // DEBUG_TRACES 510 }
else if (
const Alpha_vertex_handle* vertex = CGAL::object_cast<Alpha_vertex_handle>(&object_iterator)) {
513 std::clog <<
"from vertex - Point coordinates (" << (*vertex)->point() <<
")" << std::endl;
514 #endif // DEBUG_TRACES 515 vertex_list.push_back((*vertex));
518 Simplex_tree_vector_vertex the_simplex;
519 for (
auto the_alpha_shape_vertex : vertex_list) {
520 auto the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
521 if (the_map_iterator == map_cgal_simplex_tree.end()) {
523 Complex_vertex_handle vertex = map_cgal_simplex_tree.size();
525 std::clog <<
"Point (" << the_alpha_shape_vertex->point() <<
") not found - insert new vertex id " << vertex
527 #endif // DEBUG_TRACES 528 the_simplex.push_back(vertex);
529 map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex);
532 Complex_vertex_handle vertex = the_map_iterator->second;
534 std::clog <<
"Point (" << the_alpha_shape_vertex->point() <<
") found as vertex id " << vertex << std::endl;
535 #endif // DEBUG_TRACES 536 the_simplex.push_back(vertex);
540 Filtration_value filtr = Value_from_iterator<Complexity>::perform(alpha_value_iterator);
543 std::clog <<
"filtration = " << filtr << std::endl;
544 #endif // DEBUG_TRACES 545 complex.
insert_simplex(the_simplex, static_cast<Filtration_value>(filtr));
546 GUDHI_CHECK(alpha_value_iterator != alpha_values.end(),
"CGAL provided more simplices than values");
547 ++alpha_value_iterator;
551 std::clog <<
"vertices \t" << count_vertices << std::endl;
552 std::clog <<
"edges \t\t" << count_edges << std::endl;
553 std::clog <<
"facets \t\t" << count_facets << std::endl;
554 std::clog <<
"cells \t\t" << count_cells << std::endl;
555 #endif // DEBUG_TRACES 572 if (map_cgal_simplex_tree.size() != cgal_vertex_iterator_vector.size()) {
573 cgal_vertex_iterator_vector.resize(map_cgal_simplex_tree.size());
574 for (
auto map_iterator : map_cgal_simplex_tree) {
575 cgal_vertex_iterator_vector[map_iterator.second] = map_iterator.first;
579 auto cgal_vertex_iterator = cgal_vertex_iterator_vector.at(vertex);
580 return cgal_vertex_iterator->point();
585 std::unique_ptr<Alpha_shape_3> alpha_shape_3_ptr_;
588 std::unordered_map<Alpha_vertex_handle, std::size_t> map_cgal_simplex_tree;
590 std::vector<Alpha_vertex_handle> cgal_vertex_iterator_vector;
597 #endif // ALPHA_COMPLEX_3D_H_ bool create_complex(SimplicialComplexForAlpha3d &complex, Filtration_value max_alpha_square=std::numeric_limits< Filtration_value >::infinity())
Inserts all Delaunay triangulation into the simplicial complex. It also computes the filtration value...
Definition: Alpha_complex_3d.h:438
const Point_3 & get_point(std::size_t vertex)
get_point returns the point corresponding to the vertex given as parameter.
Definition: Alpha_complex_3d.h:571
std::size_t num_vertices()
Alpha_complex_3d(const InputPointRange &points, FT x_min, FT y_min, FT z_min, FT x_max, FT y_max, FT z_max)
Alpha_complex constructor from a list of points and an iso-cuboid coordinates.
Definition: Alpha_complex_3d.h:333
Alpha_complex_3d(const InputPointRange &points, WeightRange weights, FT x_min, FT y_min, FT z_min, FT x_max, FT y_max, FT z_max)
Alpha_complex constructor from a list of points, associated weights and an iso-cuboid coordinates...
Definition: Alpha_complex_3d.h:381
Alpha complex data structure for 3d specific case.
Definition: Alpha_complex_3d.h:121
typename Alpha_shape_3::Point Point_3
Alpha_complex_3d::Point_3 type is either a Alpha_complex_3d::Bare_point_3 (Weighted = false) or a Alp...
Definition: Alpha_complex_3d.h:251
Alpha_complex_3d(const InputPointRange &points, WeightRange weights)
Alpha_complex constructor from a list of points and associated weights.
Definition: Alpha_complex_3d.h:298
Definition: SimplicialComplexForAlpha.h:14
typename Triangulation_3< Kernel, Tds, Weighted, Periodic >::Weighted_point_3 Weighted_point_3
Gives public access to the Weighted_point_3 type. A Weighted point can be constructed as follows: ...
Definition: Alpha_complex_3d.h:246
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:20
void prune_above_filtration(Filtration_value filtration)
The concept SimplicialComplexForAlpha3d describes the requirements for a type to implement a simplici...
Definition: SimplicialComplexForAlpha3d.h:21
CGAL::Alpha_shape_3< Dt > Alpha_shape_3
The CGAL 3D Alpha Shapes type.
Definition: Alpha_complex_3d.h:218
Alpha_complex_3d(const InputPointRange &points)
Alpha_complex constructor from a list of points.
Definition: Alpha_complex_3d.h:276
void make_filtration_non_decreasing()
unspecified Filtration_value
Definition: SimplicialComplexForAlpha3d.h:25
typename Alpha_shape_3::FT FT
The alpha values type. Must be compatible with double.
Definition: Alpha_complex_3d.h:222
complexity
Alpha complex complexity template parameter possible values.
Definition: Alpha_complex_options.h:23
void insert_simplex(std::vector< Vertex_handle > const &vertex_range, Filtration_value filtration)
Inserts a simplex from a given simplex (represented by a vector of Vertex_handle) in the simplicial c...
unspecified Vertex_handle
Definition: SimplicialComplexForAlpha3d.h:23
typename Kernel::Point_3 Bare_point_3
Gives public access to the Bare_point_3 (bare aka. unweighed) type. Here is a Bare_point_3 constructo...
Definition: Alpha_complex_3d.h:233