diff options
-rw-r--r-- | src/Alpha_complex/doc/Intro_alpha_complex.h | 15 | ||||
-rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_complex.h | 117 | ||||
-rw-r--r-- | src/Alpha_complex/test/CMakeLists.txt | 4 | ||||
-rw-r--r-- | src/Alpha_complex/test/Delaunay_complex_unit_test.cpp | 72 |
4 files changed, 146 insertions, 62 deletions
diff --git a/src/Alpha_complex/doc/Intro_alpha_complex.h b/src/Alpha_complex/doc/Intro_alpha_complex.h index a8b1a106..60da7169 100644 --- a/src/Alpha_complex/doc/Intro_alpha_complex.h +++ b/src/Alpha_complex/doc/Intro_alpha_complex.h @@ -46,16 +46,17 @@ namespace alpha_complex { * \cite cgal:s-gkd-19b from CGAL as template parameter. * * \remark - * - When an \f$\alpha\f$-complex is constructed with an infinite value of \f$ \alpha^2 \f$, the complex is a Delaunay - * complex (with special filtration values). + * - When the simplicial complex is constructed with an infinite value of \f$ \alpha^2 \f$, the complex is a Delaunay + * complex with special filtration values. The Delaunay complex without filtration values is also available by passing + * `default_filtration_value=true` to `Alpha_complex::create_complex`. * - For people only interested in the topology of the \ref alpha_complex (for instance persistence), * \ref alpha_complex is equivalent to the \ref cech_complex and much smaller if you do not bound the radii. * \ref cech_complex can still make sense in higher dimension precisely because you can bound the radii. - * - Using the default `CGAL::Epeck_d` makes the construction safe. If you pass exact=true to create_complex, the - * filtration values are the exact ones converted to the filtration value type of the simplicial complex. This can be - * very slow. If you pass exact=false (the default), the filtration values are only guaranteed to have a small - * multiplicative error compared to the exact value, see <code><a class="el" target="_blank" - * href="https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html"> + * - Using the default `CGAL::Epeck_d` makes the construction safe. If you pass `exact=true` to + * `Alpha_complex::create_complex`, the filtration values are the exact ones converted to the filtration value type of + * the simplicial complex. This can be very slow. If you pass `exact=false` (the default), the filtration values are + * only guaranteed to have a small multiplicative error compared to the exact value, see <code> + * <a class="el" target="_blank" href="https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html"> * CGAL::Lazy_exact_nt<NT>::set_relative_precision_of_to_double</a></code> for details. A drawback, when computing * persistence, is that an empty exact interval [10^12,10^12] may become a non-empty approximate interval * [10^12,10^12+10^6]. Using `CGAL::Epick_d` makes the computations slightly faster, and the combinatorics are still diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 0839ae6c..65c969d2 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -248,16 +248,20 @@ class Alpha_complex { public: /** \brief Inserts all Delaunay triangulation into the simplicial complex. - * It also computes the filtration values accordingly to the \ref createcomplexalgorithm + * It also computes the filtration values accordingly to the \ref createcomplexalgorithm if default_filtration_value + * is not set. * * \tparam SimplicialComplexForAlpha must meet `SimplicialComplexForAlpha` concept. * * @param[in] complex SimplicialComplexForAlpha to be created. * @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$, and there is very - * little point using anything else since it does not save time. + * little point using anything else since it does not save time. Useless if `default_filtration_value` is set to + * `true`. * @param[in] exact Exact filtration values computation. Not exact if `Kernel` is not <a target="_blank" * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>. - * + * @param[in] default_filtration_value Set this value to `true` if filtration values are not needed to be computed. + * Default value is `false` (which means compute the filtration values). + * * @return true if creation succeeds, false otherwise. * * @pre Delaunay triangulation must be already constructed with dimension strictly greater than 0. @@ -269,7 +273,8 @@ class Alpha_complex { typename Filtration_value = typename SimplicialComplexForAlpha::Filtration_value> bool create_complex(SimplicialComplexForAlpha& complex, Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(), - bool exact = false) { + 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; @@ -316,62 +321,64 @@ class Alpha_complex { } // -------------------------------------------------------------------------------------------- - // -------------------------------------------------------------------------------------------- - // Will be re-used many times - Vector_of_CGAL_points pointVector; - // ### For i : d -> 0 - for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) { - // ### Foreach Sigma of dim i - for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) { - int f_simplex_dim = complex.dimension(f_simplex); - if (decr_dim == f_simplex_dim) { - pointVector.clear(); -#ifdef DEBUG_TRACES - std::clog << "Sigma of dim " << decr_dim << " is"; -#endif // DEBUG_TRACES - for (auto vertex : complex.simplex_vertex_range(f_simplex)) { - pointVector.push_back(get_point(vertex)); -#ifdef DEBUG_TRACES - std::clog << " " << vertex; -#endif // DEBUG_TRACES - } -#ifdef DEBUG_TRACES - std::clog << std::endl; -#endif // DEBUG_TRACES - // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma) - if (std::isnan(complex.filtration(f_simplex))) { - Filtration_value alpha_complex_filtration = 0.0; - // No need to compute squared_radius on a single point - alpha is 0.0 - if (f_simplex_dim > 0) { - // squared_radius function initialization - Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object(); - - CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv; - auto sqrad = squared_radius(pointVector.begin(), pointVector.end()); -#if CGAL_VERSION_NR >= 1050000000 - if(exact) CGAL::exact(sqrad); -#endif - alpha_complex_filtration = cv(sqrad); + if (!default_filtration_value) { + // -------------------------------------------------------------------------------------------- + // Will be re-used many times + Vector_of_CGAL_points pointVector; + // ### For i : d -> 0 + for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) { + // ### Foreach Sigma of dim i + for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) { + int f_simplex_dim = complex.dimension(f_simplex); + if (decr_dim == f_simplex_dim) { + pointVector.clear(); + #ifdef DEBUG_TRACES + std::clog << "Sigma of dim " << decr_dim << " is"; + #endif // DEBUG_TRACES + for (auto vertex : complex.simplex_vertex_range(f_simplex)) { + pointVector.push_back(get_point(vertex)); + #ifdef DEBUG_TRACES + std::clog << " " << vertex; + #endif // DEBUG_TRACES } - complex.assign_filtration(f_simplex, alpha_complex_filtration); -#ifdef DEBUG_TRACES - std::clog << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl; -#endif // DEBUG_TRACES + #ifdef DEBUG_TRACES + std::clog << std::endl; + #endif // DEBUG_TRACES + // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma) + if (std::isnan(complex.filtration(f_simplex))) { + Filtration_value alpha_complex_filtration = 0.0; + // No need to compute squared_radius on a single point - alpha is 0.0 + if (f_simplex_dim > 0) { + // squared_radius function initialization + Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object(); + + CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv; + auto sqrad = squared_radius(pointVector.begin(), pointVector.end()); + #if CGAL_VERSION_NR >= 1050000000 + if(exact) CGAL::exact(sqrad); + #endif + alpha_complex_filtration = cv(sqrad); + } + complex.assign_filtration(f_simplex, alpha_complex_filtration); + #ifdef DEBUG_TRACES + std::clog << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl; + #endif // DEBUG_TRACES + } + // No need to propagate further, unweighted points all have value 0 + if (decr_dim > 1) + propagate_alpha_filtration(complex, f_simplex); } - // No need to propagate further, unweighted points all have value 0 - if (decr_dim > 1) - propagate_alpha_filtration(complex, f_simplex); } } + // -------------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------------- + // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension + complex.make_filtration_non_decreasing(); + // Remove all simplices that have a filtration value greater than max_alpha_square + complex.prune_above_filtration(max_alpha_square); + // -------------------------------------------------------------------------------------------- } - // -------------------------------------------------------------------------------------------- - - // -------------------------------------------------------------------------------------------- - // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension - complex.make_filtration_non_decreasing(); - // Remove all simplices that have a filtration value greater than max_alpha_square - complex.prune_above_filtration(max_alpha_square); - // -------------------------------------------------------------------------------------------- return true; } diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt index 0476c6d4..fe4b23e4 100644 --- a/src/Alpha_complex/test/CMakeLists.txt +++ b/src/Alpha_complex/test/CMakeLists.txt @@ -8,11 +8,15 @@ 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 ( 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(Delaunay_complex_test_unit ${TBB_LIBRARIES}) endif() gudhi_add_boost_test(Alpha_complex_test_unit) + gudhi_add_boost_test(Delaunay_complex_test_unit) add_executable ( Alpha_complex_3d_test_unit Alpha_complex_3d_unit_test.cpp ) target_link_libraries(Alpha_complex_3d_test_unit ${CGAL_LIBRARY}) diff --git a/src/Alpha_complex/test/Delaunay_complex_unit_test.cpp b/src/Alpha_complex/test/Delaunay_complex_unit_test.cpp new file mode 100644 index 00000000..71164705 --- /dev/null +++ b/src/Alpha_complex/test/Delaunay_complex_unit_test.cpp @@ -0,0 +1,72 @@ +/* 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 "delaunay_complex" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <CGAL/Epick_d.h> +#include <CGAL/Epeck_d.h> + +#include <vector> +#include <limits> // NaN +#include <cmath> + +#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> +#include <gudhi/random_point_generators.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<5> > 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<5> > 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) { + std::cout << "*****************************************************************************************************"; + using Point = typename TestedKernel::Point_d; + std::vector<Point> points; + // 50 points on a 4-sphere + points = Gudhi::generate_points_on_sphere_d<TestedKernel>(10, 5, 1.); + + Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex(points); + + // Alpha complex + Gudhi::Simplex_tree<> stree_from_alpha_complex; + BOOST_CHECK(alpha_complex.create_complex(stree_from_alpha_complex)); + stree_from_alpha_complex.initialize_filtration(); + + // Delaunay complex + Gudhi::Simplex_tree<> stree_from_delaunay_complex; + BOOST_CHECK(alpha_complex.create_complex(stree_from_delaunay_complex, 0., false, true)); + + // Check all the simplices from alpha complex are in the Delaunay complex + for (auto f_simplex : stree_from_alpha_complex.filtration_simplex_range()) { + std::vector<Gudhi::Simplex_tree<>::Vertex_handle> simplex; + for (Gudhi::Simplex_tree<>::Vertex_handle vertex : stree_from_alpha_complex.simplex_vertex_range(f_simplex)) { + std::cout << "(" << vertex << ")"; + simplex.push_back(vertex); + } + std::cout << std::endl; + Gudhi::Simplex_tree<>::Simplex_handle sh = stree_from_delaunay_complex.find(simplex); + BOOST_CHECK(std::isnan(stree_from_delaunay_complex.filtration(sh))); + BOOST_CHECK(sh != stree_from_delaunay_complex.null_simplex()); + } +} |