diff options
author | ROUVREAU Vincent <vincent.rouvreau@inria.fr> | 2020-09-04 17:23:16 +0200 |
---|---|---|
committer | ROUVREAU Vincent <vincent.rouvreau@inria.fr> | 2020-09-04 17:23:16 +0200 |
commit | db938dbd74460e7a0fd705be8628984052f71dc0 (patch) | |
tree | 5d3e7cb24f0bcb35229145efc0093bd105bed6c8 /src/Alpha_complex/include/gudhi/Alpha_complex.h | |
parent | 644334d6bf42f8dcc750fe1b1b5ea0eebdbba167 (diff) |
Alpha_kernel_d, its tests and Alpha complex to use it
Diffstat (limited to 'src/Alpha_complex/include/gudhi/Alpha_complex.h')
-rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_complex.h | 106 |
1 files changed, 57 insertions, 49 deletions
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index ba91998d..8e9fe773 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -12,6 +12,7 @@ #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> @@ -20,6 +21,7 @@ #include <math.h> // isnan, fmax #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> @@ -91,49 +93,55 @@ 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 { public: + /** \brief Geometric traits class that provides the geometric types and predicates needed by the triangulations.*/ + using Geom_traits = typename std::conditional<Weighted, CGAL::Regular_triangulation_traits_adapter<Kernel>, + Kernel>::type; // 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; + using TDS = CGAL::Triangulation_data_structure<typename Geom_traits::Dimension, + CGAL::Triangulation_vertex<Geom_traits, std::ptrdiff_t>, + CGAL::Triangulation_full_cell<Geom_traits> >; + /** \brief A (Weighted or not) Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/ + using Triangulation = typename std::conditional<Weighted, CGAL::Regular_triangulation<Kernel, TDS>, + CGAL::Delaunay_triangulation<Kernel, TDS>>::type; + + using A_kernel_d = Alpha_kernel_d<Kernel, Weighted>; + + // Numeric type of coordinates in the kernel + using FT = typename A_kernel_d::FT; + + /** \brief If Weighted, the weighted point is cached (point + weight [= squared radius]), + * else a pair of point and squared radius is cached. + */ + using Sphere = typename A_kernel_d::Sphere; /** \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 Point_d = typename std::conditional<Weighted, typename A_kernel_d::Weighted_point_d, + typename A_kernel_d::Point_d>::type; 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; - // Vertex_iterator type from CGAL. - typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator; + using CGAL_vertex_iterator = typename Triangulation::Vertex_iterator; // size_type type from CGAL. - typedef typename Delaunay_triangulation::size_type size_type; + using size_type = typename Triangulation::size_type; // Structure to switch from simplex tree vertex handle to CGAL vertex iterator. - typedef typename std::vector< CGAL_vertex_iterator > Vector_vertex_iterator; - - // Numeric type of coordinates in the kernel - typedef typename Kernel::FT FT; + using Vector_vertex_iterator = typename 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_; + Triangulation* triangulation_; /** \brief Kernel for triangulation_ functions access.*/ - Kernel kernel_; + A_kernel_d kernel_; + /** \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. @@ -160,10 +168,10 @@ class Alpha_complex { * * 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 + * @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. + * std::end return input iterators on a Kernel::Point_d or Kernel::Weighted_point_d. */ template<typename InputPointRange > Alpha_complex(const InputPointRange& points) @@ -206,11 +214,8 @@ class Alpha_complex { 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_ = new Triangulation(kernel_.get_dimension(*first)); std::vector<Point_d> point_cloud(first, last); @@ -218,18 +223,22 @@ class Alpha_complex { 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; + using Point_property_map = boost::iterator_property_map<typename std::vector<Point_d>::iterator, + CGAL::Identity_property_map<std::ptrdiff_t>>; + 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(); + } else { + std::cout << "NULLPTR" << std::endl; + } } // -------------------------------------------------------------------------------------------- // structure to retrieve CGAL points from vertex handle - one vertex handle per point. @@ -270,9 +279,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 +289,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 +329,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"; @@ -416,8 +423,8 @@ 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; + using Vertex_handle = typename SimplicialComplexForAlpha::Vertex_handle; // ### Foreach Tau face of Sigma for (auto f_boundary : complex.boundary_simplex_range(f_simplex)) { @@ -450,7 +457,8 @@ class Alpha_complex { 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_.get_squared_distance(kernel_.get_circumcenter(cache), get_point_(extra)) >= + kernel_.get_squared_radius(cache); #ifdef DEBUG_TRACES std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << extra << std::endl; #endif // DEBUG_TRACES |