diff options
author | ROUVREAU Vincent <vincent.rouvreau@inria.fr> | 2020-06-03 14:56:43 +0200 |
---|---|---|
committer | ROUVREAU Vincent <vincent.rouvreau@inria.fr> | 2020-06-03 14:56:43 +0200 |
commit | e41edaed6e4e77438a2ab9da69862dc311602758 (patch) | |
tree | c16ef385d051b9a9fbbd28861f2398aad1d24911 /src/Collapse | |
parent | 08ffcedc11119eba2f5a0a6f22329df197c86107 (diff) |
Add an example that tests the persistence is conserved when using edge collapse, versus a classical rips
Diffstat (limited to 'src/Collapse')
3 files changed, 170 insertions, 3 deletions
diff --git a/src/Collapse/example/CMakeLists.txt b/src/Collapse/example/CMakeLists.txt index 6cf3bf07..96354489 100644 --- a/src/Collapse/example/CMakeLists.txt +++ b/src/Collapse/example/CMakeLists.txt @@ -8,3 +8,16 @@ if (TBB_FOUND) endif() add_test(NAME Edge_collapse_example_basic COMMAND $<TARGET_FILE:Edge_collapse_example_basic>) + +# Point cloud +add_executable ( Edge_collapse_conserve_persistence edge_collapse_conserve_persistence.cpp ) + +if (TBB_FOUND) + target_link_libraries(Edge_collapse_conserve_persistence ${TBB_LIBRARIES}) +endif() + +add_test(NAME Edge_collapse_conserve_persistence_1 COMMAND $<TARGET_FILE:Edge_collapse_conserve_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "1.") + + add_test(NAME Edge_collapse_conserve_persistence_2 COMMAND $<TARGET_FILE:Edge_collapse_conserve_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2.") diff --git a/src/Collapse/example/edge_collapse_conserve_persistence.cpp b/src/Collapse/example/edge_collapse_conserve_persistence.cpp new file mode 100644 index 00000000..47accb0c --- /dev/null +++ b/src/Collapse/example/edge_collapse_conserve_persistence.cpp @@ -0,0 +1,157 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2020 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <gudhi/Flag_complex_sparse_matrix.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Points_off_io.h> +#include <gudhi/graph_simplicial_complex.h> + +#include<utility> // for std::pair +#include<vector> +#include<tuple> + +// Types definition + +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Vertex_handle = Simplex_tree::Vertex_handle; +using Point = std::vector<Filtration_value>; +using Vector_of_points = std::vector<Point>; + +using Flag_complex_sparse_matrix = Gudhi::collapse::Flag_complex_sparse_matrix<Vertex_handle, Filtration_value>; +using Proximity_graph = Flag_complex_sparse_matrix::Proximity_graph; + +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; + +using Persistence_pair = std::tuple<int, Filtration_value, Filtration_value>; +/* + * Compare two intervals by dimension, then by length. + */ +struct cmp_intervals_by_dim_then_length { + explicit cmp_intervals_by_dim_then_length(Simplex_tree * 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))); + } + Simplex_tree* sc_; +}; + +std::vector<Persistence_pair> get_persistence_pairs(Simplex_tree& st, int ambient_dim) { + std::vector<Persistence_pair> ppairs; + st.expansion(ambient_dim); + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(st); + // initializes the coefficient field for homology - must be a prime number + int p = 11; + pcoh.init_coefficients(p); + + // Default min_interval_length = 0. + pcoh.compute_persistent_cohomology(); + // Custom sort and output persistence + cmp_intervals_by_dim_then_length cmp(&st); + auto persistent_pairs = pcoh.get_persistent_pairs(); + std::sort(std::begin(persistent_pairs), std::end(persistent_pairs), cmp); + for (auto pair : persistent_pairs) { + ppairs.push_back({st.dimension(get<0>(pair)), + st.filtration(get<0>(pair)), + st.filtration(get<1>(pair)) }); + } + return ppairs; +} + +int main(int argc, char* argv[]) { + if (argc != 3) { + std::cerr << "This program requires an OFF file and minimal threshold value as parameter\n"; + std::cerr << "For instance: ./Edge_collapse_conserve_persistence ../../data/points/tore3D_300.off 1.\n"; + exit(-1); // ----- >> + } + + std::string off_file_points {argv[1]}; + double threshold {atof(argv[2])}; + + Gudhi::Points_off_reader<Point> off_reader(off_file_points); + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << off_file_points << "\n"; + exit(-1); // ----- >> + } + + Vector_of_points point_vector = off_reader.get_point_cloud(); + if (point_vector.size() <= 0) { + std::cerr << "Empty point cloud." << std::endl; + exit(-1); // ----- >> + } + + Proximity_graph proximity_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(), + threshold, + Gudhi::Euclidean_distance()); + + if (num_edges(proximity_graph) <= 0) { + std::cerr << "Total number of egdes are zero." << std::endl; + exit(-1); + } + + // ***** Simplex tree from a flag complex built after collapse ***** + Flag_complex_sparse_matrix mat_filt_edge_coll(proximity_graph); + + Simplex_tree stree_from_collapse; + mat_filt_edge_coll.filtered_edge_collapse( + [&stree_from_collapse](const std::vector<Vertex_handle>& edge, Filtration_value filtration) { + // insert the 2 vertices with a 0. filtration value just like a Rips + stree_from_collapse.insert_simplex({edge[0]}, 0.); + stree_from_collapse.insert_simplex({edge[1]}, 0.); + // insert the edge + stree_from_collapse.insert_simplex(edge, filtration); + }); + + // ***** Simplex tree from the complete flag complex ***** + Simplex_tree stree_wo_collapse; + stree_wo_collapse.insert_graph(proximity_graph); + + int ambient_dim = point_vector[0].size(); + + std::vector<Persistence_pair> ppairs_from_collapse = get_persistence_pairs(stree_from_collapse, ambient_dim); + std::vector<Persistence_pair> ppairs_wo_collapse = get_persistence_pairs(stree_wo_collapse, ambient_dim); + + if (ppairs_wo_collapse.size() != ppairs_from_collapse.size()) { + std::cerr << "Number of persistence pairs with collapse is " << ppairs_from_collapse.size() << std::endl; + std::cerr << "Number of persistence pairs without collapse is " << ppairs_wo_collapse.size() << std::endl; + exit(-1); + } + + int return_value = 0; + auto ppwoc_ptr = ppairs_wo_collapse.begin(); + for (auto ppfc: ppairs_from_collapse) { + if (ppfc != *ppwoc_ptr) { + return_value++; + std::cerr << "Without collapse: " + << std::get<0>(*ppwoc_ptr) << " " << std::get<1>(*ppwoc_ptr) << " " << std::get<2>(*ppwoc_ptr) + << " - With collapse: " + << std::get<0>(ppfc) << " " << std::get<1>(ppfc) << " " << std::get<2>(ppfc) << std::endl; + } + std::clog << " ** Without collapse: " + << std::get<0>(*ppwoc_ptr) << " " << std::get<1>(*ppwoc_ptr) << " " << std::get<2>(*ppwoc_ptr) + << " - With collapse: " + << std::get<0>(ppfc) << " " << std::get<1>(ppfc) << " " << std::get<2>(ppfc) << std::endl; + ppwoc_ptr++; + } + return return_value; +} diff --git a/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp b/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp index b9130d4c..067b29e3 100644 --- a/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp +++ b/src/Collapse/utilities/point_cloud_edge_collapse_rips_persistence.cpp @@ -33,7 +33,6 @@ using Proximity_graph = Flag_complex_sparse_matrix::Proximity_graph; using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; -using Distance_matrix = std::vector<std::vector<Filtration_value>>; 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); @@ -54,8 +53,6 @@ int main(int argc, char* argv[]) { std::cout << min_persistence << ", " << threshold << ", " << dim_max << ", " << off_file_points << ", " << filediag << std::endl; - Distance_matrix sparse_distances; - Gudhi::Points_off_reader<Point> off_reader(off_file_points); if (!off_reader.is_valid()) { std::cerr << "Unable to read file " << off_file_points << "\n"; |