diff options
Diffstat (limited to 'src/Alpha_complex/include/gudhi')
| -rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_complex.h | 432 | ||||
| -rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_complex_3d.h | 579 | ||||
| -rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_complex_options.h | 33 |
3 files changed, 1044 insertions, 0 deletions
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h new file mode 100644 index 00000000..8919cdb9 --- /dev/null +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -0,0 +1,432 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL and Eigen3 + * - YYYY/MM Author: Description of the modification + */ + +#ifndef ALPHA_COMPLEX_H_ +#define ALPHA_COMPLEX_H_ + +#include <gudhi/Debug_utils.h> +// to construct Alpha_complex from a OFF file of points +#include <gudhi/Points_off_io.h> + +#include <stdlib.h> +#include <math.h> // isnan, fmax + +#include <CGAL/Delaunay_triangulation.h> +#include <CGAL/Epick_d.h> +#include <CGAL/Spatial_sort_traits_adapter_d.h> +#include <CGAL/property_map.h> // for CGAL::Identity_property_map +#include <CGAL/NT_converter.h> +#include <CGAL/version.h> // for CGAL_VERSION_NR + +#include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST + +#include <iostream> +#include <vector> +#include <string> +#include <limits> // NaN +#include <map> +#include <utility> // std::pair +#include <stdexcept> +#include <numeric> // for std::iota + +// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 +#if CGAL_VERSION_NR < 1041101000 +# error Alpha_complex_3d is only available for CGAL >= 4.11 +#endif + +#if !EIGEN_VERSION_AT_LEAST(3,1,0) +# error Alpha_complex_3d is only available for Eigen3 >= 3.1.0 installed with CGAL +#endif + +namespace Gudhi { + +namespace alpha_complex { + +/** + * \class Alpha_complex Alpha_complex.h gudhi/Alpha_complex.h + * \brief Alpha complex data structure. + * + * \ingroup alpha_complex + * + * \details + * The data structure is constructing a CGAL Delaunay triangulation (for more informations on CGAL Delaunay + * triangulation, please refer to the corresponding chapter in page http://doc.cgal.org/latest/Triangulation/) from a + * range of points or from an OFF file (cf. Points_off_reader). + * + * Please refer to \ref alpha_complex for examples. + * + * The complex is a template class requiring an Epick_d <a target="_blank" + * href="http://doc.cgal.org/latest/Kernel_d/index.html#Chapter_dD_Geometry_Kernel">dD Geometry Kernel</a> + * \cite cgal:s-gkd-15b from CGAL as template, default value is <a target="_blank" + * href="http://doc.cgal.org/latest/Kernel_d/classCGAL_1_1Epick__d.html">CGAL::Epick_d</a> + * < <a target="_blank" href="http://doc.cgal.org/latest/Kernel_23/classCGAL_1_1Dynamic__dimension__tag.html"> + * CGAL::Dynamic_dimension_tag </a> > + * + * \remark When Alpha_complex is constructed with an infinite value of alpha, the complex is a Delaunay complex. + * + */ +template<class Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>> +class Alpha_complex { + public: + // Add an int in TDS to save point index in the structure + typedef CGAL::Triangulation_data_structure<typename Kernel::Dimension, + CGAL::Triangulation_vertex<Kernel, std::ptrdiff_t>, + CGAL::Triangulation_full_cell<Kernel> > TDS; + /** \brief A Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/ + typedef CGAL::Delaunay_triangulation<Kernel, TDS> Delaunay_triangulation; + + /** \brief A point in Euclidean space.*/ + typedef typename Kernel::Point_d Point_d; + /** \brief Geometric traits class that provides the geometric types and predicates needed by Delaunay + * triangulations.*/ + typedef Kernel Geom_traits; + + private: + typedef typename Kernel::Compute_squared_radius_d Squared_Radius; + typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel; + typedef typename Kernel::Point_dimension_d Point_Dimension; + + // Type required to compute squared radius, or side of bounded sphere on a vector of points. + typedef typename std::vector<Point_d> Vector_of_CGAL_points; + + // Vertex_iterator type from CGAL. + typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator; + + // size_type type from CGAL. + typedef typename Delaunay_triangulation::size_type size_type; + + // Map type to switch from simplex tree vertex handle to CGAL vertex iterator. + typedef typename std::map< std::size_t, CGAL_vertex_iterator > Vector_vertex_iterator; + + private: + /** \brief Vertex iterator vector to switch from simplex tree vertex handle to CGAL vertex iterator. + * Vertex handles are inserted sequentially, starting at 0.*/ + Vector_vertex_iterator vertex_handle_to_iterator_; + /** \brief Pointer on the CGAL Delaunay triangulation.*/ + Delaunay_triangulation* triangulation_; + /** \brief Kernel for triangulation_ functions access.*/ + Kernel kernel_; + + public: + /** \brief Alpha_complex constructor from an OFF file name. + * + * Uses the Points_off_reader to construct the Delaunay triangulation required to initialize + * the Alpha_complex. + * + * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous. + * + * @param[in] off_file_name OFF file [path and] name. + */ + Alpha_complex(const std::string& off_file_name) + : triangulation_(nullptr) { + Gudhi::Points_off_reader<Point_d> off_reader(off_file_name); + if (!off_reader.is_valid()) { + std::cerr << "Alpha_complex - Unable to read file " << off_file_name << "\n"; + exit(-1); // ----- >> + } + + init_from_range(off_reader.get_point_cloud()); + } + + /** \brief Alpha_complex constructor from a list of points. + * + * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous. + * + * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d + * + * The type InputPointRange must be a range for which std::begin and + * std::end return input iterators on a Kernel::Point_d. + */ + template<typename InputPointRange > + Alpha_complex(const InputPointRange& points) + : triangulation_(nullptr) { + init_from_range(points); + } + + /** \brief Alpha_complex destructor deletes the Delaunay triangulation. + */ + ~Alpha_complex() { + delete triangulation_; + } + + // Forbid copy/move constructor/assignment operator + Alpha_complex(const Alpha_complex& other) = delete; + Alpha_complex& operator= (const Alpha_complex& other) = delete; + Alpha_complex (Alpha_complex&& other) = delete; + Alpha_complex& operator= (Alpha_complex&& other) = delete; + + /** \brief get_point returns the point corresponding to the vertex given as parameter. + * + * @param[in] vertex Vertex handle of the point to retrieve. + * @return The point found. + * @exception std::out_of_range In case vertex is not found (cf. std::vector::at). + */ + const Point_d& get_point(std::size_t vertex) const { + return vertex_handle_to_iterator_.at(vertex)->point(); + } + + /** \brief number_of_vertices returns the number of vertices (same as the number of points). + * + * @return The number of vertices. + */ + std::size_t number_of_vertices() const { + return vertex_handle_to_iterator_.size(); + } + + private: + template<typename InputPointRange > + void init_from_range(const InputPointRange& points) { + auto first = std::begin(points); + auto last = std::end(points); + + if (first != last) { + // point_dimension function initialization + Point_Dimension point_dimension = kernel_.point_dimension_d_object(); + + // Delaunay triangulation is point dimension. + triangulation_ = new Delaunay_triangulation(point_dimension(*first)); + + std::vector<Point_d> point_cloud(first, last); + + // Creates a vector {0, 1, ..., N-1} + std::vector<std::ptrdiff_t> indices(boost::counting_iterator<std::ptrdiff_t>(0), + boost::counting_iterator<std::ptrdiff_t>(point_cloud.size())); + + typedef boost::iterator_property_map<typename std::vector<Point_d>::iterator, + CGAL::Identity_property_map<std::ptrdiff_t>> Point_property_map; + typedef CGAL::Spatial_sort_traits_adapter_d<Kernel, Point_property_map> Search_traits_d; + + CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud))); + + typename Delaunay_triangulation::Full_cell_handle hint; + for (auto index : indices) { + typename Delaunay_triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint); + // Save index value as data to retrieve it after insertion + pos->data() = index; + hint = pos->full_cell(); + } + // -------------------------------------------------------------------------------------------- + // double map to retrieve simplex tree vertex handles from CGAL vertex iterator and vice versa + // Loop on triangulation vertices list + for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) { + if (!triangulation_->is_infinite(*vit)) { +#ifdef DEBUG_TRACES + std::cout << "Vertex insertion - " << vit->data() << " -> " << vit->point() << std::endl; +#endif // DEBUG_TRACES + vertex_handle_to_iterator_.emplace(vit->data(), vit); + } + } + // -------------------------------------------------------------------------------------------- + } + } + + public: + /** \brief Inserts all Delaunay triangulation into the simplicial complex. + * It also computes the filtration values accordingly to the \ref createcomplexalgorithm + * + * \tparam SimplicialComplexForAlpha must meet `SimplicialComplexForAlpha` concept. + * + * @param[in] complex SimplicialComplexForAlpha to be created. + * @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$, and there is very + * little point using anything else since it does not save time. + * + * @return true if creation succeeds, false otherwise. + * + * @pre Delaunay triangulation must be already constructed with dimension strictly greater than 0. + * @pre The simplicial complex must be empty (no vertices) + * + * Initialization can be launched once. + */ + template <typename SimplicialComplexForAlpha, + typename Filtration_value = typename SimplicialComplexForAlpha::Filtration_value> + bool create_complex(SimplicialComplexForAlpha& complex, + Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity()) { + // From SimplicialComplexForAlpha type required to insert into a simplicial complex (with or without subfaces). + typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle; + typedef typename SimplicialComplexForAlpha::Simplex_handle Simplex_handle; + typedef std::vector<Vertex_handle> Vector_vertex; + + if (triangulation_ == nullptr) { + std::cerr << "Alpha_complex cannot create_complex from a NULL triangulation\n"; + return false; // ----- >> + } + if (triangulation_->maximal_dimension() < 1) { + std::cerr << "Alpha_complex cannot create_complex from a zero-dimension triangulation\n"; + return false; // ----- >> + } + if (complex.num_vertices() > 0) { + std::cerr << "Alpha_complex create_complex - complex is not empty\n"; + return false; // ----- >> + } + + // -------------------------------------------------------------------------------------------- + // Simplex_tree construction from loop on triangulation finite full cells list + if (triangulation_->number_of_vertices() > 0) { + for (auto cit = triangulation_->finite_full_cells_begin(); + cit != triangulation_->finite_full_cells_end(); + ++cit) { + Vector_vertex vertexVector; +#ifdef DEBUG_TRACES + std::cout << "Simplex_tree insertion "; +#endif // DEBUG_TRACES + for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) { + if (*vit != nullptr) { +#ifdef DEBUG_TRACES + std::cout << " " << (*vit)->data(); +#endif // DEBUG_TRACES + // Vector of vertex construction for simplex_tree structure + vertexVector.push_back((*vit)->data()); + } + } +#ifdef DEBUG_TRACES + std::cout << std::endl; +#endif // DEBUG_TRACES + // Insert each simplex and its subfaces in the simplex tree - filtration is NaN + complex.insert_simplex_and_subfaces(vertexVector, std::numeric_limits<double>::quiet_NaN()); + } + } + // -------------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------------- + // Will be re-used many times + Vector_of_CGAL_points pointVector; + // ### For i : d -> 0 + for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) { + // ### Foreach Sigma of dim i + for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) { + int f_simplex_dim = complex.dimension(f_simplex); + if (decr_dim == f_simplex_dim) { + pointVector.clear(); +#ifdef DEBUG_TRACES + std::cout << "Sigma of dim " << decr_dim << " is"; +#endif // DEBUG_TRACES + for (auto vertex : complex.simplex_vertex_range(f_simplex)) { + pointVector.push_back(get_point(vertex)); +#ifdef DEBUG_TRACES + std::cout << " " << vertex; +#endif // DEBUG_TRACES + } +#ifdef DEBUG_TRACES + std::cout << std::endl; +#endif // DEBUG_TRACES + // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma) + if (std::isnan(complex.filtration(f_simplex))) { + Filtration_value alpha_complex_filtration = 0.0; + // No need to compute squared_radius on a single point - alpha is 0.0 + if (f_simplex_dim > 0) { + // squared_radius function initialization + Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object(); + CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv; + + alpha_complex_filtration = cv(squared_radius(pointVector.begin(), pointVector.end())); + } + complex.assign_filtration(f_simplex, alpha_complex_filtration); +#ifdef DEBUG_TRACES + std::cout << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl; +#endif // DEBUG_TRACES + } + // No need to propagate further, unweighted points all have value 0 + if (decr_dim > 1) + propagate_alpha_filtration(complex, f_simplex); + } + } + } + // -------------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------------- + // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension + complex.make_filtration_non_decreasing(); + // Remove all simplices that have a filtration value greater than max_alpha_square + complex.prune_above_filtration(max_alpha_square); + // -------------------------------------------------------------------------------------------- + return true; + } + + private: + template <typename SimplicialComplexForAlpha, typename Simplex_handle> + void propagate_alpha_filtration(SimplicialComplexForAlpha& complex, Simplex_handle f_simplex) { + // From SimplicialComplexForAlpha type required to assign filtration values. + typedef typename SimplicialComplexForAlpha::Filtration_value Filtration_value; +#ifdef DEBUG_TRACES + typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle; +#endif // DEBUG_TRACES + + // ### Foreach Tau face of Sigma + for (auto f_boundary : complex.boundary_simplex_range(f_simplex)) { +#ifdef DEBUG_TRACES + std::cout << " | --------------------------------------------------\n"; + std::cout << " | Tau "; + for (auto vertex : complex.simplex_vertex_range(f_boundary)) { + std::cout << vertex << " "; + } + std::cout << "is a face of Sigma\n"; + std::cout << " | isnan(complex.filtration(Tau)=" << std::isnan(complex.filtration(f_boundary)) << std::endl; +#endif // DEBUG_TRACES + // ### If filt(Tau) is not NaN + if (!std::isnan(complex.filtration(f_boundary))) { + // ### filt(Tau) = fmin(filt(Tau), filt(Sigma)) + Filtration_value alpha_complex_filtration = fmin(complex.filtration(f_boundary), + complex.filtration(f_simplex)); + complex.assign_filtration(f_boundary, alpha_complex_filtration); +#ifdef DEBUG_TRACES + std::cout << " | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << complex.filtration(f_boundary) << std::endl; +#endif // DEBUG_TRACES + // ### Else + } else { + // insert the Tau points in a vector for is_gabriel function + Vector_of_CGAL_points pointVector; +#ifdef DEBUG_TRACES + Vertex_handle vertexForGabriel = Vertex_handle(); +#endif // DEBUG_TRACES + for (auto vertex : complex.simplex_vertex_range(f_boundary)) { + pointVector.push_back(get_point(vertex)); + } + // Retrieve the Sigma point that is not part of Tau - parameter for is_gabriel function + Point_d point_for_gabriel; + for (auto vertex : complex.simplex_vertex_range(f_simplex)) { + point_for_gabriel = get_point(vertex); + if (std::find(pointVector.begin(), pointVector.end(), point_for_gabriel) == pointVector.end()) { +#ifdef DEBUG_TRACES + // vertex is not found in Tau + vertexForGabriel = vertex; +#endif // DEBUG_TRACES + // No need to continue loop + break; + } + } + // is_gabriel function initialization + Is_Gabriel is_gabriel = kernel_.side_of_bounded_sphere_d_object(); + bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), point_for_gabriel) + != CGAL::ON_BOUNDED_SIDE; +#ifdef DEBUG_TRACES + std::cout << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << vertexForGabriel << std::endl; +#endif // DEBUG_TRACES + // ### If Tau is not Gabriel of Sigma + if (false == is_gab) { + // ### filt(Tau) = filt(Sigma) + Filtration_value alpha_complex_filtration = complex.filtration(f_simplex); + complex.assign_filtration(f_boundary, alpha_complex_filtration); +#ifdef DEBUG_TRACES + std::cout << " | filt(Tau) = filt(Sigma) = " << complex.filtration(f_boundary) << std::endl; +#endif // DEBUG_TRACES + } + } + } + } +}; + +} // namespace alpha_complex + +namespace alphacomplex = alpha_complex; + +} // namespace Gudhi + +#endif // ALPHA_COMPLEX_H_ diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h new file mode 100644 index 00000000..13ebb9c1 --- /dev/null +++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h @@ -0,0 +1,579 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2018 Inria + * + * Modification(s): + * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL and Eigen3 + * - YYYY/MM Author: Description of the modification + */ + +#ifndef ALPHA_COMPLEX_3D_H_ +#define ALPHA_COMPLEX_3D_H_ + +#include <boost/version.hpp> +#include <boost/variant.hpp> + +#include <gudhi/Debug_utils.h> +#include <gudhi/Alpha_complex_options.h> + +#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> +#include <CGAL/Exact_predicates_exact_constructions_kernel.h> +#include <CGAL/Delaunay_triangulation_3.h> +#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h> +#include <CGAL/Periodic_3_Delaunay_triangulation_3.h> +#include <CGAL/Periodic_3_regular_triangulation_traits_3.h> +#include <CGAL/Periodic_3_regular_triangulation_3.h> +#include <CGAL/Regular_triangulation_3.h> +#include <CGAL/Alpha_shape_3.h> +#include <CGAL/Alpha_shape_cell_base_3.h> +#include <CGAL/Alpha_shape_vertex_base_3.h> + +#include <CGAL/Object.h> +#include <CGAL/tuple.h> +#include <CGAL/iterator.h> +#include <CGAL/version.h> // for CGAL_VERSION_NR + +#include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST + +#include <boost/container/static_vector.hpp> + +#include <iostream> +#include <vector> +#include <unordered_map> +#include <stdexcept> +#include <cstddef> +#include <memory> // for std::unique_ptr +#include <type_traits> // for std::conditional and std::enable_if +#include <limits> // for numeric_limits<> + +// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 +#if CGAL_VERSION_NR < 1041101000 +# error Alpha_complex_3d is only available for CGAL >= 4.11 +#endif + +#if !EIGEN_VERSION_AT_LEAST(3,1,0) +# error Alpha_complex_3d is only available for Eigen3 >= 3.1.0 installed with CGAL +#endif + +namespace Gudhi { + +namespace alpha_complex { + +#ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL +thread_local +#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL + double RELATIVE_PRECISION_OF_TO_DOUBLE = 0.00001; + +// Value_from_iterator returns the filtration value from an iterator on alpha shapes values +// +// FAST SAFE EXACT +// CGAL::to_double(*iterator) CGAL::to_double(*iterator) CGAL::to_double(iterator->exact()) + +template <complexity Complexity> +struct Value_from_iterator { + template <typename Iterator> + static double perform(Iterator it) { + // Default behaviour + return CGAL::to_double(*it); + } +}; + +template <> +struct Value_from_iterator<complexity::EXACT> { + template <typename Iterator> + static double perform(Iterator it) { + return CGAL::to_double(it->exact()); + } +}; + +/** + * \class Alpha_complex_3d + * \brief Alpha complex data structure for 3d specific case. + * + * \ingroup alpha_complex + * + * \details + * The data structure is constructing a <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html">CGAL 3D Alpha + * Shapes</a> from a range of points (can be read from an OFF file, cf. Points_off_reader). + * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous. + * + * \tparam Complexity shall be `Gudhi::alpha_complex::complexity` type. Default value is + * `Gudhi::alpha_complex::complexity::SAFE`. + * + * \tparam Weighted Boolean used to set/unset the weighted version of Alpha_complex_3d. Default value is false. + * + * \tparam Periodic Boolean used to set/unset the periodic version of Alpha_complex_3d. Default value is false. + * + * For the weighted version, weights values are explained on CGAL + * <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0">Alpha shapes 3d</a> and + * <a href="https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation">Regular + * triangulation</a> documentation. + * + * For the periodic version, refer to the + * <a href="https://doc.cgal.org/latest/Periodic_3_triangulation_3/index.html">CGAL’s 3D Periodic Triangulations User + * Manual </a> for more details. + * The periodicity is defined by an iso-oriented cuboid with diagonal opposite vertices (x_min, y_min, z_min) and + * (x_max, y_max, z_max). + * + * Please refer to \ref alpha_complex for examples. + * + * \remark When Alpha_complex_3d is constructed with an infinite value of alpha (default value), the complex is a + * 3d Delaunay complex. + * + */ +template <complexity Complexity = complexity::SAFE, bool Weighted = false, bool Periodic = false> +class Alpha_complex_3d { + // Epick = Exact_predicates_inexact_constructions_kernel + // Epeck = Exact_predicates_exact_constructions_kernel + // Exact_alpha_comparison_tag = exact version of CGAL Alpha_shape_3 and of its objects (Alpha_shape_vertex_base_3 and + // Alpha_shape_cell_base_3). Not available if weighted or periodic. + // Can be CGAL::Tag_false or CGAL::Tag_true. Default is False. + // cf. https://doc.cgal.org/latest/Alpha_shapes_3/classCGAL_1_1Alpha__shape__3.html + // + // We could use Epick + CGAL::Tag_true for not weighted nor periodic, but during benchmark, we found a bug + // https://github.com/CGAL/cgal/issues/3460 + // This is the reason we only use Epick + CGAL::Tag_false, or Epeck + // + // FAST SAFE EXACT + // Epick + CGAL::Tag_false Epeck Epeck + using Predicates = typename std::conditional<(Complexity == complexity::FAST), + CGAL::Exact_predicates_inexact_constructions_kernel, + CGAL::Exact_predicates_exact_constructions_kernel>::type; + + // The other way to do a conditional type. Here there are 3 possibilities + template <typename Predicates, bool Weighted_version, bool Periodic_version> + struct Kernel_3 {}; + + template <typename Predicates, bool Is_periodic> + struct Kernel_3<Predicates, Is_periodic, false> { + using Kernel = Predicates; + }; + + template <typename Predicates> + struct Kernel_3<Predicates, false, true> { + using Kernel = CGAL::Periodic_3_Delaunay_triangulation_traits_3<Predicates>; + }; + template <typename Predicates> + struct Kernel_3<Predicates, true, true> { + using Kernel = CGAL::Periodic_3_regular_triangulation_traits_3<Predicates>; + }; + + using Kernel = typename Kernel_3<Predicates, Weighted, Periodic>::Kernel; + + using TdsVb = typename std::conditional<Periodic, CGAL::Periodic_3_triangulation_ds_vertex_base_3<>, + CGAL::Triangulation_ds_vertex_base_3<>>::type; + + using Tvb = typename std::conditional<Weighted, CGAL::Regular_triangulation_vertex_base_3<Kernel, TdsVb>, + CGAL::Triangulation_vertex_base_3<Kernel, TdsVb>>::type; + + using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel, Tvb>; + + using TdsCb = typename std::conditional<Periodic, CGAL::Periodic_3_triangulation_ds_cell_base_3<>, + CGAL::Triangulation_ds_cell_base_3<>>::type; + + using Tcb = typename std::conditional<Weighted, CGAL::Regular_triangulation_cell_base_3<Kernel, TdsCb>, + CGAL::Triangulation_cell_base_3<Kernel, TdsCb>>::type; + + using Cb = CGAL::Alpha_shape_cell_base_3<Kernel, Tcb>; + using Tds = CGAL::Triangulation_data_structure_3<Vb, Cb>; + + // The other way to do a conditional type. Here there 4 possibilities, cannot use std::conditional + template <typename Kernel, typename Tds, bool Weighted_version, bool Periodic_version> + struct Triangulation_3 {}; + + template <typename Kernel, typename Tds> + struct Triangulation_3<Kernel, Tds, false, false> { + using Dt = CGAL::Delaunay_triangulation_3<Kernel, Tds>; + using Weighted_point_3 = void; + }; + template <typename Kernel, typename Tds> + struct Triangulation_3<Kernel, Tds, true, false> { + using Dt = CGAL::Regular_triangulation_3<Kernel, Tds>; + using Weighted_point_3 = typename Dt::Weighted_point; + }; + template <typename Kernel, typename Tds> + struct Triangulation_3<Kernel, Tds, false, true> { + using Dt = CGAL::Periodic_3_Delaunay_triangulation_3<Kernel, Tds>; + using Weighted_point_3 = void; + }; + template <typename Kernel, typename Tds> + struct Triangulation_3<Kernel, Tds, true, true> { + using Dt = CGAL::Periodic_3_regular_triangulation_3<Kernel, Tds>; + using Weighted_point_3 = typename Dt::Weighted_point; + }; + + /** \brief Is either Delaunay_triangulation_3 (Weighted = false and Periodic = false), + * Regular_triangulation_3 (Weighted = true and Periodic = false), + * Periodic_3_Delaunay_triangulation_3 (Weighted = false and Periodic = true) + * or Periodic_3_regular_triangulation_3 (Weighted = true and Periodic = true). + * + * This type is required by `Gudhi::alpha_complex::Alpha_complex_3d::Alpha_shape_3`. + * */ + using Dt = typename Triangulation_3<Kernel, Tds, Weighted, Periodic>::Dt; + + public: + /** \brief The <a href="https://doc.cgal.org/latest/Alpha_shapes_3/classCGAL_1_1Alpha__shape__3.html">CGAL 3D Alpha + * Shapes</a> type. + * + * The `Gudhi::alpha_complex::Alpha_complex_3d` is a wrapper on top of this class to ease the standard, weighted + * and/or periodic build of the Alpha complex 3d.*/ + using Alpha_shape_3 = CGAL::Alpha_shape_3<Dt>; + + /** \brief The alpha values type. + * Must be compatible with double. */ + using FT = typename Alpha_shape_3::FT; + + /** \brief Gives public access to the Point_3 type. Here is a Point_3 constructor example: +\code{.cpp} +using Alpha_complex_3d = Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, false, false>; + +// x0 = 1., y0 = -1.1, z0 = -1.. +Alpha_complex_3d::Point_3 p0(1., -1.1, -1.); +\endcode + * */ + using Point_3 = typename Kernel::Point_3; + + /** \brief Gives public access to the Weighted_point_3 type. A Weighted point can be constructed as follows: +\code{.cpp} +using Weighted_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, false>; + +// x0 = 1., y0 = -1.1, z0 = -1., weight = 4. +Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point_3(1., -1.1, -1.), 4.); +\endcode + * + * Note: This type is defined to void if Alpha complex is not weighted. + * + * */ + using Weighted_point_3 = typename Triangulation_3<Kernel, Tds, Weighted, Periodic>::Weighted_point_3; + + private: + using Dispatch = + CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<CGAL::Object, FT>, + CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<CGAL::Object>>, + std::back_insert_iterator<std::vector<FT>>>>; + + using Cell_handle = typename Alpha_shape_3::Cell_handle; + using Facet = typename Alpha_shape_3::Facet; + using Edge = typename Alpha_shape_3::Edge; + using Alpha_vertex_handle = typename Alpha_shape_3::Vertex_handle; + using Vertex_list = boost::container::static_vector<Alpha_vertex_handle, 4>; + + public: + /** \brief Alpha_complex constructor from a list of points. + * + * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3` or + * `Alpha_complex_3d::Weighted_point_3`. + * + * @pre Available if Alpha_complex_3d is not Periodic. + * + * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a + * `Alpha_complex_3d::Point_3` or a `Alpha_complex_3d::Weighted_point_3`. + */ + template <typename InputPointRange> + Alpha_complex_3d(const InputPointRange& points) { + static_assert(!Periodic, "This constructor is not available for periodic versions of Alpha_complex_3d"); + + alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>( + new Alpha_shape_3(std::begin(points), std::end(points), 0, Alpha_shape_3::GENERAL)); + } + + /** \brief Alpha_complex constructor from a list of points and associated weights. + * + * @exception std::invalid_argument In debug mode, if points and weights do not have the same size. + * + * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3`. + * @param[in] weights Range of weights on points. Weights shall be in double. + * + * @pre Available if Alpha_complex_3d is Weighted and not Periodic. + * + * The type InputPointRange must be a range for which std::begin and + * std::end return input iterators on a `Alpha_complex_3d::Point_3`. + * The type WeightRange must be a range for which std::begin and + * std::end return an input iterator on a double. + */ + template <typename InputPointRange, typename WeightRange> + Alpha_complex_3d(const InputPointRange& points, WeightRange weights) { + static_assert(Weighted, "This constructor is not available for non-weighted versions of Alpha_complex_3d"); + static_assert(!Periodic, "This constructor is not available for periodic versions of Alpha_complex_3d"); + GUDHI_CHECK((weights.size() == points.size()), + std::invalid_argument("Points number in range different from weights range number")); + + std::vector<Weighted_point_3> weighted_points_3; + + std::size_t index = 0; + weighted_points_3.reserve(points.size()); + while ((index < weights.size()) && (index < points.size())) { + weighted_points_3.push_back(Weighted_point_3(points[index], weights[index])); + index++; + } + + alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>( + new Alpha_shape_3(std::begin(weighted_points_3), std::end(weighted_points_3), 0, Alpha_shape_3::GENERAL)); + } + + /** \brief Alpha_complex constructor from a list of points and an iso-cuboid coordinates. + * + * @exception std::invalid_argument In debug mode, if the size of the cuboid in every directions is not the same. + * + * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3` or + * `Alpha_complex_3d::Weighted_point_3`. + * @param[in] x_min Iso-oriented cuboid x_min. + * @param[in] y_min Iso-oriented cuboid y_min. + * @param[in] z_min Iso-oriented cuboid z_min. + * @param[in] x_max Iso-oriented cuboid x_max. + * @param[in] y_max Iso-oriented cuboid y_max. + * @param[in] z_max Iso-oriented cuboid z_max. + * + * @pre Available if Alpha_complex_3d is Periodic. + * + * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a + * `Alpha_complex_3d::Point_3` or a `Alpha_complex_3d::Weighted_point_3`. + * + * @note In weighted version, please check weights are greater than zero, and lower than 1/64*cuboid length + * squared. + */ + template <typename InputPointRange> + Alpha_complex_3d(const InputPointRange& points, FT x_min, FT y_min, + FT z_min, FT x_max, FT y_max, FT z_max) { + static_assert(Periodic, "This constructor is not available for non-periodic versions of Alpha_complex_3d"); + // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it. + GUDHI_CHECK( + (x_max - x_min == y_max - y_min) && (x_max - x_min == z_max - z_min) && (z_max - z_min == y_max - y_min), + std::invalid_argument("The size of the cuboid in every directions is not the same.")); + + // Define the periodic cube + Dt pdt(typename Kernel::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max)); + // Heuristic for inserting large point sets (if pts is reasonably large) + pdt.insert(std::begin(points), std::end(points), true); + // As pdt won't be modified anymore switch to 1-sheeted cover if possible + if (!pdt.is_triangulation_in_1_sheet()) { + throw std::invalid_argument("Unable to construct a triangulation within a single periodic domain."); + } + pdt.convert_to_1_sheeted_covering(); + + // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode + // Maybe need to set it to GENERAL mode + alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>(new Alpha_shape_3(pdt, 0, Alpha_shape_3::GENERAL)); + } + + /** \brief Alpha_complex constructor from a list of points, associated weights and an iso-cuboid coordinates. + * + * @exception std::invalid_argument In debug mode, if points and weights do not have the same size. + * @exception std::invalid_argument In debug mode, if the size of the cuboid in every directions is not the same. + * @exception std::invalid_argument In debug mode, if a weight is negative, zero, or greater than 1/64*cuboid length + * squared. + * + * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3`. + * @param[in] weights Range of weights on points. Weights shall be in double. + * @param[in] x_min Iso-oriented cuboid x_min. + * @param[in] y_min Iso-oriented cuboid y_min. + * @param[in] z_min Iso-oriented cuboid z_min. + * @param[in] x_max Iso-oriented cuboid x_max. + * @param[in] y_max Iso-oriented cuboid y_max. + * @param[in] z_max Iso-oriented cuboid z_max. + * + * @pre Available if Alpha_complex_3d is Weighted and Periodic. + * + * The type InputPointRange must be a range for which std::begin and + * std::end return input iterators on a `Alpha_complex_3d::Point_3`. + * The type WeightRange must be a range for which std::begin and + * std::end return an input iterator on a double. + * The type of x_min, y_min, z_min, x_max, y_max and z_max must be a double. + */ + template <typename InputPointRange, typename WeightRange> + 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) { + static_assert(Weighted, "This constructor is not available for non-weighted versions of Alpha_complex_3d"); + static_assert(Periodic, "This constructor is not available for non-periodic versions of Alpha_complex_3d"); + GUDHI_CHECK((weights.size() == points.size()), + std::invalid_argument("Points number in range different from weights range number")); + // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it. + GUDHI_CHECK( + (x_max - x_min == y_max - y_min) && (x_max - x_min == z_max - z_min) && (z_max - z_min == y_max - y_min), + std::invalid_argument("The size of the cuboid in every directions is not the same.")); + + std::vector<Weighted_point_3> weighted_points_3; + + std::size_t index = 0; + weighted_points_3.reserve(points.size()); + +#ifdef GUDHI_DEBUG + // Defined in GUDHI_DEBUG to avoid unused variable warning for GUDHI_CHECK + FT maximal_possible_weight = 0.015625 * (x_max - x_min) * (x_max - x_min); +#endif + + while ((index < weights.size()) && (index < points.size())) { + GUDHI_CHECK((weights[index] < maximal_possible_weight) && (weights[index] >= 0), + std::invalid_argument("Invalid weight at index " + std::to_string(index + 1) + + ". Must be positive and less than maximal possible weight = 1/64*cuboid length " + "squared, which is not an acceptable input.")); + weighted_points_3.push_back(Weighted_point_3(points[index], weights[index])); + index++; + } + + // Define the periodic cube + Dt pdt(typename Kernel::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max)); + // Heuristic for inserting large point sets (if pts is reasonably large) + pdt.insert(std::begin(weighted_points_3), std::end(weighted_points_3), true); + // As pdt won't be modified anymore switch to 1-sheeted cover if possible + if (!pdt.is_triangulation_in_1_sheet()) { + throw std::invalid_argument("Unable to construct a triangulation within a single periodic domain."); + } + pdt.convert_to_1_sheeted_covering(); + + // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode + // Maybe need to set it to GENERAL mode + alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>(new Alpha_shape_3(pdt, 0, Alpha_shape_3::GENERAL)); + } + + /** \brief Inserts all Delaunay triangulation into the simplicial complex. + * It also computes the filtration values accordingly to the \ref createcomplexalgorithm + * + * \tparam SimplicialComplexForAlpha3d must meet `SimplicialComplexForAlpha3d` concept. + * + * @param[in] complex SimplicialComplexForAlpha3d to be created. + * @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$, and there is very + * little point using anything else since it does not save time. + * + * @return true if creation succeeds, false otherwise. + * + * @pre The simplicial complex must be empty (no vertices). + * + */ + template <typename SimplicialComplexForAlpha3d, + typename Filtration_value = typename SimplicialComplexForAlpha3d::Filtration_value> + bool create_complex(SimplicialComplexForAlpha3d& complex, + Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity()) { + if (complex.num_vertices() > 0) { + std::cerr << "Alpha_complex_3d create_complex - complex is not empty\n"; + return false; // ----- >> + } + + // using Filtration_value = typename SimplicialComplexForAlpha3d::Filtration_value; + using Complex_vertex_handle = typename SimplicialComplexForAlpha3d::Vertex_handle; + using Alpha_shape_simplex_tree_map = std::unordered_map<Alpha_vertex_handle, Complex_vertex_handle>; + using Simplex_tree_vector_vertex = std::vector<Complex_vertex_handle>; + +#ifdef DEBUG_TRACES + std::size_t count_vertices = 0; + std::size_t count_edges = 0; + std::size_t count_facets = 0; + std::size_t count_cells = 0; +#endif // DEBUG_TRACES + std::vector<CGAL::Object> objects; + std::vector<FT> alpha_values; + + Dispatch dispatcher = CGAL::dispatch_output<CGAL::Object, FT>(std::back_inserter(objects), + std::back_inserter(alpha_values)); + + alpha_shape_3_ptr_->filtration_with_alpha_values(dispatcher); +#ifdef DEBUG_TRACES + std::cout << "filtration_with_alpha_values returns : " << objects.size() << " objects" << std::endl; +#endif // DEBUG_TRACES + + Alpha_shape_simplex_tree_map map_cgal_simplex_tree; + using Alpha_value_iterator = typename std::vector<FT>::const_iterator; + Alpha_value_iterator alpha_value_iterator = alpha_values.begin(); + for (auto object_iterator : objects) { + Vertex_list vertex_list; + + // Retrieve Alpha shape vertex list from object + if (const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { + for (auto i = 0; i < 4; i++) { +#ifdef DEBUG_TRACES + std::cout << "from cell[" << i << "]=" << (*cell)->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + vertex_list.push_back((*cell)->vertex(i)); + } +#ifdef DEBUG_TRACES + count_cells++; +#endif // DEBUG_TRACES + } else if (const Facet* facet = CGAL::object_cast<Facet>(&object_iterator)) { + for (auto i = 0; i < 4; i++) { + if ((*facet).second != i) { +#ifdef DEBUG_TRACES + std::cout << "from facet=[" << i << "]" << (*facet).first->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + vertex_list.push_back((*facet).first->vertex(i)); + } + } +#ifdef DEBUG_TRACES + count_facets++; +#endif // DEBUG_TRACES + } else if (const Edge* edge = CGAL::object_cast<Edge>(&object_iterator)) { + for (auto i : {(*edge).second, (*edge).third}) { +#ifdef DEBUG_TRACES + std::cout << "from edge[" << i << "]=" << (*edge).first->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + vertex_list.push_back((*edge).first->vertex(i)); + } +#ifdef DEBUG_TRACES + count_edges++; +#endif // DEBUG_TRACES + } else if (const Alpha_vertex_handle* vertex = CGAL::object_cast<Alpha_vertex_handle>(&object_iterator)) { +#ifdef DEBUG_TRACES + count_vertices++; + std::cout << "from vertex=" << (*vertex)->point() << std::endl; +#endif // DEBUG_TRACES + vertex_list.push_back((*vertex)); + } + // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex + Simplex_tree_vector_vertex the_simplex; + for (auto the_alpha_shape_vertex : vertex_list) { + auto the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex); + if (the_map_iterator == map_cgal_simplex_tree.end()) { + // alpha shape not found + Complex_vertex_handle vertex = map_cgal_simplex_tree.size(); +#ifdef DEBUG_TRACES + std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl; +#endif // DEBUG_TRACES + the_simplex.push_back(vertex); + map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex); + } else { + // alpha shape found + Complex_vertex_handle vertex = the_map_iterator->second; +#ifdef DEBUG_TRACES + std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl; +#endif // DEBUG_TRACES + the_simplex.push_back(vertex); + } + } + // Construction of the simplex_tree + Filtration_value filtr = Value_from_iterator<Complexity>::perform(alpha_value_iterator); + +#ifdef DEBUG_TRACES + std::cout << "filtration = " << filtr << std::endl; +#endif // DEBUG_TRACES + complex.insert_simplex(the_simplex, static_cast<Filtration_value>(filtr)); + GUDHI_CHECK(alpha_value_iterator != alpha_values.end(), "CGAL provided more simplices than values"); + ++alpha_value_iterator; + } + +#ifdef DEBUG_TRACES + std::cout << "vertices \t" << count_vertices << std::endl; + std::cout << "edges \t\t" << count_edges << std::endl; + std::cout << "facets \t\t" << count_facets << std::endl; + std::cout << "cells \t\t" << count_cells << std::endl; +#endif // DEBUG_TRACES + // -------------------------------------------------------------------------------------------- + // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension + complex.make_filtration_non_decreasing(); + // Remove all simplices that have a filtration value greater than max_alpha_square + complex.prune_above_filtration(max_alpha_square); + // -------------------------------------------------------------------------------------------- + return true; + } + + private: + // use of a unique_ptr on cgal Alpha_shape_3, as copy and default constructor is not available - no need to be freed + std::unique_ptr<Alpha_shape_3> alpha_shape_3_ptr_; +}; + +} // namespace alpha_complex + +} // namespace Gudhi + +#endif // ALPHA_COMPLEX_3D_H_ diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_options.h b/src/Alpha_complex/include/gudhi/Alpha_complex_options.h new file mode 100644 index 00000000..85c83672 --- /dev/null +++ b/src/Alpha_complex/include/gudhi/Alpha_complex_options.h @@ -0,0 +1,33 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2018 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef ALPHA_COMPLEX_OPTIONS_H_ +#define ALPHA_COMPLEX_OPTIONS_H_ + +namespace Gudhi { + +namespace alpha_complex { + +/** + * \brief Alpha complex complexity template parameter possible values. + * + * \ingroup alpha_complex + */ +enum class complexity : char { + FAST = 'f', ///< Fast version. + SAFE = 's', ///< Safe version. + EXACT = 'e', ///< Exact version. +}; + +} // namespace alpha_complex + +} // namespace Gudhi + +#endif // ALPHA_COMPLEX_OPTIONS_H_ |