From e8cfcc7017f9002a229996232fdb2d03a41d7ea7 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Fri, 21 Aug 2020 15:40:17 +0200 Subject: Fix compilation issues when cgal but no eigen --- src/Alpha_complex/test/CMakeLists.txt | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) (limited to 'src/Alpha_complex/test/CMakeLists.txt') diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt index fe4b23e4..f38bd096 100644 --- a/src/Alpha_complex/test/CMakeLists.txt +++ b/src/Alpha_complex/test/CMakeLists.txt @@ -18,6 +18,10 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) gudhi_add_boost_test(Alpha_complex_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 +42,4 @@ 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) -- cgit v1.2.3 From db938dbd74460e7a0fd705be8628984052f71dc0 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Fri, 4 Sep 2020 17:23:16 +0200 Subject: Alpha_kernel_d, its tests and Alpha complex to use it --- src/Alpha_complex/example/CMakeLists.txt | 3 + .../example/Weighted_alpha_complex_from_points.cpp | 54 ++++++++ src/Alpha_complex/include/gudhi/Alpha_complex.h | 106 ++++++++-------- .../include/gudhi/Alpha_complex/Alpha_kernel_d.h | 137 +++++++++++++++++++++ .../test/Alpha_kernel_d_unit_test.cpp | 132 ++++++++++++++++++++ src/Alpha_complex/test/CMakeLists.txt | 4 + 6 files changed, 387 insertions(+), 49 deletions(-) create mode 100644 src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp create mode 100644 src/Alpha_complex/include/gudhi/Alpha_complex/Alpha_kernel_d.h create mode 100644 src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp (limited to 'src/Alpha_complex/test/CMakeLists.txt') diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt index 2eecd50c..c9b937ed 100644 --- a/src/Alpha_complex/example/CMakeLists.txt +++ b/src/Alpha_complex/example/CMakeLists.txt @@ -7,10 +7,13 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) target_link_libraries(Alpha_complex_example_from_off ${CGAL_LIBRARY}) add_executable ( Alpha_complex_example_fast_from_off Fast_alpha_complex_from_off.cpp ) target_link_libraries(Alpha_complex_example_fast_from_off ${CGAL_LIBRARY}) + 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(Alpha_complex_example_from_points ${TBB_LIBRARIES}) target_link_libraries(Alpha_complex_example_from_off ${TBB_LIBRARIES}) target_link_libraries(Alpha_complex_example_fast_from_off ${TBB_LIBRARIES}) + target_link_libraries(Weighted_alpha_complex_example_from_points ${TBB_LIBRARIES}) endif() add_test(NAME Alpha_complex_example_from_points COMMAND $) 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..19a04282 --- /dev/null +++ b/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp @@ -0,0 +1,54 @@ +#include +// to construct a simplex_tree from alpha complex +#include + +#include + +#include +#include + +// Explicit dimension 2 Epeck_d kernel +using Kernel = CGAL::Epeck_d< CGAL::Dimension_tag<2> >; +using Bare_point = Kernel::Point_d; +using Weighted_point = Kernel::Weighted_point_d; +using Vector_of_points = std::vector; + +int main() { + // ---------------------------------------------------------------------------- + // Init of a list of points + // ---------------------------------------------------------------------------- + Vector_of_points points; + points.push_back(Weighted_point(Bare_point(1.0, 1.0) , 1.)); + points.push_back(Weighted_point(Bare_point(7.0, 0.0) , 1.)); + points.push_back(Weighted_point(Bare_point(4.0, 6.0) , 1.)); + points.push_back(Weighted_point(Bare_point(9.0, 6.0) , 1.)); + points.push_back(Weighted_point(Bare_point(0.0, 14.0), 1.)); + points.push_back(Weighted_point(Bare_point(2.0, 19.0), 1.)); + points.push_back(Weighted_point(Bare_point(9.0, 17.0), 1.)); + + // ---------------------------------------------------------------------------- + // Init of an alpha complex from the list of points + // ---------------------------------------------------------------------------- + Gudhi::alpha_complex::Alpha_complex 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/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index ba91998d..8e9fe773 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -12,6 +12,7 @@ #ifndef ALPHA_COMPLEX_H_ #define ALPHA_COMPLEX_H_ +#include #include // to construct Alpha_complex from a OFF file of points #include @@ -20,6 +21,7 @@ #include // isnan, fmax #include +#include // aka. Weighted Delaunay triangulation #include // For EXACT or SAFE version #include // For FAST version #include @@ -91,49 +93,55 @@ template struct Is_Epeck_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> +template, bool Weighted = false> class Alpha_complex { public: + /** \brief Geometric traits class that provides the geometric types and predicates needed by the triangulations.*/ + using Geom_traits = typename std::conditional, + Kernel>::type; // Add an int in TDS to save point index in the structure - typedef CGAL::Triangulation_data_structure, - CGAL::Triangulation_full_cell > TDS; - /** \brief A Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/ - typedef CGAL::Delaunay_triangulation Delaunay_triangulation; + using TDS = CGAL::Triangulation_data_structure, + CGAL::Triangulation_full_cell >; + /** \brief A (Weighted or not) Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/ + using Triangulation = typename std::conditional, + CGAL::Delaunay_triangulation>::type; + + using A_kernel_d = Alpha_kernel_d; + + // Numeric type of coordinates in the kernel + using FT = typename A_kernel_d::FT; + + /** \brief If Weighted, the weighted point is cached (point + weight [= squared radius]), + * else a pair of point and squared radius is cached. + */ + using Sphere = typename A_kernel_d::Sphere; /** \brief A point in Euclidean space.*/ - typedef typename Kernel::Point_d Point_d; - /** \brief Geometric traits class that provides the geometric types and predicates needed by Delaunay - * triangulations.*/ - typedef Kernel Geom_traits; + using Point_d = typename std::conditional::type; private: - typedef typename Kernel::Compute_squared_radius_d Squared_Radius; - typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel; - typedef typename Kernel::Point_dimension_d Point_Dimension; - // Vertex_iterator type from CGAL. - typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator; + using CGAL_vertex_iterator = typename Triangulation::Vertex_iterator; // size_type type from CGAL. - typedef typename Delaunay_triangulation::size_type size_type; + using size_type = typename Triangulation::size_type; // Structure to switch from simplex tree vertex handle to CGAL vertex iterator. - typedef typename std::vector< CGAL_vertex_iterator > Vector_vertex_iterator; - - // Numeric type of coordinates in the kernel - typedef typename Kernel::FT FT; + using Vector_vertex_iterator = typename std::vector< CGAL_vertex_iterator >; private: /** \brief Vertex iterator vector to switch from simplex tree vertex handle to CGAL vertex iterator. * Vertex handles are inserted sequentially, starting at 0.*/ Vector_vertex_iterator vertex_handle_to_iterator_; /** \brief Pointer on the CGAL Delaunay triangulation.*/ - Delaunay_triangulation* triangulation_; + Triangulation* triangulation_; /** \brief Kernel for triangulation_ functions access.*/ - Kernel kernel_; + A_kernel_d kernel_; + /** \brief Cache for geometric constructions: circumcenter and squared radius of a simplex.*/ - std::vector> cache_, old_cache_; + std::vector cache_, old_cache_; public: /** \brief Alpha_complex constructor from an OFF file name. @@ -160,10 +168,10 @@ class Alpha_complex { * * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous. * - * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d + * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d or Kernel::Weighted_point_d. * * The type InputPointRange must be a range for which std::begin and - * std::end return input iterators on a Kernel::Point_d. + * std::end return input iterators on a Kernel::Point_d or Kernel::Weighted_point_d. */ template Alpha_complex(const InputPointRange& points) @@ -206,11 +214,8 @@ class Alpha_complex { auto last = std::end(points); if (first != last) { - // point_dimension function initialization - Point_Dimension point_dimension = kernel_.point_dimension_d_object(); - - // Delaunay triangulation is point dimension. - triangulation_ = new Delaunay_triangulation(point_dimension(*first)); + // Delaunay triangulation init with point dimension. + triangulation_ = new Triangulation(kernel_.get_dimension(*first)); std::vector point_cloud(first, last); @@ -218,18 +223,22 @@ class Alpha_complex { std::vector indices(boost::counting_iterator(0), boost::counting_iterator(point_cloud.size())); - typedef boost::iterator_property_map::iterator, - CGAL::Identity_property_map> Point_property_map; - typedef CGAL::Spatial_sort_traits_adapter_d Search_traits_d; + using Point_property_map = boost::iterator_property_map::iterator, + CGAL::Identity_property_map>; + using Search_traits_d = CGAL::Spatial_sort_traits_adapter_d; CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud))); - typename Delaunay_triangulation::Full_cell_handle hint; + typename Triangulation::Full_cell_handle hint; for (auto index : indices) { - typename Delaunay_triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint); - // Save index value as data to retrieve it after insertion - pos->data() = index; - hint = pos->full_cell(); + typename Triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint); + if (pos != nullptr) { + // Save index value as data to retrieve it after insertion + pos->data() = index; + hint = pos->full_cell(); + } else { + std::cout << "NULLPTR" << std::endl; + } } // -------------------------------------------------------------------------------------------- // structure to retrieve CGAL points from vertex handle - one vertex handle per point. @@ -270,9 +279,7 @@ class Alpha_complex { v.clear(); for (auto vertex : cplx.simplex_vertex_range(s)) v.push_back(get_point_(vertex)); - Point_d c = kernel_.construct_circumcenter_d_object()(v.cbegin(), v.cend()); - FT r = kernel_.squared_distance_d_object()(c, v[0]); - cache_.emplace_back(std::move(c), std::move(r)); + cache_.emplace_back(kernel_.get_sphere(v.cbegin(), v.cend())); } return cache_[k]; } @@ -282,13 +289,13 @@ class Alpha_complex { auto radius(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) { auto k = cplx.key(s); if(k!=cplx.null_key()) - return old_cache_[k].second; + return kernel_.get_squared_radius(old_cache_[k]); // Using a transform_range is slower, currently. thread_local std::vector v; v.clear(); for (auto vertex : cplx.simplex_vertex_range(s)) v.push_back(get_point_(vertex)); - return kernel_.compute_squared_radius_d_object()(v.cbegin(), v.cend()); + return kernel_.get_squared_radius(v.cbegin(), v.cend()); } public: @@ -322,9 +329,9 @@ class Alpha_complex { bool exact = false, bool default_filtration_value = false) { // From SimplicialComplexForAlpha type required to insert into a simplicial complex (with or without subfaces). - typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle; - typedef typename SimplicialComplexForAlpha::Simplex_handle Simplex_handle; - typedef std::vector Vector_vertex; + using Vertex_handle = typename SimplicialComplexForAlpha::Vertex_handle; + using Simplex_handle = typename SimplicialComplexForAlpha::Simplex_handle; + using Vector_vertex = std::vector; if (triangulation_ == nullptr) { std::cerr << "Alpha_complex cannot create_complex from a NULL triangulation\n"; @@ -416,8 +423,8 @@ class Alpha_complex { template void propagate_alpha_filtration(SimplicialComplexForAlpha& complex, Simplex_handle f_simplex) { // From SimplicialComplexForAlpha type required to assign filtration values. - typedef typename SimplicialComplexForAlpha::Filtration_value Filtration_value; - typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle; + using Filtration_value = typename SimplicialComplexForAlpha::Filtration_value; + using Vertex_handle = typename SimplicialComplexForAlpha::Vertex_handle; // ### Foreach Tau face of Sigma for (auto f_boundary : complex.boundary_simplex_range(f_simplex)) { @@ -450,7 +457,8 @@ class Alpha_complex { while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; } Vertex_handle extra = *longiter; auto const& cache=get_cache(complex, f_boundary); - bool is_gab = kernel_.squared_distance_d_object()(cache.first, get_point_(extra)) >= cache.second; + bool is_gab = kernel_.get_squared_distance(kernel_.get_circumcenter(cache), get_point_(extra)) >= + kernel_.get_squared_radius(cache); #ifdef DEBUG_TRACES std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << extra << std::endl; #endif // DEBUG_TRACES 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..87604ec4 --- /dev/null +++ b/src/Alpha_complex/include/gudhi/Alpha_complex/Alpha_kernel_d.h @@ -0,0 +1,137 @@ +/* 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 // For EXACT or SAFE version +#include // For FAST version +#include // for CGAL_VERSION_NR + +#include // for EIGEN_VERSION_AT_LEAST + +#include // 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 { + +template < typename Kernel, bool Weighted = false > +class Alpha_kernel_d { +}; + +// Unweighted Kernel_d version +template < typename Kernel > +class Alpha_kernel_d { + 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; + + int get_dimension(const Point_d& p0) const { + return kernel_.point_dimension_d_object()(p0); + } + + template + 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 + 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; + } + + template + Point get_circumcenter(const Sphere& sph) const { + return sph.first; + } + + template + FT get_squared_distance(const Point& first, const Point& second) const { + return kernel_.squared_distance_d_object()(first, second); + } +}; + +// Weighted Kernel_d version +template < typename Kernel > +class Alpha_kernel_d { + 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 + Sphere get_sphere(PointIterator begin, PointIterator end) const { + return kernel_.power_center_d_object()(begin, end); + } + + template + FT get_squared_radius(PointIterator begin, PointIterator end) const { + Sphere sph = get_sphere(begin, end); + return sph.weight(); + } + + FT get_squared_radius(const Sphere& sph) const { + return sph.weight(); + } + + template + Point get_circumcenter(const Sphere& sph) const { + return sph.point(); + } + + template + FT get_squared_distance(const Point& first, const Point& second) const { + return kernel_.power_distance_d_object()(first, second); + } +}; + +} // 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/test/Alpha_kernel_d_unit_test.cpp b/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp new file mode 100644 index 00000000..70b65aad --- /dev/null +++ b/src/Alpha_complex/test/Alpha_kernel_d_unit_test.cpp @@ -0,0 +1,132 @@ +/* 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 +#include + +#include +#include +#include + +#include +#include +#include // for std::pair + +#include +#include + +// 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 list_of_kernel_variants; + +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_kernel_d_dimension, TestedKernel, list_of_kernel_variants) { + Gudhi::alpha_complex::Alpha_kernel_d kernel; + std::vector 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 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) { + using Unweighted_kernel = Gudhi::alpha_complex::Alpha_kernel_d; + // Sphere: (x-1)² + (y-1)² + z² + t² = 1 + // At least 5 points for a 3-sphere + std::vector p0 {1., 0., 0., 0.}; + std::vector p1 {0., 1., 0., 0.}; + std::vector p2 {1., 1., 1., 0.}; + std::vector p3 {1., 1., 0., 1.}; + std::vector p4 {1., 1., -1., 0.}; + + using Point_d = typename Unweighted_kernel::Point_d; + std::vector unw_pts; + unw_pts.push_back(Point_d(p0.begin(), p0.end())); + unw_pts.push_back(Point_d(p1.begin(), p1.end())); + unw_pts.push_back(Point_d(p2.begin(), p2.end())); + unw_pts.push_back(Point_d(p3.begin(), p3.end())); + unw_pts.push_back(Point_d(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; + + using Weighted_point_d = typename Weighted_kernel::Weighted_point_d; + using Bare_point_d = typename Weighted_kernel::Bare_point_d; + std::vector w_pts; + w_pts.push_back(Weighted_point_d(Bare_point_d(p0.begin(), p0.end()), 0.)); + w_pts.push_back(Weighted_point_d(Bare_point_d(p1.begin(), p1.end()), 0.)); + w_pts.push_back(Weighted_point_d(Bare_point_d(p2.begin(), p2.end()), 0.)); + w_pts.push_back(Weighted_point_d(Bare_point_d(p3.begin(), p3.end()), 0.)); + w_pts.push_back(Weighted_point_d(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 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)); +} + + +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_kernel_d_distance, TestedKernel, list_of_kernel_variants) { + using Unweighted_kernel = Gudhi::alpha_complex::Alpha_kernel_d; + + std::vector p0 {1., 0., 0., 0.}; + std::vector p1 {0., 1., 0., 0.}; + + using Point_d = typename Unweighted_kernel::Point_d; + Unweighted_kernel kernel; + auto dist_01 = kernel.get_squared_distance(Point_d(p0.begin(), p0.end()), Point_d(p1.begin(), p1.end())); + std::clog << "Distance is " << dist_01 << std::endl; + + using Weighted_kernel = Gudhi::alpha_complex::Alpha_kernel_d; + + using Weighted_point_d = typename Weighted_kernel::Weighted_point_d; + using Bare_point_d = typename Weighted_kernel::Bare_point_d; + std::vector w_pts; + + Weighted_kernel w_kernel; + auto w_dist_01 = w_kernel.get_squared_distance(Weighted_point_d(Bare_point_d(p0.begin(), p0.end()), 0.), + Weighted_point_d(Bare_point_d(p1.begin(), p1.end()), 0.)); + std::clog << "Distance is " << w_dist_01 << std::endl; + + CGAL::NT_converter cast_to_double; + // The results shall be the same with weights = 0. + GUDHI_TEST_FLOAT_EQUALITY_CHECK(cast_to_double(dist_01), cast_to_double(w_dist_01)); +} diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt index fe4b23e4..bd61c0d8 100644 --- a/src/Alpha_complex/test/CMakeLists.txt +++ b/src/Alpha_complex/test/CMakeLists.txt @@ -10,13 +10,17 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) 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}) + 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_complex_test_unit ${TBB_LIBRARIES}) target_link_libraries(Delaunay_complex_test_unit ${TBB_LIBRARIES}) + target_link_libraries(Alpha_kernel_d_test_unit ${TBB_LIBRARIES}) endif() gudhi_add_boost_test(Alpha_complex_test_unit) gudhi_add_boost_test(Delaunay_complex_test_unit) + gudhi_add_boost_test(Alpha_kernel_d_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}) -- cgit v1.2.3 From 171ddab9b7a50f0303d7201fa547dbfb445f9698 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Fri, 11 Sep 2020 11:55:55 +0200 Subject: Requires CGAL >= 5.1. Fix is_gabriel computation --- src/Alpha_complex/example/CMakeLists.txt | 14 +++++++++----- .../example/Weighted_alpha_complex_from_points.cpp | 14 +++++++------- src/Alpha_complex/include/gudhi/Alpha_complex.h | 8 ++++++-- .../include/gudhi/Alpha_complex/Alpha_kernel_d.h | 17 ++++------------- src/Alpha_complex/test/CMakeLists.txt | 13 +++++++++---- 5 files changed, 35 insertions(+), 31 deletions(-) (limited to 'src/Alpha_complex/test/CMakeLists.txt') diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt index 17dc896c..1fc2330a 100644 --- a/src/Alpha_complex/example/CMakeLists.txt +++ b/src/Alpha_complex/example/CMakeLists.txt @@ -7,13 +7,10 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) target_link_libraries(Alpha_complex_example_from_off ${CGAL_LIBRARY}) add_executable ( Alpha_complex_example_fast_from_off Fast_alpha_complex_from_off.cpp ) target_link_libraries(Alpha_complex_example_fast_from_off ${CGAL_LIBRARY}) - 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(Alpha_complex_example_from_points ${TBB_LIBRARIES}) target_link_libraries(Alpha_complex_example_from_off ${TBB_LIBRARIES}) target_link_libraries(Alpha_complex_example_fast_from_off ${TBB_LIBRARIES}) - target_link_libraries(Weighted_alpha_complex_example_from_points ${TBB_LIBRARIES}) endif() add_test(NAME Alpha_complex_example_from_points COMMAND $) @@ -28,8 +25,6 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) add_test(NAME Alpha_complex_example_fast_from_off_32 COMMAND $ "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "32.0" "${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_32.txt") - add_test(NAME Weighted_alpha_complex_example_from_points COMMAND $) - 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}/) @@ -69,3 +64,12 @@ if (NOT CGAL_VERSION VERSION_LESS 4.11.0) COMMAND $) 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 $) +endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) diff --git a/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp b/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp index 05858084..19a04282 100644 --- a/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp +++ b/src/Alpha_complex/example/Weighted_alpha_complex_from_points.cpp @@ -18,13 +18,13 @@ int main() { // Init of a list of points // ---------------------------------------------------------------------------- Vector_of_points points; - points.push_back(Weighted_point(Bare_point(1.0, 1.0) , 0.)); - points.push_back(Weighted_point(Bare_point(7.0, 0.0) , 0.)); - points.push_back(Weighted_point(Bare_point(4.0, 6.0) , 0.)); - points.push_back(Weighted_point(Bare_point(9.0, 6.0) , 0.)); - points.push_back(Weighted_point(Bare_point(0.0, 14.0), 0.)); - points.push_back(Weighted_point(Bare_point(2.0, 19.0), 0.)); - points.push_back(Weighted_point(Bare_point(9.0, 17.0), 0.)); + points.push_back(Weighted_point(Bare_point(1.0, 1.0) , 1.)); + points.push_back(Weighted_point(Bare_point(7.0, 0.0) , 1.)); + points.push_back(Weighted_point(Bare_point(4.0, 6.0) , 1.)); + points.push_back(Weighted_point(Bare_point(9.0, 6.0) , 1.)); + points.push_back(Weighted_point(Bare_point(0.0, 14.0), 1.)); + points.push_back(Weighted_point(Bare_point(2.0, 19.0), 1.)); + points.push_back(Weighted_point(Bare_point(9.0, 17.0), 1.)); // ---------------------------------------------------------------------------- // Init of an alpha complex from the list of points diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 8e9fe773..4316d9bc 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -210,6 +210,11 @@ 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); @@ -457,8 +462,7 @@ class Alpha_complex { while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; } Vertex_handle extra = *longiter; auto const& cache=get_cache(complex, f_boundary); - bool is_gab = kernel_.get_squared_distance(kernel_.get_circumcenter(cache), get_point_(extra)) >= - kernel_.get_squared_radius(cache); + bool is_gab = kernel_.is_gabriel(cache, get_point_(extra)); #ifdef DEBUG_TRACES std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << extra << std::endl; #endif // DEBUG_TRACES 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 index a3e3845a..b64e4f59 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex/Alpha_kernel_d.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex/Alpha_kernel_d.h @@ -71,13 +71,8 @@ class Alpha_kernel_d { return sph.second; } - auto get_circumcenter(const Sphere& sph) const { - return sph.first; - } - - template - FT get_squared_distance(const Point& first, const Point& second) const { - return kernel_.squared_distance_d_object()(first, second); + bool is_gabriel(const Sphere& circumcenter, const Point_d& point) { + return kernel_.squared_distance_d_object()(circumcenter.first, point) >= circumcenter.second; } }; @@ -116,12 +111,8 @@ class Alpha_kernel_d { return sph.weight(); } - auto get_circumcenter(const Sphere& sph) const { - return sph; - } - - FT get_squared_distance(const Weighted_point_d& first, const Weighted_point_d& second) const { - return kernel_.power_distance_d_object()(first, second); + bool is_gabriel(const Sphere& circumcenter, const Weighted_point_d& point) { + return kernel_.power_distance_d_object()(circumcenter, point) >= 0; } }; diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt index 71e4ea7c..837d2948 100644 --- a/src/Alpha_complex/test/CMakeLists.txt +++ b/src/Alpha_complex/test/CMakeLists.txt @@ -10,17 +10,13 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) 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}) - 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_complex_test_unit ${TBB_LIBRARIES}) target_link_libraries(Delaunay_complex_test_unit ${TBB_LIBRARIES}) - target_link_libraries(Alpha_kernel_d_test_unit ${TBB_LIBRARIES}) endif() gudhi_add_boost_test(Alpha_complex_test_unit) gudhi_add_boost_test(Delaunay_complex_test_unit) - gudhi_add_boost_test(Alpha_kernel_d_test_unit) endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) @@ -47,3 +43,12 @@ if (NOT CGAL_VERSION VERSION_LESS 4.11.0) gudhi_add_boost_test(Weighted_periodic_alpha_complex_3d_test_unit) 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) +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) \ No newline at end of file -- cgit v1.2.3 From 54e0cec31e556c88d9f8107394ad972883fdbbdf Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Fri, 11 Sep 2020 17:41:39 +0200 Subject: Add weighted ctor and its test --- src/Alpha_complex/include/gudhi/Alpha_complex.h | 26 +++++++ src/Alpha_complex/test/CMakeLists.txt | 8 +++ .../test/Weighted_alpha_complex_unit_test.cpp | 82 ++++++++++++++++++++++ 3 files changed, 116 insertions(+) create mode 100644 src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp (limited to 'src/Alpha_complex/test/CMakeLists.txt') diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 4316d9bc..00231e1c 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -31,6 +31,10 @@ #include // for EIGEN_VERSION_AT_LEAST +#include +#include +#include + #include #include #include @@ -179,6 +183,28 @@ class Alpha_complex { init_from_range(points); } + /** \brief Alpha_complex constructor from a list of points and weights. + * + * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous. + * + * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d or Kernel::Weighted_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. + */ + template + Alpha_complex(const InputPointRange& points, WeightRange weights) { + static_assert(Weighted, "This constructor is not available for non-weighted versions of Alpha_complex_3d"); + // 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); + } + /** \brief Alpha_complex destructor deletes the Delaunay triangulation. */ ~Alpha_complex() { diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt index 837d2948..db5d840f 100644 --- a/src/Alpha_complex/test/CMakeLists.txt +++ b/src/Alpha_complex/test/CMakeLists.txt @@ -51,4 +51,12 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) 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) + endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) \ 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..b4fc76de --- /dev/null +++ b/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp @@ -0,0 +1,82 @@ +/* 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 +#include + +#include +#include + +#include // float comparison +#include +#include + +#include +#include +#include + +// 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 list_of_kernel_variants; + +BOOST_AUTO_TEST_CASE_TEMPLATE(Zero_weighted_alpha_complex, Kernel, list_of_kernel_variants) { + // Random points construction + using Point_d = typename Kernel::Point_d; + std::vector points; + std::uniform_real_distribution rd_pts(-10., 10.); + std::random_device rand_dev; + std::mt19937 rand_engine(rand_dev()); + for (int idx = 0; idx < 20; idx++) { + std::vector point {rd_pts(rand_engine), rd_pts(rand_engine), rd_pts(rand_engine), rd_pts(rand_engine)}; + points.emplace_back(Point_d(point.begin(), point.end())); + } + + // Alpha complex from points + Gudhi::alpha_complex::Alpha_complex alpha_complex_from_points(points); + Gudhi::Simplex_tree<> simplex; + BOOST_CHECK(alpha_complex_from_points.create_complex(simplex)); + 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::clog << std::endl; + } + + // Alpha complex from zero weighted points + std::vector weights(20, 0.); + Gudhi::alpha_complex::Alpha_complex 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)); + + 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::clog << std::endl; + } + + BOOST_CHECK(zw_simplex == simplex); +} \ No newline at end of file -- cgit v1.2.3