diff options
Diffstat (limited to 'src/Alpha_complex/include/gudhi/Alpha_complex.h')
| -rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_complex.h | 227 |
1 files changed, 141 insertions, 86 deletions
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index ba91998d..a7372f19 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -12,14 +12,17 @@ #ifndef ALPHA_COMPLEX_H_ #define ALPHA_COMPLEX_H_ +#include <gudhi/Alpha_complex/Alpha_kernel_d.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 <cmath> // isnan, fmax +#include <memory> // for std::unique_ptr +#include <cstddef> // for std::size_t #include <CGAL/Delaunay_triangulation.h> +#include <CGAL/Regular_triangulation.h> // aka. Weighted Delaunay triangulation #include <CGAL/Epeck_d.h> // For EXACT or SAFE version #include <CGAL/Epick_d.h> // For FAST version #include <CGAL/Spatial_sort_traits_adapter_d.h> @@ -29,6 +32,10 @@ #include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST +#include <boost/range/size.hpp> +#include <boost/range/combine.hpp> +#include <boost/range/adaptor/transformed.hpp> + #include <iostream> #include <vector> #include <string> @@ -37,6 +44,7 @@ #include <utility> // std::pair #include <stdexcept> #include <numeric> // for std::iota +#include <algorithm> // for std::sort // Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 #if CGAL_VERSION_NR < 1041101000 @@ -61,7 +69,7 @@ template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool val * \ingroup alpha_complex * * \details - * The data structure is constructing a CGAL Delaunay triangulation (for more informations on CGAL Delaunay + * The data structure is constructing a CGAL Delaunay triangulation (for more information 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). * @@ -91,49 +99,61 @@ template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool val * guarantee that the output is a valid filtration (faces have a filtration value no larger than their cofaces). * - For performances reasons, it is advised to use `Alpha_complex` with \ref cgal ≥ 5.0.0. */ -template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>> +template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>, bool Weighted = false> class Alpha_complex { + private: + // Vertex_handle internal type (required by triangulation_ and vertices_). + using Internal_vertex_handle = std::ptrdiff_t; + public: + /** \brief Geometric traits class that provides the geometric types and predicates needed by the triangulations.*/ + using Geom_traits = std::conditional_t<Weighted, CGAL::Regular_triangulation_traits_adapter<Kernel>, Kernel>; + // 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; + using TDS = CGAL::Triangulation_data_structure<typename Geom_traits::Dimension, + CGAL::Triangulation_vertex<Geom_traits, Internal_vertex_handle>, + CGAL::Triangulation_full_cell<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; + /** \brief A (Weighted or not) Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/ + using Triangulation = std::conditional_t<Weighted, CGAL::Regular_triangulation<Kernel, TDS>, + CGAL::Delaunay_triangulation<Kernel, TDS>>; - // Vertex_iterator type from CGAL. - typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator; + /** \brief CGAL kernel container for computations in function of the weighted or not characteristics.*/ + using A_kernel_d = Alpha_kernel_d<Kernel, Weighted>; - // size_type type from CGAL. - typedef typename Delaunay_triangulation::size_type size_type; + // Numeric type of coordinates in the kernel + using FT = typename A_kernel_d::FT; - // Structure to switch from simplex tree vertex handle to CGAL vertex iterator. - typedef typename std::vector< CGAL_vertex_iterator > Vector_vertex_iterator; + /** \brief Sphere is a std::pair<Kernel::Point_d, Kernel::FT> (aka. circurmcenter and squared radius). + * If Weighted, Sphere is a Kernel::Weighted_point_d (aka. circurmcenter and the weight value is the squared radius). + */ + using Sphere = typename A_kernel_d::Sphere; - // Numeric type of coordinates in the kernel - typedef typename Kernel::FT FT; + /** \brief A point, or a weighted point in Euclidean space.*/ + using Point_d = typename Geom_traits::Point_d; + + private: + // Vertex_iterator type from CGAL. + using CGAL_vertex_iterator = typename Triangulation::Vertex_iterator; + + // Structure to switch from simplex tree vertex handle to CGAL vertex iterator. + using Vector_vertex_iterator = std::vector< CGAL_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_; + std::unique_ptr<Triangulation> triangulation_; /** \brief Kernel for triangulation_ functions access.*/ - Kernel kernel_; + A_kernel_d kernel_; + /** \brief Vertices to be inserted first by the create_complex method to avoid quadratic complexity. + * It isn't just [0, n) if some points have multiplicity (only one copy appears in the complex). + */ + std::vector<Internal_vertex_handle> vertices_; + /** \brief Cache for geometric constructions: circumcenter and squared radius of a simplex.*/ - std::vector<std::pair<Point_d, FT>> cache_, old_cache_; + std::vector<Sphere> cache_, old_cache_; public: /** \brief Alpha_complex constructor from an OFF file name. @@ -145,8 +165,7 @@ class Alpha_complex { * * @param[in] off_file_name OFF file [path and] name. */ - Alpha_complex(const std::string& off_file_name) - : triangulation_(nullptr) { + Alpha_complex(const std::string& off_file_name) { 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"; @@ -158,23 +177,40 @@ class Alpha_complex { /** \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. + * The vertices may be not contiguous as some points may be discarded in the triangulation (duplicate points, + * weighted hidden point, ...). * - * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d + * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d or Kernel::Weighted_point_d. * - * The type InputPointRange must be a range for which std::begin and - * std::end return input iterators on a 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 or Kernel::Weighted_point_d. */ template<typename InputPointRange > - Alpha_complex(const InputPointRange& points) - : triangulation_(nullptr) { + Alpha_complex(const InputPointRange& points) { init_from_range(points); } - /** \brief Alpha_complex destructor deletes the Delaunay triangulation. + /** \brief Alpha_complex constructor from a list of points and weights. + * + * The vertices may be not contiguous as some points may be discarded in the triangulation (duplicate points, + * weighted hidden point, ...). + * + * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d. + * + * @param[in] weights Range of points weights. Weights must be in Kernel::FT. + * + * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a + * Kernel::Point_d. */ - ~Alpha_complex() { - delete triangulation_; + template <typename InputPointRange, typename WeightRange> + Alpha_complex(const InputPointRange& points, WeightRange weights) { + static_assert(Weighted, "This constructor is not available for non-weighted versions of Alpha_complex"); + // FIXME: this test is only valid if we have a forward range + GUDHI_CHECK(boost::size(weights) == boost::size(points), + std::invalid_argument("Points number in range different from weights range number")); + auto weighted_points = boost::range::combine(points, weights) + | boost::adaptors::transformed([](auto const&t){return Point_d(boost::get<0>(t), boost::get<1>(t));}); + init_from_range(weighted_points); } // Forbid copy/move constructor/assignment operator @@ -183,6 +219,15 @@ class Alpha_complex { Alpha_complex (Alpha_complex&& other) = delete; Alpha_complex& operator= (Alpha_complex&& other) = delete; + /** \brief Returns the number of finite vertices in the triangulation. + */ + std::size_t num_vertices() const { + if (triangulation_ == nullptr) + return 0; + else + return triangulation_->number_of_vertices(); + } + /** \brief get_point returns the point corresponding to the vertex given as parameter. * * @param[in] vertex Vertex handle of the point to retrieve. @@ -202,39 +247,46 @@ class Alpha_complex { << std::endl; #endif +#if CGAL_VERSION_NR < 1050101000 + // Make compilation fail if weighted and CGAL < 5.1 + static_assert(!Weighted, "Weighted Alpha_complex is only available for CGAL >= 5.1"); +#endif + 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)); + // Delaunay triangulation init with point dimension. + triangulation_ = std::make_unique<Triangulation>(kernel_.get_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())); + std::vector<Internal_vertex_handle> indices(boost::counting_iterator<Internal_vertex_handle>(0), + boost::counting_iterator<Internal_vertex_handle>(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; + using Point_property_map = boost::iterator_property_map<typename std::vector<Point_d>::iterator, + CGAL::Identity_property_map<Internal_vertex_handle>>; + using Search_traits_d = CGAL::Spatial_sort_traits_adapter_d<Geom_traits, Point_property_map>; CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud))); - typename Delaunay_triangulation::Full_cell_handle hint; + typename 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(); + typename Triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint); + if (pos != nullptr) { + // Save index value as data to retrieve it after insertion + pos->data() = index; + hint = pos->full_cell(); + } } // -------------------------------------------------------------------------------------------- // structure to retrieve CGAL points from vertex handle - one vertex handle per point. // Needs to be constructed before as vertex handles arrives in no particular order. vertex_handle_to_iterator_.resize(point_cloud.size()); + // List of sorted unique vertices in the triangulation. We take advantage of the existing loop to construct it + // Vertices list avoids quadratic complexity with the Simplex_tree. We should not fill it up with Toplex_map e.g. + vertices_.reserve(triangulation_->number_of_vertices()); // Loop on triangulation vertices list for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) { if (!triangulation_->is_infinite(*vit)) { @@ -242,8 +294,10 @@ class Alpha_complex { std::clog << "Vertex insertion - " << vit->data() << " -> " << vit->point() << std::endl; #endif // DEBUG_TRACES vertex_handle_to_iterator_[vit->data()] = vit; + vertices_.push_back(vit->data()); } } + std::sort(vertices_.begin(), vertices_.end()); // -------------------------------------------------------------------------------------------- } } @@ -270,9 +324,7 @@ class Alpha_complex { v.clear(); for (auto vertex : cplx.simplex_vertex_range(s)) v.push_back(get_point_(vertex)); - Point_d c = kernel_.construct_circumcenter_d_object()(v.cbegin(), v.cend()); - FT r = kernel_.squared_distance_d_object()(c, v[0]); - cache_.emplace_back(std::move(c), std::move(r)); + cache_.emplace_back(kernel_.get_sphere(v.cbegin(), v.cend())); } return cache_[k]; } @@ -282,13 +334,13 @@ class Alpha_complex { auto radius(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) { auto k = cplx.key(s); if(k!=cplx.null_key()) - return old_cache_[k].second; + return kernel_.get_squared_radius(old_cache_[k]); // Using a transform_range is slower, currently. thread_local std::vector<Point_d> v; v.clear(); for (auto vertex : cplx.simplex_vertex_range(s)) v.push_back(get_point_(vertex)); - return kernel_.compute_squared_radius_d_object()(v.cbegin(), v.cend()); + return kernel_.get_squared_radius(v.cbegin(), v.cend()); } public: @@ -322,9 +374,9 @@ class Alpha_complex { bool exact = false, bool default_filtration_value = false) { // 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; + using Vertex_handle = typename SimplicialComplexForAlpha::Vertex_handle; + using Simplex_handle = typename SimplicialComplexForAlpha::Simplex_handle; + using Vector_vertex = std::vector<Vertex_handle>; if (triangulation_ == nullptr) { std::cerr << "Alpha_complex cannot create_complex from a NULL triangulation\n"; @@ -341,13 +393,22 @@ class Alpha_complex { // -------------------------------------------------------------------------------------------- // Simplex_tree construction from loop on triangulation finite full cells list - if (triangulation_->number_of_vertices() > 0) { + if (num_vertices() > 0) { + std::vector<Vertex_handle> one_vertex(1); + for (auto vertex : vertices_) { +#ifdef DEBUG_TRACES + std::clog << "SimplicialComplex insertion " << vertex << std::endl; +#endif // DEBUG_TRACES + one_vertex[0] = vertex; + complex.insert_simplex_and_subfaces(one_vertex, std::numeric_limits<double>::quiet_NaN()); + } + for (auto cit = triangulation_->finite_full_cells_begin(); cit != triangulation_->finite_full_cells_end(); ++cit) { Vector_vertex vertexVector; #ifdef DEBUG_TRACES - std::clog << "Simplex_tree insertion "; + std::clog << "SimplicialComplex insertion "; #endif // DEBUG_TRACES for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) { if (*vit != nullptr) { @@ -368,6 +429,7 @@ class Alpha_complex { // -------------------------------------------------------------------------------------------- if (!default_filtration_value) { + CGAL::NT_converter<FT, Filtration_value> cgal_converter; // -------------------------------------------------------------------------------------------- // ### For i : d -> 0 for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) { @@ -378,14 +440,13 @@ class Alpha_complex { // ### 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) { + // No need to compute squared_radius on a non-weighted single point - alpha is 0.0 + if (Weighted || f_simplex_dim > 0) { auto const& sqrad = radius(complex, f_simplex); #if CGAL_VERSION_NR >= 1050000000 if(exact) CGAL::exact(sqrad); #endif - CGAL::NT_converter<FT, Filtration_value> cv; - alpha_complex_filtration = cv(sqrad); + alpha_complex_filtration = cgal_converter(sqrad); } complex.assign_filtration(f_simplex, alpha_complex_filtration); #ifdef DEBUG_TRACES @@ -393,7 +454,7 @@ class Alpha_complex { #endif // DEBUG_TRACES } // No need to propagate further, unweighted points all have value 0 - if (decr_dim > 1) + if (decr_dim > !Weighted) propagate_alpha_filtration(complex, f_simplex); } } @@ -403,8 +464,10 @@ class Alpha_complex { // -------------------------------------------------------------------------------------------- // -------------------------------------------------------------------------------------------- - // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension - complex.make_filtration_non_decreasing(); + if (!exact) + // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension + // Only in not exact version, cf. https://github.com/GUDHI/gudhi-devel/issues/57 + 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); // -------------------------------------------------------------------------------------------- @@ -416,11 +479,11 @@ class Alpha_complex { 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; - typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle; + using Filtration_value = typename SimplicialComplexForAlpha::Filtration_value; // ### Foreach Tau face of Sigma - for (auto f_boundary : complex.boundary_simplex_range(f_simplex)) { + for (auto face_opposite_vertex : complex.boundary_opposite_vertex_simplex_range(f_simplex)) { + auto f_boundary = face_opposite_vertex.first; #ifdef DEBUG_TRACES std::clog << " | --------------------------------------------------\n"; std::clog << " | Tau "; @@ -441,18 +504,10 @@ class Alpha_complex { #endif // DEBUG_TRACES // ### Else } else { - // Find which vertex of f_simplex is missing in f_boundary. We could actually write a variant of boundary_simplex_range that gives pairs (f_boundary, vertex). We rely on the fact that simplex_vertex_range is sorted. - auto longlist = complex.simplex_vertex_range(f_simplex); - auto shortlist = complex.simplex_vertex_range(f_boundary); - auto longiter = std::begin(longlist); - auto shortiter = std::begin(shortlist); - auto enditer = std::end(shortlist); - while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; } - Vertex_handle extra = *longiter; auto const& cache=get_cache(complex, f_boundary); - bool is_gab = kernel_.squared_distance_d_object()(cache.first, get_point_(extra)) >= cache.second; + bool is_gab = kernel_.is_gabriel(cache, get_point_(face_opposite_vertex.second)); #ifdef DEBUG_TRACES - std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << extra << std::endl; + std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << face_opposite_vertex.second << std::endl; #endif // DEBUG_TRACES // ### If Tau is not Gabriel of Sigma if (false == is_gab) { |