From 0586a149b5bb3a4b65b63b2ab7d3ecdd9682ee1b Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Tue, 20 Feb 2018 16:03:52 +0000 Subject: tests and utils fix git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/cechcomplex_vincent@3253 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 5786a8a7e4b16750f29fac99ca61926158542cfd --- src/Cech_complex/test/test_cech_complex.cpp | 274 ++++++++++++++++++++++++++++ 1 file changed, 274 insertions(+) create mode 100644 src/Cech_complex/test/test_cech_complex.cpp (limited to 'src/Cech_complex/test/test_cech_complex.cpp') diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp new file mode 100644 index 00000000..aa42d322 --- /dev/null +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -0,0 +1,274 @@ +/* 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 . + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "cech_complex" +#include + +#include // float comparison +#include +#include +#include +#include // std::max + +#include +// to construct Cech_complex from a OFF file of points +#include +#include +#include +#include + +#include + +// Type definitions +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Point = std::vector; +using Point_cloud = std::vector; +using Points_off_reader = Gudhi::Points_off_reader; +using Cech_complex = Gudhi::cech_complex::Cech_complex; + +using Point_iterator = Point_cloud::const_iterator; +using Coordinate_iterator = Point::const_iterator; +using Min_sphere = Miniball::Miniball>; + +BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { + // ---------------------------------------------------------------------------- + // + // Init of a Cech complex from a OFF 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; + + 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()); + + std::size_t i = 0; + for (; i < point_cloud.size(); i++) { + BOOST_CHECK(point_cloud[i] == *(cech_complex_from_file.point_iterator(i))); + } +#ifdef GUDHI_DEBUG + BOOST_CHECK_THROW (cech_complex_from_file.point_iterator(i+1), std::out_of_range); +#endif // GUDHI_DEBUG + + const int DIMENSION_1 = 1; + Simplex_tree st; + cech_complex_from_file.create_complex(st, DIMENSION_1); + std::cout << "st.dimension()=" << st.dimension() << std::endl; + BOOST_CHECK(st.dimension() == DIMENSION_1); + + const int NUMBER_OF_VERTICES = 7; + std::cout << "st.num_vertices()=" << st.num_vertices() << std::endl; + BOOST_CHECK(st.num_vertices() == NUMBER_OF_VERTICES); + + std::cout << "st.num_simplices()=" << st.num_simplices() << std::endl; + BOOST_CHECK(st.num_simplices() == 18); + + // Check filtration values of vertices is 0.0 + for (auto f_simplex : st.skeleton_simplex_range(0)) { + BOOST_CHECK(st.filtration(f_simplex) == 0.0); + } + + // Check filtration values of edges + for (auto f_simplex : st.skeleton_simplex_range(DIMENSION_1)) { + if (DIMENSION_1 == st.dimension(f_simplex)) { + std::vector vp; + std::cout << "vertex = ("; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + 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)) << + " - 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))); + } + } + + const int DIMENSION_2 = 2; + Simplex_tree st2; + cech_complex_from_file.create_complex(st2, DIMENSION_2); + std::cout << "st2.dimension()=" << st2.dimension() << std::endl; + BOOST_CHECK(st2.dimension() == DIMENSION_2); + + std::cout << "st2.num_vertices()=" << st2.num_vertices() << std::endl; + BOOST_CHECK(st2.num_vertices() == NUMBER_OF_VERTICES); + + std::cout << "st2.num_simplices()=" << st2.num_simplices() << std::endl; + BOOST_CHECK(st2.num_simplices() == 23); + + Point_cloud points012; + for (std::size_t vertex = 0; vertex <= 2; vertex++) { + 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()); + + 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; + + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(ms012.squared_radius())); + + Point_cloud points456; + for (std::size_t vertex = 4; vertex <= 6; vertex++) { + 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()); + + 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; + + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f456, std::sqrt(ms456.squared_radius())); + + const int DIMENSION_3 = 3; + Simplex_tree st3; + cech_complex_from_file.create_complex(st3, DIMENSION_3); + std::cout << "st3.dimension()=" << st3.dimension() << std::endl; + BOOST_CHECK(st3.dimension() == DIMENSION_3); + + std::cout << "st3.num_vertices()=" << st3.num_vertices() << std::endl; + BOOST_CHECK(st3.num_vertices() == NUMBER_OF_VERTICES); + + std::cout << "st3.num_simplices()=" << st3.num_simplices() << std::endl; + BOOST_CHECK(st3.num_simplices() == 24); + + Point_cloud points0123; + for (std::size_t vertex = 0; vertex <= 3; vertex++) { + 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()); + + 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; + + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f0123, std::sqrt(ms0123.squared_radius())); + + + + Point_cloud points01; + for (std::size_t vertex = 0; vertex <= 1; vertex++) { + 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()); + + Simplex_tree::Filtration_value f01 = st2.filtration(st2.find({0, 1})); + std::cout << "f01= " << f01 << " | ms01_radius= " << std::sqrt(ms01.squared_radius()) << std::endl; + +} + +BOOST_AUTO_TEST_CASE(Cech_complex_from_points) { + // ---------------------------------------------------------------------------- + // Init of a list of points + // ---------------------------------------------------------------------------- + Point_cloud points; + std::vector coords = { 0.0, 0.0, 0.0, 1.0 }; + points.push_back(Point(coords.begin(), coords.end())); + coords = { 0.0, 0.0, 1.0, 0.0 }; + points.push_back(Point(coords.begin(), coords.end())); + coords = { 0.0, 1.0, 0.0, 0.0 }; + points.push_back(Point(coords.begin(), coords.end())); + coords = { 1.0, 0.0, 0.0, 0.0 }; + points.push_back(Point(coords.begin(), coords.end())); + + // ---------------------------------------------------------------------------- + // Init of a Cech complex from the list of points + // ---------------------------------------------------------------------------- + Cech_complex cech_complex_from_points(points, 2.0, Gudhi::Euclidean_distance()); + + std::cout << "========== cech_complex_from_points ==========" << std::endl; + Simplex_tree st; + const int DIMENSION = 3; + cech_complex_from_points.create_complex(st, DIMENSION); + + // Another way to check num_simplices + std::cout << "Iterator on Cech complex simplices in the filtration order, with [filtration value]:" << std::endl; + int num_simplices = 0; + for (auto f_simplex : st.filtration_simplex_range()) { + num_simplices++; + std::cout << " ( "; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " << "[" << st.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + BOOST_CHECK(num_simplices == 15); + std::cout << "st.num_simplices()=" << st.num_simplices() << std::endl; + BOOST_CHECK(st.num_simplices() == 15); + + std::cout << "st.dimension()=" << st.dimension() << std::endl; + BOOST_CHECK(st.dimension() == DIMENSION); + std::cout << "st.num_vertices()=" << st.num_vertices() << std::endl; + BOOST_CHECK(st.num_vertices() == 4); + + for (auto f_simplex : st.filtration_simplex_range()) { + std::cout << "dimension(" << st.dimension(f_simplex) << ") - f = " << st.filtration(f_simplex) << std::endl; + switch (st.dimension(f_simplex)) { + case 0: + 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); + break; + case 2: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.816497, .00001); + break; + case 3: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.866025, .00001); + break; + default: + BOOST_CHECK(false); // Shall not happen + break; + } + } +} + +#ifdef GUDHI_DEBUG +BOOST_AUTO_TEST_CASE(Cech_create_complex_throw) { + // ---------------------------------------------------------------------------- + // + // Init of a Cech complex from a OFF 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; + + Gudhi::Points_off_reader off_reader(off_file_name); + Cech_complex cech_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + Simplex_tree stree; + std::vector simplex = {0, 1, 2}; + stree.insert_simplex_and_subfaces(simplex); + std::cout << "Check exception throw in debug mode" << std::endl; + // throw excpt because stree is not empty + BOOST_CHECK_THROW (cech_complex_from_file.create_complex(stree, 1), std::invalid_argument); +} +#endif -- cgit v1.2.3 From 3cd1e01f0b0d4fdb46f49ec640c389374ca2fe70 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Thu, 22 Feb 2018 23:16:55 +0000 Subject: Fix Cech with radius distance Add a meta generation script for off_file_generator git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/cechcomplex_vincent@3256 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: fb13baa124ddc97c0dc61835ab0c72595d666155 --- scripts/metagen.sh | 15 ++ src/Cech_complex/benchmark/CMakeLists.txt | 12 ++ .../benchmark/cech_complex_benchmark.cpp | 153 +++++++++++++++++++++ src/Cech_complex/doc/cech_one_skeleton.png | Bin 47651 -> 12070 bytes .../example/cech_complex_example_from_points.cpp | 43 ++---- .../example/cech_complex_step_by_step.cpp | 32 ++--- src/Cech_complex/include/gudhi/Cech_complex.h | 33 ++--- .../include/gudhi/Cech_complex_blocker.h | 19 +-- src/Cech_complex/test/test_cech_complex.cpp | 38 ++--- src/Cech_complex/utilities/CMakeLists.txt | 2 +- src/Cech_complex/utilities/cech_persistence.cpp | 14 +- src/common/include/gudhi/distance_functions.h | 33 +++++ 12 files changed, 284 insertions(+), 110 deletions(-) create mode 100755 scripts/metagen.sh create mode 100644 src/Cech_complex/benchmark/CMakeLists.txt create mode 100644 src/Cech_complex/benchmark/cech_complex_benchmark.cpp (limited to 'src/Cech_complex/test/test_cech_complex.cpp') 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 . + */ + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "boost/filesystem.hpp" // includes all needed Boost.Filesystem declarations + +#include +#include + + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Point = std::vector; +using Point_cloud = std::vector; +using Points_off_reader = Gudhi::Points_off_reader; +using Proximity_graph = Gudhi::Proximity_graph; +using Rips_complex = Gudhi::rips_complex::Rips_complex; +using Cech_complex = Gudhi::cech_complex::Cech_complex; + + +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::type>::value_type + operator()(const Point& p1, const Point& p2) const { + // Type def + using Point_cloud = std::vector; + using Point_iterator = typename Point_cloud::const_iterator; + using Coordinate_iterator = typename Point::const_iterator; + using Min_sphere = typename Miniball::Miniball>; + + 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(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(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(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 index 1028770e..ffa9c329 100644 Binary files a/src/Cech_complex/doc/cech_one_skeleton.png and b/src/Cech_complex/doc/cech_one_skeleton.png differ 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 . - */ - #include #include #include @@ -37,23 +15,22 @@ int main() { using Cech_complex = Gudhi::cech_complex::Cech_complex; 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_cloud) + Cech_blocker(Simplex_tree& simplex_tree, Filtration_value max_radius, const std::vector& 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_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(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(&threshold)->default_value(std::numeric_limits::infinity()), + ("max-radius,r", + po::value(&max_radius)->default_value(std::numeric_limits::infinity()), "Maximal length of an edge for the Rips complex construction.") ("cpx-dimension,d", po::value(&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 // for Gudhi::Squared_radius #include // for Gudhi::Proximity_graph #include // for GUDHI_CHECK #include // for Gudhi::cech_complex::Cech_blocker #include -#include // for std::size_t #include // 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 - 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(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(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 // Cech_blocker is using a pointer on Gudhi::cech_complex::Cech_complex - -#include +#include // Cech_blocker is using a pointer on Gudhi::cech_complex::Cech_complex +#include // for Gudhi::Squared_radius #include #include @@ -58,9 +57,6 @@ class Cech_blocker { private: using Point = std::vector; using Point_cloud = std::vector; - using Point_iterator = Point_cloud::const_iterator; - using Coordinate_iterator = Point::const_iterator; - using Min_sphere = Miniball::Miniball>; using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; using Filtration_value = typename SimplicialComplexForCech::Filtration_value; using Cech_complex = Gudhi::cech_complex::Cech_complex; @@ -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 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 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 $ - "${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; 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(&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(&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(&threshold)->default_value(std::numeric_limits::infinity()), + "max-radius,r", + po::value(&max_radius)->default_value(std::numeric_limits::infinity()), "Maximal length of an edge for the Cech complex construction.")( "cpx-dimension,d", po::value(&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 +#include + #include #include // 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::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::type, + typename Point= typename std::iterator_traits::value_type, + typename Coordinate_iterator = typename boost::range_const_iterator::type, + typename Coordinate = typename std::iterator_traits::value_type> + Coordinate + operator()(const Point_cloud& point_cloud) const { + using Min_sphere = Miniball::Miniball>; + + //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_ -- cgit v1.2.3 From b3a64294af818c977804c4b67a317782d872e2b5 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Mon, 5 Mar 2018 13:38:57 +0000 Subject: Fix doc and tests git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/cechcomplex_vincent@3262 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: c76eb981f5960938221aa2498cb87b0707391733 --- .../benchmark/cech_complex_benchmark.cpp | 10 +- src/Cech_complex/doc/Intro_cech_complex.h | 52 ++++--- .../doc/cech_complex_representation.ipe | 160 +++++++++++---------- .../doc/cech_complex_representation.png | Bin 15677 -> 54399 bytes src/Cech_complex/doc/cech_one_skeleton.ipe | 154 +++++++++----------- src/Cech_complex/doc/cech_one_skeleton.png | Bin 12070 -> 29354 bytes .../example/cech_complex_example_from_points.cpp | 22 +-- .../cech_complex_example_from_points_for_doc.txt | 45 ++++-- src/Cech_complex/include/gudhi/Cech_complex.h | 4 +- .../include/gudhi/Cech_complex_blocker.h | 11 +- src/Cech_complex/test/test_cech_complex.cpp | 12 +- src/Doxyfile | 3 +- src/common/include/gudhi/distance_functions.h | 8 +- 13 files changed, 255 insertions(+), 226 deletions(-) (limited to 'src/Cech_complex/test/test_cech_complex.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 71c88982..83ef9dca 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -93,12 +93,12 @@ int main(int argc, char * argv[]) { Radius_distance()); std::cout << radius_clock << std::endl; - Gudhi::Clock squared_radius_clock("Gudhi::Radius_distance()"); + Gudhi::Clock common_radius_clock("Gudhi::Radius_distance()"); // Compute the proximity graph of the points Proximity_graph sq_radius_prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), threshold, Gudhi::Radius_distance()); - std::cout << squared_radius_clock << std::endl; + std::cout << common_radius_clock << std::endl; boost::filesystem::path full_path(boost::filesystem::current_path()); @@ -116,6 +116,7 @@ int main(int argc, char * argv[]) { if ( itr->path().extension() == ".off" ) // see below { Points_off_reader off_reader(itr->path().string()); + Point p0 = off_reader.get_point_cloud()[0]; for (Filtration_value radius = 0.1; radius < 0.4; radius += 0.1) { std::cout << itr->path().stem() << ";"; @@ -123,7 +124,7 @@ int main(int argc, char * argv[]) { 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); + rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); // ------------------------------------------ // Display information about the Rips complex // ------------------------------------------ @@ -133,7 +134,7 @@ int main(int argc, char * argv[]) { 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); + cech_complex_from_points.create_complex(cech_stree, p0.size() - 1); // ------------------------------------------ // Display information about the Cech complex // ------------------------------------------ @@ -141,6 +142,7 @@ int main(int argc, char * argv[]) { std::cout << cech_sec << ";"; std::cout << cech_sec / rips_sec << ";"; + assert(rips_stree.num_simplices() >= cech_stree.num_simplices()); std::cout << rips_stree.num_simplices() << ";"; std::cout << cech_stree.num_simplices() << ";" << std::endl; } diff --git a/src/Cech_complex/doc/Intro_cech_complex.h b/src/Cech_complex/doc/Intro_cech_complex.h index f2052763..8b6c7851 100644 --- a/src/Cech_complex/doc/Intro_cech_complex.h +++ b/src/Cech_complex/doc/Intro_cech_complex.h @@ -29,7 +29,7 @@ namespace cech_complex { /** \defgroup cech_complex Cech complex * - * \author Clément Maria, Pawel Dlotko, Vincent Rouvreau + * \author Vincent Rouvreau * * @{ * @@ -40,40 +40,54 @@ namespace cech_complex { * proximity graph that allows to construct a * simplicial complex * from it. - * The input can be a point cloud with a given distance function. + * The input shall be a point cloud in an Euclidean space. * - * The filtration value of each edge is computed from a user-given distance function. + * The filtration value of each edge of the `Gudhi::Proximity_graph` is computed from `Gudhi::Radius_distance` function. * - * All edges that have a filtration value strictly greater than a given threshold value are not inserted into - * the complex. + * All edges that have a filtration value strictly greater than a user given maximal radius value, \f$max\_radius\f$, + * are not inserted into the complex. * - * When creating a simplicial complex from this proximity graph, Cech inserts the proximity graph into the data - * structure, and then expands the simplicial complex when required. - * * Vertex name correspond to the index of the point in the given range (aka. the point cloud). * - * \image html "cech_complex_representation.png" "Cech complex proximity graph representation" - * - * On this example, as edges (4,5), (4,6) and (5,6) are in the complex, simplex (4,5,6) is added with the filtration - * value set with \f$max(filtration(4,5), filtration(4,6), filtration(5,6))\f$. - * And so on for simplex (0,1,2,3). + * \image html "cech_one_skeleton.png" "Cech complex proximity graph representation" * + * When creating a simplicial complex from this proximity graph, Cech inserts the proximity graph into the simplicial + * complex data structure, and then expands the simplicial complex when required. + * + * On this example, as edges \f$(x,y)\f$, \f$(y,z)\f$ and \f$(z,y)\f$ are in the complex, the minimal ball radius + * containing the points \f$(x,y,z)\f$ is computed. + * + * \f$(x,y,z)\f$ is inserted to the simplicial complex with the filtration value set with + * \f$mini\_ball\_radius(x,y,z))\f$ iff \f$mini\_ball\_radius(x,y,z)) \leq max\_radius\f$. + * + * And so on for higher dimensions. + * + * \image html "cech_complex_representation.png" "Cech complex expansion" + * + * The minimal ball radius computation is insured by + * + * the miniball software (V3.0) - Smallest Enclosing Balls of Points - and distributed with GUDHI. + * + * Please refer to + * + * the miniball software design description for more information about this computation. + * * If the Cech_complex interfaces are not detailed enough for your need, please refer to - * - * cech_persistence_step_by_step.cpp example, where the graph construction over the Simplex_tree is more detailed. + * + * cech_complex_step_by_step.cpp example, where the graph construction over the Simplex_tree is more detailed. * - * \section cechpointsdistance Point cloud and distance function + * \section cechpointsdistance Point cloud * - * \subsection cechpointscloudexample Example from a point cloud and a distance function + * \subsection cechpointscloudexample Example from a point cloud * - * This example builds the proximity graph from the given points, threshold value, and distance function. + * This example builds the proximity graph from the given points, and maximal radius values. * Then it creates a `Simplex_tree` with it. * * Then, it is asked to display information about the simplicial complex. * * \include Cech_complex/cech_complex_example_from_points.cpp * - * When launching (Cech maximal distance between 2 points is 7.1, is expanded until dimension 2): + * When launching (Cech maximal distance between 2 points is 1., is expanded until dimension 2): * * \code $> ./Cech_complex_example_from_points * \endcode diff --git a/src/Cech_complex/doc/cech_complex_representation.ipe b/src/Cech_complex/doc/cech_complex_representation.ipe index 7f6028f4..c64d7596 100644 --- a/src/Cech_complex/doc/cech_complex_representation.ipe +++ b/src/Cech_complex/doc/cech_complex_representation.ipe @@ -1,7 +1,7 @@ - + @@ -232,95 +232,99 @@ h - -109.771 601.912 m -159.595 601.797 l -140.058 541.915 l + +48 640 m +80 672 l +48 672 l h - -79.8776 552.169 m -109.756 601.699 l -139.812 542.209 l +Cech complex +0 +1 +2 +3 +4 +5 +6 + + + + + + + + + + + + + +32 0 0 32 304 672 e + + +304 672 m +336 672 l + +Maximal radius +7 +8 +9 + +112 576 m +144 608 l + + +144 672 m +144 608 l +200 640 l h - -69.8453 682.419 m -159.925 712.208 l -90.12 732.039 l + +48 576 m +112 576 l +80 544 l h -Rips complex -0 -1 -2 -3 -4 -5 -6 - -60 710 m -40 660 l - - -40 660 m -130 690 l - - -130 690 m -60 710 l - - -40 660 m -80 580 l - - -80 580 m -130 580 l -130 580 l - - -130 580 m -110 520 l - - -110 520 m -50 530 l -50 530 l -50 530 l + + +80 672 m +144 672 l +112 728 l +h - -50 530 m -80 580 l + + + + +48 576 m +48 640 l +32 608 l +h - -130 580 m -130 690 l + + + + + + + +32 0 0 32 80 576 e - - - - - - -150.038 609.9 m -179.929 549.727 l + +22.6274 0 0 22.6274 64 656 e - - - -158.7 593.269 m -81.4925 544.805 l + +37.1429 0 0 37.1429 112 690.857 e - -256.324 639.958 m -370.055 639.958 l + +37.1429 0 0 37.1429 162.857 640 e - -56.8567 0 0 56.8567 313.217 639.756 e + +10 + +32 0 0 32 48 608 e - - -Rips threshold + + diff --git a/src/Cech_complex/doc/cech_complex_representation.png b/src/Cech_complex/doc/cech_complex_representation.png index e901d92e..4d103a56 100644 Binary files a/src/Cech_complex/doc/cech_complex_representation.png and b/src/Cech_complex/doc/cech_complex_representation.png differ diff --git a/src/Cech_complex/doc/cech_one_skeleton.ipe b/src/Cech_complex/doc/cech_one_skeleton.ipe index 3a35970c..345e6d7b 100644 --- a/src/Cech_complex/doc/cech_one_skeleton.ipe +++ b/src/Cech_complex/doc/cech_one_skeleton.ipe @@ -1,7 +1,7 @@ - + @@ -232,95 +232,83 @@ h - -109.771 601.912 m -159.595 601.797 l -140.058 541.915 l +Proximity graph +0 +1 + +304 672 m +336 672 l + +2 +3 +4 +5 +6 + + + + + + + + + + + +32 0 0 32 304 672 e + +Maximal radius +7 +8 +9 + +112 576 m +144 608 l + + +144 672 m +144 608 l +200 640 l h - -79.8776 552.169 m -109.756 601.699 l -139.812 542.209 l + +48 640 m +80 672 l +48 672 l h - -69.8453 682.419 m -159.925 712.208 l -90.12 732.039 l + +48 576 m +112 576 l +80 544 l h -One skeleton graph -0 -1 -2 -3 -4 -5 -6 - -60 710 m -40 660 l - - -40 660 m -130 690 l - - -130 690 m -60 710 l - - -40 660 m -80 580 l - - -80 580 m -130 580 l -130 580 l - - -130 580 m -110 520 l - - -110 520 m -50 530 l -50 530 l -50 530 l - - -50 530 m -80 580 l - - -130 580 m -130 690 l - - - - - - - -150.038 609.9 m -179.929 549.727 l - - - - -158.7 593.269 m -81.4925 544.805 l - - -256.324 639.958 m -370.055 639.958 l + + +80 672 m +144 672 l +112 728 l +h - -56.8567 0 0 56.8567 313.217 639.756 e + + +48 576 m +48 640 l +32 608 l +h - - -Rips threshold + + + + + + + + + + +10 + + diff --git a/src/Cech_complex/doc/cech_one_skeleton.png b/src/Cech_complex/doc/cech_one_skeleton.png index ffa9c329..807e0936 100644 Binary files a/src/Cech_complex/doc/cech_one_skeleton.png and b/src/Cech_complex/doc/cech_one_skeleton.png differ 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 97327e69..3b889d56 100644 --- a/src/Cech_complex/example/cech_complex_example_from_points.cpp +++ b/src/Cech_complex/example/cech_complex_example_from_points.cpp @@ -15,22 +15,26 @@ int main() { using Cech_complex = Gudhi::cech_complex::Cech_complex; Point_cloud points; - 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.}); + 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 // ---------------------------------------------------------------------------- // Init of a Cech complex from points // ---------------------------------------------------------------------------- - // 5. is a magic number to force one blocker, and one non-blocker - Filtration_value max_radius = 12.; + Filtration_value max_radius = 1.; Cech_complex cech_complex_from_points(points, max_radius); Simplex_tree stree; - cech_complex_from_points.create_complex(stree, -1); + cech_complex_from_points.create_complex(stree, 2); // ---------------------------------------------------------------------------- // Display information about the one skeleton Cech complex // ---------------------------------------------------------------------------- diff --git a/src/Cech_complex/example/cech_complex_example_from_points_for_doc.txt b/src/Cech_complex/example/cech_complex_example_from_points_for_doc.txt index 684e120b..be0afc76 100644 --- a/src/Cech_complex/example/cech_complex_example_from_points_for_doc.txt +++ b/src/Cech_complex/example/cech_complex_example_from_points_for_doc.txt @@ -1,16 +1,31 @@ -Cech complex is of dimension 2 - 14 simplices - 7 vertices. Iterator on Cech complex simplices in the filtration order, with [filtration value]: - ( 0 ) -> [0] - ( 1 ) -> [0] - ( 2 ) -> [0] - ( 3 ) -> [0] - ( 4 ) -> [0] - ( 5 ) -> [0] - ( 6 ) -> [0] - ( 3 2 ) -> [5] - ( 5 4 ) -> [5.38516] - ( 2 0 ) -> [5.83095] - ( 1 0 ) -> [6.08276] - ( 3 1 ) -> [6.32456] - ( 2 1 ) -> [6.7082] - ( 3 2 1 ) -> [7.07107] + ( 0 ) -> [0] + ( 1 ) -> [0] + ( 2 ) -> [0] + ( 3 ) -> [0] + ( 4 ) -> [0] + ( 5 ) -> [0] + ( 6 ) -> [0] + ( 7 ) -> [0] + ( 8 ) -> [0] + ( 9 ) -> [0] + ( 10 ) -> [0] + ( 9 4 ) -> [0.5] + ( 9 6 ) -> [0.5] + ( 10 1 ) -> [0.559017] + ( 10 4 ) -> [0.559017] + ( 1 0 ) -> [0.707107] + ( 2 0 ) -> [0.707107] + ( 3 2 ) -> [0.707107] + ( 6 4 ) -> [0.707107] + ( 9 6 4 ) -> [0.707107] + ( 2 1 ) -> [1] + ( 2 1 0 ) -> [1] + ( 4 1 ) -> [1] + ( 5 3 ) -> [1] + ( 7 3 ) -> [1] + ( 7 5 ) -> [1] + ( 7 6 ) -> [1] + ( 8 6 ) -> [1] + ( 8 7 ) -> [1] + ( 10 4 1 ) -> [1] diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index a50ed9fa..612c73c3 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -23,7 +23,7 @@ #ifndef CECH_COMPLEX_H_ #define CECH_COMPLEX_H_ -#include // for Gudhi::Squared_radius +#include // for Gudhi::Radius_distance #include // for Gudhi::Proximity_graph #include // for GUDHI_CHECK #include // for Gudhi::cech_complex::Cech_blocker @@ -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 max_radius. Edge length is computed from `Gudhi::Squared_radius` distance function. + * to a given max_radius. Edge length is computed from `Gudhi::Radius_distance` distance function. * * \tparam SimplicialComplexForProximityGraph furnishes `Vertex_handle` and `Filtration_value` type definition required * by `Gudhi::Proximity_graph`. diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index d718b56e..5ba17c51 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -24,7 +24,7 @@ #define CECH_COMPLEX_BLOCKER_H_ #include // Cech_blocker is using a pointer on Gudhi::cech_complex::Cech_complex -#include // for Gudhi::Squared_radius +#include // for Gudhi::Radius_distance #include #include @@ -75,12 +75,13 @@ class Cech_blocker { std::cout << "#(" << vertex << ")#"; #endif // DEBUG_TRACES } - Filtration_value squared_radius = Gudhi::Radius_distance()(points); + Filtration_value radius = Gudhi::Radius_distance()(points); #ifdef DEBUG_TRACES - std::cout << "squared_radius = " << squared_radius << " - " << (squared_radius > cc_ptr_->max_radius()) << std::endl; + if (radius > cc_ptr_->max_radius()) + std::cout << "radius > max_radius => expansion is blocked\n"; #endif // DEBUG_TRACES - simplicial_complex_.assign_filtration(sh, squared_radius); - return (squared_radius > cc_ptr_->max_radius()); + simplicial_complex_.assign_filtration(sh, radius); + return (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 626f1d82..eae8778c 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -86,7 +86,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { BOOST_CHECK(st.num_vertices() == NUMBER_OF_VERTICES); std::cout << "st.num_simplices()=" << st.num_simplices() << std::endl; - BOOST_CHECK(st.num_simplices() == 18); + BOOST_CHECK(st.num_simplices() == 28); // Check filtration values of vertices is 0.0 for (auto f_simplex : st.skeleton_simplex_range(0)) { @@ -119,7 +119,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { BOOST_CHECK(st2.num_vertices() == NUMBER_OF_VERTICES); std::cout << "st2.num_simplices()=" << st2.num_simplices() << std::endl; - BOOST_CHECK(st2.num_simplices() == 23); + BOOST_CHECK(st2.num_simplices() == 63); Point_cloud points012; for (std::size_t vertex = 0; vertex <= 2; vertex++) { @@ -156,7 +156,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { BOOST_CHECK(st3.num_vertices() == NUMBER_OF_VERTICES); std::cout << "st3.num_simplices()=" << st3.num_simplices() << std::endl; - BOOST_CHECK(st3.num_simplices() == 24); + BOOST_CHECK(st3.num_simplices() == 98); Point_cloud points0123; for (std::size_t vertex = 0; vertex <= 3; vertex++) { @@ -236,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), 0.5); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.707107, .00001); break; case 2: - GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.666667, .00001); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.816497, .00001); break; case 3: - GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.75); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.866025, .00001); break; default: BOOST_CHECK(false); // Shall not happen diff --git a/src/Doxyfile b/src/Doxyfile index f1981e2e..52de65b0 100644 --- a/src/Doxyfile +++ b/src/Doxyfile @@ -843,7 +843,8 @@ EXAMPLE_RECURSIVE = NO IMAGE_PATH = doc/Skeleton_blocker/ \ doc/Alpha_complex/ \ doc/common/ \ - doc/Contraction/ \ + doc/Cech_complex/ \ + doc/Contraction/ \ doc/Simplex_tree/ \ doc/Persistent_cohomology/ \ doc/Witness_complex/ \ diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index 3ce51ad1..7749b98e 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -90,11 +90,11 @@ class Radius_distance { operator()(const Point_cloud& point_cloud) const { using Min_sphere = Miniball::Miniball>; - //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()); + Min_sphere ms(point_cloud.begin()->end() - point_cloud.begin()->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; + std::cout << "Radius_distance = " << std::sqrt(ms.squared_radius()) << " | nb points = " + << point_cloud.end() - point_cloud.begin() << " | dimension = " + << point_cloud.begin()->end() - point_cloud.begin()->begin() << std::endl; #endif // DEBUG_TRACES return std::sqrt(ms.squared_radius()); -- cgit v1.2.3 From 2a9f6eb628979a9634c647be93e575e3177c15da Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Mon, 5 Mar 2018 15:43:13 +0000 Subject: Fix documentation git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/cechcomplex_vincent@3264 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: af6241dfe42d988b542224b9dae7364f3f22b6fd --- src/Cech_complex/doc/Intro_cech_complex.h | 4 +- .../doc/cech_complex_representation.ipe | 4 +- .../doc/cech_complex_representation.png | Bin 54399 -> 39938 bytes src/Cech_complex/doc/cech_one_skeleton.ipe | 4 +- src/Cech_complex/doc/cech_one_skeleton.png | Bin 29354 -> 24662 bytes src/Cech_complex/test/test_cech_complex.cpp | 124 ++++++++++----------- src/common/doc/main_page.h | 18 ++- 7 files changed, 85 insertions(+), 69 deletions(-) (limited to 'src/Cech_complex/test/test_cech_complex.cpp') diff --git a/src/Cech_complex/doc/Intro_cech_complex.h b/src/Cech_complex/doc/Intro_cech_complex.h index 8b6c7851..ec0d35e2 100644 --- a/src/Cech_complex/doc/Intro_cech_complex.h +++ b/src/Cech_complex/doc/Intro_cech_complex.h @@ -67,11 +67,13 @@ namespace cech_complex { * The minimal ball radius computation is insured by * * the miniball software (V3.0) - Smallest Enclosing Balls of Points - and distributed with GUDHI. - * * Please refer to * * the miniball software design description for more information about this computation. * + * This radius computation is the reason why the Cech_complex is taking much more time to be computed than the + * \ref rips_complex but it offers more topological guarantees. + * * If the Cech_complex interfaces are not detailed enough for your need, please refer to * * cech_complex_step_by_step.cpp example, where the graph construction over the Simplex_tree is more detailed. diff --git a/src/Cech_complex/doc/cech_complex_representation.ipe b/src/Cech_complex/doc/cech_complex_representation.ipe index c64d7596..377745a3 100644 --- a/src/Cech_complex/doc/cech_complex_representation.ipe +++ b/src/Cech_complex/doc/cech_complex_representation.ipe @@ -1,7 +1,7 @@ - + @@ -267,7 +267,7 @@ h Maximal radius 7 -8 +8 9 112 576 m diff --git a/src/Cech_complex/doc/cech_complex_representation.png b/src/Cech_complex/doc/cech_complex_representation.png index 4d103a56..d0eb85a5 100644 Binary files a/src/Cech_complex/doc/cech_complex_representation.png and b/src/Cech_complex/doc/cech_complex_representation.png differ diff --git a/src/Cech_complex/doc/cech_one_skeleton.ipe b/src/Cech_complex/doc/cech_one_skeleton.ipe index 345e6d7b..ed66e132 100644 --- a/src/Cech_complex/doc/cech_one_skeleton.ipe +++ b/src/Cech_complex/doc/cech_one_skeleton.ipe @@ -1,7 +1,7 @@ - + @@ -259,7 +259,7 @@ h Maximal radius 7 -8 +8 9 112 576 m diff --git a/src/Cech_complex/doc/cech_one_skeleton.png b/src/Cech_complex/doc/cech_one_skeleton.png index 807e0936..cc636616 100644 Binary files a/src/Cech_complex/doc/cech_one_skeleton.png and b/src/Cech_complex/doc/cech_one_skeleton.png differ diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index eae8778c..8cbfe431 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -42,7 +42,7 @@ // Type definitions using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; -using Point = std::vector; +using Point = std::vector; using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; using Cech_complex = Gudhi::cech_complex::Cech_complex; @@ -51,42 +51,52 @@ using Point_iterator = Point_cloud::const_iterator; using Coordinate_iterator = Point::const_iterator; using Min_sphere = Miniball::Miniball>; -BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { +BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { // ---------------------------------------------------------------------------- // - // Init of a Cech complex from a OFF file + // Init of a Cech complex from a point cloud // // ---------------------------------------------------------------------------- - std::string off_file_name("alphacomplexdoc.off"); - double max_radius = 12.0; - std::cout << "========== OFF FILE NAME = " << off_file_name << " - Cech max_radius=" << + 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 + + Filtration_value max_radius = 1.0; + std::cout << "========== NUMBER OF POINTS = " << points.size() << " - 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, max_radius); + Cech_complex cech_complex_for_doc(points, max_radius); - GUDHI_TEST_FLOAT_EQUALITY_CHECK(cech_complex_from_file.max_radius(), max_radius); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(cech_complex_for_doc.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))); + for (; i < points.size(); i++) { + BOOST_CHECK(points[i] == *(cech_complex_for_doc.point_iterator(i))); } #ifdef GUDHI_DEBUG - BOOST_CHECK_THROW (cech_complex_from_file.point_iterator(i+1), std::out_of_range); + BOOST_CHECK_THROW (cech_complex_for_doc.point_iterator(i+1), std::out_of_range); #endif // GUDHI_DEBUG const int DIMENSION_1 = 1; Simplex_tree st; - cech_complex_from_file.create_complex(st, DIMENSION_1); + cech_complex_for_doc.create_complex(st, DIMENSION_1); std::cout << "st.dimension()=" << st.dimension() << std::endl; BOOST_CHECK(st.dimension() == DIMENSION_1); - const int NUMBER_OF_VERTICES = 7; + const int NUMBER_OF_VERTICES = 11; std::cout << "st.num_vertices()=" << st.num_vertices() << std::endl; BOOST_CHECK(st.num_vertices() == NUMBER_OF_VERTICES); std::cout << "st.num_simplices()=" << st.num_simplices() << std::endl; - BOOST_CHECK(st.num_simplices() == 28); + BOOST_CHECK(st.num_simplices() == 27); // Check filtration values of vertices is 0.0 for (auto f_simplex : st.skeleton_simplex_range(0)) { @@ -100,7 +110,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { std::cout << "vertex = ("; for (auto vertex : st.simplex_vertex_range(f_simplex)) { std::cout << vertex << ","; - vp.push_back(off_reader.get_point_cloud().at(vertex)); + vp.push_back(points.at(vertex)); } std::cout << ") - distance =" << Gudhi::Radius_distance()(vp.at(0), vp.at(1)) << " - filtration =" << st.filtration(f_simplex) << std::endl; @@ -110,8 +120,13 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { } const int DIMENSION_2 = 2; + +#ifdef GUDHI_DEBUG + BOOST_CHECK_THROW (cech_complex_for_doc.create_complex(st, DIMENSION_2), std::invalid_argument); +#endif + Simplex_tree st2; - cech_complex_from_file.create_complex(st2, DIMENSION_2); + cech_complex_for_doc.create_complex(st2, DIMENSION_2); std::cout << "st2.dimension()=" << st2.dimension() << std::endl; BOOST_CHECK(st2.dimension() == DIMENSION_2); @@ -119,14 +134,14 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { BOOST_CHECK(st2.num_vertices() == NUMBER_OF_VERTICES); std::cout << "st2.num_simplices()=" << st2.num_simplices() << std::endl; - BOOST_CHECK(st2.num_simplices() == 63); + BOOST_CHECK(st2.num_simplices() == 30); Point_cloud points012; for (std::size_t vertex = 0; vertex <= 2; vertex++) { - points012.push_back(Point(cech_complex_from_file.point_iterator(vertex)->begin(), - cech_complex_from_file.point_iterator(vertex)->end())); + points012.push_back(Point(cech_complex_for_doc.point_iterator(vertex)->begin(), + cech_complex_for_doc.point_iterator(vertex)->end())); } - std::size_t dimension = point_cloud[0].end() - point_cloud[0].begin(); + std::size_t dimension = points[0].end() - points[0].begin(); Min_sphere ms012(dimension, points012.begin(),points012.end()); Simplex_tree::Filtration_value f012 = st2.filtration(st2.find({0, 1, 2})); @@ -134,53 +149,36 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_file) { GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(ms012.squared_radius())); - Point_cloud points456; - for (std::size_t vertex = 4; vertex <= 6; vertex++) { - points456.push_back(Point(cech_complex_from_file.point_iterator(vertex)->begin(), - cech_complex_from_file.point_iterator(vertex)->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; + Point_cloud points1410; + points1410.push_back(Point(cech_complex_for_doc.point_iterator(1)->begin(), + cech_complex_for_doc.point_iterator(1)->end())); + points1410.push_back(Point(cech_complex_for_doc.point_iterator(4)->begin(), + cech_complex_for_doc.point_iterator(4)->end())); + points1410.push_back(Point(cech_complex_for_doc.point_iterator(10)->begin(), + cech_complex_for_doc.point_iterator(10)->end())); + Min_sphere ms1410(dimension, points1410.begin(),points1410.end()); - GUDHI_TEST_FLOAT_EQUALITY_CHECK(f456, std::sqrt(ms456.squared_radius())); + Simplex_tree::Filtration_value f1410 = st2.filtration(st2.find({1, 4, 10})); + std::cout << "f1410= " << f1410 << " | ms1410_radius= " << std::sqrt(ms1410.squared_radius()) << std::endl; - const int DIMENSION_3 = 3; - Simplex_tree st3; - cech_complex_from_file.create_complex(st3, DIMENSION_3); - std::cout << "st3.dimension()=" << st3.dimension() << std::endl; - BOOST_CHECK(st3.dimension() == DIMENSION_3); - - std::cout << "st3.num_vertices()=" << st3.num_vertices() << std::endl; - BOOST_CHECK(st3.num_vertices() == NUMBER_OF_VERTICES); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, std::sqrt(ms1410.squared_radius())); - std::cout << "st3.num_simplices()=" << st3.num_simplices() << std::endl; - BOOST_CHECK(st3.num_simplices() == 98); - - Point_cloud points0123; - for (std::size_t vertex = 0; vertex <= 3; vertex++) { - points0123.push_back(Point(cech_complex_from_file.point_iterator(vertex)->begin(), - cech_complex_from_file.point_iterator(vertex)->end())); - } - Min_sphere ms0123(dimension, points0123.begin(),points0123.end()); + Point_cloud points469; + points469.push_back(Point(cech_complex_for_doc.point_iterator(4)->begin(), + cech_complex_for_doc.point_iterator(4)->end())); + points469.push_back(Point(cech_complex_for_doc.point_iterator(6)->begin(), + cech_complex_for_doc.point_iterator(6)->end())); + points469.push_back(Point(cech_complex_for_doc.point_iterator(9)->begin(), + cech_complex_for_doc.point_iterator(9)->end())); + Min_sphere ms469(dimension, points469.begin(),points469.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; + Simplex_tree::Filtration_value f469 = st2.filtration(st2.find({4, 6, 9})); + std::cout << "f469= " << f469 << " | ms469_radius= " << std::sqrt(ms469.squared_radius()) << std::endl; - GUDHI_TEST_FLOAT_EQUALITY_CHECK(f0123, std::sqrt(ms0123.squared_radius())); - - - - Point_cloud points01; - for (std::size_t vertex = 0; vertex <= 1; vertex++) { - points01.push_back(Point(cech_complex_from_file.point_iterator(vertex)->begin(), - cech_complex_from_file.point_iterator(vertex)->end())); - } - Min_sphere ms01(dimension, points01.begin(),points01.end()); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f469, std::sqrt(ms469.squared_radius())); - Simplex_tree::Filtration_value f01 = st2.filtration(st2.find({0, 1})); - std::cout << "f01= " << f01 << " | ms01_radius= " << std::sqrt(ms01.squared_radius()) << std::endl; + BOOST_CHECK((st2.find({6, 7, 8}) == st2.null_simplex())); + BOOST_CHECK((st2.find({3, 5, 7}) == st2.null_simplex())); } diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h index b3e9ea03..3851dde7 100644 --- a/src/common/doc/main_page.h +++ b/src/common/doc/main_page.h @@ -41,6 +41,22 @@ User manual: \ref alpha_complex - Reference manual: Gudhi::alpha_complex::Alpha_complex + + \subsection CechComplexDataStructure Cech complex + \image html "cech_complex_representation.png" "Cech complex representation" + + + + +
+ Author: Vincent Rouvreau
+ Introduced in: GUDHI 2.2.0
+ Copyright: GPL v3
+ 		+ The Cech complex is a proximity graph that allows to construct a simplicial complex from it. + It guarantees the simplices inserted are inside a ball of a given radius.
+ User manual: \ref cech_complex - Reference manual: Gudhi::cech_complex::Cech_complex +
\subsection CubicalComplexDataStructure Cubical complex \image html "Cubical_complex_representation.png" "Cubical complex representation" @@ -57,6 +73,7 @@ User manual: \ref cubical_complex - Reference manual: Gudhi::cubical_complex::Bitmap_cubical_complex + \subsection RipsComplexDataStructure Rips complex \image html "rips_complex_representation.png" "Rips complex representation" @@ -74,7 +91,6 @@ User manual: \ref rips_complex - Reference manual: Gudhi::rips_complex::Rips_complex -
\subsection SimplexTreeDataStructure Simplex tree \image html "Simplex_tree_representation.png" "Simplex tree representation" -- cgit v1.2.3 From 0d0ca116e7fef77cc950b7e85380495661311d91 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Thu, 22 Mar 2018 08:36:00 +0000 Subject: Move Miniball in gudhi git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/cechcomplex_vincent@3302 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 736e5c33a9593e4dfb5b19ddc745ab9852d71b40 --- src/CMakeLists.txt | 1 - .../benchmark/cech_complex_benchmark.cpp | 5 +- src/Cech_complex/example/CMakeLists.txt | 2 + .../example/cech_complex_step_by_step.cpp | 4 +- .../include/Miniball/Miniball.COPYRIGHT | 4 - src/Cech_complex/include/Miniball/Miniball.README | 23 - src/Cech_complex/include/Miniball/Miniball.hpp | 520 -------------------- src/Cech_complex/include/gudhi/Miniball.COPYRIGHT | 4 + src/Cech_complex/include/gudhi/Miniball.README | 23 + src/Cech_complex/include/gudhi/Miniball.hpp | 523 +++++++++++++++++++++ src/Cech_complex/test/test_cech_complex.cpp | 5 +- src/common/include/gudhi/distance_functions.h | 2 +- 12 files changed, 559 insertions(+), 557 deletions(-) delete mode 100644 src/Cech_complex/include/Miniball/Miniball.COPYRIGHT delete mode 100644 src/Cech_complex/include/Miniball/Miniball.README delete mode 100644 src/Cech_complex/include/Miniball/Miniball.hpp create mode 100644 src/Cech_complex/include/gudhi/Miniball.COPYRIGHT create mode 100644 src/Cech_complex/include/gudhi/Miniball.README create mode 100644 src/Cech_complex/include/gudhi/Miniball.hpp (limited to 'src/Cech_complex/test/test_cech_complex.cpp') diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ff7acafc..7cccb19f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -105,5 +105,4 @@ install(FILES # install the include file on "make install" install(DIRECTORY include/gudhi DESTINATION include) -install(DIRECTORY include/Miniball DESTINATION include) #--------------------------------------------------------------------------------------- diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 83ef9dca..6ba52419 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -27,8 +27,7 @@ #include #include #include - -#include +#include #include "boost/filesystem.hpp" // includes all needed Boost.Filesystem declarations @@ -57,7 +56,7 @@ class Radius_distance { using Point_cloud = std::vector; using Point_iterator = typename Point_cloud::const_iterator; using Coordinate_iterator = typename Point::const_iterator; - using Min_sphere = typename Miniball::Miniball>; + using Min_sphere = typename Gudhi::Miniball::Miniball>; Point_cloud point_cloud; point_cloud.push_back(p1); diff --git a/src/Cech_complex/example/CMakeLists.txt b/src/Cech_complex/example/CMakeLists.txt index ac32ff95..ab391215 100644 --- a/src/Cech_complex/example/CMakeLists.txt +++ b/src/Cech_complex/example/CMakeLists.txt @@ -6,6 +6,8 @@ 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 $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-d" "3") add_executable ( Cech_complex_example_from_points cech_complex_example_from_points.cpp) if (TBB_FOUND) 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 8705a3e5..0e7c4bbd 100644 --- a/src/Cech_complex/example/cech_complex_step_by_step.cpp +++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp @@ -25,7 +25,7 @@ #include #include -#include +#include #include @@ -55,7 +55,7 @@ class Cech_blocker { using Point_cloud = std::vector; using Point_iterator = Point_cloud::const_iterator; using Coordinate_iterator = Point::const_iterator; - using Min_sphere = Miniball::Miniball >; + using Min_sphere = Gudhi::Miniball::Miniball >; public: bool operator()(Simplex_handle sh) { std::vector points; diff --git a/src/Cech_complex/include/Miniball/Miniball.COPYRIGHT b/src/Cech_complex/include/Miniball/Miniball.COPYRIGHT deleted file mode 100644 index dbe4c553..00000000 --- a/src/Cech_complex/include/Miniball/Miniball.COPYRIGHT +++ /dev/null @@ -1,4 +0,0 @@ -The miniball software is available under the GNU General Public License (GPLv3 - https://www.gnu.org/copyleft/gpl.html). -If your intended use is not compliant with this license, please buy a commercial license (EUR 500 - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/license.html). -You need a license if the software that you develop using Miniball V3.0 is not open source. - diff --git a/src/Cech_complex/include/Miniball/Miniball.README b/src/Cech_complex/include/Miniball/Miniball.README deleted file mode 100644 index 86a96f08..00000000 --- a/src/Cech_complex/include/Miniball/Miniball.README +++ /dev/null @@ -1,23 +0,0 @@ -https://people.inf.ethz.ch/gaertner/subdir/software/miniball.html - -Smallest Enclosing Balls of Points - Fast and Robust in C++. -(high-quality software for smallest enclosing balls of balls is available in the computational geometry algorithms library CGAL) - - -This is the miniball software (V3.0) for computing smallest enclosing balls of points in arbitrary dimensions. It consists of a C++ header file Miniball.hpp (around 500 lines of code) and two example programs miniball_example.cpp and miniball_example_containers.cpp that demonstrate the usage. The first example stores the coordinates of the input points in a two-dimensional array, the second example uses a list of vectors to show how generic containers can be used. - -Credits: Aditya Gupta and Alexandros Konstantinakis-Karmis have significantly contributed to this version of the software. - -Changes - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/changes.txt - from previous versions. - -The theory - https://people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf - behind the miniball software (Proc. 7th Annual European Symposium on Algorithms (ESA), Lecture Notes in Computer Science 1643, Springer-Verlag, pp.325-338, 1999). - -Main Features: - - Very fast in low dimensions. 1 million points in 5-space are processed within 0.05 seconds on any recent machine. - - High numerical stability. Almost all input degeneracies (cospherical points, multiple points, points very close together) are routinely handled. - - Easily integrates into your code. You can freely choose the coordinate type of your points and the container to store the points. If you still need to adapt the code, the header is small and readable and contains documentation for all major methods. - - diff --git a/src/Cech_complex/include/Miniball/Miniball.hpp b/src/Cech_complex/include/Miniball/Miniball.hpp deleted file mode 100644 index a42d62a7..00000000 --- a/src/Cech_complex/include/Miniball/Miniball.hpp +++ /dev/null @@ -1,520 +0,0 @@ -// Copright (C) 1999-2013, Bernd Gaertner -// $Rev: 3581 $ -// -// 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 . -// -// Contact: -// -------- -// Bernd Gaertner -// Institute of Theoretical Computer Science -// ETH Zuerich -// CAB G31.1 -// CH-8092 Zuerich, Switzerland -// http://www.inf.ethz.ch/personal/gaertner - -#ifndef MINIBALL_HPP_ -#define MINIBALL_HPP_ - -#include -#include -#include -#include -#include - -namespace Miniball { - - // Global Functions - // ================ - template - inline NT mb_sqr (NT r) {return r*r;} - - // Functors - // ======== - - // functor to map a point iterator to the corresponding coordinate iterator; - // generic version for points whose coordinate containers have begin() - template < typename Pit_, typename Cit_ > - struct CoordAccessor { - typedef Pit_ Pit; - typedef Cit_ Cit; - inline Cit operator() (Pit it) const { return (*it).begin(); } - }; - - // partial specialization for points whose coordinate containers are arrays - template < typename Pit_, typename Cit_ > - struct CoordAccessor { - typedef Pit_ Pit; - typedef Cit_* Cit; - inline Cit operator() (Pit it) const { return *it; } - }; - - // Class Declaration - // ================= - - template - class Miniball { - private: - // types - // The iterator type to go through the input points - typedef typename CoordAccessor::Pit Pit; - // The iterator type to go through the coordinates of a single point. - typedef typename CoordAccessor::Cit Cit; - // The coordinate type - typedef typename std::iterator_traits::value_type NT; - // The iterator to go through the support points - typedef typename std::list::iterator Sit; - - // data members... - const int d; // dimension - Pit points_begin; - Pit points_end; - CoordAccessor coord_accessor; - double time; - const NT nt0; // NT(0) - - //...for the algorithms - std::list L; - Sit support_end; - int fsize; // number of forced points - int ssize; // number of support points - - // ...for the ball updates - NT* current_c; - NT current_sqr_r; - NT** c; - NT* sqr_r; - - // helper arrays - NT* q0; - NT* z; - NT* f; - NT** v; - NT** a; - - public: - // The iterator type to go through the support points - typedef typename std::list::const_iterator SupportPointIterator; - - // PRE: [begin, end) is a nonempty range - // POST: computes the smallest enclosing ball of the points in the range - // [begin, end); the functor a maps a point iterator to an iterator - // through the d coordinates of the point - Miniball (int d_, Pit begin, Pit end, CoordAccessor ca = CoordAccessor()); - - // POST: returns a pointer to the first element of an array that holds - // the d coordinates of the center of the computed ball - const NT* center () const; - - // POST: returns the squared radius of the computed ball - NT squared_radius () const; - - // POST: returns the number of support points of the computed ball; - // the support points form a minimal set with the same smallest - // enclosing ball as the input set; in particular, the support - // points are on the boundary of the computed ball, and their - // number is at most d+1 - int nr_support_points () const; - - // POST: returns an iterator to the first support point - SupportPointIterator support_points_begin () const; - - // POST: returns a past-the-end iterator for the range of support points - SupportPointIterator support_points_end () const; - - // POST: returns the maximum excess of any input point w.r.t. the computed - // ball, divided by the squared radius of the computed ball. The - // excess of a point is the difference between its squared distance - // from the center and the squared radius; Ideally, the return value - // is 0. subopt is set to the absolute value of the most negative - // coefficient in the affine combination of the support points that - // yields the center. Ideally, this is a convex combination, and there - // is no negative coefficient in which case subopt is set to 0. - NT relative_error (NT& subopt) const; - - // POST: return true if the relative error is at most tol, and the - // suboptimality is 0; the default tolerance is 10 times the - // coordinate type's machine epsilon - bool is_valid (NT tol = NT(10) * std::numeric_limits::epsilon()) const; - - // POST: returns the time in seconds taken by the constructor call for - // computing the smallest enclosing ball - double get_time() const; - - // POST: deletes dynamically allocated arrays - ~Miniball(); - - private: - void mtf_mb (Sit n); - void mtf_move_to_front (Sit j); - void pivot_mb (Pit n); - void pivot_move_to_front (Pit j); - NT excess (Pit pit) const; - void pop (); - bool push (Pit pit); - NT suboptimality () const; - void create_arrays(); - void delete_arrays(); - }; - - // Class Definition - // ================ - template - Miniball::Miniball (int d_, Pit begin, Pit end, - CoordAccessor ca) - : d (d_), - points_begin (begin), - points_end (end), - coord_accessor (ca), - time (clock()), - nt0 (NT(0)), - L(), - support_end (L.begin()), - fsize(0), - ssize(0), - current_c (NULL), - current_sqr_r (NT(-1)), - c (NULL), - sqr_r (NULL), - q0 (NULL), - z (NULL), - f (NULL), - v (NULL), - a (NULL) - { - assert (points_begin != points_end); - create_arrays(); - - // set initial center - for (int j=0; j - Miniball::~Miniball() - { - delete_arrays(); - } - - template - void Miniball::create_arrays() - { - c = new NT*[d+1]; - v = new NT*[d+1]; - a = new NT*[d+1]; - for (int i=0; i - void Miniball::delete_arrays() - { - delete[] f; - delete[] z; - delete[] q0; - delete[] sqr_r; - for (int i=0; i - const typename Miniball::NT* - Miniball::center () const - { - return current_c; - } - - template - typename Miniball::NT - Miniball::squared_radius () const - { - return current_sqr_r; - } - - template - int Miniball::nr_support_points () const - { - assert (ssize < d+2); - return ssize; - } - - template - typename Miniball::SupportPointIterator - Miniball::support_points_begin () const - { - return L.begin(); - } - - template - typename Miniball::SupportPointIterator - Miniball::support_points_end () const - { - return support_end; - } - - template - typename Miniball::NT - Miniball::relative_error (NT& subopt) const - { - NT e, max_e = nt0; - // compute maximum absolute excess of support points - for (SupportPointIterator it = support_points_begin(); - it != support_points_end(); ++it) { - e = excess (*it); - if (e < nt0) e = -e; - if (e > max_e) { - max_e = e; - } - } - // compute maximum excess of any point - for (Pit i = points_begin; i != points_end; ++i) - if ((e = excess (i)) > max_e) - max_e = e; - - subopt = suboptimality(); - assert (current_sqr_r > nt0 || max_e == nt0); - return (current_sqr_r == nt0 ? nt0 : max_e / current_sqr_r); - } - - template - bool Miniball::is_valid (NT tol) const - { - NT suboptimality; - return ( (relative_error (suboptimality) <= tol) && (suboptimality == 0) ); - } - - template - double Miniball::get_time() const - { - return time; - } - - template - void Miniball::mtf_mb (Sit n) - { - // Algorithm 1: mtf_mb (L_{n-1}, B), where L_{n-1} = [L.begin, n) - // B: the set of forced points, defining the current ball - // S: the superset of support points computed by the algorithm - // -------------------------------------------------------------- - // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, - // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf - - // PRE: B = S - assert (fsize == ssize); - - support_end = L.begin(); - if ((fsize) == d+1) return; - - // incremental construction - for (Sit i = L.begin(); i != n;) - { - // INV: (support_end - L.begin() == |S|-|B|) - assert (std::distance (L.begin(), support_end) == ssize - fsize); - - Sit j = i++; - if (excess(*j) > nt0) - if (push(*j)) { // B := B + p_i - mtf_mb (j); // mtf_mb (L_{i-1}, B + p_i) - pop(); // B := B - p_i - mtf_move_to_front(j); - } - } - // POST: the range [L.begin(), support_end) stores the set S\B - } - - template - void Miniball::mtf_move_to_front (Sit j) - { - if (support_end == j) - support_end++; - L.splice (L.begin(), L, j); - } - - template - void Miniball::pivot_mb (Pit n) - { - // Algorithm 2: pivot_mb (L_{n-1}), where L_{n-1} = [L.begin, n) - // -------------------------------------------------------------- - // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, - // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf - NT old_sqr_r; - const NT* c; - Pit pivot, k; - NT e, max_e, sqr_r; - Cit p; - do { - old_sqr_r = current_sqr_r; - sqr_r = current_sqr_r; - - pivot = points_begin; - max_e = nt0; - for (k = points_begin; k != n; ++k) { - p = coord_accessor(k); - e = -sqr_r; - c = current_c; - for (int j=0; j(*p++-*c++); - if (e > max_e) { - max_e = e; - pivot = k; - } - } - - if (max_e > nt0) { - // check if the pivot is already contained in the support set - if (std::find(L.begin(), support_end, pivot) == support_end) { - assert (fsize == 0); - if (push (pivot)) { - mtf_mb(support_end); - pop(); - pivot_move_to_front(pivot); - } - } - } - } while (old_sqr_r < current_sqr_r); - } - - template - void Miniball::pivot_move_to_front (Pit j) - { - L.push_front(j); - if (std::distance(L.begin(), support_end) == d+2) - support_end--; - } - - template - inline typename Miniball::NT - Miniball::excess (Pit pit) const - { - Cit p = coord_accessor(pit); - NT e = -current_sqr_r; - NT* c = current_c; - for (int k=0; k(*p++-*c++); - } - return e; - } - - template - void Miniball::pop () - { - --fsize; - } - - template - bool Miniball::push (Pit pit) - { - int i, j; - NT eps = mb_sqr(std::numeric_limits::epsilon()); - - Cit cit = coord_accessor(pit); - Cit p = cit; - - if (fsize==0) { - for (i=0; i(v[fsize][j]); - z[fsize]*=2; - - // reject push if z_fsize too small - if (z[fsize](*p++-c[fsize-1][i]); - f[fsize]=e/z[fsize]; - - for (i=0; i - typename Miniball::NT - Miniball::suboptimality () const - { - NT* l = new NT[d+1]; - NT min_l = nt0; - l[0] = NT(1); - for (int i=ssize-1; i>0; --i) { - l[i] = f[i]; - for (int k=ssize-1; k>i; --k) - l[i]-=a[k][i]*l[k]; - if (l[i] < min_l) min_l = l[i]; - l[0] -= l[i]; - } - if (l[0] < min_l) min_l = l[0]; - delete[] l; - if (min_l < nt0) - return -min_l; - return nt0; - } - -} // end Namespace Miniball - -#endif // MINIBALL_HPP_ diff --git a/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT b/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT new file mode 100644 index 00000000..dbe4c553 --- /dev/null +++ b/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT @@ -0,0 +1,4 @@ +The miniball software is available under the GNU General Public License (GPLv3 - https://www.gnu.org/copyleft/gpl.html). +If your intended use is not compliant with this license, please buy a commercial license (EUR 500 - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/license.html). +You need a license if the software that you develop using Miniball V3.0 is not open source. + diff --git a/src/Cech_complex/include/gudhi/Miniball.README b/src/Cech_complex/include/gudhi/Miniball.README new file mode 100644 index 00000000..86a96f08 --- /dev/null +++ b/src/Cech_complex/include/gudhi/Miniball.README @@ -0,0 +1,23 @@ +https://people.inf.ethz.ch/gaertner/subdir/software/miniball.html + +Smallest Enclosing Balls of Points - Fast and Robust in C++. +(high-quality software for smallest enclosing balls of balls is available in the computational geometry algorithms library CGAL) + + +This is the miniball software (V3.0) for computing smallest enclosing balls of points in arbitrary dimensions. It consists of a C++ header file Miniball.hpp (around 500 lines of code) and two example programs miniball_example.cpp and miniball_example_containers.cpp that demonstrate the usage. The first example stores the coordinates of the input points in a two-dimensional array, the second example uses a list of vectors to show how generic containers can be used. + +Credits: Aditya Gupta and Alexandros Konstantinakis-Karmis have significantly contributed to this version of the software. + +Changes - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/changes.txt - from previous versions. + +The theory - https://people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf - behind the miniball software (Proc. 7th Annual European Symposium on Algorithms (ESA), Lecture Notes in Computer Science 1643, Springer-Verlag, pp.325-338, 1999). + +Main Features: + + Very fast in low dimensions. 1 million points in 5-space are processed within 0.05 seconds on any recent machine. + + High numerical stability. Almost all input degeneracies (cospherical points, multiple points, points very close together) are routinely handled. + + Easily integrates into your code. You can freely choose the coordinate type of your points and the container to store the points. If you still need to adapt the code, the header is small and readable and contains documentation for all major methods. + + diff --git a/src/Cech_complex/include/gudhi/Miniball.hpp b/src/Cech_complex/include/gudhi/Miniball.hpp new file mode 100644 index 00000000..ce6cbb5b --- /dev/null +++ b/src/Cech_complex/include/gudhi/Miniball.hpp @@ -0,0 +1,523 @@ +// Copright (C) 1999-2013, Bernd Gaertner +// $Rev: 3581 $ +// +// 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 . +// +// Contact: +// -------- +// Bernd Gaertner +// Institute of Theoretical Computer Science +// ETH Zuerich +// CAB G31.1 +// CH-8092 Zuerich, Switzerland +// http://www.inf.ethz.ch/personal/gaertner + +#ifndef MINIBALL_HPP_ +#define MINIBALL_HPP_ + +#include +#include +#include +#include +#include + +namespace Gudhi { + +namespace Miniball { + + // Global Functions + // ================ + template + inline NT mb_sqr (NT r) {return r*r;} + + // Functors + // ======== + + // functor to map a point iterator to the corresponding coordinate iterator; + // generic version for points whose coordinate containers have begin() + template < typename Pit_, typename Cit_ > + struct CoordAccessor { + typedef Pit_ Pit; + typedef Cit_ Cit; + inline Cit operator() (Pit it) const { return (*it).begin(); } + }; + + // partial specialization for points whose coordinate containers are arrays + template < typename Pit_, typename Cit_ > + struct CoordAccessor { + typedef Pit_ Pit; + typedef Cit_* Cit; + inline Cit operator() (Pit it) const { return *it; } + }; + + // Class Declaration + // ================= + + template + class Miniball { + private: + // types + // The iterator type to go through the input points + typedef typename CoordAccessor::Pit Pit; + // The iterator type to go through the coordinates of a single point. + typedef typename CoordAccessor::Cit Cit; + // The coordinate type + typedef typename std::iterator_traits::value_type NT; + // The iterator to go through the support points + typedef typename std::list::iterator Sit; + + // data members... + const int d; // dimension + Pit points_begin; + Pit points_end; + CoordAccessor coord_accessor; + double time; + const NT nt0; // NT(0) + + //...for the algorithms + std::list L; + Sit support_end; + int fsize; // number of forced points + int ssize; // number of support points + + // ...for the ball updates + NT* current_c; + NT current_sqr_r; + NT** c; + NT* sqr_r; + + // helper arrays + NT* q0; + NT* z; + NT* f; + NT** v; + NT** a; + + public: + // The iterator type to go through the support points + typedef typename std::list::const_iterator SupportPointIterator; + + // PRE: [begin, end) is a nonempty range + // POST: computes the smallest enclosing ball of the points in the range + // [begin, end); the functor a maps a point iterator to an iterator + // through the d coordinates of the point + Miniball (int d_, Pit begin, Pit end, CoordAccessor ca = CoordAccessor()); + + // POST: returns a pointer to the first element of an array that holds + // the d coordinates of the center of the computed ball + const NT* center () const; + + // POST: returns the squared radius of the computed ball + NT squared_radius () const; + + // POST: returns the number of support points of the computed ball; + // the support points form a minimal set with the same smallest + // enclosing ball as the input set; in particular, the support + // points are on the boundary of the computed ball, and their + // number is at most d+1 + int nr_support_points () const; + + // POST: returns an iterator to the first support point + SupportPointIterator support_points_begin () const; + + // POST: returns a past-the-end iterator for the range of support points + SupportPointIterator support_points_end () const; + + // POST: returns the maximum excess of any input point w.r.t. the computed + // ball, divided by the squared radius of the computed ball. The + // excess of a point is the difference between its squared distance + // from the center and the squared radius; Ideally, the return value + // is 0. subopt is set to the absolute value of the most negative + // coefficient in the affine combination of the support points that + // yields the center. Ideally, this is a convex combination, and there + // is no negative coefficient in which case subopt is set to 0. + NT relative_error (NT& subopt) const; + + // POST: return true if the relative error is at most tol, and the + // suboptimality is 0; the default tolerance is 10 times the + // coordinate type's machine epsilon + bool is_valid (NT tol = NT(10) * std::numeric_limits::epsilon()) const; + + // POST: returns the time in seconds taken by the constructor call for + // computing the smallest enclosing ball + double get_time() const; + + // POST: deletes dynamically allocated arrays + ~Miniball(); + + private: + void mtf_mb (Sit n); + void mtf_move_to_front (Sit j); + void pivot_mb (Pit n); + void pivot_move_to_front (Pit j); + NT excess (Pit pit) const; + void pop (); + bool push (Pit pit); + NT suboptimality () const; + void create_arrays(); + void delete_arrays(); + }; + + // Class Definition + // ================ + template + Miniball::Miniball (int d_, Pit begin, Pit end, + CoordAccessor ca) + : d (d_), + points_begin (begin), + points_end (end), + coord_accessor (ca), + time (clock()), + nt0 (NT(0)), + L(), + support_end (L.begin()), + fsize(0), + ssize(0), + current_c (NULL), + current_sqr_r (NT(-1)), + c (NULL), + sqr_r (NULL), + q0 (NULL), + z (NULL), + f (NULL), + v (NULL), + a (NULL) + { + assert (points_begin != points_end); + create_arrays(); + + // set initial center + for (int j=0; j + Miniball::~Miniball() + { + delete_arrays(); + } + + template + void Miniball::create_arrays() + { + c = new NT*[d+1]; + v = new NT*[d+1]; + a = new NT*[d+1]; + for (int i=0; i + void Miniball::delete_arrays() + { + delete[] f; + delete[] z; + delete[] q0; + delete[] sqr_r; + for (int i=0; i + const typename Miniball::NT* + Miniball::center () const + { + return current_c; + } + + template + typename Miniball::NT + Miniball::squared_radius () const + { + return current_sqr_r; + } + + template + int Miniball::nr_support_points () const + { + assert (ssize < d+2); + return ssize; + } + + template + typename Miniball::SupportPointIterator + Miniball::support_points_begin () const + { + return L.begin(); + } + + template + typename Miniball::SupportPointIterator + Miniball::support_points_end () const + { + return support_end; + } + + template + typename Miniball::NT + Miniball::relative_error (NT& subopt) const + { + NT e, max_e = nt0; + // compute maximum absolute excess of support points + for (SupportPointIterator it = support_points_begin(); + it != support_points_end(); ++it) { + e = excess (*it); + if (e < nt0) e = -e; + if (e > max_e) { + max_e = e; + } + } + // compute maximum excess of any point + for (Pit i = points_begin; i != points_end; ++i) + if ((e = excess (i)) > max_e) + max_e = e; + + subopt = suboptimality(); + assert (current_sqr_r > nt0 || max_e == nt0); + return (current_sqr_r == nt0 ? nt0 : max_e / current_sqr_r); + } + + template + bool Miniball::is_valid (NT tol) const + { + NT suboptimality; + return ( (relative_error (suboptimality) <= tol) && (suboptimality == 0) ); + } + + template + double Miniball::get_time() const + { + return time; + } + + template + void Miniball::mtf_mb (Sit n) + { + // Algorithm 1: mtf_mb (L_{n-1}, B), where L_{n-1} = [L.begin, n) + // B: the set of forced points, defining the current ball + // S: the superset of support points computed by the algorithm + // -------------------------------------------------------------- + // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, + // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf + + // PRE: B = S + assert (fsize == ssize); + + support_end = L.begin(); + if ((fsize) == d+1) return; + + // incremental construction + for (Sit i = L.begin(); i != n;) + { + // INV: (support_end - L.begin() == |S|-|B|) + assert (std::distance (L.begin(), support_end) == ssize - fsize); + + Sit j = i++; + if (excess(*j) > nt0) + if (push(*j)) { // B := B + p_i + mtf_mb (j); // mtf_mb (L_{i-1}, B + p_i) + pop(); // B := B - p_i + mtf_move_to_front(j); + } + } + // POST: the range [L.begin(), support_end) stores the set S\B + } + + template + void Miniball::mtf_move_to_front (Sit j) + { + if (support_end == j) + support_end++; + L.splice (L.begin(), L, j); + } + + template + void Miniball::pivot_mb (Pit n) + { + // Algorithm 2: pivot_mb (L_{n-1}), where L_{n-1} = [L.begin, n) + // -------------------------------------------------------------- + // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, + // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf + NT old_sqr_r; + const NT* c; + Pit pivot, k; + NT e, max_e, sqr_r; + Cit p; + do { + old_sqr_r = current_sqr_r; + sqr_r = current_sqr_r; + + pivot = points_begin; + max_e = nt0; + for (k = points_begin; k != n; ++k) { + p = coord_accessor(k); + e = -sqr_r; + c = current_c; + for (int j=0; j(*p++-*c++); + if (e > max_e) { + max_e = e; + pivot = k; + } + } + + if (max_e > nt0) { + // check if the pivot is already contained in the support set + if (std::find(L.begin(), support_end, pivot) == support_end) { + assert (fsize == 0); + if (push (pivot)) { + mtf_mb(support_end); + pop(); + pivot_move_to_front(pivot); + } + } + } + } while (old_sqr_r < current_sqr_r); + } + + template + void Miniball::pivot_move_to_front (Pit j) + { + L.push_front(j); + if (std::distance(L.begin(), support_end) == d+2) + support_end--; + } + + template + inline typename Miniball::NT + Miniball::excess (Pit pit) const + { + Cit p = coord_accessor(pit); + NT e = -current_sqr_r; + NT* c = current_c; + for (int k=0; k(*p++-*c++); + } + return e; + } + + template + void Miniball::pop () + { + --fsize; + } + + template + bool Miniball::push (Pit pit) + { + int i, j; + NT eps = mb_sqr(std::numeric_limits::epsilon()); + + Cit cit = coord_accessor(pit); + Cit p = cit; + + if (fsize==0) { + for (i=0; i(v[fsize][j]); + z[fsize]*=2; + + // reject push if z_fsize too small + if (z[fsize](*p++-c[fsize-1][i]); + f[fsize]=e/z[fsize]; + + for (i=0; i + typename Miniball::NT + Miniball::suboptimality () const + { + NT* l = new NT[d+1]; + NT min_l = nt0; + l[0] = NT(1); + for (int i=ssize-1; i>0; --i) { + l[i] = f[i]; + for (int k=ssize-1; k>i; --k) + l[i]-=a[k][i]*l[k]; + if (l[i] < min_l) min_l = l[i]; + l[0] -= l[i]; + } + if (l[0] < min_l) min_l = l[0]; + delete[] l; + if (min_l < nt0) + return -min_l; + return nt0; + } +} // namespace Miniball + +} // namespace Gudhi + +#endif // MINIBALL_HPP_ diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 8cbfe431..4aa85057 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -36,8 +36,7 @@ #include #include #include - -#include +#include // Type definitions using Simplex_tree = Gudhi::Simplex_tree<>; @@ -49,7 +48,7 @@ using Cech_complex = Gudhi::cech_complex::Cech_complex>; +using Min_sphere = Gudhi::Miniball::Miniball>; BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { // ---------------------------------------------------------------------------- diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index 7749b98e..93956a69 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -25,7 +25,7 @@ #include -#include +#include #include -- cgit v1.2.3 From 1f3716292673a56413d3501b4b98b54416d193ed Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Thu, 22 Mar 2018 10:10:08 +0000 Subject: code review : Rename Radius_distance Minimal_enclosing_ball_radius git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/cechcomplex_vincent@3304 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 204aa7a7098773b2d290e35a12a1c92449743d7d --- .../benchmark/cech_complex_benchmark.cpp | 32 ++++++++++++---------- .../example/cech_complex_step_by_step.cpp | 4 +-- src/Cech_complex/include/gudhi/Cech_complex.h | 11 ++++---- .../include/gudhi/Cech_complex_blocker.h | 4 +-- src/Cech_complex/test/test_cech_complex.cpp | 4 +-- src/common/include/gudhi/distance_functions.h | 8 +++--- 6 files changed, 34 insertions(+), 29 deletions(-) (limited to 'src/Cech_complex/test/test_cech_complex.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 6ba52419..2fa255ed 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -46,7 +46,7 @@ using Rips_complex = Gudhi::rips_complex::Rips_complex; using Cech_complex = Gudhi::cech_complex::Cech_complex; -class Radius_distance { +class Minimal_enclosing_ball_radius { public: // boost::range_value is not SFINAE-friendly so we cannot use it in the return type template< typename Point > @@ -80,24 +80,26 @@ int main(int argc, char * argv[]) { Gudhi::Clock euclidean_clock("Gudhi::Euclidean_distance"); // Compute the proximity graph of the points Proximity_graph euclidean_prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), - threshold, - Gudhi::Euclidean_distance()); + threshold, + Gudhi::Euclidean_distance()); std::cout << euclidean_clock << std::endl; - Gudhi::Clock radius_clock("Radius_distance"); + Gudhi::Clock miniball_clock("Minimal_enclosing_ball_radius"); // Compute the proximity graph of the points - Proximity_graph radius_prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), - threshold, - Radius_distance()); - std::cout << radius_clock << std::endl; + Proximity_graph miniball_prox_graph = + Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), + threshold, + Minimal_enclosing_ball_radius()); + std::cout << miniball_clock << std::endl; - Gudhi::Clock common_radius_clock("Gudhi::Radius_distance()"); + Gudhi::Clock common_miniball_clock("Gudhi::Minimal_enclosing_ball_radius()"); // Compute the proximity graph of the points - Proximity_graph sq_radius_prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), - threshold, - Gudhi::Radius_distance()); - std::cout << common_radius_clock << std::endl; + Proximity_graph common_miniball_prox_graph = + Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), + threshold, + Gudhi::Minimal_enclosing_ball_radius()); + std::cout << common_miniball_clock << std::endl; boost::filesystem::path full_path(boost::filesystem::current_path()); @@ -121,7 +123,9 @@ int main(int argc, char * argv[]) { 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()); + Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), + radius, + Gudhi::Minimal_enclosing_ball_radius()); Simplex_tree rips_stree; rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); // ------------------------------------------ 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 0e7c4bbd..760b53dc 100644 --- a/src/Cech_complex/example/cech_complex_step_by_step.cpp +++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp @@ -65,7 +65,7 @@ class Cech_blocker { std::cout << "#(" << vertex << ")#"; #endif // DEBUG_TRACES } - Filtration_value radius = Gudhi::Radius_distance()(points); + Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); #ifdef DEBUG_TRACES std::cout << "radius = " << radius << " - " << (radius > max_radius_) << std::endl; #endif // DEBUG_TRACES @@ -105,7 +105,7 @@ int main(int argc, char * argv[]) { // Compute the proximity graph of the points Proximity_graph prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), max_radius, - Gudhi::Radius_distance()); + Gudhi::Minimal_enclosing_ball_radius()); // Construct the Rips complex in a Simplex Tree Simplex_tree st; diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 612c73c3..8b1a9221 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -23,7 +23,7 @@ #ifndef CECH_COMPLEX_H_ #define CECH_COMPLEX_H_ -#include // for Gudhi::Radius_distance +#include // for Gudhi::Minimal_enclosing_ball_radius #include // for Gudhi::Proximity_graph #include // for GUDHI_CHECK #include // for Gudhi::cech_complex::Cech_blocker @@ -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 max_radius. Edge length is computed from `Gudhi::Radius_distance` distance function. + * to a given max_radius. Edge length is computed from `Gudhi::Minimal_enclosing_ball_radius` distance function. * * \tparam SimplicialComplexForProximityGraph furnishes `Vertex_handle` and `Filtration_value` type definition required * by `Gudhi::Proximity_graph`. @@ -71,9 +71,10 @@ class Cech_complex { Cech_complex(const ForwardPointRange& points, Filtration_value max_radius) : max_radius_(max_radius), point_cloud_(points) { - cech_skeleton_graph_ = Gudhi::compute_proximity_graph(point_cloud_, - max_radius_, - Gudhi::Radius_distance()); + cech_skeleton_graph_ = + Gudhi::compute_proximity_graph(point_cloud_, + max_radius_, + Gudhi::Minimal_enclosing_ball_radius()); } /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 5ba17c51..c082815d 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -24,7 +24,7 @@ #define CECH_COMPLEX_BLOCKER_H_ #include // Cech_blocker is using a pointer on Gudhi::cech_complex::Cech_complex -#include // for Gudhi::Radius_distance +#include // for Gudhi::Minimal_enclosing_ball_radius #include #include @@ -75,7 +75,7 @@ class Cech_blocker { std::cout << "#(" << vertex << ")#"; #endif // DEBUG_TRACES } - Filtration_value radius = Gudhi::Radius_distance()(points); + Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); #ifdef DEBUG_TRACES if (radius > cc_ptr_->max_radius()) std::cout << "radius > max_radius => expansion is blocked\n"; diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 4aa85057..8658729b 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -111,10 +111,10 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { std::cout << vertex << ","; vp.push_back(points.at(vertex)); } - std::cout << ") - distance =" << Gudhi::Radius_distance()(vp.at(0), vp.at(1)) << + std::cout << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius()(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::Radius_distance()(vp.at(0), vp.at(1))); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1))); } } diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index 93956a69..429a5fa7 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -70,9 +70,9 @@ 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 { +/** @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: // boost::range_value is not SFINAE-friendly so we cannot use it in the return type template< typename Point > @@ -92,7 +92,7 @@ class Radius_distance { Min_sphere ms(point_cloud.begin()->end() - point_cloud.begin()->begin(), point_cloud.begin(),point_cloud.end()); #ifdef DEBUG_TRACES - std::cout << "Radius_distance = " << std::sqrt(ms.squared_radius()) << " | nb points = " + std::cout << "Minimal_enclosing_ball_radius = " << std::sqrt(ms.squared_radius()) << " | nb points = " << point_cloud.end() - point_cloud.begin() << " | dimension = " << point_cloud.begin()->end() - point_cloud.begin()->begin() << std::endl; #endif // DEBUG_TRACES -- cgit v1.2.3 From 51ce9b513116f5fed2b4dc109f0b52595a2cd538 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Tue, 10 Apr 2018 14:53:20 +0000 Subject: Code review : Cech_blocker was hardcoding double replace point_iterator function in a get_point const function that returns a Point InputPointRange description Cech blocker is now templated with Cech complex, which is no more included. Deep copy of the point cloud on Cech complex ctor git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/cechcomplex_vincent@3365 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 93e3cf3fe61290b88a0714b9e55ad80e01a34f36 --- .../example/cech_complex_example_from_points.cpp | 1 - src/Cech_complex/include/gudhi/Cech_complex.h | 39 ++++++++++++++-------- .../include/gudhi/Cech_complex_blocker.h | 15 +++------ src/Cech_complex/test/test_cech_complex.cpp | 26 +++++---------- 4 files changed, 38 insertions(+), 43 deletions(-) (limited to 'src/Cech_complex/test/test_cech_complex.cpp') 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 3b889d56..2b36e33c 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,5 @@ #include #include -#include #include #include diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 52f03d6b..abad0c21 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -54,11 +54,20 @@ namespace cech_complex { template class Cech_complex { private: + // Required by compute_proximity_graph using Vertex_handle = typename SimplicialComplexForProximityGraph::Vertex_handle; using Filtration_value = typename SimplicialComplexForProximityGraph::Filtration_value; using Proximity_graph = Gudhi::Proximity_graph; - using Point_iterator = typename boost::range_const_iterator::type; - using Point= typename std::iterator_traits::value_type; + + // Retrieve Coordinate type from InputPointRange + using Point_from_range_iterator = typename boost::range_const_iterator::type; + using Point_from_range = typename std::iterator_traits::value_type; + using Coordinate_iterator = typename boost::range_const_iterator::type; + using Coordinate= typename std::iterator_traits::value_type; + + public: + // Point and Point_cloud type definition + using Point = std::vector; using Point_cloud = std::vector; public: @@ -67,13 +76,20 @@ class Cech_complex { * @param[in] points Range of points. * @param[in] max_radius Maximal radius value. * - * \tparam InputPointRange 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. + * \tparam InputPointRange must be a range of Point. Point must be a range of copyable Cartesian coordinates. * */ Cech_complex(const InputPointRange& points, Filtration_value max_radius) - : max_radius_(max_radius), - point_cloud_(std::begin(points), std::end(points)) { + : max_radius_(max_radius) { + // Point cloud deep copy + auto points_begin_itr = std::begin(points); + auto points_end_itr = std::end(points); + + point_cloud_.reserve(points_end_itr - points_begin_itr); + for (auto point_itr = points_begin_itr; point_itr < points_end_itr; point_itr++) { + point_cloud_.push_back(Point(std::begin(*point_itr), std::end(*point_itr))); + } + cech_skeleton_graph_ = Gudhi::compute_proximity_graph(point_cloud_, max_radius_, @@ -97,7 +113,7 @@ class Cech_complex { complex.insert_graph(cech_skeleton_graph_); // expand the graph until dimension dim_max complex.expansion_with_blockers(dim_max, - Cech_blocker(&complex, this)); + Cech_blocker(&complex, this)); } /** @return max_radius value given at construction. */ @@ -106,13 +122,10 @@ class Cech_complex { } /** @param[in] vertex Point position in the range. - * @return A const iterator on the point. - * @exception std::out_of_range In debug mode, if point position in the range is out. + * @return The point. */ - typename InputPointRange::const_iterator point_iterator(std::size_t vertex) const { - GUDHI_CHECK((std::begin(point_cloud_) + vertex) < std::end(point_cloud_), - std::out_of_range("Cech_complex::point - simplicial complex is not empty")); - return (std::begin(point_cloud_) + vertex); + const Point& get_point(Vertex_handle vertex) const{ + return point_cloud_[vertex]; } private: diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index ab56c10d..2ecef9cf 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -23,7 +23,6 @@ #ifndef CECH_COMPLEX_BLOCKER_H_ #define CECH_COMPLEX_BLOCKER_H_ -#include // Cech_blocker is using a pointer on Gudhi::cech_complex::Cech_complex #include // for Gudhi::Minimal_enclosing_ball_radius #include @@ -34,10 +33,6 @@ namespace Gudhi { namespace cech_complex { -// Just declaring Cech_complex class because used and not yet defined. -template -class Cech_complex; - /** \internal * \class Cech_blocker * \brief Cech complex blocker. @@ -52,14 +47,13 @@ class Cech_complex; * * \tparam InputPointRange is required by the pointer on Chech_complex for type definition. */ -template +template class Cech_blocker { private: - using Point = std::vector; - using Point_cloud = std::vector; + using Point_cloud = typename Cech_complex::Point_cloud; + using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; using Filtration_value = typename SimplicialComplexForCech::Filtration_value; - using Cech_complex = Gudhi::cech_complex::Cech_complex; public: /** \internal \brief Cech complex blocker operator() - the oracle - assigns the filtration value from the simplex @@ -69,8 +63,7 @@ class Cech_blocker { bool operator()(Simplex_handle sh) { Point_cloud points; for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { - points.push_back(Point(cc_ptr_->point_iterator(vertex)->begin(), - cc_ptr_->point_iterator(vertex)->end())); + points.push_back(cc_ptr_->get_point(vertex)); #ifdef DEBUG_TRACES std::cout << "#(" << vertex << ")#"; #endif // DEBUG_TRACES diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 8658729b..5ca25db4 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -78,11 +78,8 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { GUDHI_TEST_FLOAT_EQUALITY_CHECK(cech_complex_for_doc.max_radius(), max_radius); std::size_t i = 0; for (; i < points.size(); i++) { - BOOST_CHECK(points[i] == *(cech_complex_for_doc.point_iterator(i))); + BOOST_CHECK(points[i] == cech_complex_for_doc.get_point(i)); } -#ifdef GUDHI_DEBUG - BOOST_CHECK_THROW (cech_complex_for_doc.point_iterator(i+1), std::out_of_range); -#endif // GUDHI_DEBUG const int DIMENSION_1 = 1; Simplex_tree st; @@ -137,8 +134,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { Point_cloud points012; for (std::size_t vertex = 0; vertex <= 2; vertex++) { - points012.push_back(Point(cech_complex_for_doc.point_iterator(vertex)->begin(), - cech_complex_for_doc.point_iterator(vertex)->end())); + 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()); @@ -149,12 +145,9 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(ms012.squared_radius())); Point_cloud points1410; - points1410.push_back(Point(cech_complex_for_doc.point_iterator(1)->begin(), - cech_complex_for_doc.point_iterator(1)->end())); - points1410.push_back(Point(cech_complex_for_doc.point_iterator(4)->begin(), - cech_complex_for_doc.point_iterator(4)->end())); - points1410.push_back(Point(cech_complex_for_doc.point_iterator(10)->begin(), - cech_complex_for_doc.point_iterator(10)->end())); + 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})); @@ -163,12 +156,9 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, std::sqrt(ms1410.squared_radius())); Point_cloud points469; - points469.push_back(Point(cech_complex_for_doc.point_iterator(4)->begin(), - cech_complex_for_doc.point_iterator(4)->end())); - points469.push_back(Point(cech_complex_for_doc.point_iterator(6)->begin(), - cech_complex_for_doc.point_iterator(6)->end())); - points469.push_back(Point(cech_complex_for_doc.point_iterator(9)->begin(), - cech_complex_for_doc.point_iterator(9)->end())); + 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})); -- cgit v1.2.3 From cfe5c6b05435cb7d8cbd1d615e0c402f4a5e1674 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Tue, 29 May 2018 19:49:30 +0000 Subject: Clang format files git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/cechcomplex_vincent@3488 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 7d1da847f4c39a191afd5d8e56fe34ee79ade495 --- .../benchmark/cech_complex_benchmark.cpp | 58 ++++++++------------ .../example/cech_complex_example_from_points.cpp | 33 ++++++------ .../example/cech_complex_step_by_step.cpp | 51 +++++++----------- src/Cech_complex/include/gudhi/Cech_complex.h | 30 ++++------- .../include/gudhi/Cech_complex_blocker.h | 9 ++-- src/Cech_complex/test/test_cech_complex.cpp | 63 +++++++++++----------- 6 files changed, 104 insertions(+), 140 deletions(-) (limited to 'src/Cech_complex/test/test_cech_complex.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 2fa255ed..86314930 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -29,12 +29,11 @@ #include #include -#include "boost/filesystem.hpp" // includes all needed Boost.Filesystem declarations +#include "boost/filesystem.hpp" // includes all needed Boost.Filesystem declarations #include #include - // Types definition using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; @@ -45,32 +44,31 @@ using Proximity_graph = Gudhi::Proximity_graph; using Rips_complex = Gudhi::rips_complex::Rips_complex; using Cech_complex = Gudhi::cech_complex::Cech_complex; - class Minimal_enclosing_ball_radius { public: // boost::range_value is not SFINAE-friendly so we cannot use it in the return type - template< typename Point > - typename std::iterator_traits::type>::value_type - operator()(const Point& p1, const Point& p2) const { + template + typename std::iterator_traits::type>::value_type operator()( + const Point& p1, const Point& p2) const { // Type def using Point_cloud = std::vector; using Point_iterator = typename Point_cloud::const_iterator; using Coordinate_iterator = typename Point::const_iterator; - using Min_sphere = typename Gudhi::Miniball::Miniball>; + using Min_sphere = + typename Gudhi::Miniball::Miniball>; 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()); + 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[]) { +int main(int argc, char* argv[]) { std::string off_file_points = "tore3D_1307.off"; Filtration_value threshold = 1e20; @@ -79,42 +77,32 @@ int main(int argc, char * argv[]) { Gudhi::Clock euclidean_clock("Gudhi::Euclidean_distance"); // Compute the proximity graph of the points - Proximity_graph euclidean_prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), - threshold, - Gudhi::Euclidean_distance()); + Proximity_graph euclidean_prox_graph = Gudhi::compute_proximity_graph( + off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); std::cout << euclidean_clock << std::endl; Gudhi::Clock miniball_clock("Minimal_enclosing_ball_radius"); // Compute the proximity graph of the points - Proximity_graph miniball_prox_graph = - Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), - threshold, - Minimal_enclosing_ball_radius()); + Proximity_graph miniball_prox_graph = Gudhi::compute_proximity_graph( + off_reader.get_point_cloud(), threshold, Minimal_enclosing_ball_radius()); std::cout << miniball_clock << std::endl; Gudhi::Clock common_miniball_clock("Gudhi::Minimal_enclosing_ball_radius()"); // Compute the proximity graph of the points - Proximity_graph common_miniball_prox_graph = - Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), - threshold, - Gudhi::Minimal_enclosing_ball_radius()); + Proximity_graph common_miniball_prox_graph = Gudhi::compute_proximity_graph( + off_reader.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius()); std::cout << common_miniball_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 + 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()); Point p0 = off_reader.get_point_cloud()[0]; @@ -123,8 +111,7 @@ int main(int argc, char * argv[]) { 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, + Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), radius, Gudhi::Minimal_enclosing_ball_radius()); Simplex_tree rips_stree; rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); @@ -153,6 +140,5 @@ int main(int argc, char * argv[]) { } } - return 0; } 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 2b36e33c..3cc5a4df 100644 --- a/src/Cech_complex/example/cech_complex_example_from_points.cpp +++ b/src/Cech_complex/example/cech_complex_example_from_points.cpp @@ -14,17 +14,17 @@ int main() { using Cech_complex = Gudhi::cech_complex::Cech_complex; 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 + 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 // ---------------------------------------------------------------------------- // Init of a Cech complex from points @@ -37,18 +37,17 @@ int main() { // ---------------------------------------------------------------------------- // Display information about the one skeleton Cech complex // ---------------------------------------------------------------------------- - std::cout << "Cech complex is of dimension " << stree.dimension() << - " - " << stree.num_simplices() << " simplices - " << - stree.num_vertices() << " vertices." << std::endl; + std::cout << "Cech complex is of dimension " << stree.dimension() << " - " << stree.num_simplices() << " simplices - " + << stree.num_vertices() << " vertices." << std::endl; - std::cout << "Iterator on Cech complex simplices in the filtration order, with [filtration value]:" << - std::endl; + std::cout << "Iterator on Cech complex simplices in the filtration order, with [filtration value]:" << std::endl; for (auto f_simplex : stree.filtration_simplex_range()) { std::cout << " ( "; for (auto vertex : stree.simplex_vertex_range(f_simplex)) { std::cout << vertex << " "; } - std::cout << ") -> " << "[" << stree.filtration(f_simplex) << "] "; + std::cout << ") -> " + << "[" << stree.filtration(f_simplex) << "] "; std::cout << std::endl; } return 0; 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 760b53dc..d2dc8b65 100644 --- a/src/Cech_complex/example/cech_complex_step_by_step.cpp +++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp @@ -31,7 +31,7 @@ #include #include -#include // infinity +#include // infinity #include // for pair #include @@ -55,7 +55,8 @@ class Cech_blocker { using Point_cloud = std::vector; using Point_iterator = Point_cloud::const_iterator; using Coordinate_iterator = Point::const_iterator; - using Min_sphere = Gudhi::Miniball::Miniball >; + using Min_sphere = Gudhi::Miniball::Miniball>; + public: bool operator()(Simplex_handle sh) { std::vector points; @@ -73,11 +74,10 @@ class Cech_blocker { return (radius > max_radius_); } Cech_blocker(Simplex_tree& simplex_tree, Filtration_value max_radius, const std::vector& point_cloud) - : simplex_tree_(simplex_tree), - max_radius_(max_radius), - point_cloud_(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_; @@ -85,14 +85,9 @@ class Cech_blocker { int dimension_; }; +void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& max_radius, int& dim_max); -void program_options(int argc, char * argv[] - , std::string & off_file_points - , Filtration_value & max_radius - , int & dim_max); - - -int main(int argc, char * argv[]) { +int main(int argc, char* argv[]) { std::string off_file_points; Filtration_value max_radius; int dim_max; @@ -103,8 +98,7 @@ int main(int argc, char * argv[]) { Points_off_reader off_reader(off_file_points); // Compute the proximity graph of the points - Proximity_graph prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), - max_radius, + Proximity_graph prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), max_radius, Gudhi::Minimal_enclosing_ball_radius()); // Construct the Rips complex in a Simplex Tree @@ -125,7 +119,8 @@ int main(int argc, char * argv[]) { 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()) { - std::cout << " " << "[" << st.filtration(f_simplex) << "] "; + std::cout << " " + << "[" << st.filtration(f_simplex) << "] "; for (auto vertex : st.simplex_vertex_range(f_simplex)) { std::cout << static_cast(vertex) << " "; } @@ -136,24 +131,19 @@ int main(int argc, char * argv[]) { return 0; } -void program_options(int argc, char * argv[] - , std::string & off_file_points - , Filtration_value & max_radius - , int & dim_max) { +void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& max_radius, int& dim_max) { namespace po = boost::program_options; po::options_description hidden("Hidden options"); - hidden.add_options() - ("input-file", po::value(&off_file_points), - "Name of an OFF file containing a point set.\n"); + hidden.add_options()("input-file", po::value(&off_file_points), + "Name of an OFF file containing a point set.\n"); po::options_description visible("Allowed options", 100); - visible.add_options() - ("help,h", "produce help message") - ("max-radius,r", - po::value(&max_radius)->default_value(std::numeric_limits::infinity()), - "Maximal length of an edge for the Rips complex construction.") - ("cpx-dimension,d", po::value(&dim_max)->default_value(1), - "Maximal dimension of the Rips complex we want to compute."); + visible.add_options()("help,h", "produce help message")( + "max-radius,r", + po::value(&max_radius)->default_value(std::numeric_limits::infinity()), + "Maximal length of an edge for the Rips complex construction.")( + "cpx-dimension,d", po::value(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute."); po::positional_options_description pos; pos.add("input-file", 1); @@ -162,8 +152,7 @@ void program_options(int argc, char * argv[] all.add(visible).add(hidden); po::variables_map vm; - po::store(po::command_line_parser(argc, argv). - options(all).positional(pos).run(), vm); + po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm); po::notify(vm); if (vm.count("help") || !vm.count("input-file")) { diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index def46c79..f4fb4288 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -38,9 +38,9 @@ namespace cech_complex { /** * \class Cech_complex * \brief Cech complex data structure. - * + * * \ingroup cech_complex - * + * * \details * The data structure is a proximity graph, containing edges when the edge length is less or equal * to a given max_radius. Edge length is computed from `Gudhi::Minimal_enclosing_ball_radius` distance function. @@ -51,7 +51,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 +template class Cech_complex { private: // Required by compute_proximity_graph @@ -63,7 +63,7 @@ class Cech_complex { using Point_from_range_iterator = typename boost::range_const_iterator::type; using Point_from_range = typename std::iterator_traits::value_type; using Coordinate_iterator = typename boost::range_const_iterator::type; - using Coordinate= typename std::iterator_traits::value_type; + using Coordinate = typename std::iterator_traits::value_type; public: // Point and Point_cloud type definition @@ -79,17 +79,13 @@ class Cech_complex { * \tparam ForwardPointRange must be a range of Point. Point must be a range of copyable Cartesian coordinates. * */ - Cech_complex(const ForwardPointRange& points, Filtration_value max_radius) - : max_radius_(max_radius) { + 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)); + for (auto&& point : points) point_cloud_.emplace_back(std::begin(point), std::end(point)); - cech_skeleton_graph_ = - Gudhi::compute_proximity_graph(point_cloud_, - max_radius_, - Gudhi::Minimal_enclosing_ball_radius()); + cech_skeleton_graph_ = Gudhi::compute_proximity_graph( + point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); } /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal @@ -100,7 +96,7 @@ class Cech_complex { * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0. * */ - template + template 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")); @@ -113,16 +109,12 @@ class Cech_complex { } /** @return max_radius value given at construction. */ - Filtration_value max_radius() const { - return max_radius_; - } + Filtration_value max_radius() const { return max_radius_; } /** @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& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; } private: Proximity_graph cech_skeleton_graph_; diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 697d2246..b0d347b1 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -70,18 +70,15 @@ class Cech_blocker { } Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); #ifdef DEBUG_TRACES - if (radius > cc_ptr_->max_radius()) - std::cout << "radius > max_radius => expansion is blocked\n"; + if (radius > cc_ptr_->max_radius()) std::cout << "radius > max_radius => expansion is blocked\n"; #endif // DEBUG_TRACES sc_ptr_->assign_filtration(sh, radius); return (radius > cc_ptr_->max_radius()); } /** \internal \brief Čech complex blocker constructor. */ - Cech_blocker(SimplicialComplexForCech* sc_ptr, Cech_complex* cc_ptr) - : sc_ptr_(sc_ptr), - cc_ptr_(cc_ptr) { - } + Cech_blocker(SimplicialComplexForCech* sc_ptr, Cech_complex* cc_ptr) : sc_ptr_(sc_ptr), cc_ptr_(cc_ptr) {} + private: SimplicialComplexForCech* sc_ptr_; Cech_complex* cc_ptr_; diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 5b35735b..9039169c 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -28,7 +28,7 @@ #include #include #include -#include // std::max +#include // std::max #include // to construct Cech_complex from a OFF file of points @@ -57,21 +57,21 @@ 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 + 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 Filtration_value max_radius = 1.0; - std::cout << "========== NUMBER OF POINTS = " << points.size() << " - Cech max_radius = " << - max_radius << "==========" << std::endl; + std::cout << "========== NUMBER OF POINTS = " << points.size() << " - Cech max_radius = " << max_radius + << "==========" << std::endl; Cech_complex cech_complex_for_doc(points, max_radius); @@ -108,24 +108,25 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { std::cout << vertex << ","; vp.push_back(points.at(vertex)); } - std::cout << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1)) << - " - filtration =" << st.filtration(f_simplex) << std::endl; + std::cout << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius()(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_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), + Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1))); } } const int DIMENSION_2 = 2; #ifdef GUDHI_DEBUG - BOOST_CHECK_THROW (cech_complex_for_doc.create_complex(st, DIMENSION_2), std::invalid_argument); + BOOST_CHECK_THROW(cech_complex_for_doc.create_complex(st, DIMENSION_2), std::invalid_argument); #endif Simplex_tree st2; cech_complex_for_doc.create_complex(st2, DIMENSION_2); std::cout << "st2.dimension()=" << st2.dimension() << std::endl; BOOST_CHECK(st2.dimension() == DIMENSION_2); - + std::cout << "st2.num_vertices()=" << st2.num_vertices() << std::endl; BOOST_CHECK(st2.num_vertices() == NUMBER_OF_VERTICES); @@ -137,7 +138,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { 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()); + 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; @@ -148,7 +149,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { 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()); + Min_sphere ms1410(dimension, points1410.begin(), points1410.end()); Simplex_tree::Filtration_value f1410 = st2.filtration(st2.find({1, 4, 10})); std::cout << "f1410= " << f1410 << " | ms1410_radius= " << std::sqrt(ms1410.squared_radius()) << std::endl; @@ -159,7 +160,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { 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()); + Min_sphere ms469(dimension, points469.begin(), points469.end()); Simplex_tree::Filtration_value f469 = st2.filtration(st2.find({4, 6, 9})); std::cout << "f469= " << f469 << " | ms469_radius= " << std::sqrt(ms469.squared_radius()) << std::endl; @@ -168,7 +169,6 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { BOOST_CHECK((st2.find({6, 7, 8}) == st2.null_simplex())); BOOST_CHECK((st2.find({3, 5, 7}) == st2.null_simplex())); - } BOOST_AUTO_TEST_CASE(Cech_complex_from_points) { @@ -176,13 +176,13 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_points) { // Init of a list of points // ---------------------------------------------------------------------------- Point_cloud points; - std::vector coords = { 0.0, 0.0, 0.0, 1.0 }; + std::vector coords = {0.0, 0.0, 0.0, 1.0}; points.push_back(Point(coords.begin(), coords.end())); - coords = { 0.0, 0.0, 1.0, 0.0 }; + coords = {0.0, 0.0, 1.0, 0.0}; points.push_back(Point(coords.begin(), coords.end())); - coords = { 0.0, 1.0, 0.0, 0.0 }; + coords = {0.0, 1.0, 0.0, 0.0}; points.push_back(Point(coords.begin(), coords.end())); - coords = { 1.0, 0.0, 0.0, 0.0 }; + coords = {1.0, 0.0, 0.0, 0.0}; points.push_back(Point(coords.begin(), coords.end())); // ---------------------------------------------------------------------------- @@ -204,7 +204,8 @@ BOOST_AUTO_TEST_CASE(Cech_complex_from_points) { for (auto vertex : st.simplex_vertex_range(f_simplex)) { std::cout << vertex << " "; } - std::cout << ") -> " << "[" << st.filtration(f_simplex) << "] "; + std::cout << ") -> " + << "[" << st.filtration(f_simplex) << "] "; std::cout << std::endl; } BOOST_CHECK(num_simplices == 15); @@ -247,8 +248,8 @@ BOOST_AUTO_TEST_CASE(Cech_create_complex_throw) { // ---------------------------------------------------------------------------- std::string off_file_name("alphacomplexdoc.off"); double max_radius = 12.0; - std::cout << "========== OFF FILE NAME = " << off_file_name << " - Cech max_radius=" << - max_radius << "==========" << std::endl; + std::cout << "========== OFF FILE NAME = " << off_file_name << " - Cech max_radius=" << max_radius + << "==========" << std::endl; Gudhi::Points_off_reader off_reader(off_file_name); Cech_complex cech_complex_from_file(off_reader.get_point_cloud(), max_radius); @@ -258,6 +259,6 @@ BOOST_AUTO_TEST_CASE(Cech_create_complex_throw) { stree.insert_simplex_and_subfaces(simplex); std::cout << "Check exception throw in debug mode" << std::endl; // throw excpt because stree is not empty - BOOST_CHECK_THROW (cech_complex_from_file.create_complex(stree, 1), std::invalid_argument); + BOOST_CHECK_THROW(cech_complex_from_file.create_complex(stree, 1), std::invalid_argument); } #endif -- cgit v1.2.3