diff options
-rw-r--r-- | src/Cech_complex/benchmark/CMakeLists.txt | 18 | ||||
-rw-r--r-- | src/Cech_complex/benchmark/cech_complex_benchmark.cpp | 27 | ||||
-rw-r--r-- | src/Cech_complex/concept/SimplicialComplexForCech.h | 4 | ||||
-rw-r--r-- | src/Cech_complex/example/CMakeLists.txt | 26 | ||||
-rw-r--r-- | src/Cech_complex/example/cech_complex_example_from_points.cpp | 47 | ||||
-rw-r--r-- | src/Cech_complex/example/cech_complex_step_by_step.cpp | 28 | ||||
-rw-r--r-- | src/Cech_complex/include/gudhi/Cech_complex.h | 45 | ||||
-rw-r--r-- | src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h | 107 | ||||
-rw-r--r-- | src/Cech_complex/include/gudhi/Cech_complex_blocker.h | 98 | ||||
-rw-r--r-- | src/Cech_complex/test/CMakeLists.txt | 19 | ||||
-rw-r--r-- | src/Cech_complex/test/test_cech_complex.cpp | 73 | ||||
-rw-r--r-- | src/Cech_complex/utilities/CMakeLists.txt | 24 | ||||
-rw-r--r-- | src/Cech_complex/utilities/cech_persistence.cpp | 10 | ||||
-rw-r--r-- | src/common/include/gudhi/distance_functions.h | 49 |
14 files changed, 382 insertions, 193 deletions
diff --git a/src/Cech_complex/benchmark/CMakeLists.txt b/src/Cech_complex/benchmark/CMakeLists.txt index bc54c0f3..a6b3d70b 100644 --- a/src/Cech_complex/benchmark/CMakeLists.txt +++ b/src/Cech_complex/benchmark/CMakeLists.txt @@ -1,13 +1,15 @@ project(Cech_complex_benchmark) -# Do not forget to copy test files in current binary dir -file(COPY "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.0.1) + # Do not forget to copy test files in current binary dir + file(COPY "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) -if(TARGET Boost::filesystem) - add_executable(cech_complex_benchmark cech_complex_benchmark.cpp) - target_link_libraries(cech_complex_benchmark Boost::filesystem) + if(TARGET Boost::filesystem) + add_executable(cech_complex_benchmark cech_complex_benchmark.cpp) + target_link_libraries(cech_complex_benchmark Boost::filesystem) - if (TBB_FOUND) - target_link_libraries(cech_complex_benchmark ${TBB_LIBRARIES}) + if (TBB_FOUND) + target_link_libraries(cech_complex_benchmark ${TBB_LIBRARIES}) + endif() endif() -endif()
\ No newline at end of file +endif() diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 2e4adce4..94c5fa4f 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -17,6 +17,8 @@ #include <gudhi/Simplex_tree.h> #include <gudhi/Miniball.hpp> +#include <CGAL/Epick_d.h> + #include "boost/filesystem.hpp" // includes all needed Boost.Filesystem declarations #include <string> @@ -30,7 +32,11 @@ using Point_cloud = std::vector<Point>; using Points_off_reader = Gudhi::Points_off_reader<Point>; using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>; using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; -using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>; +using Kernel = CGAL::Epick_d<CGAL::Dimension_tag<3>>; +using Point_cgal = typename Kernel::Point_d; +using Point_cloud_cgal = std::vector<Point_cgal>; +using Points_off_reader_cgal = Gudhi::Points_off_reader<Point_cgal>; +using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud_cgal, Kernel, Simplex_tree>; class Minimal_enclosing_ball_radius { public: @@ -62,6 +68,7 @@ int main(int argc, char* argv[]) { // Extract the points from the file filepoints Points_off_reader off_reader(off_file_points); + Points_off_reader_cgal off_reader_cgal(off_file_points); Gudhi::Clock euclidean_clock("Gudhi::Euclidean_distance"); // Compute the proximity graph of the points @@ -76,16 +83,16 @@ int main(int argc, char* argv[]) { off_reader.get_point_cloud(), threshold, Minimal_enclosing_ball_radius()); std::clog << miniball_clock << std::endl; - Gudhi::Clock common_miniball_clock("Gudhi::Minimal_enclosing_ball_radius()"); + Gudhi::Clock cgal_miniball_clock("Gudhi::Minimal_enclosing_ball_radius_cgal()"); // Compute the proximity graph of the points - Proximity_graph common_miniball_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>( - off_reader.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius()); - std::clog << common_miniball_clock << std::endl; + Proximity_graph cgal_miniball_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>( + off_reader_cgal.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius<Kernel>()); + std::clog << cgal_miniball_clock << std::endl; boost::filesystem::path full_path(boost::filesystem::current_path()); std::clog << "Current path is : " << full_path << std::endl; - std::clog << "File name;Radius;Rips time;Cech time; Ratio Rips/Cech time;Rips nb simplices;Cech nb simplices;" + std::clog << "File name; Radius; Rips time; Cech time; Ratio Rips/Cech time; Rips nb simplices; Cech nb simplices;" << std::endl; boost::filesystem::directory_iterator end_itr; // default construction yields past-the-end for (boost::filesystem::directory_iterator itr(boost::filesystem::current_path()); itr != end_itr; ++itr) { @@ -93,14 +100,16 @@ int main(int argc, char* argv[]) { if (itr->path().extension() == ".off") // see below { Points_off_reader off_reader(itr->path().string()); + Points_off_reader_cgal off_reader_cgal(itr->path().string()); + Point p0 = off_reader.get_point_cloud()[0]; for (Filtration_value radius = 0.1; radius < 0.4; radius += 0.1) { std::clog << itr->path().stem() << ";"; std::clog << radius << ";"; Gudhi::Clock rips_clock("Rips computation"); - Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), radius, - Gudhi::Minimal_enclosing_ball_radius()); + Rips_complex rips_complex_from_points(off_reader_cgal.get_point_cloud(), radius, + Gudhi::Minimal_enclosing_ball_radius<Kernel>()); Simplex_tree rips_stree; rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); // ------------------------------------------ @@ -110,7 +119,7 @@ int main(int argc, char* argv[]) { std::clog << rips_sec << ";"; Gudhi::Clock cech_clock("Cech computation"); - Cech_complex cech_complex_from_points(off_reader.get_point_cloud(), radius); + Cech_complex cech_complex_from_points(off_reader_cgal.get_point_cloud(), radius); Simplex_tree cech_stree; cech_complex_from_points.create_complex(cech_stree, p0.size() - 1); // ------------------------------------------ diff --git a/src/Cech_complex/concept/SimplicialComplexForCech.h b/src/Cech_complex/concept/SimplicialComplexForCech.h index 00c7df3a..6202fe92 100644 --- a/src/Cech_complex/concept/SimplicialComplexForCech.h +++ b/src/Cech_complex/concept/SimplicialComplexForCech.h @@ -47,8 +47,8 @@ struct SimplicialComplexForCech { }; -} // namespace alpha_complex +} // namespace cech_complex } // namespace Gudhi -#endif // CONCEPT_ALPHA_COMPLEX_SIMPLICIAL_COMPLEX_FOR_ALPHA_H_ +#endif // CONCEPT_CECH_COMPLEX_SIMPLICIAL_COMPLEX_FOR_CECH_H_ diff --git a/src/Cech_complex/example/CMakeLists.txt b/src/Cech_complex/example/CMakeLists.txt index 1b08c7cb..4d11ace2 100644 --- a/src/Cech_complex/example/CMakeLists.txt +++ b/src/Cech_complex/example/CMakeLists.txt @@ -1,17 +1,19 @@ project(Cech_complex_examples) -if (TARGET Boost::program_options) - add_executable ( Cech_complex_example_step_by_step cech_complex_step_by_step.cpp ) - target_link_libraries(Cech_complex_example_step_by_step Boost::program_options) - if (TBB_FOUND) - target_link_libraries(Cech_complex_example_step_by_step ${TBB_LIBRARIES}) +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.0.1) + if (TARGET Boost::program_options) + add_executable ( Cech_complex_example_step_by_step cech_complex_step_by_step.cpp ) + target_link_libraries(Cech_complex_example_step_by_step Boost::program_options) + if (TBB_FOUND) + target_link_libraries(Cech_complex_example_step_by_step ${TBB_LIBRARIES}) + endif() + add_test(NAME Cech_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:Cech_complex_example_step_by_step> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-d" "3") endif() - add_test(NAME Cech_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:Cech_complex_example_step_by_step> - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-d" "3") -endif() -add_executable ( Cech_complex_example_from_points cech_complex_example_from_points.cpp) -if (TBB_FOUND) - target_link_libraries(Cech_complex_example_from_points ${TBB_LIBRARIES}) + add_executable ( Cech_complex_example_from_points cech_complex_example_from_points.cpp) + if (TBB_FOUND) + target_link_libraries(Cech_complex_example_from_points ${TBB_LIBRARIES}) + endif() + add_test(NAME Cech_complex_example_from_points COMMAND $<TARGET_FILE:Cech_complex_example_from_points>) endif() -add_test(NAME Cech_complex_example_from_points COMMAND $<TARGET_FILE:Cech_complex_example_from_points>) diff --git a/src/Cech_complex/example/cech_complex_example_from_points.cpp b/src/Cech_complex/example/cech_complex_example_from_points.cpp index 1a1f708c..78861951 100644 --- a/src/Cech_complex/example/cech_complex_example_from_points.cpp +++ b/src/Cech_complex/example/cech_complex_example_from_points.cpp @@ -1,6 +1,8 @@ #include <gudhi/Cech_complex.h> #include <gudhi/Simplex_tree.h> +#include <CGAL/Epeck_d.h> // For EXACT or SAFE version + #include <iostream> #include <string> #include <vector> @@ -8,32 +10,47 @@ int main() { // Type definitions - using Point_cloud = std::vector<std::array<double, 2>>; using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; using Filtration_value = Simplex_tree::Filtration_value; - using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>; + using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>; + using FT = typename Kernel::FT; + using Point = typename Kernel::Point_d; + using Point_cloud = std::vector<Point>; + using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel, Simplex_tree>; Point_cloud points; - points.push_back({1., 0.}); // 0 - points.push_back({0., 1.}); // 1 - points.push_back({2., 1.}); // 2 - points.push_back({3., 2.}); // 3 - points.push_back({0., 3.}); // 4 - points.push_back({3. + std::sqrt(3.), 3.}); // 5 - points.push_back({1., 4.}); // 6 - points.push_back({3., 4.}); // 7 - points.push_back({2., 4. + std::sqrt(3.)}); // 8 - points.push_back({0., 4.}); // 9 - points.push_back({-0.5, 2.}); // 10 + + std::vector<FT> point0({1., 0.}); + points.emplace_back(point0.begin(), point0.end()); + std::vector<FT> point1({0., 1.}); + points.emplace_back(point1.begin(), point1.end()); + std::vector<FT> point2({2., 1.}); + points.emplace_back(point2.begin(), point2.end()); + std::vector<FT> point3({3., 2.}); + points.emplace_back(point3.begin(), point3.end()); + std::vector<FT> point4({0., 3.}); + points.emplace_back(point4.begin(), point4.end()); + std::vector<FT> point5({3. + std::sqrt(3.), 3.}); + points.emplace_back(point5.begin(), point5.end()); + std::vector<FT> point6({1., 4.}); + points.emplace_back(point6.begin(), point6.end()); + std::vector<FT> point7({3., 4.}); + points.emplace_back(point7.begin(), point7.end()); + std::vector<FT> point8({2., 4. + std::sqrt(3.)}); + points.emplace_back(point8.begin(), point8.end()); + std::vector<FT> point9({0., 4.}); + points.emplace_back(point9.begin(), point9.end()); + std::vector<FT> point10({-0.5, 2.}); + points.emplace_back(point10.begin(), point10.end()); // ---------------------------------------------------------------------------- // Init of a Cech complex from points // ---------------------------------------------------------------------------- - Filtration_value max_radius = 1.; + Filtration_value max_radius = 100.; //100.; Cech_complex cech_complex_from_points(points, max_radius); Simplex_tree stree; - cech_complex_from_points.create_complex(stree, 2); + cech_complex_from_points.create_complex(stree, 6); //6 // ---------------------------------------------------------------------------- // Display information about the one skeleton Cech complex // ---------------------------------------------------------------------------- diff --git a/src/Cech_complex/example/cech_complex_step_by_step.cpp b/src/Cech_complex/example/cech_complex_step_by_step.cpp index f59f0293..c8dd1585 100644 --- a/src/Cech_complex/example/cech_complex_step_by_step.cpp +++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp @@ -9,11 +9,11 @@ */ #include <gudhi/graph_simplicial_complex.h> -#include <gudhi/distance_functions.h> +#include <gudhi/Cech_complex/Cech_kernel.h> #include <gudhi/Simplex_tree.h> #include <gudhi/Points_off_io.h> -#include <gudhi/Miniball.hpp> +#include <CGAL/Epeck_d.h> #include <boost/program_options.hpp> @@ -24,9 +24,9 @@ #include <map> // ---------------------------------------------------------------------------- -// rips_persistence_step_by_step is an example of each step that is required to -// build a Rips over a Simplex_tree. Please refer to rips_persistence to see -// how to do the same thing with the Rips_complex wrapper for less detailed +// cech_complex_step_by_step is an example of each step that is required to +// build a Cech over a Simplex_tree. Please refer to cech_complex_example_from_points to see +// how to do the same thing with the Cech complex wrapper for less detailed // steps. // ---------------------------------------------------------------------------- @@ -34,16 +34,14 @@ using Simplex_tree = Gudhi::Simplex_tree<>; using Simplex_handle = Simplex_tree::Simplex_handle; using Filtration_value = Simplex_tree::Filtration_value; -using Point = std::vector<double>; +using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>; +using Point = typename Kernel::Point_d; using Points_off_reader = Gudhi::Points_off_reader<Point>; using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>; class Cech_blocker { private: using Point_cloud = std::vector<Point>; - using Point_iterator = Point_cloud::const_iterator; - using Coordinate_iterator = Point::const_iterator; - using Min_sphere = Gudhi::Miniball::Miniball<Gudhi::Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>; public: bool operator()(Simplex_handle sh) { @@ -54,7 +52,7 @@ class Cech_blocker { std::clog << "#(" << vertex << ")#"; #endif // DEBUG_TRACES } - Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); + Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius<Kernel>()(points); #ifdef DEBUG_TRACES std::clog << "radius = " << radius << " - " << (radius > max_radius_) << std::endl; #endif // DEBUG_TRACES @@ -63,14 +61,12 @@ class Cech_blocker { } Cech_blocker(Simplex_tree& simplex_tree, Filtration_value max_radius, const std::vector<Point>& point_cloud) : simplex_tree_(simplex_tree), max_radius_(max_radius), point_cloud_(point_cloud) { - dimension_ = point_cloud_[0].size(); } private: Simplex_tree simplex_tree_; Filtration_value max_radius_; std::vector<Point> point_cloud_; - int dimension_; }; void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& max_radius, int& dim_max); @@ -87,9 +83,9 @@ int main(int argc, char* argv[]) { // Compute the proximity graph of the points Proximity_graph prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(), max_radius, - Gudhi::Minimal_enclosing_ball_radius()); + Gudhi::Minimal_enclosing_ball_radius<Kernel>()); - // Construct the Rips complex in a Simplex Tree + // Construct the Cech complex in a Simplex Tree Simplex_tree st; // insert the proximity graph in the simplex tree st.insert_graph(prox_graph); @@ -129,9 +125,9 @@ void program_options(int argc, char* argv[], std::string& off_file_points, Filtr visible.add_options()("help,h", "produce help message")( "max-radius,r", po::value<Filtration_value>(&max_radius)->default_value(std::numeric_limits<Filtration_value>::infinity()), - "Maximal length of an edge for the Rips complex construction.")( + "Maximal length of an edge for the Cech complex construction.")( "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), - "Maximal dimension of the Rips complex we want to compute."); + "Maximal dimension of the Cech complex we want to compute."); po::positional_options_description pos; pos.add("input-file", 1); diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index b0871e10..0031d861 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -11,14 +11,13 @@ #ifndef CECH_COMPLEX_H_ #define CECH_COMPLEX_H_ -#include <gudhi/distance_functions.h> // for Gudhi::Minimal_enclosing_ball_radius +#include <gudhi/Cech_complex/Cech_kernel.h> // for Gudhi::Minimal_enclosing_ball_radius #include <gudhi/graph_simplicial_complex.h> // for Gudhi::Proximity_graph #include <gudhi/Debug_utils.h> // for GUDHI_CHECK #include <gudhi/Cech_complex_blocker.h> // for Gudhi::cech_complex::Cech_blocker #include <iostream> #include <stdexcept> // for exception management -#include <vector> namespace Gudhi { @@ -40,7 +39,7 @@ namespace cech_complex { * \tparam ForwardPointRange must be a range for which `std::begin()` and `std::end()` methods return input * iterators on a point. `std::begin()` and `std::end()` methods are also required for a point. */ -template <typename SimplicialComplexForProximityGraph, typename ForwardPointRange> +template <typename SimplicialComplexForProximityGraph, typename ForwardPointRange, typename Kernel, typename SimplicialComplexForCechComplex> class Cech_complex { private: // Required by compute_proximity_graph @@ -48,16 +47,17 @@ class Cech_complex { using Filtration_value = typename SimplicialComplexForProximityGraph::Filtration_value; using Proximity_graph = Gudhi::Proximity_graph<SimplicialComplexForProximityGraph>; - // Retrieve Coordinate type from ForwardPointRange - using Point_from_range_iterator = typename boost::range_const_iterator<ForwardPointRange>::type; - using Point_from_range = typename std::iterator_traits<Point_from_range_iterator>::value_type; - using Coordinate_iterator = typename boost::range_const_iterator<Point_from_range>::type; - using Coordinate = typename std::iterator_traits<Coordinate_iterator>::value_type; - public: - // Point and Point_cloud type definition - using Point = std::vector<Coordinate>; - using Point_cloud = std::vector<Point>; + + using cech_blocker = Cech_blocker<SimplicialComplexForCechComplex, Cech_complex, Kernel>; + + using Point_d = typename cech_blocker::Point_d; + using Point_cloud = std::vector<Point_d>; + + // Numeric type of coordinates in the kernel + using FT = typename cech_blocker::FT; + // Sphere is a pair of point and squared radius. + using Sphere = typename std::pair<Point_d, FT>; public: /** \brief Cech_complex constructor from a list of points. @@ -70,11 +70,14 @@ class Cech_complex { */ Cech_complex(const ForwardPointRange& points, Filtration_value max_radius) : max_radius_(max_radius) { // Point cloud deep copy - point_cloud_.reserve(boost::size(points)); - for (auto&& point : points) point_cloud_.emplace_back(std::begin(point), std::end(point)); + +// point_cloud_.reserve(boost::size(points)); +// for (auto&& point : points) point_cloud_.emplace_back(point.cartesian_begin(), point.cartesian_end()); + + point_cloud_.assign(points.begin(), points.end()); cech_skeleton_graph_ = Gudhi::compute_proximity_graph<SimplicialComplexForProximityGraph>( - point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); + point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius<Kernel>()); } /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal @@ -85,7 +88,6 @@ class Cech_complex { * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0. * */ - template <typename SimplicialComplexForCechComplex> void create_complex(SimplicialComplexForCechComplex& complex, int dim_max) { GUDHI_CHECK(complex.num_vertices() == 0, std::invalid_argument("Cech_complex::create_complex - simplicial complex is not empty")); @@ -93,8 +95,7 @@ class Cech_complex { // insert the proximity graph in the simplicial complex complex.insert_graph(cech_skeleton_graph_); // expand the graph until dimension dim_max - complex.expansion_with_blockers(dim_max, - Cech_blocker<SimplicialComplexForCechComplex, Cech_complex>(&complex, this)); + complex.expansion_with_blockers(dim_max, cech_blocker(&complex, this)); } /** @return max_radius value given at construction. */ @@ -103,12 +104,18 @@ class Cech_complex { /** @param[in] vertex Point position in the range. * @return The point. */ - const Point& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; } + const Point_d& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; } + + /** + * @return Vector of cached spheres. + */ + std::vector<Sphere> & get_cache() { return cache_; } private: Proximity_graph cech_skeleton_graph_; Filtration_value max_radius_; Point_cloud point_cloud_; + std::vector<Sphere> cache_; }; } // namespace cech_complex diff --git a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h new file mode 100644 index 00000000..89012206 --- /dev/null +++ b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h @@ -0,0 +1,107 @@ +/* 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): Hind Montassif + * + * Copyright (C) 2021 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CECH_KERNEL_H_ +#define CECH_KERNEL_H_ + +#include <CGAL/Epeck_d.h> // for #include <CGAL/NewKernel_d/KernelD_converter.h> + +#include <cmath> // for std::sqrt +#include <vector> + +namespace Gudhi { + +// namespace cech_complex { + +/** @brief Compute the radius of the minimal enclosing ball between Points given by a range of coordinates. + * The points are assumed to have the same dimension. */ +template<typename Kernel> +class Minimal_enclosing_ball_radius { + private: + Kernel kernel_; + public: + using Point = typename Kernel::Point_d; + using Point_cloud = typename std::vector<Point>; + + /** \brief Enclosing ball radius from two points using CGAL. + * + * @param[in] point_1 + * @param[in] point_2 + * @return Enclosing ball radius for the two points. + * \tparam Point must be a Kernel::Point_d from CGAL. + * + */ + double operator()(const Point& point_1, const Point& point_2) const { + return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; + } + + + /** \brief Enclosing ball radius from a point cloud using CGAL. + * + * @param[in] point_cloud The points. + * @return Enclosing ball radius for the points. + * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. + * + */ + double operator()(const Point_cloud& point_cloud) const { + return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); + } + +}; + +/** + * \class Cech_kernel + * \brief Cech complex kernel container. + * + * \details + * The Cech complex kernel container stores CGAL Kernel and dispatch basic computations. + */ + +// template < typename Kernel > +// class Cech_kernel<Kernel> { +// private: +// // Kernel for functions access. +// Kernel kernel_; +// public: +// 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; +// } +// }; + + +//} // namespace cech_complex + +// namespace cechcomplex = cech_complex; + +} // namespace Gudhi + +#endif // CECH_KERNEL_H_ diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 31b9aab5..f7f86534 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -11,10 +11,11 @@ #ifndef CECH_COMPLEX_BLOCKER_H_ #define CECH_COMPLEX_BLOCKER_H_ -#include <gudhi/distance_functions.h> // for Gudhi::Minimal_enclosing_ball_radius +#include <CGAL/NT_converter.h> // for casting from FT to double #include <iostream> #include <vector> +#include <set> #include <cmath> // for std::sqrt namespace Gudhi { @@ -35,28 +36,104 @@ namespace cech_complex { * * \tparam Chech_complex is required by the blocker. */ -template <typename SimplicialComplexForCech, typename Cech_complex> +template <typename SimplicialComplexForCech, typename Cech_complex, typename Kernel> class Cech_blocker { private: - using Point_cloud = typename Cech_complex::Point_cloud; using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; using Filtration_value = typename SimplicialComplexForCech::Filtration_value; public: + + 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>; + + template<class PointIterator> + FT get_squared_radius(PointIterator begin, PointIterator end) const { + return kernel_.compute_squared_radius_d_object()(begin, end); + } + + 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)); + } + + /** \internal \brief Čech complex blocker operator() - the oracle - assigns the filtration value from the simplex * radius and returns if the simplex expansion must be blocked. * \param[in] sh The Simplex_handle. * \return true if the simplex radius is greater than the Cech_complex max_radius*/ bool operator()(Simplex_handle sh) { - Point_cloud points; - for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { - points.push_back(cc_ptr_->get_point(vertex)); -#ifdef DEBUG_TRACES - std::clog << "#(" << vertex << ")#"; -#endif // DEBUG_TRACES + using Point_cloud = std::vector<Point_d>; + CGAL::NT_converter<FT, double> cast_to_double; + Filtration_value radius = 0.; + + // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices + Sphere min_enclos_ball; + CGAL::NT_converter<double, FT> cast_to_FT; + min_enclos_ball.second = cast_to_FT(std::numeric_limits<double>::max()); + Point_cloud face_points; + for (auto face : sc_ptr_->boundary_simplex_range(sh)) { + // Find which vertex of sh is missing in face. We rely on the fact that simplex_vertex_range is sorted. + auto longlist = sc_ptr_->simplex_vertex_range(sh); + auto shortlist = sc_ptr_->simplex_vertex_range(face); + + auto longiter = std::begin(longlist); + auto shortiter = std::begin(shortlist); + auto enditer = std::end(shortlist); + while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; } + auto extra = *longiter; // Vertex_handle + + for (auto vertex : sc_ptr_->simplex_vertex_range(face)) { + face_points.push_back(cc_ptr_->get_point(vertex)); + #ifdef DEBUG_TRACES + std::clog << "#(" << vertex << ")#"; + #endif // DEBUG_TRACES + } + Sphere sph; + auto k = sc_ptr_->key(face); + if(k != sc_ptr_->null_key()) { + sph = cc_ptr_->get_cache().at(k); + } + else { + sph = get_sphere(face_points.cbegin(), face_points.cend()); + } + face_points.clear(); + + if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) <= sph.second) { + radius = std::sqrt(cast_to_double(sph.second)); + #ifdef DEBUG_TRACES + std::clog << "circumcenter: " << sph.first << ", radius: " << radius << std::endl; + #endif // DEBUG_TRACES + if (cast_to_double(sph.second) < cast_to_double(min_enclos_ball.second)) + min_enclos_ball = sph; + } + } + // Get the minimal radius of all faces enclosing balls if exists + if(cast_to_double(min_enclos_ball.second) != std::numeric_limits<double>::max()) { + radius = std::sqrt(cast_to_double(min_enclos_ball.second)); + + sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); + cc_ptr_->get_cache().push_back(min_enclos_ball); } - Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); + + if (radius == 0.) { // Spheres of each face don't contain the whole simplex + Point_cloud points; + for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { + points.push_back(cc_ptr_->get_point(vertex)); + } + Sphere sph = get_sphere(points.cbegin(), points.cend()); + radius = std::sqrt(cast_to_double(sph.second)); + + sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); + cc_ptr_->get_cache().push_back(sph); + } + #ifdef DEBUG_TRACES if (radius > cc_ptr_->max_radius()) std::clog << "radius > max_radius => expansion is blocked\n"; #endif // DEBUG_TRACES @@ -70,6 +147,7 @@ class Cech_blocker { private: SimplicialComplexForCech* sc_ptr_; Cech_complex* cc_ptr_; + Kernel kernel_; }; } // namespace cech_complex diff --git a/src/Cech_complex/test/CMakeLists.txt b/src/Cech_complex/test/CMakeLists.txt index e6a2a18f..2d736f27 100644 --- a/src/Cech_complex/test/CMakeLists.txt +++ b/src/Cech_complex/test/CMakeLists.txt @@ -1,11 +1,14 @@ -include(GUDHI_boost_test) +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.0.1) + include(GUDHI_boost_test) -add_executable ( Cech_complex_test_unit test_cech_complex.cpp ) -if (TBB_FOUND) - target_link_libraries(Cech_complex_test_unit ${TBB_LIBRARIES}) -endif() + add_executable ( Cech_complex_test_unit test_cech_complex.cpp ) + if (TBB_FOUND) + target_link_libraries(Cech_complex_test_unit ${TBB_LIBRARIES}) + endif() + + # Do not forget to copy test files in current binary dir + file(COPY "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) -# Do not forget to copy test files in current binary dir -file(COPY "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + gudhi_add_boost_test(Cech_complex_test_unit) -gudhi_add_boost_test(Cech_complex_test_unit) +endif() diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 6e00d7b5..ca7a9778 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -22,21 +22,21 @@ // to construct Cech_complex from a OFF file of points #include <gudhi/Points_off_io.h> #include <gudhi/Simplex_tree.h> -#include <gudhi/distance_functions.h> +#include <gudhi/Cech_complex/Cech_kernel.h> #include <gudhi/Unitary_tests_utils.h> -#include <gudhi/Miniball.hpp> + +#include <CGAL/Epeck_d.h> // For EXACT or SAFE version // Type definitions using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; -using Point = std::vector<Filtration_value>; +using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>; +using FT = typename Kernel::FT; +using Point = typename Kernel::Point_d; + using Point_cloud = std::vector<Point>; using Points_off_reader = Gudhi::Points_off_reader<Point>; -using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>; - -using Point_iterator = Point_cloud::const_iterator; -using Coordinate_iterator = Point::const_iterator; -using Min_sphere = Gudhi::Miniball::Miniball<Gudhi::Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>; +using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel, Simplex_tree>; BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { // ---------------------------------------------------------------------------- @@ -45,17 +45,29 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { // // ---------------------------------------------------------------------------- Point_cloud points; - points.push_back({1., 0.}); // 0 - points.push_back({0., 1.}); // 1 - points.push_back({2., 1.}); // 2 - points.push_back({3., 2.}); // 3 - points.push_back({0., 3.}); // 4 - points.push_back({3. + std::sqrt(3.), 3.}); // 5 - points.push_back({1., 4.}); // 6 - points.push_back({3., 4.}); // 7 - points.push_back({2., 4. + std::sqrt(3.)}); // 8 - points.push_back({0., 4.}); // 9 - points.push_back({-0.5, 2.}); // 10 + + std::vector<FT> point0({1., 0.}); + points.emplace_back(point0.begin(), point0.end()); + std::vector<FT> point1({0., 1.}); + points.emplace_back(point1.begin(), point1.end()); + std::vector<FT> point2({2., 1.}); + points.emplace_back(point2.begin(), point2.end()); + std::vector<FT> point3({3., 2.}); + points.emplace_back(point3.begin(), point3.end()); + std::vector<FT> point4({0., 3.}); + points.emplace_back(point4.begin(), point4.end()); + std::vector<FT> point5({3. + std::sqrt(3.), 3.}); + points.emplace_back(point5.begin(), point5.end()); + std::vector<FT> point6({1., 4.}); + points.emplace_back(point6.begin(), point6.end()); + std::vector<FT> point7({3., 4.}); + points.emplace_back(point7.begin(), point7.end()); + std::vector<FT> point8({2., 4. + std::sqrt(3.)}); + points.emplace_back(point8.begin(), point8.end()); + std::vector<FT> point9({0., 4.}); + points.emplace_back(point9.begin(), point9.end()); + std::vector<FT> point10({-0.5, 2.}); + points.emplace_back(point10.begin(), point10.end()); Filtration_value max_radius = 1.0; std::clog << "========== NUMBER OF POINTS = " << points.size() << " - Cech max_radius = " << max_radius @@ -96,11 +108,11 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { std::clog << vertex << ","; vp.push_back(points.at(vertex)); } - std::clog << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1)) + std::clog << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius<Kernel>()(vp.at(0), vp.at(1)) << " - filtration =" << st.filtration(f_simplex) << std::endl; BOOST_CHECK(vp.size() == 2); GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), - Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1))); + Gudhi::Minimal_enclosing_ball_radius<Kernel>()(vp.at(0), vp.at(1))); } } @@ -125,35 +137,34 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { for (std::size_t vertex = 0; vertex <= 2; vertex++) { points012.push_back(cech_complex_for_doc.get_point(vertex)); } - std::size_t dimension = points[0].end() - points[0].begin(); - Min_sphere ms012(dimension, points012.begin(), points012.end()); + Kernel kern; Simplex_tree::Filtration_value f012 = st2.filtration(st2.find({0, 1, 2})); - std::clog << "f012= " << f012 << " | ms012_radius= " << std::sqrt(ms012.squared_radius()) << std::endl; + std::clog << "f012= " << f012 << std::endl; - GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(ms012.squared_radius())); + CGAL::NT_converter<FT, double> cast_to_double; + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points012.begin(), points012.end())))); Point_cloud points1410; points1410.push_back(cech_complex_for_doc.get_point(1)); points1410.push_back(cech_complex_for_doc.get_point(4)); points1410.push_back(cech_complex_for_doc.get_point(10)); - Min_sphere ms1410(dimension, points1410.begin(), points1410.end()); Simplex_tree::Filtration_value f1410 = st2.filtration(st2.find({1, 4, 10})); - std::clog << "f1410= " << f1410 << " | ms1410_radius= " << std::sqrt(ms1410.squared_radius()) << std::endl; + std::clog << "f1410= " << f1410 << std::endl; - GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, std::sqrt(ms1410.squared_radius())); + // In this case, the computed sphere using CGAL kernel does not match the minimal enclosing ball; the filtration value check is therefore done against a hardcoded value + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, 1.); Point_cloud points469; points469.push_back(cech_complex_for_doc.get_point(4)); points469.push_back(cech_complex_for_doc.get_point(6)); points469.push_back(cech_complex_for_doc.get_point(9)); - Min_sphere ms469(dimension, points469.begin(), points469.end()); Simplex_tree::Filtration_value f469 = st2.filtration(st2.find({4, 6, 9})); - std::clog << "f469= " << f469 << " | ms469_radius= " << std::sqrt(ms469.squared_radius()) << std::endl; + std::clog << "f469= " << f469 << std::endl; - GUDHI_TEST_FLOAT_EQUALITY_CHECK(f469, std::sqrt(ms469.squared_radius())); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f469, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points469.begin(), points469.end())))); BOOST_CHECK((st2.find({6, 7, 8}) == st2.null_simplex())); BOOST_CHECK((st2.find({3, 5, 7}) == st2.null_simplex())); diff --git a/src/Cech_complex/utilities/CMakeLists.txt b/src/Cech_complex/utilities/CMakeLists.txt index b183c8d8..e80a698e 100644 --- a/src/Cech_complex/utilities/CMakeLists.txt +++ b/src/Cech_complex/utilities/CMakeLists.txt @@ -1,15 +1,17 @@ -project(Cech_complex_utilities) +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.0.1) + project(Cech_complex_utilities) -if (TARGET Boost::program_options) - add_executable(cech_persistence cech_persistence.cpp) - target_link_libraries(cech_persistence Boost::program_options) + if (TARGET Boost::program_options) + add_executable(cech_persistence cech_persistence.cpp) + target_link_libraries(cech_persistence Boost::program_options) - if (TBB_FOUND) - target_link_libraries(cech_persistence ${TBB_LIBRARIES}) - endif() + if (TBB_FOUND) + target_link_libraries(cech_persistence ${TBB_LIBRARIES}) + endif() - add_test(NAME Cech_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:cech_persistence> - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") + add_test(NAME Cech_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:cech_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") - install(TARGETS cech_persistence DESTINATION bin) -endif()
\ No newline at end of file + install(TARGETS cech_persistence DESTINATION bin) + endif() +endif() diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp index daea08e2..ccf63e3e 100644 --- a/src/Cech_complex/utilities/cech_persistence.cpp +++ b/src/Cech_complex/utilities/cech_persistence.cpp @@ -9,13 +9,15 @@ */ #include <gudhi/Cech_complex.h> -#include <gudhi/distance_functions.h> +#include <gudhi/Cech_complex/Cech_kernel.h> #include <gudhi/Simplex_tree.h> #include <gudhi/Persistent_cohomology.h> #include <gudhi/Points_off_io.h> #include <boost/program_options.hpp> +#include <CGAL/Epeck_d.h> // For EXACT or SAFE version + #include <string> #include <vector> #include <limits> // infinity @@ -23,10 +25,12 @@ // Types definition using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; using Filtration_value = Simplex_tree::Filtration_value; -using Point = std::vector<double>; + +using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>; +using Point = typename Kernel::Point_d; using Point_cloud = std::vector<Point>; using Points_off_reader = Gudhi::Points_off_reader<Point>; -using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>; +using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud, Kernel, Simplex_tree>; using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index 9bbc62b7..5e5a1e31 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -13,8 +13,6 @@ #include <gudhi/Debug_utils.h> -#include <gudhi/Miniball.hpp> - #include <boost/range/metafunctions.hpp> #include <boost/range/size.hpp> @@ -59,53 +57,6 @@ class Euclidean_distance { } }; -/** @brief Compute the radius of the minimal enclosing ball between Points given by a range of coordinates. - * The points are assumed to have the same dimension. */ -class Minimal_enclosing_ball_radius { - public: - /** \brief Minimal_enclosing_ball_radius from two points. - * - * @param[in] point_1 First point. - * @param[in] point_2 second point. - * @return The minimal enclosing ball radius for the two points (aka. Euclidean distance / 2.). - * - * \tparam Point must be a range of Cartesian coordinates. - * - */ - template< typename Point > - typename std::iterator_traits<typename boost::range_iterator<Point>::type>::value_type - operator()(const Point& point_1, const Point& point_2) const { - return Euclidean_distance()(point_1, point_2) / 2.; - } - /** \brief Minimal_enclosing_ball_radius from a point cloud. - * - * @param[in] point_cloud The points. - * @return The minimal enclosing ball radius for the points. - * - * \tparam Point_cloud must be a range of points with Cartesian coordinates. - * Point_cloud is a range over a range of Coordinate. - * - */ - template< typename Point_cloud, - typename Point_iterator = typename boost::range_const_iterator<Point_cloud>::type, - typename Point = typename std::iterator_traits<Point_iterator>::value_type, - typename Coordinate_iterator = typename boost::range_const_iterator<Point>::type, - typename Coordinate = typename std::iterator_traits<Coordinate_iterator>::value_type> - Coordinate - operator()(const Point_cloud& point_cloud) const { - using Min_sphere = Miniball::Miniball<Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>; - - Min_sphere ms(boost::size(*point_cloud.begin()), point_cloud.begin(), point_cloud.end()); -#ifdef DEBUG_TRACES - std::clog << "Minimal_enclosing_ball_radius = " << std::sqrt(ms.squared_radius()) << " | nb points = " - << boost::size(point_cloud) << " | dimension = " - << boost::size(*point_cloud.begin()) << std::endl; -#endif // DEBUG_TRACES - - return std::sqrt(ms.squared_radius()); - } -}; - } // namespace Gudhi #endif // DISTANCE_FUNCTIONS_H_ |