diff options
Diffstat (limited to 'src/Persistent_cohomology/example')
4 files changed, 175 insertions, 31 deletions
diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index 1f32da17..c37b4229 100644 --- a/src/Persistent_cohomology/example/CMakeLists.txt +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -21,11 +21,11 @@ add_executable(persistence_from_file persistence_from_file.cpp) target_link_libraries(persistence_from_file ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) if (TBB_FOUND) - target_link_libraries(plain_homology ${TBB_RELEASE_LIBRARY}) - target_link_libraries(persistence_from_simple_simplex_tree ${TBB_RELEASE_LIBRARY}) - target_link_libraries(rips_persistence ${TBB_RELEASE_LIBRARY}) - target_link_libraries(rips_persistence_via_boundary_matrix ${TBB_RELEASE_LIBRARY}) - target_link_libraries(persistence_from_file ${TBB_RELEASE_LIBRARY}) + target_link_libraries(plain_homology ${TBB_LIBRARIES}) + target_link_libraries(persistence_from_simple_simplex_tree ${TBB_LIBRARIES}) + target_link_libraries(rips_persistence ${TBB_LIBRARIES}) + target_link_libraries(rips_persistence_via_boundary_matrix ${TBB_LIBRARIES}) + target_link_libraries(persistence_from_file ${TBB_LIBRARIES}) endif() add_test(plain_homology ${CMAKE_CURRENT_BINARY_DIR}/plain_homology) @@ -35,28 +35,34 @@ add_test(rips_persistence_via_boundary_matrix_3 ${CMAKE_CURRENT_BINARY_DIR}/rips add_test(persistence_from_file_3_2_0 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file ${CMAKE_SOURCE_DIR}/data/points/bunny_5000.st -p 2 -m 0) add_test(persistence_from_file_3_3_100 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file ${CMAKE_SOURCE_DIR}/data/points/bunny_5000.st -p 3 -m 100) -if(GMPXX_FOUND AND GMP_FOUND) - message("GMPXX_LIBRARIES = ${GMPXX_LIBRARIES}") +if(GMP_FOUND) message("GMP_LIBRARIES = ${GMP_LIBRARIES}") + if(GMPXX_FOUND) + message("GMPXX_LIBRARIES = ${GMPXX_LIBRARIES}") + add_executable(rips_multifield_persistence rips_multifield_persistence.cpp ) target_link_libraries(rips_multifield_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) add_executable ( performance_rips_persistence performance_rips_persistence.cpp ) target_link_libraries(performance_rips_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) if (TBB_FOUND) - target_link_libraries(rips_multifield_persistence ${TBB_RELEASE_LIBRARY}) - target_link_libraries(performance_rips_persistence ${TBB_RELEASE_LIBRARY}) - endif() + target_link_libraries(rips_multifield_persistence ${TBB_LIBRARIES}) + target_link_libraries(performance_rips_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) add_test(rips_multifield_persistence_2_71 ${CMAKE_CURRENT_BINARY_DIR}/rips_multifield_persistence ${CMAKE_SOURCE_DIR}/data/points/Kl.txt -r 0.2 -d 3 -p 2 -q 71 -m 100) + endif(GMPXX_FOUND) +else() + # message(WARNING "GMP not found.") +endif(GMP_FOUND) if(CGAL_FOUND) add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp) - target_link_libraries(alpha_complex_3d_persistence ${Boost_SYSTEM_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES} ${CGAL_LIBRARY}) + target_link_libraries(alpha_complex_3d_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) if (TBB_FOUND) - target_link_libraries(alpha_complex_3d_persistence ${TBB_RELEASE_LIBRARY}) - endif() + target_link_libraries(alpha_complex_3d_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) add_test(alpha_complex_3d_persistence_2_0_5 ${CMAKE_CURRENT_BINARY_DIR}/alpha_complex_3d_persistence ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off 2 0.45) @@ -72,23 +78,27 @@ if(GMPXX_FOUND AND GMP_FOUND) target_link_libraries(alpha_complex_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) add_executable(periodic_alpha_complex_3d_persistence periodic_alpha_complex_3d_persistence.cpp) - target_link_libraries(periodic_alpha_complex_3d_persistence ${Boost_SYSTEM_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES} ${CGAL_LIBRARY}) + target_link_libraries(periodic_alpha_complex_3d_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + + add_executable(custom_persistence_sort custom_persistence_sort.cpp) + target_link_libraries(custom_persistence_sort ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) if (TBB_FOUND) - target_link_libraries(alpha_complex_persistence ${TBB_RELEASE_LIBRARY}) - target_link_libraries(periodic_alpha_complex_3d_persistence ${TBB_RELEASE_LIBRARY}) - endif() + target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES}) + target_link_libraries(periodic_alpha_complex_3d_persistence ${TBB_LIBRARIES}) + target_link_libraries(custom_persistence_sort ${TBB_LIBRARIES}) + endif(TBB_FOUND) add_test(alpha_complex_persistence_2_0_45 ${CMAKE_CURRENT_BINARY_DIR}/alpha_complex_persistence ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -m 0.45 -p 2) - add_test(periodic_alpha_complex_3d_persistence_2_0 ${CMAKE_CURRENT_BINARY_DIR}/periodic_alpha_complex_3d_persistence ${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off 2 0) + add_test(periodic_alpha_complex_3d_persistence_2_0 ${CMAKE_CURRENT_BINARY_DIR}/periodic_alpha_complex_3d_persistence ${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off ${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt 2 0) + add_test(custom_persistence_sort ${CMAKE_CURRENT_BINARY_DIR}/custom_persistence_sort) else() message(WARNING "Eigen3 not found. Version 3.1.0 is required for Alpha shapes feature.") - endif() + endif(EIGEN3_FOUND) else() message(WARNING "CGAL version: ${CGAL_VERSION} is too old to compile Alpha shapes feature. Version 4.6.0 is required.") endif () else() # message(WARNING "CGAL not found.") - endif() + endif(CGAL_FOUND) -endif() diff --git a/src/Persistent_cohomology/example/custom_persistence_sort.cpp b/src/Persistent_cohomology/example/custom_persistence_sort.cpp new file mode 100644 index 00000000..d591a0d0 --- /dev/null +++ b/src/Persistent_cohomology/example/custom_persistence_sort.cpp @@ -0,0 +1,116 @@ +/* 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) 2014 INRIA Saclay (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <CGAL/Epick_d.h> +#include <CGAL/point_generators_d.h> +#include <CGAL/algorithm.h> +#include <CGAL/assertions.h> + +#include <gudhi/Alpha_complex.h> +#include <gudhi/Persistent_cohomology.h> + +#include <iostream> +#include <iterator> +#include <vector> +#include <fstream> // for std::ofstream +#include <algorithm> // for std::sort + + +using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<3> >; +using Point = Kernel::Point_d; +using Alpha_complex = Gudhi::alphacomplex::Alpha_complex<Kernel>; + +std::vector<Point> random_points() { + // Instanciate a random point generator + CGAL::Random rng(0); + + // Generate "points_number" random points in a vector + std::vector<Point> points; + + // Generates 1000 random 3D points on a sphere of radius 4.0 + CGAL::Random_points_on_sphere_d<Point> rand_outside(3, 4.0, rng); + CGAL::cpp11::copy_n(rand_outside, 1000, std::back_inserter(points)); + // Generates 2000 random 3D points in a sphere of radius 3.0 + CGAL::Random_points_in_ball_d<Point> rand_inside(3, 3.0, rng); + CGAL::cpp11::copy_n(rand_inside, 2000, std::back_inserter(points)); + + return points; +} + +/* + * Compare two intervals by dimension, then by length. + */ +struct cmp_intervals_by_dim_then_length { + explicit cmp_intervals_by_dim_then_length(Alpha_complex * sc) + : sc_(sc) { } + + template<typename Persistent_interval> + bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) { + if (sc_->dimension(get < 0 > (p1)) == sc_->dimension(get < 0 > (p2))) + return (sc_->filtration(get < 1 > (p1)) - sc_->filtration(get < 0 > (p1)) + > sc_->filtration(get < 1 > (p2)) - sc_->filtration(get < 0 > (p2))); + else + return (sc_->dimension(get < 0 > (p1)) > sc_->dimension(get < 0 > (p2))); + } + Alpha_complex* sc_; +}; + +int main(int argc, char **argv) { + std::vector<Point> points = random_points(); + + // Alpha complex persistence computation from generated points + Alpha_complex alpha_complex_from_points(points, 0.6); + + using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology< Alpha_complex, + Gudhi::persistent_cohomology::Field_Zp >; + Persistent_cohomology pcoh(alpha_complex_from_points); + + // initializes the coefficient field for homology - Z/3Z + pcoh.init_coefficients(3); + pcoh.compute_persistent_cohomology(0.2); + + // Custom sort and output persistence + cmp_intervals_by_dim_then_length cmp(&alpha_complex_from_points); + auto persistent_pairs = pcoh.get_persistent_pairs(); + std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); + for (auto pair : persistent_pairs) { + std::cout << alpha_complex_from_points.dimension(get<0>(pair)) << " " + << alpha_complex_from_points.filtration(get<0>(pair)) << " " + << alpha_complex_from_points.filtration(get<1>(pair)) << std::endl; + } + + // Persistent Betti numbers + std::cout << "The persistent Betti numbers in interval [0.40, 0.41] are : "; + for (int dim = 0; dim < alpha_complex_from_points.dimension(); dim++) + std::cout << "b" << dim << " = " << pcoh.persistent_betti_number(dim, 0.40, 0.41) << " ; "; + std::cout << std::endl; + + // Betti numbers + std::vector<int> betti_numbers = pcoh.betti_numbers(); + std::cout << "The Betti numbers are : "; + for (std::size_t i = 0; i < betti_numbers.size(); i++) + std::cout << "b" << i << " = " << betti_numbers[i] << " ; "; + std::cout << std::endl; + + return 0; +} + diff --git a/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp index b8e1097d..a199fea1 100644 --- a/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp +++ b/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp @@ -125,28 +125,28 @@ Vertex_list from(const Alpha_shape_3::Vertex_handle& vh) { void usage(char * const progName) { std::cerr << "Usage: " << progName << - " path_to_file_graph coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n"; + " path_to_file_graph path_to_iso_cuboid_3_file coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n"; exit(-1); } int main(int argc, char * const argv[]) { // program args management - if (argc != 4) { + if (argc != 5) { std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; usage(argv[0]); } int coeff_field_characteristic = 0; - int returnedScanValue = sscanf(argv[2], "%d", &coeff_field_characteristic); + int returnedScanValue = sscanf(argv[3], "%d", &coeff_field_characteristic); if ((returnedScanValue == EOF) || (coeff_field_characteristic <= 0)) { - std::cerr << "Error: " << argv[2] << " is not correct\n"; + std::cerr << "Error: " << argv[3] << " is not correct\n"; usage(argv[0]); } Filtration_value min_persistence = 0.0; - returnedScanValue = sscanf(argv[3], "%lf", &min_persistence); + returnedScanValue = sscanf(argv[4], "%lf", &min_persistence); if ((returnedScanValue == EOF) || (min_persistence < -1.0)) { - std::cerr << "Error: " << argv[3] << " is not correct\n"; + std::cerr << "Error: " << argv[4] << " is not correct\n"; usage(argv[0]); } @@ -160,11 +160,21 @@ int main(int argc, char * const argv[]) { usage(argv[0]); } + // Read iso_cuboid_3 information from file + std::ifstream iso_cuboid_str(argv[2]); + double x_min, y_min, z_min, x_max, y_max, z_max; + if (iso_cuboid_str.good()) { + iso_cuboid_str >> x_min >> y_min >> z_min >> x_max >> y_max >> z_max; + } else { + std::cerr << "Unable to read file " << argv[2] << std::endl; + usage(argv[0]); + } + // Retrieve the triangulation std::vector<Point_3> lp = off_reader.get_point_cloud(); // Define the periodic cube - P3DT3 pdt(PK::Iso_cuboid_3(0, 0, 0, 1, 1, 1)); + P3DT3 pdt(PK::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max)); // Heuristic for inserting large point sets (if pts is reasonably large) pdt.insert(lp.begin(), lp.end(), true); // As pdt won't be modified anymore switch to 1-sheeted cover if possible diff --git a/src/Persistent_cohomology/example/plain_homology.cpp b/src/Persistent_cohomology/example/plain_homology.cpp index 0a692717..2156ffd9 100644 --- a/src/Persistent_cohomology/example/plain_homology.cpp +++ b/src/Persistent_cohomology/example/plain_homology.cpp @@ -24,6 +24,7 @@ #include <gudhi/Persistent_cohomology.h> #include <iostream> +#include <vector> using namespace Gudhi; @@ -76,9 +77,16 @@ int main() { // Print the result. The format is, on each line: 2 dim 0 inf // where 2 represents the field, dim the dimension of the feature. - // 2 0 0 inf - // 2 0 0 inf - // 2 1 0 inf + // 2 0 0 inf + // 2 0 0 inf + // 2 1 0 inf // means that in Z/2Z-homology, the Betti numbers are b0=2 and b1=1. pcoh.output_diagram(); + + // Print the Betti numbers are b0=2 and b1=1. + std::cout << std::endl; + std::cout << "The Betti numbers are : "; + for (int i = 0; i < st.dimension(); i++) + std::cout << "b" << i << " = " << pcoh.betti_number(i) << " ; "; + std::cout << std::endl; } |