diff options
Diffstat (limited to 'src/Alpha_complex/include/gudhi')
| -rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_complex.h | 75 | ||||
| -rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_complex_3d.h | 81 |
2 files changed, 104 insertions, 52 deletions
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 8919cdb9..f2a05e95 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -20,11 +20,12 @@ #include <math.h> // isnan, fmax #include <CGAL/Delaunay_triangulation.h> -#include <CGAL/Epick_d.h> +#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> #include <CGAL/property_map.h> // for CGAL::Identity_property_map -#include <CGAL/NT_converter.h> #include <CGAL/version.h> // for CGAL_VERSION_NR +#include <CGAL/NT_converter.h> #include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST @@ -39,17 +40,20 @@ // 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 +# error Alpha_complex 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 +# error Alpha_complex is only available for Eigen3 >= 3.1.0 installed with CGAL #endif namespace Gudhi { namespace alpha_complex { +template<typename D> struct Is_Epeck_D { static const bool value = false; }; +template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool value = true; }; + /** * \class Alpha_complex Alpha_complex.h gudhi/Alpha_complex.h * \brief Alpha complex data structure. @@ -63,17 +67,31 @@ namespace alpha_complex { * * Please refer to \ref alpha_complex for examples. * - * The complex is a template class requiring an Epick_d <a target="_blank" + * The complex is a template class requiring an <a target="_blank" + * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>, + * or an <a target="_blank" + * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epick__d.html">CGAL::Epick_d</a> <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> + * \cite cgal:s-gkd-19b from CGAL as template, default value is <a target="_blank" + * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_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. - * + * \remark + * - When Alpha_complex is constructed with an infinite value of alpha, the complex is a Delaunay complex. + * - Using the default `CGAL::Epeck_d` makes the construction safe. If you pass exact=true to create_complex, the + * filtration values are the exact ones converted to the filtration value type of the simplicial complex. This can be + * very slow. If you pass exact=false (the default), the filtration values are only guaranteed to have a small + * multiplicative error compared to the exact value, see <code><a class="el" target="_blank" + * href="https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html"> + * CGAL::Lazy_exact_nt<NT>::set_relative_precision_of_to_double</a></code> for details. A drawback, when computing + * persistence, is that an empty exact interval [10^12,10^12] may become a non-empty approximate interval + * [10^12,10^12+10^6]. Using `CGAL::Epick_d` makes the computations slightly faster, and the combinatorics are still + * exact, but the computation of filtration values can exceptionally be arbitrarily bad. In all cases, we still + * 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::Epick_d<CGAL::Dynamic_dimension_tag>> +template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>> class Alpha_complex { public: // Add an int in TDS to save point index in the structure @@ -103,8 +121,8 @@ class Alpha_complex { // 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; + // Structure to switch from simplex tree vertex handle to CGAL vertex iterator. + typedef typename std::vector< CGAL_vertex_iterator > Vector_vertex_iterator; private: /** \brief Vertex iterator vector to switch from simplex tree vertex handle to CGAL vertex iterator. @@ -173,17 +191,15 @@ class Alpha_complex { 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) { + #if CGAL_VERSION_NR < 1050000000 + if (Is_Epeck_D<Kernel>::value) + std::cerr << "It is strongly advised to use a CGAL version >= 5.0 with Epeck_d Kernel for performance reasons." + << std::endl; + #endif + auto first = std::begin(points); auto last = std::end(points); @@ -214,14 +230,16 @@ class Alpha_complex { hint = pos->full_cell(); } // -------------------------------------------------------------------------------------------- - // double map to retrieve simplex tree vertex handles from CGAL vertex iterator and vice versa + // 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()); // 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); + vertex_handle_to_iterator_[vit->data()] = vit; } } // -------------------------------------------------------------------------------------------- @@ -237,6 +255,8 @@ class Alpha_complex { * @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. + * @param[in] exact Exact filtration values computation. Not exact if `Kernel` is not <a target="_blank" + * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>. * * @return true if creation succeeds, false otherwise. * @@ -248,7 +268,8 @@ class Alpha_complex { 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()) { + Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(), + bool exact = 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; @@ -324,9 +345,13 @@ class Alpha_complex { 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())); + CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv; + auto sqrad = squared_radius(pointVector.begin(), pointVector.end()); +#if CGAL_VERSION_NR >= 1050000000 + if(exact) CGAL::exact(sqrad); +#endif + alpha_complex_filtration = cv(sqrad); } complex.assign_filtration(f_simplex, alpha_complex_filtration); #ifdef DEBUG_TRACES diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h index 13ebb9c1..7f96c94c 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h @@ -43,7 +43,7 @@ #include <vector> #include <unordered_map> #include <stdexcept> -#include <cstddef> +#include <cstddef> // for std::size_t #include <memory> // for std::unique_ptr #include <type_traits> // for std::conditional and std::enable_if #include <limits> // for numeric_limits<> @@ -97,7 +97,7 @@ struct Value_from_iterator<complexity::EXACT> { * \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. + * Duplicate points are inserted once in the Alpha_complex. * * \tparam Complexity shall be `Gudhi::alpha_complex::complexity` type. Default value is * `Gudhi::alpha_complex::complexity::SAFE`. @@ -225,23 +225,23 @@ class Alpha_complex_3d { * 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: + /** \brief Gives public access to the Bare_point_3 (bare aka. unweighed) type. + * Here is a Bare_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.); +Alpha_complex_3d::Bare_point_3 p0(1., -1.1, -1.); \endcode * */ - using Point_3 = typename Kernel::Point_3; + using Bare_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>; +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.); +Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_point_3(1., -1.1, -1.), 4.); \endcode * * Note: This type is defined to void if Alpha complex is not weighted. @@ -249,6 +249,11 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point * */ using Weighted_point_3 = typename Triangulation_3<Kernel, Tds, Weighted, Periodic>::Weighted_point_3; + /** \brief `Alpha_complex_3d::Point_3` type is either a `Alpha_complex_3d::Bare_point_3` (Weighted = false) or a + * `Alpha_complex_3d::Weighted_point_3` (Weighted = true). + */ + using Point_3 = typename Alpha_shape_3::Point; + private: using Dispatch = CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<CGAL::Object, FT>, @@ -264,13 +269,12 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point 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`. + * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::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`. + * `Alpha_complex_3d::Point_3`. */ template <typename InputPointRange> Alpha_complex_3d(const InputPointRange& points) { @@ -284,13 +288,13 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point * * @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] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Bare_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`. + * std::end return input iterators on a `Alpha_complex_3d::Bare_point_3`. * The type WeightRange must be a range for which std::begin and * std::end return an input iterator on a double. */ @@ -318,8 +322,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point * * @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] points Range of points to triangulate. Points must be in `Alpha_complex_3d::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. @@ -330,7 +333,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point * @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`. + * `Alpha_complex_3d::Point_3`. * * @note In weighted version, please check weights are greater than zero, and lower than 1/64*cuboid length * squared. @@ -366,7 +369,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point * @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] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Bare_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. @@ -378,7 +381,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point * @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`. + * std::end return input iterators on a `Alpha_complex_3d::Bare_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. @@ -452,9 +455,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point 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 @@ -474,7 +475,6 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point 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) { @@ -484,7 +484,8 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point 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; + std::cout << "from cell[" << i << "] - Point coordinates (" << (*cell)->vertex(i)->point() << ")" + << std::endl; #endif // DEBUG_TRACES vertex_list.push_back((*cell)->vertex(i)); } @@ -495,7 +496,8 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point 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; + std::cout << "from facet=[" << i << "] - Point coordinates (" << (*facet).first->vertex(i)->point() << ")" + << std::endl; #endif // DEBUG_TRACES vertex_list.push_back((*facet).first->vertex(i)); } @@ -506,7 +508,8 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point } 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; + std::cout << "from edge[" << i << "] - Point coordinates (" << (*edge).first->vertex(i)->point() << ")" + << std::endl; #endif // DEBUG_TRACES vertex_list.push_back((*edge).first->vertex(i)); } @@ -516,7 +519,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point } 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; + std::cout << "from vertex - Point coordinates (" << (*vertex)->point() << ")" << std::endl; #endif // DEBUG_TRACES vertex_list.push_back((*vertex)); } @@ -528,7 +531,8 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point // 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; + std::cout << "Point (" << the_alpha_shape_vertex->point() << ") not found - insert new vertex id " << vertex + << std::endl; #endif // DEBUG_TRACES the_simplex.push_back(vertex); map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex); @@ -536,7 +540,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point // 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; + std::cout << "Point (" << the_alpha_shape_vertex->point() << ") found as vertex id " << vertex << std::endl; #endif // DEBUG_TRACES the_simplex.push_back(vertex); } @@ -567,9 +571,32 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point return true; } + /** \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_3& get_point(std::size_t vertex) { + if (map_cgal_simplex_tree.size() != cgal_vertex_iterator_vector.size()) { + cgal_vertex_iterator_vector.resize(map_cgal_simplex_tree.size()); + for (auto map_iterator : map_cgal_simplex_tree) { + cgal_vertex_iterator_vector[map_iterator.second] = map_iterator.first; + } + + } + auto cgal_vertex_iterator = cgal_vertex_iterator_vector.at(vertex); + return cgal_vertex_iterator->point(); + } + 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_; + + // Map type to switch from CGAL vertex iterator to simplex tree vertex handle. + std::unordered_map<Alpha_vertex_handle, std::size_t> map_cgal_simplex_tree; + // Vector type to switch from simplex tree vertex handle to CGAL vertex iterator. + std::vector<Alpha_vertex_handle> cgal_vertex_iterator_vector; }; } // namespace alpha_complex |