diff options
Diffstat (limited to 'src/Alpha_complex')
19 files changed, 1134 insertions, 337 deletions
diff --git a/src/Alpha_complex/doc/Intro_alpha_complex.h b/src/Alpha_complex/doc/Intro_alpha_complex.h index 60da7169..41e5e16d 100644 --- a/src/Alpha_complex/doc/Intro_alpha_complex.h +++ b/src/Alpha_complex/doc/Intro_alpha_complex.h @@ -22,6 +22,18 @@ namespace alpha_complex { * * @{ * +<div class="toc"> +Table of Contents +<ul> +<li class="level1"><a href="#definition">Definition</a></li> +<li class="level1"><a href="#pointsexample">Example from points</a></li> +<li class="level1"><a href="#createcomplexalgorithm">Create complex algorithm</a></li> +<li class="level1"><a href="#weightedversion">Weighted specific version</a></li> +<li class="level1"><a href="#offexample">Example from OFF file</a></li> +<li class="level1"><a href="#weighted3dexample">3d specific version</a></li> +</ul> +</div> + * \section definition Definition * * Alpha_complex is a <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a> @@ -62,7 +74,7 @@ namespace alpha_complex { * [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. + * - For performances reasons, it is advised to use \ref eigen ≥ 3.3.5 and \ref cgal ≥ 5.2.0. * * \section pointsexample Example from points * @@ -71,7 +83,7 @@ namespace alpha_complex { * * Then, it is asked to display information about the simplicial complex. * - * \include Alpha_complex/Alpha_complex_from_points.cpp + * \include Alpha_complex_from_points.cpp * * When launching: * @@ -80,7 +92,7 @@ namespace alpha_complex { * * the program output is: * - * \include Alpha_complex/alphaoffreader_for_doc_60.txt + * \include alphaoffreader_for_doc_60.txt * * \section createcomplexalgorithm Create complex algorithm * @@ -95,6 +107,7 @@ namespace alpha_complex { * \subsection filtrationcomputation Filtration value computation algorithm * <br> * \f$ + * \begin{array}{l} * \textbf{for } \text{i : dimension } \rightarrow 0 \textbf{ do}\\ * \quad \textbf{for all } \sigma \text{ of dimension i}\\ * \quad\quad \textbf{if } \text{filtration(} \sigma ) \text{ is NaN} \textbf{ then}\\ @@ -115,6 +128,7 @@ namespace alpha_complex { * \textbf{end for}\\ * \text{make_filtration_non_decreasing()}\\ * \text{prune_above_filtration()}\\ + * \end{array} * \f$ * * \subsubsection dimension2 Dimension 2 @@ -140,12 +154,38 @@ namespace alpha_complex { * not quite define a proper filtration (i.e. non-decreasing with respect to inclusion). * We fix that up by calling `SimplicialComplexForAlpha::make_filtration_non_decreasing()`. * + * \note This is not the case in `exact` version, this is the reason why it is not called in this case. + * * \subsubsection pruneabove Prune above given filtration value * * The simplex tree is pruned from the given maximum \f$ \alpha^2 \f$ value (cf. * `SimplicialComplexForAlpha::prune_above_filtration()`). * In the following example, the value is given by the user as argument of the program. * + * \section weightedversion Weighted specific version + * <b>Requires:</b> \ref eigen ≥ 3.1.0 and \ref cgal ≥ 5.1.0. + * + * A weighted version for Alpha complex is available (cf. Alpha_complex). It is like a usual Alpha complex, but based + * on a <a href="https://doc.cgal.org/latest/Triangulation/index.html#TriangulationSecRT">CGAL regular triangulation</a> instead + * of Delaunay. + * + * This example builds the CGAL weighted alpha shapes from a small molecule, and initializes the alpha complex with + * it. This example is taken from <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#AlphaShape_3DExampleforWeightedAlphaShapes">CGAL 3d + * weighted alpha shapes</a>. + * + * Then, it is asked to display information about the alpha complex. + * + * \include Weighted_alpha_complex_from_points.cpp + * + * When launching: + * + * \code $> ./Weighted_alpha_complex_example_from_points + * \endcode + * + * the program output is: + * + * \include weightedalpha3dfrompoints_for_doc.txt + * * * \section offexample Example from OFF file * @@ -154,7 +194,7 @@ namespace alpha_complex { * * Then, it is asked to display information about the alpha complex. * - * \include Alpha_complex/Alpha_complex_from_off.cpp + * \include Alpha_complex_from_off.cpp * * When launching: * @@ -163,10 +203,10 @@ namespace alpha_complex { * * the program output is: * - * \include Alpha_complex/alphaoffreader_for_doc_32.txt + * \include alphaoffreader_for_doc_32.txt * * - * \section weighted3dexample 3d specific example + * \section weighted3dexample 3d specific version * * A specific module for Alpha complex is available in 3d (cf. Alpha_complex_3d) and allows to construct standard, * weighted, periodic or weighted and periodic versions of alpha complexes. Alpha values computation can be @@ -174,21 +214,14 @@ namespace alpha_complex { * Gudhi::alpha_complex::complexity::EXACT. * * This example builds the CGAL 3d weighted alpha shapes from a small molecule, and initializes the alpha complex with - * it. This example is taken from <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title13">CGAL 3d + * it. This example is taken from <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#AlphaShape_3DExampleforWeightedAlphaShapes">CGAL 3d * weighted alpha shapes</a>. * * Then, it is asked to display information about the alpha complex. * - * \include Alpha_complex/Weighted_alpha_complex_3d_from_points.cpp - * - * When launching: - * - * \code $> ./Alpha_complex_example_weighted_3d_from_points - * \endcode - * - * the program output is: + * \include Weighted_alpha_complex_3d_from_points.cpp * - * \include Alpha_complex/weightedalpha3dfrompoints_for_doc.txt + * The results will be the same as in \ref weightedversion . * */ /** @} */ // end defgroup alpha_complex diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt index 2eecd50c..1fc2330a 100644 --- a/src/Alpha_complex/example/CMakeLists.txt +++ b/src/Alpha_complex/example/CMakeLists.txt @@ -25,7 +25,7 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) add_test(NAME Alpha_complex_example_fast_from_off_32 COMMAND $<TARGET_FILE:Alpha_complex_example_fast_from_off> "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "32.0" "${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_32.txt") -if (DIFF_PATH) + if (DIFF_PATH) # Do not forget to copy test results files in current binary dir file(COPY "alphaoffreader_for_doc_32.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) file(COPY "alphaoffreader_for_doc_60.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) @@ -43,8 +43,10 @@ if (DIFF_PATH) add_test(Alpha_complex_example_fast_from_off_32_diff_files ${DIFF_PATH} ${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt) set_tests_properties(Alpha_complex_example_fast_from_off_32_diff_files PROPERTIES DEPENDS Alpha_complex_example_fast_from_off_32) -endif() + endif() +endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) +if (NOT CGAL_VERSION VERSION_LESS 4.11.0) add_executable ( Alpha_complex_example_weighted_3d_from_points Weighted_alpha_complex_3d_from_points.cpp ) target_link_libraries(Alpha_complex_example_weighted_3d_from_points ${CGAL_LIBRARY}) if (TBB_FOUND) @@ -61,4 +63,13 @@ endif() add_test(NAME Alpha_complex_example_3d_from_points COMMAND $<TARGET_FILE:Alpha_complex_example_3d_from_points>) -endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) +endif(NOT CGAL_VERSION VERSION_LESS 4.11.0) + +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) + add_executable ( Weighted_alpha_complex_example_from_points Weighted_alpha_complex_from_points.cpp ) + target_link_libraries(Weighted_alpha_complex_example_from_points ${CGAL_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Weighted_alpha_complex_example_from_points ${TBB_LIBRARIES}) + endif() + add_test(NAME Weighted_alpha_complex_example_from_points COMMAND $<TARGET_FILE:Weighted_alpha_complex_example_from_points>) +endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) diff --git a/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp b/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp index c044194e..ee12d418 100644 --- a/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp +++ b/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp @@ -7,7 +7,7 @@ #include <vector> #include <limits> // for numeric limits -// Complexity = FAST, weighted = true, periodic = false +// Complexity = SAFE, weighted = true, periodic = false using Weighted_alpha_complex_3d = Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, false>; using Bare_point = Weighted_alpha_complex_3d::Bare_point_3; @@ -18,11 +18,11 @@ int main(int argc, char **argv) { // Init of a list of points and weights from a small molecule // ---------------------------------------------------------------------------- std::vector<Weighted_point> weighted_points; - weighted_points.push_back(Weighted_point(Bare_point(1, -1, -1), 4.)); - weighted_points.push_back(Weighted_point(Bare_point(-1, 1, -1), 4.)); - weighted_points.push_back(Weighted_point(Bare_point(-1, -1, 1), 4.)); - weighted_points.push_back(Weighted_point(Bare_point(1, 1, 1), 4.)); - weighted_points.push_back(Weighted_point(Bare_point(2, 2, 2), 1.)); + weighted_points.emplace_back(Bare_point(1, -1, -1), 4.); + weighted_points.emplace_back(Bare_point(-1, 1, -1), 4.); + weighted_points.emplace_back(Bare_point(-1, -1, 1), 4.); + weighted_points.emplace_back(Bare_point(1, 1, 1), 4.); + weighted_points.emplace_back(Bare_point(2, 2, 2), 1.); // ---------------------------------------------------------------------------- // Init of an alpha complex from the list of points @@ -34,10 +34,10 @@ int main(int argc, char **argv) { // ---------------------------------------------------------------------------- // Display information about the alpha complex // ---------------------------------------------------------------------------- - std::clog << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() + std::clog << "Weighted alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; - std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + std::clog << "Iterator on weighted alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; for (auto f_simplex : simplex.filtration_simplex_range()) { std::clog << " ( "; for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { diff --git a/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp b/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp new file mode 100644 index 00000000..d1f3e436 --- /dev/null +++ b/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp @@ -0,0 +1,52 @@ +#include <gudhi/Alpha_complex.h> +// to construct a simplex_tree from alpha complex +#include <gudhi/Simplex_tree.h> + +#include <CGAL/Epeck_d.h> + +#include <iostream> +#include <vector> + +// Explicit dimension 2 Epeck_d kernel +using Kernel = CGAL::Epeck_d< CGAL::Dimension_tag<3> >; +using Bare_point = Kernel::Point_d; +using Weighted_point = Kernel::Weighted_point_d; +using Vector_of_points = std::vector<Weighted_point>; + +int main() { + // ---------------------------------------------------------------------------- + // Init of a list of points and weights from a small molecule + // ---------------------------------------------------------------------------- + Vector_of_points points; + points.emplace_back(Bare_point(1, -1, -1), 4.); + points.emplace_back(Bare_point(-1, 1, -1), 4.); + points.emplace_back(Bare_point(-1, -1, 1), 4.); + points.emplace_back(Bare_point(1, 1, 1), 4.); + points.emplace_back(Bare_point(2, 2, 2), 1.); + + // ---------------------------------------------------------------------------- + // Init of an alpha complex from the list of points + // ---------------------------------------------------------------------------- + Gudhi::alpha_complex::Alpha_complex<Kernel, true> alpha_complex_from_weighted_points(points); + + Gudhi::Simplex_tree<> simplex; + if (alpha_complex_from_weighted_points.create_complex(simplex)) { + // ---------------------------------------------------------------------------- + // Display information about the alpha complex + // ---------------------------------------------------------------------------- + std::clog << "Weighted alpha complex is of dimension " << simplex.dimension() << + " - " << simplex.num_simplices() << " simplices - " << + simplex.num_vertices() << " vertices." << std::endl; + + std::clog << "Iterator on weighted alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + for (auto f_simplex : simplex.filtration_simplex_range()) { + std::clog << " ( "; + for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { + std::clog << vertex << " "; + } + std::clog << ") -> " << "[" << simplex.filtration(f_simplex) << "] "; + std::clog << std::endl; + } + } + return 0; +} diff --git a/src/Alpha_complex/example/weightedalpha3dfrompoints_for_doc.txt b/src/Alpha_complex/example/weightedalpha3dfrompoints_for_doc.txt index 7a09998d..f0695f1a 100644 --- a/src/Alpha_complex/example/weightedalpha3dfrompoints_for_doc.txt +++ b/src/Alpha_complex/example/weightedalpha3dfrompoints_for_doc.txt @@ -1,5 +1,5 @@ -Alpha complex is of dimension 3 - 29 simplices - 5 vertices. -Iterator on alpha complex simplices in the filtration order, with [filtration value]: +Weighted alpha complex is of dimension 3 - 29 simplices - 5 vertices. +Iterator on weighted alpha complex simplices in the filtration order, with [filtration value]: ( 0 ) -> [-4] ( 1 ) -> [-4] ( 2 ) -> [-4] 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) { diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex/Alpha_kernel_d.h b/src/Alpha_complex/include/gudhi/Alpha_complex/Alpha_kernel_d.h new file mode 100644 index 00000000..28d6d7a9 --- /dev/null +++ b/src/Alpha_complex/include/gudhi/Alpha_complex/Alpha_kernel_d.h @@ -0,0 +1,141 @@ +/* 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) 2020 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef ALPHA_COMPLEX_ALPHA_KERNEL_D_H_ +#define ALPHA_COMPLEX_ALPHA_KERNEL_D_H_ + +#include <CGAL/version.h> // for CGAL_VERSION_NR + +#include <Eigen/Core> // for EIGEN_VERSION_AT_LEAST + +#include <utility> // for std::make_pair + +// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 +#if CGAL_VERSION_NR < 1041101000 +# error Alpha_complex is only available for CGAL >= 4.11 +#endif + +#if !EIGEN_VERSION_AT_LEAST(3,1,0) +# error Alpha_complex is only available for Eigen3 >= 3.1.0 installed with CGAL +#endif + +namespace Gudhi { + +namespace alpha_complex { + +/** + * \class Alpha_kernel_d + * \brief Alpha complex kernel container. + * + * \details + * The Alpha complex kernel container stores CGAL Kernel and dispatch basics computations in function of the weighted + * or not version of the Alpha complex. + */ +template < typename Kernel, bool Weighted = false > +class Alpha_kernel_d { +}; + +// Unweighted Kernel_d version +template < typename Kernel > +class Alpha_kernel_d<Kernel, false> { + private: + // Kernel for functions access. + Kernel kernel_; + public: + // Fake type for compilation to succeed (cf. std::conditional in Alpha_complex.h) + using Weighted_point_d = void; + using Point_d = typename Kernel::Point_d; + // Numeric type of coordinates in the kernel + using FT = typename Kernel::FT; + // Sphere is a pair of point and squared radius. + using Sphere = typename std::pair<Point_d, FT>; + + int get_dimension(const Point_d& p0) const { + return kernel_.point_dimension_d_object()(p0); + } + + template<class PointIterator> + Sphere get_sphere(PointIterator begin, PointIterator end) const { + Point_d c = kernel_.construct_circumcenter_d_object()(begin, end); + FT r = kernel_.squared_distance_d_object()(c, *begin); + return std::make_pair(std::move(c), std::move(r)); + } + + template<class PointIterator> + FT get_squared_radius(PointIterator begin, PointIterator end) const { + return kernel_.compute_squared_radius_d_object()(begin, end); + } + + FT get_squared_radius(const Sphere& sph) const { + return sph.second; + } + + bool is_gabriel(const Sphere& circumcenter, const Point_d& point) { + return kernel_.squared_distance_d_object()(circumcenter.first, point) >= circumcenter.second; + } +}; + +// Weighted Kernel_d version +template < typename Kernel > +class Alpha_kernel_d<Kernel, true> { + private: + // Kernel for functions access. + Kernel kernel_; + + public: + // Fake type for compilation to succeed (cf. std::conditional in Alpha_complex.h) + using Point_d = void; + using Weighted_point_d = typename Kernel::Weighted_point_d; + using Bare_point_d = typename Kernel::Point_d; + // Numeric type of coordinates in the kernel + using FT = typename Kernel::FT; + // Sphere is a weighted point (point + weight [= squared radius]). + using Sphere = Weighted_point_d; + + int get_dimension(const Weighted_point_d& p0) const { + return kernel_.point_dimension_d_object()(p0.point()); + } + + template<class PointIterator> + Sphere get_sphere(PointIterator begin, PointIterator end) const { + // power_center_d_object has been renamed between CGAL 5.1 and 5.2 +#if CGAL_VERSION_NR < 1050200000 + return kernel_.power_center_d_object()(begin, end); +#else + return kernel_.construct_power_sphere_d_object()(begin, end); +#endif + } + + template<class PointIterator> + FT get_squared_radius(PointIterator begin, PointIterator end) const { + return kernel_.compute_squared_radius_smallest_orthogonal_sphere_d_object()(begin, end); + } + + FT get_squared_radius(const Sphere& sph) const { + return sph.weight(); + } + + bool is_gabriel(const Sphere& circumcenter, const Weighted_point_d& point) { + // power_center_d_object has been renamed between CGAL 5.1 and 5.2 +#if CGAL_VERSION_NR < 1050200000 + return kernel_.power_distance_d_object()(circumcenter, point) >= 0; +#else + return kernel_.compute_power_product_d_object()(circumcenter, point) >= 0; +#endif + } +}; + +} // namespace alpha_complex + +namespace alphacomplex = alpha_complex; + +} // namespace Gudhi + +#endif // ALPHA_COMPLEX_ALPHA_KERNEL_D_H_
\ No newline at end of file diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h index c19ebb79..562ef139 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h @@ -12,8 +12,10 @@ #ifndef ALPHA_COMPLEX_3D_H_ #define ALPHA_COMPLEX_3D_H_ -#include <boost/version.hpp> #include <boost/variant.hpp> +#include <boost/range/size.hpp> +#include <boost/range/combine.hpp> +#include <boost/range/adaptor/transformed.hpp> #include <gudhi/Debug_utils.h> #include <gudhi/Alpha_complex_options.h> @@ -35,8 +37,6 @@ #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> @@ -53,16 +53,10 @@ # 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 { -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 @@ -104,7 +98,7 @@ struct Value_from_iterator<complexity::EXACT> { * \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/Alpha_shapes_3/index.html#Alpha_shapes_3Definitions">Alpha shapes 3d</a> and * <a href="https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation">Regular * triangulation</a> documentation. * @@ -157,8 +151,10 @@ class Alpha_complex_3d { using Kernel = CGAL::Periodic_3_regular_triangulation_traits_3<Predicates>; }; + public: using Kernel = typename Kernel_3<Predicates, Weighted, Periodic>::Kernel; + private: using TdsVb = typename std::conditional<Periodic, CGAL::Periodic_3_triangulation_ds_vertex_base_3<>, CGAL::Triangulation_ds_vertex_base_3<>>::type; @@ -277,8 +273,8 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ 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)); + alpha_shape_3_ptr_ = std::make_unique<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. @@ -299,20 +295,15 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ 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()), + // 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")); - std::vector<Weighted_point_3> weighted_points_3; + auto weighted_points_3 = boost::range::combine(points, weights) + | boost::adaptors::transformed([](auto const&t){return Weighted_point_3(boost::get<0>(t), boost::get<1>(t));}); - 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)); + alpha_shape_3_ptr_ = std::make_unique<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. @@ -356,7 +347,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ // 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)); + alpha_shape_3_ptr_ = std::make_unique<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. @@ -388,31 +379,27 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ 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()), + // 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")); // 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++; - } + auto weighted_points_3 = boost::range::combine(points, weights) + | boost::adaptors::transformed([=](auto const&t){ + auto w = boost::get<1>(t); + GUDHI_CHECK((w < maximal_possible_weight) && (w >= 0), + std::invalid_argument("Invalid weight " + std::to_string(w) + + ". Must be non-negative and less than maximal possible weight = 1/64*cuboid length squared.")); + return Weighted_point_3(boost::get<0>(t), w); + }); // Define the periodic cube Dt pdt(typename Kernel::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max)); @@ -426,7 +413,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ // 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)); + alpha_shape_3_ptr_ = std::make_unique<Alpha_shape_3>(pdt, 0, Alpha_shape_3::GENERAL); } /** \brief Inserts all Delaunay triangulation into the simplicial complex. @@ -471,6 +458,10 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ #ifdef DEBUG_TRACES std::clog << "filtration_with_alpha_values returns : " << objects.size() << " objects" << std::endl; #endif // DEBUG_TRACES + if (objects.size() == 0) { + std::cerr << "Alpha_complex_3d create_complex - no triangulation as points are on a 2d plane\n"; + return false; // ----- >> + } using Alpha_value_iterator = typename std::vector<FT>::const_iterator; Alpha_value_iterator alpha_value_iterator = alpha_values.begin(); @@ -560,8 +551,10 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ std::clog << "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(); + if (Complexity == complexity::FAST) + // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension + // Only in FAST 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); // -------------------------------------------------------------------------------------------- diff --git a/src/Alpha_complex/test/Alpha_complex_dim3_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_dim3_unit_test.cpp new file mode 100644 index 00000000..e7c261f1 --- /dev/null +++ b/src/Alpha_complex/test/Alpha_complex_dim3_unit_test.cpp @@ -0,0 +1,117 @@ +/* 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): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "alpha_complex_dim3" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <CGAL/Epick_d.h> +#include <CGAL/Epeck_d.h> + +#include <stdexcept> // std::out_of_range +#include <string> +#include <vector> + +#include <gudhi/Alpha_complex.h> +#include <gudhi/Simplex_tree.h> + +// Use dynamic_dimension_tag for the user to be able to set dimension +typedef CGAL::Epeck_d< CGAL::Dynamic_dimension_tag > Exact_kernel_d; +// Use static dimension_tag for the user not to be able to set dimension +typedef CGAL::Epeck_d< CGAL::Dimension_tag<3> > Exact_kernel_s; +// Use dynamic_dimension_tag for the user to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Inexact_kernel_d; +// Use static dimension_tag for the user not to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dimension_tag<3> > Inexact_kernel_s; +// The triangulation uses the default instantiation of the TriangulationDataStructure template parameter + +typedef boost::mpl::list<Exact_kernel_d, Exact_kernel_s, Inexact_kernel_d, Inexact_kernel_s> list_of_kernel_variants; + +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of_kernel_variants) { + // ---------------------------------------------------------------------------- + // + // Init of an alpha-complex from a OFF file + // + // ---------------------------------------------------------------------------- + std::string off_file_name("alphacomplexdoc.off"); + double max_alpha_square_value = 60.0; + std::clog << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << + max_alpha_square_value << "==========" << std::endl; + + Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_file(off_file_name); + + Gudhi::Simplex_tree<> simplex_tree_60; + BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_60, max_alpha_square_value)); + + std::clog << "alpha_complex_from_file.num_vertices()=" << alpha_complex_from_file.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_file.num_vertices() == 7); + + std::clog << "simplex_tree_60.dimension()=" << simplex_tree_60.dimension() << std::endl; + BOOST_CHECK(simplex_tree_60.dimension() == 2); + + std::clog << "simplex_tree_60.num_vertices()=" << simplex_tree_60.num_vertices() << std::endl; + BOOST_CHECK(simplex_tree_60.num_vertices() == 7); + + std::clog << "simplex_tree_60.num_simplices()=" << simplex_tree_60.num_simplices() << std::endl; + BOOST_CHECK(simplex_tree_60.num_simplices() == 25); + + max_alpha_square_value = 59.0; + std::clog << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << + max_alpha_square_value << "==========" << std::endl; + + Gudhi::Simplex_tree<> simplex_tree_59; + BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_59, max_alpha_square_value)); + + std::clog << "alpha_complex_from_file.num_vertices()=" << alpha_complex_from_file.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_file.num_vertices() == 7); + + std::clog << "simplex_tree_59.dimension()=" << simplex_tree_59.dimension() << std::endl; + BOOST_CHECK(simplex_tree_59.dimension() == 2); + + std::clog << "simplex_tree_59.num_vertices()=" << simplex_tree_59.num_vertices() << std::endl; + BOOST_CHECK(simplex_tree_59.num_vertices() == 7); + + std::clog << "simplex_tree_59.num_simplices()=" << simplex_tree_59.num_simplices() << std::endl; + BOOST_CHECK(simplex_tree_59.num_simplices() == 23); +} + + +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_empty_points, TestedKernel, list_of_kernel_variants) { + std::clog << "========== Alpha_complex_from_empty_points ==========" << std::endl; + + // ---------------------------------------------------------------------------- + // Init of an empty list of points + // ---------------------------------------------------------------------------- + std::vector<typename TestedKernel::Point_d> points; + + // ---------------------------------------------------------------------------- + // Init of an alpha complex from the list of points + // ---------------------------------------------------------------------------- + Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_points(points); + + std::clog << "alpha_complex_from_points.num_vertices()=" << alpha_complex_from_points.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_points.num_vertices() == points.size()); + + // Test to the limit + BOOST_CHECK_THROW (alpha_complex_from_points.get_point(0), std::out_of_range); + + Gudhi::Simplex_tree<> simplex_tree; + BOOST_CHECK(!alpha_complex_from_points.create_complex(simplex_tree)); + + std::clog << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl; + BOOST_CHECK(simplex_tree.num_simplices() == 0); + + std::clog << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl; + BOOST_CHECK(simplex_tree.dimension() == -1); + + std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; + BOOST_CHECK(simplex_tree.num_vertices() == points.size()); +} diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp index 4b37e4bd..b474917f 100644 --- a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp +++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp @@ -13,75 +13,17 @@ #include <boost/test/unit_test.hpp> #include <boost/mpl/list.hpp> -#include <CGAL/Delaunay_triangulation.h> #include <CGAL/Epick_d.h> #include <CGAL/Epeck_d.h> -#include <cmath> // float comparison -#include <limits> +#include <stdexcept> // std::out_of_range #include <string> #include <vector> #include <gudhi/Alpha_complex.h> -// to construct a simplex_tree from Delaunay_triangulation -#include <gudhi/graph_simplicial_complex.h> #include <gudhi/Simplex_tree.h> #include <gudhi/Unitary_tests_utils.h> -// Use dynamic_dimension_tag for the user to be able to set dimension -typedef CGAL::Epeck_d< CGAL::Dynamic_dimension_tag > Exact_kernel_d; -// Use static dimension_tag for the user not to be able to set dimension -typedef CGAL::Epeck_d< CGAL::Dimension_tag<3> > Exact_kernel_s; -// Use dynamic_dimension_tag for the user to be able to set dimension -typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Inexact_kernel_d; -// Use static dimension_tag for the user not to be able to set dimension -typedef CGAL::Epick_d< CGAL::Dimension_tag<3> > Inexact_kernel_s; -// The triangulation uses the default instantiation of the TriangulationDataStructure template parameter - -typedef boost::mpl::list<Exact_kernel_d, Exact_kernel_s, Inexact_kernel_d, Inexact_kernel_s> list_of_kernel_variants; - -BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of_kernel_variants) { - // ---------------------------------------------------------------------------- - // - // Init of an alpha-complex from a OFF file - // - // ---------------------------------------------------------------------------- - std::string off_file_name("alphacomplexdoc.off"); - double max_alpha_square_value = 60.0; - std::clog << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << - max_alpha_square_value << "==========" << std::endl; - - Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_file(off_file_name); - - Gudhi::Simplex_tree<> simplex_tree_60; - BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_60, max_alpha_square_value)); - - std::clog << "simplex_tree_60.dimension()=" << simplex_tree_60.dimension() << std::endl; - BOOST_CHECK(simplex_tree_60.dimension() == 2); - - std::clog << "simplex_tree_60.num_vertices()=" << simplex_tree_60.num_vertices() << std::endl; - BOOST_CHECK(simplex_tree_60.num_vertices() == 7); - - std::clog << "simplex_tree_60.num_simplices()=" << simplex_tree_60.num_simplices() << std::endl; - BOOST_CHECK(simplex_tree_60.num_simplices() == 25); - - max_alpha_square_value = 59.0; - std::clog << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << - max_alpha_square_value << "==========" << std::endl; - - Gudhi::Simplex_tree<> simplex_tree_59; - BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_59, max_alpha_square_value)); - - std::clog << "simplex_tree_59.dimension()=" << simplex_tree_59.dimension() << std::endl; - BOOST_CHECK(simplex_tree_59.dimension() == 2); - - std::clog << "simplex_tree_59.num_vertices()=" << simplex_tree_59.num_vertices() << std::endl; - BOOST_CHECK(simplex_tree_59.num_vertices() == 7); - - std::clog << "simplex_tree_59.num_simplices()=" << simplex_tree_59.num_simplices() << std::endl; - BOOST_CHECK(simplex_tree_59.num_simplices() == 23); -} - // Use static dimension_tag for the user not to be able to set dimension typedef CGAL::Epeck_d< CGAL::Dimension_tag<4> > Kernel_4; typedef Kernel_4::Point_d Point_4; @@ -120,6 +62,9 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { Gudhi::Simplex_tree<> simplex_tree; BOOST_CHECK(alpha_complex_from_points.create_complex(simplex_tree)); + std::clog << "alpha_complex_from_points.num_vertices()=" << alpha_complex_from_points.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_points.num_vertices() == points.size()); + // Another way to check num_simplices std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; int num_simplices = 0; @@ -227,34 +172,6 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { } -BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_empty_points, TestedKernel, list_of_kernel_variants) { - std::clog << "========== Alpha_complex_from_empty_points ==========" << std::endl; - - // ---------------------------------------------------------------------------- - // Init of an empty list of points - // ---------------------------------------------------------------------------- - std::vector<typename TestedKernel::Point_d> points; - - // ---------------------------------------------------------------------------- - // Init of an alpha complex from the list of points - // ---------------------------------------------------------------------------- - Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_points(points); - - // Test to the limit - BOOST_CHECK_THROW (alpha_complex_from_points.get_point(0), std::out_of_range); - - Gudhi::Simplex_tree<> simplex_tree; - BOOST_CHECK(!alpha_complex_from_points.create_complex(simplex_tree)); - - std::clog << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl; - BOOST_CHECK(simplex_tree.num_simplices() == 0); - - std::clog << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl; - BOOST_CHECK(simplex_tree.dimension() == -1); - - std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; - BOOST_CHECK(simplex_tree.num_vertices() == points.size()); -} using Inexact_kernel_2 = CGAL::Epick_d< CGAL::Dimension_tag<2> >; using Exact_kernel_2 = CGAL::Epeck_d< CGAL::Dimension_tag<2> >; @@ -291,6 +208,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_with_duplicated_points, TestedKernel std::clog << "create_complex" << std::endl; BOOST_CHECK(alpha_complex_from_points.create_complex(simplex_tree)); + std::clog << "alpha_complex_from_points.num_vertices()=" << alpha_complex_from_points.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_points.num_vertices() < points.size()); + std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; BOOST_CHECK(simplex_tree.num_vertices() < points.size()); diff --git a/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp b/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp new file mode 100644 index 00000000..6da4c084 --- /dev/null +++ b/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp @@ -0,0 +1,109 @@ +/* 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) 2020 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "alpha_kernel_d" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <CGAL/Epick_d.h> +#include <CGAL/Epeck_d.h> +#include <CGAL/NT_converter.h> + +#include <iostream> +#include <vector> +#include <utility> // for std::pair + +#include <gudhi/Alpha_complex/Alpha_kernel_d.h> +#include <gudhi/Unitary_tests_utils.h> + +// Use dynamic_dimension_tag for the user to be able to set dimension +typedef CGAL::Epeck_d< CGAL::Dynamic_dimension_tag > Exact_kernel_d; +// Use static dimension_tag for the user not to be able to set dimension +typedef CGAL::Epeck_d< CGAL::Dimension_tag<4> > Exact_kernel_s; +// Use dynamic_dimension_tag for the user to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Inexact_kernel_d; +// Use static dimension_tag for the user not to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dimension_tag<4> > Inexact_kernel_s; +// The triangulation uses the default instantiation of the TriangulationDataStructure template parameter + +typedef boost::mpl::list<Exact_kernel_d, Exact_kernel_s, Inexact_kernel_d, Inexact_kernel_s> list_of_kernel_variants; + +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_kernel_d_dimension, TestedKernel, list_of_kernel_variants) { + // Test for a point (weighted or not) in 4d, that the dimension is 4. + + Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, false> kernel; + std::vector<double> p0 {0., 1., 2., 3.}; + typename TestedKernel::Point_d p0_d(p0.begin(), p0.end()); + + std::clog << "Dimension is " << kernel.get_dimension(p0_d) << std::endl; + BOOST_CHECK(kernel.get_dimension(p0_d) == 4); + + Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, true> w_kernel; + typename TestedKernel::Weighted_point_d w_p0_d(p0_d, 10.); + + std::clog << "Dimension is " << w_kernel.get_dimension(w_p0_d) << std::endl; + BOOST_CHECK(w_kernel.get_dimension(w_p0_d) == 4); +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_kernel_d_sphere, TestedKernel, list_of_kernel_variants) { + // Test with 5 points on a 3-sphere, that get_sphere returns the same center and squared radius + // for dD unweighted and for dD weighted with all weights at 0. + + using Unweighted_kernel = Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, false>; + // Sphere: (x-1)² + (y-1)² + z² + t² = 1 + // At least 5 points for a 3-sphere + std::vector<double> p0 {1., 0., 0., 0.}; + std::vector<double> p1 {0., 1., 0., 0.}; + std::vector<double> p2 {1., 1., 1., 0.}; + std::vector<double> p3 {1., 1., 0., 1.}; + std::vector<double> p4 {1., 1., -1., 0.}; + + using Point_d = typename Unweighted_kernel::Point_d; + std::vector<Point_d> unw_pts; + unw_pts.emplace_back(p0.begin(), p0.end()); + unw_pts.emplace_back(p1.begin(), p1.end()); + unw_pts.emplace_back(p2.begin(), p2.end()); + unw_pts.emplace_back(p3.begin(), p3.end()); + unw_pts.emplace_back(p4.begin(), p4.end()); + + Unweighted_kernel kernel; + auto unw_sphere = kernel.get_sphere(unw_pts.cbegin(), unw_pts.cend()); + + std::clog << "Center is " << unw_sphere.first << " - squared radius is " << unw_sphere.second << std::endl; + + using Weighted_kernel = Gudhi::alpha_complex::Alpha_kernel_d<TestedKernel, true>; + + using Weighted_point_d = typename Weighted_kernel::Weighted_point_d; + using Bare_point_d = typename Weighted_kernel::Bare_point_d; + std::vector<Weighted_point_d> w_pts; + w_pts.emplace_back(Bare_point_d(p0.begin(), p0.end()), 0.); + w_pts.emplace_back(Bare_point_d(p1.begin(), p1.end()), 0.); + w_pts.emplace_back(Bare_point_d(p2.begin(), p2.end()), 0.); + w_pts.emplace_back(Bare_point_d(p3.begin(), p3.end()), 0.); + w_pts.emplace_back(Bare_point_d(p4.begin(), p4.end()), 0.); + + Weighted_kernel w_kernel; + auto w_sphere = w_kernel.get_sphere(w_pts.cbegin(), w_pts.cend()); + + std::clog << "Center is " << w_sphere.point() << " - squared radius is " << w_sphere.weight() << std::endl; + + CGAL::NT_converter<typename Weighted_kernel::FT, double> cast_to_double; + // The results shall be the same with weights = 0. + GUDHI_TEST_FLOAT_EQUALITY_CHECK(cast_to_double(unw_sphere.second), cast_to_double(w_sphere.weight())); + BOOST_CHECK(unw_sphere.first == w_sphere.point()); + + auto unw_sq_rd = kernel.get_squared_radius(unw_pts.cbegin(), unw_pts.cend()); + std::clog << "Squared radius is " << unw_sq_rd << std::endl; + GUDHI_TEST_FLOAT_EQUALITY_CHECK(cast_to_double(unw_sphere.second), cast_to_double(unw_sq_rd)); + auto w_sq_rd = w_kernel.get_squared_radius(w_pts.cbegin(), w_pts.cend()); + std::clog << "Squared radius is " << w_sq_rd << std::endl; + GUDHI_TEST_FLOAT_EQUALITY_CHECK(cast_to_double(w_sphere.weight()), cast_to_double(w_sq_rd)); +} diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt index fe4b23e4..dd2c235f 100644 --- a/src/Alpha_complex/test/CMakeLists.txt +++ b/src/Alpha_complex/test/CMakeLists.txt @@ -8,16 +8,24 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) add_executable ( Alpha_complex_test_unit Alpha_complex_unit_test.cpp ) target_link_libraries(Alpha_complex_test_unit ${CGAL_LIBRARY}) + add_executable ( Alpha_complex_dim3_test_unit Alpha_complex_dim3_unit_test.cpp ) + target_link_libraries(Alpha_complex_dim3_test_unit ${CGAL_LIBRARY}) add_executable ( Delaunay_complex_test_unit Delaunay_complex_unit_test.cpp ) target_link_libraries(Delaunay_complex_test_unit ${CGAL_LIBRARY}) if (TBB_FOUND) target_link_libraries(Alpha_complex_test_unit ${TBB_LIBRARIES}) + target_link_libraries(Alpha_complex_dim3_test_unit ${TBB_LIBRARIES}) target_link_libraries(Delaunay_complex_test_unit ${TBB_LIBRARIES}) endif() gudhi_add_boost_test(Alpha_complex_test_unit) + gudhi_add_boost_test(Alpha_complex_dim3_test_unit) gudhi_add_boost_test(Delaunay_complex_test_unit) +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + +if (NOT CGAL_VERSION VERSION_LESS 4.11.0) + add_executable ( Alpha_complex_3d_test_unit Alpha_complex_3d_unit_test.cpp ) target_link_libraries(Alpha_complex_3d_test_unit ${CGAL_LIBRARY}) add_executable ( Weighted_alpha_complex_3d_test_unit Weighted_alpha_complex_3d_unit_test.cpp ) @@ -38,4 +46,35 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) gudhi_add_boost_test(Periodic_alpha_complex_3d_test_unit) gudhi_add_boost_test(Weighted_periodic_alpha_complex_3d_test_unit) -endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) +endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) + +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) + add_executable ( Alpha_kernel_d_test_unit Alpha_kernel_d_unit_test.cpp ) + target_link_libraries(Alpha_kernel_d_test_unit ${CGAL_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Alpha_kernel_d_test_unit ${TBB_LIBRARIES}) + endif() + gudhi_add_boost_test(Alpha_kernel_d_test_unit) + + add_executable ( Weighted_alpha_complex_test_unit Weighted_alpha_complex_unit_test.cpp ) + target_link_libraries(Weighted_alpha_complex_test_unit ${CGAL_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Weighted_alpha_complex_test_unit ${TBB_LIBRARIES}) + endif() + gudhi_add_boost_test(Weighted_alpha_complex_test_unit) + + add_executable ( Weighted_alpha_complex_non_visible_points_test_unit Weighted_alpha_complex_non_visible_points_unit_test.cpp ) + target_link_libraries(Weighted_alpha_complex_non_visible_points_test_unit ${CGAL_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Weighted_alpha_complex_non_visible_points_test_unit ${TBB_LIBRARIES}) + endif() + gudhi_add_boost_test(Weighted_alpha_complex_non_visible_points_test_unit) + + add_executable ( Zero_weighted_alpha_complex_test_unit Zero_weighted_alpha_complex_unit_test.cpp ) + target_link_libraries(Zero_weighted_alpha_complex_test_unit ${CGAL_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Zero_weighted_alpha_complex_test_unit ${TBB_LIBRARIES}) + endif() + gudhi_add_boost_test(Zero_weighted_alpha_complex_test_unit) + +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) diff --git a/src/Alpha_complex/test/Weighted_alpha_complex_non_visible_points_unit_test.cpp b/src/Alpha_complex/test/Weighted_alpha_complex_non_visible_points_unit_test.cpp new file mode 100644 index 00000000..dd83c1da --- /dev/null +++ b/src/Alpha_complex/test/Weighted_alpha_complex_non_visible_points_unit_test.cpp @@ -0,0 +1,60 @@ +/* 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) 2020 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "weighted_alpha_complex_non_visible_points" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <CGAL/Epick_d.h> +#include <CGAL/Epeck_d.h> + +#include <vector> + +#include <gudhi/Alpha_complex.h> +#include <gudhi/Simplex_tree.h> + + +using list_of_1d_kernel_variants = boost::mpl::list<CGAL::Epeck_d< CGAL::Dynamic_dimension_tag >, + CGAL::Epeck_d< CGAL::Dimension_tag<1>>, + CGAL::Epick_d< CGAL::Dynamic_dimension_tag >, + CGAL::Epick_d< CGAL::Dimension_tag<1>> + >; + +BOOST_AUTO_TEST_CASE_TEMPLATE(Weighted_alpha_complex_non_visible_points, Kernel, list_of_1d_kernel_variants) { + // check that for 2 closed weighted 1-d points, one with a high weight to hide the second one with a small weight, + // that the point with a small weight has the same high filtration value than the edge formed by the 2 points + using Point_d = typename Kernel::Point_d; + std::vector<Point_d> points; + std::vector<double> p1 {0.}; + points.emplace_back(p1.begin(), p1.end()); + // closed enough points + std::vector<double> p2 {0.1}; + points.emplace_back(p2.begin(), p2.end()); + std::vector<typename Kernel::FT> weights {100., 0.01}; + + Gudhi::alpha_complex::Alpha_complex<Kernel, true> alpha_complex(points, weights); + Gudhi::Simplex_tree<> stree; + BOOST_CHECK(alpha_complex.create_complex(stree)); + + std::clog << "Iterator on weighted alpha complex simplices in the filtration order, with [filtration value]:" + << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + std::clog << " ( "; + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + std::clog << vertex << " "; + } + std::clog << ") -> " << "[" << stree.filtration(f_simplex) << "] " << std::endl; + } + + BOOST_CHECK(stree.filtration(stree.find({0})) == -100.); + BOOST_CHECK(stree.filtration(stree.find({1})) == stree.filtration(stree.find({0, 1}))); + BOOST_CHECK(stree.filtration(stree.find({1})) > 100000); +}
\ No newline at end of file diff --git a/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp new file mode 100644 index 00000000..875704ee --- /dev/null +++ b/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp @@ -0,0 +1,127 @@ +/* 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) 2020 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "weighted_alpha_complex" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <CGAL/Epeck_d.h> + +#include <vector> +#include <random> +#include <array> +#include <cmath> // for std::fabs + +#include <gudhi/Alpha_complex.h> +#include <gudhi/Alpha_complex_3d.h> +#include <gudhi/Simplex_tree.h> + +BOOST_AUTO_TEST_CASE(Weighted_alpha_complex_3d_comparison) { + // check that for random weighted 3d points in safe mode the 3D and dD codes give the same result with some tolerance + + // Random points construction + using Kernel_dD = CGAL::Epeck_d< CGAL::Dimension_tag<3> >; + using Bare_point_d = typename Kernel_dD::Point_d; + using Weighted_point_d = typename Kernel_dD::Weighted_point_d; + std::vector<Weighted_point_d> w_points_d; + + using Exact_weighted_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, true, false>; + using Bare_point_3 = typename Exact_weighted_alpha_complex_3d::Bare_point_3; + using Weighted_point_3 = typename Exact_weighted_alpha_complex_3d::Weighted_point_3; + std::vector<Weighted_point_3> w_points_3; + + std::uniform_real_distribution<double> rd_pts(-10., 10.); + std::uniform_real_distribution<double> rd_wghts(-0.5, 0.5); + std::random_device rand_dev; + std::mt19937 rand_engine(rand_dev()); + for (int idx = 0; idx < 20; idx++) { + std::vector<double> point {rd_pts(rand_engine), rd_pts(rand_engine), rd_pts(rand_engine)}; + double weight = rd_wghts(rand_engine); + w_points_d.emplace_back(Bare_point_d(point.begin(), point.end()), weight); + w_points_3.emplace_back(Bare_point_3(point[0], point[1], point[2]), weight); + } + + // Structures necessary for comparison + using Points = std::vector<std::array<double,3>>; + using Points_and_filtrations = std::map<Points, double>; + Points_and_filtrations pts_fltr_dD; + Points_and_filtrations pts_fltr_3d; + + // Weighted alpha complex for dD version + Gudhi::alpha_complex::Alpha_complex<Kernel_dD, true> alpha_complex_dD_from_weighted_points(w_points_d); + Gudhi::Simplex_tree<> w_simplex_d; + BOOST_CHECK(alpha_complex_dD_from_weighted_points.create_complex(w_simplex_d)); + + std::clog << "Iterator on weighted alpha complex dD simplices in the filtration order, with [filtration value]:" + << std::endl; + for (auto f_simplex : w_simplex_d.filtration_simplex_range()) { + Points points; + for (auto vertex : w_simplex_d.simplex_vertex_range(f_simplex)) { + CGAL::NT_converter<Kernel_dD::RT, double> cgal_converter; + Bare_point_d pt = alpha_complex_dD_from_weighted_points.get_point(vertex).point(); + points.push_back({cgal_converter(pt[0]), cgal_converter(pt[1]), cgal_converter(pt[2])}); + } + std::clog << " ( "; + std::sort (points.begin(), points.end()); + for (auto point : points) { + std::clog << point[0] << " " << point[1] << " " << point[2] << " | "; + } + std::clog << ") -> " << "[" << w_simplex_d.filtration(f_simplex) << "] "; + std::clog << std::endl; + pts_fltr_dD[points] = w_simplex_d.filtration(f_simplex); + } + + // Weighted alpha complex for 3D version + Exact_weighted_alpha_complex_3d alpha_complex_3D_from_weighted_points(w_points_3); + Gudhi::Simplex_tree<> w_simplex_3; + BOOST_CHECK(alpha_complex_3D_from_weighted_points.create_complex(w_simplex_3)); + + std::clog << "Iterator on weighted alpha complex 3D simplices in the filtration order, with [filtration value]:" + << std::endl; + for (auto f_simplex : w_simplex_3.filtration_simplex_range()) { + Points points; + for (auto vertex : w_simplex_3.simplex_vertex_range(f_simplex)) { + Bare_point_3 pt = alpha_complex_3D_from_weighted_points.get_point(vertex).point(); + CGAL::NT_converter<Exact_weighted_alpha_complex_3d::Kernel::RT, double> cgal_converter; + points.push_back({cgal_converter(pt[0]), cgal_converter(pt[1]), cgal_converter(pt[2])}); + } + std::clog << " ( "; + std::sort (points.begin(), points.end()); + for (auto point : points) { + std::clog << point[0] << " " << point[1] << " " << point[2] << " | "; + } + std::clog << ") -> " << "[" << w_simplex_3.filtration(f_simplex) << "] " << std::endl; + pts_fltr_3d[points] = w_simplex_d.filtration(f_simplex); + } + + // Compares structures + auto d3_itr = pts_fltr_3d.begin(); + auto dD_itr = pts_fltr_dD.begin(); + for (; d3_itr != pts_fltr_3d.end() && dD_itr != pts_fltr_dD.end(); ++d3_itr) { + if (d3_itr->first != dD_itr->first) { + for(auto point : d3_itr->first) + std::clog << point[0] << " " << point[1] << " " << point[2] << " | "; + std::clog << " versus "; + for(auto point : dD_itr->first) + std::clog << point[0] << " " << point[1] << " " << point[2] << " | "; + std::clog << std::endl; + BOOST_CHECK(false); + } + // In safe mode, relative error is less than 1e-5 (can be changed with set_relative_precision_of_to_double) + if (std::fabs(d3_itr->second - dD_itr->second) > 1e-5 * (std::fabs(d3_itr->second) + std::fabs(dD_itr->second))) { + std::clog << d3_itr->second << " versus " << dD_itr->second << " diff " + << std::fabs(d3_itr->second - dD_itr->second) << std::endl; + BOOST_CHECK(false); + } + ++dD_itr; + } +}
\ No newline at end of file diff --git a/src/Alpha_complex/test/Zero_weighted_alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Zero_weighted_alpha_complex_unit_test.cpp new file mode 100644 index 00000000..b7df07c7 --- /dev/null +++ b/src/Alpha_complex/test/Zero_weighted_alpha_complex_unit_test.cpp @@ -0,0 +1,77 @@ +/* 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) 2020 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "zero_weighted_alpha_complex" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <CGAL/Epeck_d.h> + +#include <vector> +#include <random> +#include <cmath> // for std::fabs + +#include <gudhi/Alpha_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Unitary_tests_utils.h> + +using list_of_exact_kernel_variants = boost::mpl::list<CGAL::Epeck_d< CGAL::Dynamic_dimension_tag >, + CGAL::Epeck_d< CGAL::Dimension_tag<4> > + > ; + +BOOST_AUTO_TEST_CASE_TEMPLATE(Zero_weighted_alpha_complex, Kernel, list_of_exact_kernel_variants) { + // Check that in exact mode for static dimension 4 the code for dD unweighted and for dD weighted with all weights + // 0 give exactly the same simplex tree (simplices and filtration values). + + // Random points construction + using Point_d = typename Kernel::Point_d; + std::vector<Point_d> points; + std::uniform_real_distribution<double> rd_pts(-10., 10.); + std::random_device rand_dev; + std::mt19937 rand_engine(rand_dev()); + for (int idx = 0; idx < 20; idx++) { + std::vector<double> point {rd_pts(rand_engine), rd_pts(rand_engine), rd_pts(rand_engine), rd_pts(rand_engine)}; + points.emplace_back(point.begin(), point.end()); + } + + // Alpha complex from points + Gudhi::alpha_complex::Alpha_complex<Kernel, false> alpha_complex_from_points(points); + Gudhi::Simplex_tree<> simplex; + Gudhi::Simplex_tree<>::Filtration_value infty = std::numeric_limits<Gudhi::Simplex_tree<>::Filtration_value>::infinity(); + BOOST_CHECK(alpha_complex_from_points.create_complex(simplex, infty, true)); + std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" + << std::endl; + for (auto f_simplex : simplex.filtration_simplex_range()) { + std::clog << " ( "; + for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { + std::clog << vertex << " "; + } + std::clog << ") -> " << "[" << simplex.filtration(f_simplex) << "] " << std::endl; + } + + // Alpha complex from zero weighted points + std::vector<typename Kernel::FT> weights(20, 0.); + Gudhi::alpha_complex::Alpha_complex<Kernel, true> alpha_complex_from_zero_weighted_points(points, weights); + Gudhi::Simplex_tree<> zw_simplex; + BOOST_CHECK(alpha_complex_from_zero_weighted_points.create_complex(zw_simplex, infty, true)); + + std::clog << "Iterator on zero weighted alpha complex simplices in the filtration order, with [filtration value]:" + << std::endl; + for (auto f_simplex : zw_simplex.filtration_simplex_range()) { + std::clog << " ( "; + for (auto vertex : zw_simplex.simplex_vertex_range(f_simplex)) { + std::clog << vertex << " "; + } + std::clog << ") -> " << "[" << zw_simplex.filtration(f_simplex) << "] " << std::endl; + } + + BOOST_CHECK(zw_simplex == simplex); +}
\ No newline at end of file diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt index 2ffbdde0..303bd0a6 100644 --- a/src/Alpha_complex/utilities/CMakeLists.txt +++ b/src/Alpha_complex/utilities/CMakeLists.txt @@ -1,78 +1,82 @@ project(Alpha_complex_utilities) -if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) - add_executable (alpha_complex_persistence alpha_complex_persistence.cpp) - target_link_libraries(alpha_complex_persistence ${CGAL_LIBRARY} Boost::program_options) +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) + if (TARGET Boost::program_options) + add_executable (alpha_complex_persistence alpha_complex_persistence.cpp) + target_link_libraries(alpha_complex_persistence ${CGAL_LIBRARY} Boost::program_options) - if (TBB_FOUND) - target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES}) - endif(TBB_FOUND) - add_test(NAME Alpha_complex_utilities_safe_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence> - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "safe.pers") - add_test(NAME Alpha_complex_utilities_fast_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence> - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "fast.pers" "-f") - add_test(NAME Alpha_complex_utilities_exact_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence> - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "exact.pers" "-e") - if (DIFF_PATH) - add_test(Alpha_complex_utilities_diff_exact_alpha_complex ${DIFF_PATH} - "exact.pers" "safe.pers") - set_tests_properties(Alpha_complex_utilities_diff_exact_alpha_complex PROPERTIES DEPENDS - "Alpha_complex_utilities_exact_alpha_complex_persistence;Alpha_complex_utilities_safe_alpha_complex_persistence") - - add_test(Alpha_complex_utilities_diff_fast_alpha_complex ${DIFF_PATH} - "fast.pers" "safe.pers") - set_tests_properties(Alpha_complex_utilities_diff_fast_alpha_complex PROPERTIES DEPENDS - "Alpha_complex_utilities_fast_alpha_complex_persistence;Alpha_complex_utilities_safe_alpha_complex_persistence") + if (TBB_FOUND) + target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) + add_test(NAME Alpha_complex_utilities_safe_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "safe.pers") + add_test(NAME Alpha_complex_utilities_fast_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "fast.pers" "-f") + add_test(NAME Alpha_complex_utilities_exact_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "exact.pers" "-e") + if (DIFF_PATH) + add_test(Alpha_complex_utilities_diff_exact_alpha_complex ${DIFF_PATH} + "exact.pers" "safe.pers") + set_tests_properties(Alpha_complex_utilities_diff_exact_alpha_complex PROPERTIES DEPENDS + "Alpha_complex_utilities_exact_alpha_complex_persistence;Alpha_complex_utilities_safe_alpha_complex_persistence") - endif() - - install(TARGETS alpha_complex_persistence DESTINATION bin) + add_test(Alpha_complex_utilities_diff_fast_alpha_complex ${DIFF_PATH} + "fast.pers" "safe.pers") + set_tests_properties(Alpha_complex_utilities_diff_fast_alpha_complex PROPERTIES DEPENDS + "Alpha_complex_utilities_fast_alpha_complex_persistence;Alpha_complex_utilities_safe_alpha_complex_persistence") + endif() - add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp) - target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY} Boost::program_options) - if (TBB_FOUND) - target_link_libraries(alpha_complex_3d_persistence ${TBB_LIBRARIES}) - endif(TBB_FOUND) + install(TARGETS alpha_complex_persistence DESTINATION bin) + endif() +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) - add_test(NAME Alpha_complex_utilities_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" - "-p" "2" "-m" "0.45" "-o" "safe_3d.pers") +if (NOT CGAL_VERSION VERSION_LESS 4.11.0) + if (TARGET Boost::program_options) + add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp) + target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY} Boost::program_options) + if (TBB_FOUND) + target_link_libraries(alpha_complex_3d_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) - add_test(NAME Alpha_complex_utilities_exact_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" - "-p" "2" "-m" "0.45" "-o" "exact_3d.pers" "-e") + add_test(NAME Alpha_complex_utilities_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" + "-p" "2" "-m" "0.45" "-o" "safe_3d.pers") - add_test(NAME Alpha_complex_utilities_fast_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" - "-p" "2" "-m" "0.45" "-o" "fast_3d.pers" "-f") + add_test(NAME Alpha_complex_utilities_exact_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" + "-p" "2" "-m" "0.45" "-o" "exact_3d.pers" "-e") - if (DIFF_PATH) - add_test(Alpha_complex_utilities_diff_exact_alpha_complex_3d ${DIFF_PATH} - "exact_3d.pers" "safe_3d.pers") - set_tests_properties(Alpha_complex_utilities_diff_exact_alpha_complex_3d PROPERTIES DEPENDS - "Alpha_complex_utilities_exact_alpha_complex_3d;Alpha_complex_utilities_alpha_complex_3d") - add_test(Alpha_complex_utilities_diff_fast_alpha_complex_3d ${DIFF_PATH} - "fast_3d.pers" "safe_3d.pers") - set_tests_properties(Alpha_complex_utilities_diff_fast_alpha_complex_3d PROPERTIES DEPENDS - "Alpha_complex_utilities_fast_alpha_complex_3d;Alpha_complex_utilities_alpha_complex_3d") - endif() + add_test(NAME Alpha_complex_utilities_fast_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" + "-p" "2" "-m" "0.45" "-o" "fast_3d.pers" "-f") - add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> - "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" - "-c" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" - "-p" "2" "-m" "0") + if (DIFF_PATH) + add_test(Alpha_complex_utilities_diff_exact_alpha_complex_3d ${DIFF_PATH} + "exact_3d.pers" "safe_3d.pers") + set_tests_properties(Alpha_complex_utilities_diff_exact_alpha_complex_3d PROPERTIES DEPENDS + "Alpha_complex_utilities_exact_alpha_complex_3d;Alpha_complex_utilities_alpha_complex_3d") + add_test(Alpha_complex_utilities_diff_fast_alpha_complex_3d ${DIFF_PATH} + "fast_3d.pers" "safe_3d.pers") + set_tests_properties(Alpha_complex_utilities_diff_fast_alpha_complex_3d PROPERTIES DEPENDS + "Alpha_complex_utilities_fast_alpha_complex_3d;Alpha_complex_utilities_alpha_complex_3d") + endif() - add_test(NAME Alpha_complex_utilities_weighted_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> - "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" - "-w" "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.weights" - "-p" "2" "-m" "0") + add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" + "-c" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" + "-p" "2" "-m" "0") - add_test(NAME Alpha_complex_utilities_weighted_periodic_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> - "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" - "-w" "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.weights" - "-c" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" - "-p" "2" "-m" "0" "-e") + add_test(NAME Alpha_complex_utilities_weighted_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" + "-w" "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.weights" + "-p" "2" "-m" "0") - install(TARGETS alpha_complex_3d_persistence DESTINATION bin) - -endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + add_test(NAME Alpha_complex_utilities_weighted_periodic_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" + "-w" "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.weights" + "-c" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" + "-p" "2" "-m" "0" "-e") + + install(TARGETS alpha_complex_3d_persistence DESTINATION bin) + endif() +endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp index 91899040..e65d8c6f 100644 --- a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp @@ -263,7 +263,7 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool "cuboid-file,c", po::value<std::string>(&cuboid_file), "Name of file describing the periodic domain. Format is:\n min_hx min_hy min_hz\n max_hx max_hy max_hz")( "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog")( + "Name of file in which the persistence diagram is written. Default print in standard output")( "max-alpha-square-value,r", po::value<Filtration_value>(&alpha_square_max_value) ->default_value(std::numeric_limits<Filtration_value>::infinity()), diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp index 7c898dfd..29edbd8e 100644 --- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -11,24 +11,88 @@ #include <boost/program_options.hpp> #include <CGAL/Epick_d.h> +#include <CGAL/Epeck_d.h> #include <gudhi/Alpha_complex.h> #include <gudhi/Persistent_cohomology.h> // to construct a simplex_tree from alpha complex #include <gudhi/Simplex_tree.h> +#include <gudhi/Points_off_io.h> #include <iostream> #include <string> #include <limits> // for numeric_limits +#include <vector> +#include <fstream> using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast, - std::string &output_file_diag, Filtration_value &alpha_square_max_value, + std::string &weight_file, std::string &output_file_diag, Filtration_value &alpha_square_max_value, int &coeff_field_characteristic, Filtration_value &min_persistence); +template<class Point_d> +std::vector<Point_d> read_off(const std::string &off_file_points) { + Gudhi::Points_off_reader<Point_d> off_reader(off_file_points); + if (!off_reader.is_valid()) { + std::cerr << "Alpha_complex - Unable to read file " << off_file_points << "\n"; + exit(-1); // ----- >> + } + return off_reader.get_point_cloud(); +} + +std::vector<double> read_weight_file(const std::string &weight_file) { + std::vector<double> weights; + // Read weights information from file + std::ifstream weights_ifstr(weight_file); + if (weights_ifstr.good()) { + double weight = 0.0; + // Attempt read the weight in a double format, return false if it fails + while (weights_ifstr >> weight) { + weights.push_back(weight); + } + } else { + std::cerr << "Unable to read weights file " << weight_file << std::endl; + exit(-1); + } + return weights; +} + +template<class Kernel> +Simplex_tree create_simplex_tree(const std::string &off_file_points, const std::string &weight_file, + bool exact_version, Filtration_value alpha_square_max_value) { + Simplex_tree stree; + auto points = read_off<typename Kernel::Point_d>(off_file_points); + + if (weight_file != std::string()) { + std::vector<double> weights = read_weight_file(weight_file); + if (points.size() != weights.size()) { + std::cerr << "Alpha_complex - Inconsistency between number of points (" << points.size() + << ") and number of weights (" << weights.size() << ")" << "\n"; + exit(-1); // ----- >> + } + // Init of an alpha complex from an OFF file + Gudhi::alpha_complex::Alpha_complex<Kernel, true> alpha_complex_from_file(points, weights); + + if (!alpha_complex_from_file.create_complex(stree, alpha_square_max_value, exact_version)) { + std::cerr << "Alpha complex simplicial complex creation failed." << std::endl; + exit(-1); + } + } else { + // Init of an alpha complex from an OFF file + Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(points); + + if (!alpha_complex_from_file.create_complex(stree, alpha_square_max_value, exact_version)) { + std::cerr << "Alpha complex simplicial complex creation failed." << std::endl; + exit(-1); + } + } + return stree; +} + int main(int argc, char **argv) { + std::string weight_file; std::string off_file_points; std::string output_file_diag; bool exact_version = false; @@ -37,48 +101,34 @@ int main(int argc, char **argv) { int coeff_field_characteristic; Filtration_value min_persistence; - program_options(argc, argv, off_file_points, exact_version, fast_version, output_file_diag, alpha_square_max_value, - coeff_field_characteristic, min_persistence); + program_options(argc, argv, off_file_points, exact_version, fast_version, weight_file, output_file_diag, + alpha_square_max_value, coeff_field_characteristic, min_persistence); if ((exact_version) && (fast_version)) { std::cerr << "You cannot set the exact and the fast version." << std::endl; exit(-1); } - Simplex_tree simplex; + Simplex_tree stree; if (fast_version) { // WARNING : CGAL::Epick_d is fast but not safe (unlike CGAL::Epeck_d) // (i.e. when the points are on a grid) using Fast_kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; - - // Init of an alpha complex from an OFF file - Gudhi::alpha_complex::Alpha_complex<Fast_kernel> alpha_complex_from_file(off_file_points); - - if (!alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) { - std::cerr << "Fast Alpha complex simplicial complex creation failed." << std::endl; - exit(-1); - } + stree = create_simplex_tree<Fast_kernel>(off_file_points, weight_file, exact_version, alpha_square_max_value); } else { using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>; - - // Init of an alpha complex from an OFF file - Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_points); - - if (!alpha_complex_from_file.create_complex(simplex, alpha_square_max_value, exact_version)) { - std::cerr << "Alpha complex simplicial complex creation failed." << std::endl; - exit(-1); - } + stree = create_simplex_tree<Kernel>(off_file_points, weight_file, exact_version, alpha_square_max_value); } // ---------------------------------------------------------------------------- // Display information about the alpha complex // ---------------------------------------------------------------------------- - std::clog << "Simplicial complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() - << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; + std::clog << "Simplicial complex is of dimension " << stree.dimension() << " - " << stree.num_simplices() + << " simplices - " << stree.num_vertices() << " vertices." << std::endl; - std::clog << "Simplex_tree dim: " << simplex.dimension() << std::endl; + std::clog << "Simplex_tree dim: " << stree.dimension() << std::endl; // Compute the persistence diagram of the complex Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh( - simplex); + stree); // initializes the coefficient field for homology pcoh.init_coefficients(coeff_field_characteristic); @@ -97,7 +147,7 @@ int main(int argc, char **argv) { } void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast, - std::string &output_file_diag, Filtration_value &alpha_square_max_value, + std::string &weight_file, std::string &output_file_diag, Filtration_value &alpha_square_max_value, int &coeff_field_characteristic, Filtration_value &min_persistence) { namespace po = boost::program_options; po::options_description hidden("Hidden options"); @@ -110,8 +160,10 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool "To activate exact version of Alpha complex (default is false, not available if fast is set)")( "fast,f", po::bool_switch(&fast), "To activate fast version of Alpha complex (default is false, not available if exact is set)")( + "weight-file,w", po::value<std::string>(&weight_file)->default_value(std::string()), + "Name of file containing a point weights. Format is one weight per line:\n W1\n ...\n Wn ")( "output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::clog")( + "Name of file in which the persistence diagram is written. Default print in standard output")( "max-alpha-square-value,r", po::value<Filtration_value>(&alpha_square_max_value) ->default_value(std::numeric_limits<Filtration_value>::infinity()), "Maximal alpha square value for the Alpha complex construction.")( diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md index 527598a9..1e3b8fab 100644 --- a/src/Alpha_complex/utilities/alphacomplex.md +++ b/src/Alpha_complex/utilities/alphacomplex.md @@ -46,6 +46,9 @@ for the Alpha complex construction. coefficient field Z/pZ for computing homology.
* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature
to be recorded. Enter a negative value to see zero length intervals.
+* `-w [ --weight-file ]` is the path to the file containing the weights of the
+points (one value per line).
+Default version is not weighted.
* `-e [ --exact ]` for the exact computation version.
* `-f [ --fast ]` for the fast computation version.
@@ -58,6 +61,10 @@ to be recorded. Enter a negative value to see zero length intervals. N.B.:
* Filtration values are alpha square values.
+* Weights values are explained on CGAL
+[dD Triangulations](https://doc.cgal.org/latest/Triangulation/index.html)
+and
+[Regular triangulation](https://doc.cgal.org/latest/Triangulation/index.html#TriangulationSecRT) documentation.
## alpha_complex_3d_persistence ##
@@ -124,6 +131,6 @@ N.B.: * `alpha_complex_3d_persistence` only accepts OFF files in dimension 3.
* Filtration values are alpha square values.
* Weights values are explained on CGAL
-[Alpha shape](https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0)
+[Alpha shape](https://doc.cgal.org/latest/Alpha_shapes_3/index.html#Alpha_shapes_3Definitions)
and
[Regular triangulation](https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation) documentation.
|