diff options
| -rwxr-xr-x | scripts/metagen.sh | 15 | ||||
| -rw-r--r-- | src/Cech_complex/benchmark/CMakeLists.txt | 12 | ||||
| -rw-r--r-- | src/Cech_complex/benchmark/cech_complex_benchmark.cpp | 153 | ||||
| -rw-r--r-- | src/Cech_complex/doc/cech_one_skeleton.png | bin | 47651 -> 12070 bytes | |||
| -rw-r--r-- | src/Cech_complex/example/cech_complex_example_from_points.cpp | 43 | ||||
| -rw-r--r-- | src/Cech_complex/example/cech_complex_step_by_step.cpp | 32 | ||||
| -rw-r--r-- | src/Cech_complex/include/gudhi/Cech_complex.h | 33 | ||||
| -rw-r--r-- | src/Cech_complex/include/gudhi/Cech_complex_blocker.h | 19 | ||||
| -rw-r--r-- | src/Cech_complex/test/test_cech_complex.cpp | 38 | ||||
| -rw-r--r-- | src/Cech_complex/utilities/CMakeLists.txt | 2 | ||||
| -rw-r--r-- | src/Cech_complex/utilities/cech_persistence.cpp | 14 | ||||
| -rw-r--r-- | src/common/include/gudhi/distance_functions.h | 33 |
12 files changed, 284 insertions, 110 deletions
diff --git a/scripts/metagen.sh b/scripts/metagen.sh new file mode 100755 index 00000000..4483d24e --- /dev/null +++ b/scripts/metagen.sh @@ -0,0 +1,15 @@ +#!/bin/bash +sep="_" +for geom in "sphere" "klein" "torus" +do + for number in 10 100 1000 + do + for dim in {3..5} + do + echo "./off_file_from_shape_generator on $geom $geom$sep$number$sep$dim.off $number $dim" + ./off_file_from_shape_generator on $geom $geom$sep$number$sep$dim.off $number $dim + done + done +done + +#./off_file_from_shape_generator in|on sphere|cube off_file_name points_number[integer > 0] dimension[integer > 1] radius[double > 0.0 | default = 1.0] diff --git a/src/Cech_complex/benchmark/CMakeLists.txt b/src/Cech_complex/benchmark/CMakeLists.txt new file mode 100644 index 00000000..2a65865b --- /dev/null +++ b/src/Cech_complex/benchmark/CMakeLists.txt @@ -0,0 +1,12 @@ +cmake_minimum_required(VERSION 2.6) +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}/) + +add_executable(cech_complex_benchmark cech_complex_benchmark.cpp) +target_link_libraries(cech_complex_benchmark ${Boost_FILESYSTEM_LIBRARY}) + +if (TBB_FOUND) + target_link_libraries(cech_complex_benchmark ${TBB_LIBRARIES}) +endif() diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp new file mode 100644 index 00000000..71c88982 --- /dev/null +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -0,0 +1,153 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2018 Inria + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Points_off_io.h> +#include <gudhi/distance_functions.h> +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Clock.h> +#include <gudhi/Rips_complex.h> +#include <gudhi/Cech_complex.h> +#include <gudhi/Simplex_tree.h> + +#include <Miniball/Miniball.hpp> + +#include "boost/filesystem.hpp" // includes all needed Boost.Filesystem declarations + +#include <string> +#include <vector> + + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Point = std::vector<Filtration_value>; +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>; + + +class Radius_distance { + public: + // boost::range_value is not SFINAE-friendly so we cannot use it in the return type + template< typename Point > + typename std::iterator_traits<typename boost::range_iterator<Point>::type>::value_type + operator()(const Point& p1, const Point& p2) const { + // Type def + using Point_cloud = std::vector<Point>; + using Point_iterator = typename Point_cloud::const_iterator; + using Coordinate_iterator = typename Point::const_iterator; + using Min_sphere = typename Miniball::Miniball<Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>; + + Point_cloud point_cloud; + point_cloud.push_back(p1); + point_cloud.push_back(p2); + + GUDHI_CHECK((p1.end()-p1.begin()) != (p2.end()-p2.begin()), "inconsistent point dimensions"); + Min_sphere min_sphere(p1.end()-p1.begin(), point_cloud.begin(),point_cloud.end()); + + return std::sqrt(min_sphere.squared_radius()); + } +}; + + +int main(int argc, char * argv[]) { + std::string off_file_points = "tore3D_1307.off"; + Filtration_value threshold = 1e20; + + // Extract the points from the file filepoints + Points_off_reader off_reader(off_file_points); + + Gudhi::Clock euclidean_clock("Gudhi::Euclidean_distance"); + // Compute the proximity graph of the points + Proximity_graph euclidean_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(), + threshold, + Gudhi::Euclidean_distance()); + + std::cout << euclidean_clock << std::endl; + + Gudhi::Clock radius_clock("Radius_distance"); + // Compute the proximity graph of the points + Proximity_graph radius_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(), + threshold, + Radius_distance()); + std::cout << radius_clock << std::endl; + + Gudhi::Clock squared_radius_clock("Gudhi::Radius_distance()"); + // Compute the proximity graph of the points + Proximity_graph sq_radius_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(), + threshold, + Gudhi::Radius_distance()); + std::cout << squared_radius_clock << std::endl; + + + boost::filesystem::path full_path(boost::filesystem::current_path()); + std::cout << "Current path is : " << full_path << std::endl; + + + std::cout << "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 ) + { + if ( ! boost::filesystem::is_directory(itr->status()) ) + { + if ( itr->path().extension() == ".off" ) // see below + { + Points_off_reader off_reader(itr->path().string()); + + for (Filtration_value radius = 0.1; radius < 0.4; radius += 0.1) { + std::cout << itr->path().stem() << ";"; + std::cout << radius << ";"; + Gudhi::Clock rips_clock("Rips computation"); + Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), radius, Gudhi::Radius_distance()); + Simplex_tree rips_stree; + rips_complex_from_points.create_complex(rips_stree, 2); + // ------------------------------------------ + // Display information about the Rips complex + // ------------------------------------------ + double rips_sec = rips_clock.num_seconds(); + std::cout << rips_sec << ";"; + + Gudhi::Clock cech_clock("Cech computation"); + Cech_complex cech_complex_from_points(off_reader.get_point_cloud(), radius); + Simplex_tree cech_stree; + cech_complex_from_points.create_complex(cech_stree, 2); + // ------------------------------------------ + // Display information about the Cech complex + // ------------------------------------------ + double cech_sec = cech_clock.num_seconds(); + std::cout << cech_sec << ";"; + std::cout << cech_sec / rips_sec << ";"; + + std::cout << rips_stree.num_simplices() << ";"; + std::cout << cech_stree.num_simplices() << ";" << std::endl; + } + } + } + } + + + return 0; +} diff --git a/src/Cech_complex/doc/cech_one_skeleton.png b/src/Cech_complex/doc/cech_one_skeleton.png Binary files differindex 1028770e..ffa9c329 100644 --- a/src/Cech_complex/doc/cech_one_skeleton.png +++ b/src/Cech_complex/doc/cech_one_skeleton.png 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 882849c3..97327e69 100644 --- a/src/Cech_complex/example/cech_complex_example_from_points.cpp +++ b/src/Cech_complex/example/cech_complex_example_from_points.cpp @@ -1,25 +1,3 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Author(s): Vincent Rouvreau - * - * Copyright (C) 2018 Inria - * - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see <http://www.gnu.org/licenses/>. - */ - #include <gudhi/Cech_complex.h> #include <gudhi/Simplex_tree.h> #include <gudhi/distance_functions.h> @@ -37,23 +15,22 @@ int main() { using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>; Point_cloud points; - points.push_back({1.0, 1.0}); - points.push_back({7.0, 0.0}); - points.push_back({4.0, 6.0}); - points.push_back({9.0, 6.0}); - points.push_back({0.0, 14.0}); - points.push_back({2.0, 19.0}); - points.push_back({9.0, 17.0}); + points.push_back({0., 0.}); + points.push_back({0., 2.}); + points.push_back({std::sqrt(3.), 1.}); + points.push_back({1., 0.}); + points.push_back({1., 2.}); + points.push_back({1. - std::sqrt(3.), 1.}); // ---------------------------------------------------------------------------- // Init of a Cech complex from points // ---------------------------------------------------------------------------- - // 7.1 is a magic number to force one blocker, and one non-blocker - Filtration_value threshold = 7.1; - Cech_complex cech_complex_from_points(points, threshold, Gudhi::Euclidean_distance()); + // 5. is a magic number to force one blocker, and one non-blocker + Filtration_value max_radius = 12.; + 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, -1); // ---------------------------------------------------------------------------- // 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 e71086b6..8705a3e5 100644 --- a/src/Cech_complex/example/cech_complex_step_by_step.cpp +++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp @@ -65,23 +65,22 @@ class Cech_blocker { std::cout << "#(" << vertex << ")#"; #endif // DEBUG_TRACES } - Min_sphere ms(dimension_, points.begin(),points.end()); - Filtration_value radius = std::sqrt(ms.squared_radius()); + Filtration_value radius = Gudhi::Radius_distance()(points); #ifdef DEBUG_TRACES - std::cout << "radius = " << radius << " - " << (radius > threshold_) << std::endl; + std::cout << "radius = " << radius << " - " << (radius > max_radius_) << std::endl; #endif // DEBUG_TRACES simplex_tree_.assign_filtration(sh, radius); - return (radius > threshold_); + return (radius > max_radius_); } - Cech_blocker(Simplex_tree& simplex_tree, Filtration_value threshold, const std::vector<Point>& point_cloud) + Cech_blocker(Simplex_tree& simplex_tree, Filtration_value max_radius, const std::vector<Point>& point_cloud) : simplex_tree_(simplex_tree), - threshold_(threshold), + max_radius_(max_radius), point_cloud_(point_cloud) { dimension_ = point_cloud_[0].size(); } private: Simplex_tree simplex_tree_; - Filtration_value threshold_; + Filtration_value max_radius_; std::vector<Point> point_cloud_; int dimension_; }; @@ -89,31 +88,31 @@ class Cech_blocker { void program_options(int argc, char * argv[] , std::string & off_file_points - , Filtration_value & threshold + , Filtration_value & max_radius , int & dim_max); int main(int argc, char * argv[]) { std::string off_file_points; - Filtration_value threshold; + Filtration_value max_radius; int dim_max; - program_options(argc, argv, off_file_points, threshold, dim_max); + program_options(argc, argv, off_file_points, max_radius, dim_max); // Extract the points from the file filepoints Points_off_reader off_reader(off_file_points); // Compute the proximity graph of the points Proximity_graph prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(), - threshold, - Gudhi::Euclidean_distance()); + max_radius, + Gudhi::Radius_distance()); // Construct the Rips complex in a Simplex Tree Simplex_tree st; // insert the proximity graph in the simplex tree st.insert_graph(prox_graph); // expand the graph until dimension dim_max - st.expansion_with_blockers(dim_max, Cech_blocker(st, threshold, off_reader.get_point_cloud())); + st.expansion_with_blockers(dim_max, Cech_blocker(st, max_radius, off_reader.get_point_cloud())); std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; std::cout << " and has dimension " << st.dimension() << " \n"; @@ -123,7 +122,6 @@ int main(int argc, char * argv[]) { #if DEBUG_TRACES std::cout << "********************************************************************\n"; - // Display the Simplex_tree - Can not be done in the middle of 2 inserts std::cout << "* The complex contains " << st.num_simplices() << " simplices - dimension=" << st.dimension() << "\n"; std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; for (auto f_simplex : st.filtration_simplex_range()) { @@ -140,7 +138,7 @@ int main(int argc, char * argv[]) { void program_options(int argc, char * argv[] , std::string & off_file_points - , Filtration_value & threshold + , Filtration_value & max_radius , int & dim_max) { namespace po = boost::program_options; po::options_description hidden("Hidden options"); @@ -151,8 +149,8 @@ void program_options(int argc, char * argv[] po::options_description visible("Allowed options", 100); visible.add_options() ("help,h", "produce help message") - ("max-edge-length,r", - po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), + ("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.") ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), "Maximal dimension of the Rips complex we want to compute."); diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index e847c97f..a50ed9fa 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -23,12 +23,12 @@ #ifndef CECH_COMPLEX_H_ #define CECH_COMPLEX_H_ +#include <gudhi/distance_functions.h> // for Gudhi::Squared_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 <cstddef> // for std::size_t #include <stdexcept> // for exception management namespace Gudhi { @@ -43,7 +43,7 @@ namespace cech_complex { * * \details * The data structure is a proximity graph, containing edges when the edge length is less or equal - * to a given threshold. Edge length is computed from a user given point cloud with a given distance function. + * to a given max_radius. Edge length is computed from `Gudhi::Squared_radius` distance function. * * \tparam SimplicialComplexForProximityGraph furnishes `Vertex_handle` and `Filtration_value` type definition required * by `Gudhi::Proximity_graph`. @@ -62,23 +62,18 @@ class Cech_complex { /** \brief Cech_complex constructor from a list of points. * * @param[in] points Range of points. - * @param[in] threshold Rips value. - * @param[in] distance distance function that returns a `Filtration_value` from 2 given points. + * @param[in] max_radius Maximal radius value. * * \tparam ForwardPointRange must be a range for which `.size()`, `.begin()` and `.end()` methods return input * iterators on a point. `.begin()` and `.end()` methods are required for a point. * - * \tparam Distance furnishes `operator()(const Point& p1, const Point& p2)`, where - * `Point` is a point from the `ForwardPointRange`, and that returns a `Filtration_value`. */ - template<typename Distance > - Cech_complex(const ForwardPointRange& points, Filtration_value threshold, Distance distance) - : threshold_(threshold), + Cech_complex(const ForwardPointRange& points, Filtration_value max_radius) + : max_radius_(max_radius), point_cloud_(points) { - dimension_ = points.begin()->end() - points.begin()->begin(); cech_skeleton_graph_ = Gudhi::compute_proximity_graph<SimplicialComplexForProximityGraph>(point_cloud_, - threshold_, - distance); + max_radius_, + Gudhi::Radius_distance()); } /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal @@ -101,14 +96,9 @@ class Cech_complex { Cech_blocker<SimplicialComplexForCechComplex, ForwardPointRange>(complex, this)); } - /** @return Threshold value given at construction. */ - Filtration_value threshold() const { - return threshold_; - } - - /** @return Dimension value given at construction by the first point dimension. */ - std::size_t dimension() const { - return dimension_; + /** @return max_radius value given at construction. */ + Filtration_value max_radius() const { + return max_radius_; } /** @param[in] vertex Point position in the range. @@ -123,9 +113,8 @@ class Cech_complex { private: Proximity_graph cech_skeleton_graph_; - Filtration_value threshold_; + Filtration_value max_radius_; ForwardPointRange point_cloud_; - std::size_t dimension_; }; } // namespace cech_complex diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index fb52f712..d718b56e 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -23,9 +23,8 @@ #ifndef CECH_COMPLEX_BLOCKER_H_ #define CECH_COMPLEX_BLOCKER_H_ -#include <gudhi/Cech_complex.h> // Cech_blocker is using a pointer on Gudhi::cech_complex::Cech_complex - -#include <Miniball/Miniball.hpp> +#include <gudhi/Cech_complex.h> // Cech_blocker is using a pointer on Gudhi::cech_complex::Cech_complex +#include <gudhi/distance_functions.h> // for Gudhi::Squared_radius #include <iostream> #include <vector> @@ -58,9 +57,6 @@ class Cech_blocker { private: using Point = std::vector<double>; using Point_cloud = std::vector<Point>; - using Point_iterator = Point_cloud::const_iterator; - using Coordinate_iterator = Point::const_iterator; - using Min_sphere = Miniball::Miniball<Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>; using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; using Filtration_value = typename SimplicialComplexForCech::Filtration_value; using Cech_complex = Gudhi::cech_complex::Cech_complex<SimplicialComplexForCech, ForwardPointRange>; @@ -69,7 +65,7 @@ class Cech_blocker { /** \internal \brief Cech 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 threshold*/ + * \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 : simplicial_complex_.simplex_vertex_range(sh)) { @@ -79,13 +75,12 @@ class Cech_blocker { std::cout << "#(" << vertex << ")#"; #endif // DEBUG_TRACES } - Min_sphere ms(cc_ptr_->dimension(), points.begin(),points.end()); - Filtration_value diameter = 2 * std::sqrt(ms.squared_radius()); + Filtration_value squared_radius = Gudhi::Radius_distance()(points); #ifdef DEBUG_TRACES - std::cout << "diameter = " << diameter << " - " << (diameter > cc_ptr_->threshold()) << std::endl; + std::cout << "squared_radius = " << squared_radius << " - " << (squared_radius > cc_ptr_->max_radius()) << std::endl; #endif // DEBUG_TRACES - simplicial_complex_.assign_filtration(sh, diameter); - return (diameter > cc_ptr_->threshold()); + simplicial_complex_.assign_filtration(sh, squared_radius); + return (squared_radius > cc_ptr_->max_radius()); } /** \internal \brief Cech complex blocker constructor. */ diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index aa42d322..626f1d82 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -58,14 +58,15 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { // // ---------------------------------------------------------------------------- std::string off_file_name("alphacomplexdoc.off"); - double threshold = 12.0; - std::cout << "========== OFF FILE NAME = " << off_file_name << " - Cech threshold=" << - threshold << "==========" << std::endl; + double max_radius = 12.0; + std::cout << "========== OFF FILE NAME = " << off_file_name << " - Cech max_radius=" << + max_radius << "==========" << std::endl; Points_off_reader off_reader(off_file_name); Point_cloud point_cloud = off_reader.get_point_cloud(); - Cech_complex cech_complex_from_file(point_cloud, threshold, Gudhi::Euclidean_distance()); + Cech_complex cech_complex_from_file(point_cloud, max_radius); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(cech_complex_from_file.max_radius(), max_radius); std::size_t i = 0; for (; i < point_cloud.size(); i++) { BOOST_CHECK(point_cloud[i] == *(cech_complex_from_file.point_iterator(i))); @@ -101,10 +102,10 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { std::cout << vertex << ","; vp.push_back(off_reader.get_point_cloud().at(vertex)); } - std::cout << ") - distance =" << Gudhi::Euclidean_distance()(vp.at(0), vp.at(1)) << + std::cout << ") - distance =" << Gudhi::Radius_distance()(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::Euclidean_distance()(vp.at(0), vp.at(1))); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), Gudhi::Radius_distance()(vp.at(0), vp.at(1))); } } @@ -125,7 +126,8 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { points012.push_back(Point(cech_complex_from_file.point_iterator(vertex)->begin(), cech_complex_from_file.point_iterator(vertex)->end())); } - Min_sphere ms012(cech_complex_from_file.dimension(), points012.begin(),points012.end()); + std::size_t dimension = point_cloud[0].end() - point_cloud[0].begin(); + Min_sphere ms012(dimension, points012.begin(),points012.end()); Simplex_tree::Filtration_value f012 = st2.filtration(st2.find({0, 1, 2})); std::cout << "f012= " << f012 << " | ms012_radius= " << std::sqrt(ms012.squared_radius()) << std::endl; @@ -137,7 +139,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { points456.push_back(Point(cech_complex_from_file.point_iterator(vertex)->begin(), cech_complex_from_file.point_iterator(vertex)->end())); } - Min_sphere ms456(cech_complex_from_file.dimension(), points456.begin(),points456.end()); + Min_sphere ms456(dimension, points456.begin(),points456.end()); Simplex_tree::Filtration_value f456 = st2.filtration(st2.find({4, 5, 6})); std::cout << "f456= " << f456 << " | ms456_radius= " << std::sqrt(ms456.squared_radius()) << std::endl; @@ -161,7 +163,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { points0123.push_back(Point(cech_complex_from_file.point_iterator(vertex)->begin(), cech_complex_from_file.point_iterator(vertex)->end())); } - Min_sphere ms0123(cech_complex_from_file.dimension(), points0123.begin(),points0123.end()); + Min_sphere ms0123(dimension, points0123.begin(),points0123.end()); Simplex_tree::Filtration_value f0123 = st3.filtration(st3.find({0, 1, 2, 3})); std::cout << "f0123= " << f0123 << " | ms0123_radius= " << std::sqrt(ms0123.squared_radius()) << std::endl; @@ -175,7 +177,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { points01.push_back(Point(cech_complex_from_file.point_iterator(vertex)->begin(), cech_complex_from_file.point_iterator(vertex)->end())); } - Min_sphere ms01(cech_complex_from_file.dimension(), points01.begin(),points01.end()); + Min_sphere ms01(dimension, points01.begin(),points01.end()); Simplex_tree::Filtration_value f01 = st2.filtration(st2.find({0, 1})); std::cout << "f01= " << f01 << " | ms01_radius= " << std::sqrt(ms01.squared_radius()) << std::endl; @@ -199,7 +201,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_points) { // ---------------------------------------------------------------------------- // Init of a Cech complex from the list of points // ---------------------------------------------------------------------------- - Cech_complex cech_complex_from_points(points, 2.0, Gudhi::Euclidean_distance()); + Cech_complex cech_complex_from_points(points, 2.0); std::cout << "========== cech_complex_from_points ==========" << std::endl; Simplex_tree st; @@ -234,13 +236,13 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_points) { GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.0); break; case 1: - GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 1.41421, .00001); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.5); break; case 2: - GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.816497, .00001); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.666667, .00001); break; case 3: - GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.866025, .00001); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.75); break; default: BOOST_CHECK(false); // Shall not happen @@ -257,12 +259,12 @@ BOOST_AUTO_TEST_CASE(Cech_create_complex_throw) { // // ---------------------------------------------------------------------------- std::string off_file_name("alphacomplexdoc.off"); - double threshold = 12.0; - std::cout << "========== OFF FILE NAME = " << off_file_name << " - Cech threshold=" << - threshold << "==========" << std::endl; + double max_radius = 12.0; + std::cout << "========== OFF FILE NAME = " << off_file_name << " - Cech max_radius=" << + max_radius << "==========" << std::endl; Gudhi::Points_off_reader<Point> off_reader(off_file_name); - Cech_complex cech_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + Cech_complex cech_complex_from_file(off_reader.get_point_cloud(), max_radius); Simplex_tree stree; std::vector<int> simplex = {0, 1, 2}; diff --git a/src/Cech_complex/utilities/CMakeLists.txt b/src/Cech_complex/utilities/CMakeLists.txt index a4f89d2c..30b99729 100644 --- a/src/Cech_complex/utilities/CMakeLists.txt +++ b/src/Cech_complex/utilities/CMakeLists.txt @@ -9,6 +9,6 @@ if (TBB_FOUND) endif() add_test(NAME Cech_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:cech_persistence> - "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") install(TARGETS cech_persistence DESTINATION bin) diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp index e93596d4..93a200ff 100644 --- a/src/Cech_complex/utilities/cech_persistence.cpp +++ b/src/Cech_complex/utilities/cech_persistence.cpp @@ -43,20 +43,20 @@ using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, - Filtration_value& threshold, int& dim_max, int& p, Filtration_value& min_persistence); + Filtration_value& max_radius, int& dim_max, int& p, Filtration_value& min_persistence); int main(int argc, char* argv[]) { std::string off_file_points; std::string filediag; - Filtration_value threshold; + Filtration_value max_radius; int dim_max; int p; Filtration_value min_persistence; - program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence); + program_options(argc, argv, off_file_points, filediag, max_radius, dim_max, p, min_persistence); Points_off_reader off_reader(off_file_points); - Cech_complex cech_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + Cech_complex cech_complex_from_file(off_reader.get_point_cloud(), max_radius); // Construct the Cech complex in a Simplex Tree Simplex_tree simplex_tree; @@ -88,7 +88,7 @@ int main(int argc, char* argv[]) { } void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, - Filtration_value& threshold, int& dim_max, int& p, Filtration_value& min_persistence) { + Filtration_value& max_radius, int& dim_max, int& p, Filtration_value& min_persistence) { namespace po = boost::program_options; po::options_description hidden("Hidden options"); hidden.add_options()("input-file", po::value<std::string>(&off_file_points), @@ -98,8 +98,8 @@ void program_options(int argc, char* argv[], std::string& off_file_points, std:: visible.add_options()("help,h", "produce help message")( "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), "Name of file in which the persistence diagram is written. Default print in std::cout")( - "max-edge-length,r", - po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()), + "max-radius,r", + po::value<Filtration_value>(&max_radius)->default_value(std::numeric_limits<Filtration_value>::infinity()), "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 Cech complex we want to compute.")( diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index 3a5d1fd5..3ce51ad1 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -25,6 +25,8 @@ #include <gudhi/Debug_utils.h> +#include <Miniball/Miniball.hpp> + #include <boost/range/metafunctions.hpp> #include <cmath> // for std::sqrt @@ -68,6 +70,37 @@ class Euclidean_distance { } }; +/** @brief Compute the squared radius between Points given by a range of coordinates. The points are assumed to + * have the same dimension. */ +class Radius_distance { + public: + // boost::range_value is not SFINAE-friendly so we cannot use it in the return type + template< typename Point > + typename std::iterator_traits<typename boost::range_iterator<Point>::type>::value_type + operator()(const Point& p1, const Point& p2) const { + return Euclidean_distance()(p1, p2) / 2.; + } + // boost::range_value is not SFINAE-friendly so we cannot use it in the return type + 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(point_cloud.begin()->end() - point_cloud.begin()->begin(), point_cloud.begin(),point_cloud.end()); + Min_sphere ms(point_cloud.end() - point_cloud.begin(), point_cloud.begin(),point_cloud.end()); +#ifdef DEBUG_TRACES + std::cout << "Radius on " << point_cloud.end() - point_cloud.begin() << " points = " + << std::sqrt(ms.squared_radius()) << std::endl; +#endif // DEBUG_TRACES + + return std::sqrt(ms.squared_radius()); + } +}; + } // namespace Gudhi #endif // DISTANCE_FUNCTIONS_H_ |
