diff options
author | vrouvrea <vrouvrea@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2017-03-16 12:44:20 +0000 |
---|---|---|
committer | vrouvrea <vrouvrea@636b058d-ea47-450e-bf9e-a15bfbe3eedb> | 2017-03-16 12:44:20 +0000 |
commit | 19bb0e8b6925c07d68a152cf70c7b5618a1f4af0 (patch) | |
tree | a057c11e74ae88056f82bf5643e38f9138e7f3fc /src/Persistent_cohomology | |
parent | 11af9b60e8d2b007d5bedba9f599f4eeb5468fea (diff) | |
parent | 90f6c9e1ee1c60e8e5d1de361bcecd89948b2630 (diff) |
Merge last trunk modifications
git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/ST_cythonize@2194 636b058d-ea47-450e-bf9e-a15bfbe3eedb
Former-commit-id: a8fa0a9eeb6b7064cb1d3d3780d5ab58818f49ec
Diffstat (limited to 'src/Persistent_cohomology')
8 files changed, 699 insertions, 140 deletions
diff --git a/src/Persistent_cohomology/doc/COPYRIGHT b/src/Persistent_cohomology/doc/COPYRIGHT new file mode 100644 index 00000000..34345bef --- /dev/null +++ b/src/Persistent_cohomology/doc/COPYRIGHT @@ -0,0 +1,19 @@ +The files of this directory are part of the Gudhi Library. The Gudhi library +(Geometric Understanding in Higher Dimensions) is a generic C++ library for +computational topology. + +Author(s): Clément Maria + +Copyright (C) 2015 INRIA + +This program is free software: you can redistribute it and/or modify it under +the terms of the GNU General Public License as published by the Free Software +Foundation, either version 3 of the License, or (at your option) any later +version. + +This program is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along with +this program. If not, see <http://www.gnu.org/licenses/>. diff --git a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h index 40dd3f93..e17e5926 100644 --- a/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h +++ b/src/Persistent_cohomology/doc/Intro_persistent_cohomology.h @@ -172,6 +172,28 @@ Persistent_cohomology/alpha_complex_3d_persistence.cpp</a> computes the persiste 2 1 0.0934117 1.00003 2 2 0.56444 1.03938 \endcode +\li <a href="_persistent_cohomology_2exact_alpha_complex_3d_persistence_8cpp-example.html"> +Persistent_cohomology/exact_alpha_complex_3d_persistence.cpp</a> computes the persistent homology with +\f$\mathbb{Z}/2\mathbb{Z}\f$ coefficients of the alpha complex on points sampling from an OFF file. +Here, as CGAL computes the exact values, it is slower, but it is necessary when points are on a grid +for instance. +\code $> ./exact_alpha_complex_3d_persistence ../../data/points/sphere3D_pts_on_grid.off 2 0.1 \endcode +\code Simplex_tree dim: 3 +2 0 0 inf +2 2 0.0002 0.2028 \endcode + +\li <a href="_persistent_cohomology_2weighted_alpha_complex_3d_persistence_8cpp-example.html"> +Persistent_cohomology/weighted_alpha_complex_3d_persistence.cpp</a> computes the persistent homology with +\f$\mathbb{Z}/2\mathbb{Z}\f$ coefficients of the weighted alpha complex on points sampling from an OFF file +and a weights file. +\code $> ./weighted_alpha_complex_3d_persistence ../../data/points/tore3D_300.off +../../data/points/tore3D_300.weights 2 0.45 \endcode +\code Simplex_tree dim: 3 +2 -0 0 inf +2 1 0.0682162 1.0001 +2 1 0.0934117 1.00003 +2 2 0.56444 1.03938 \endcode + \li <a href="_persistent_cohomology_2alpha_complex_persistence_8cpp-example.html"> Persistent_cohomology/alpha_complex_persistence.cpp</a> computes the persistent homology with \f$\mathbb{Z}/p\mathbb{Z}\f$ coefficients of the alpha complex on points sampling from an OFF file. diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index d2a84b1e..20700d86 100644 --- a/src/Persistent_cohomology/example/CMakeLists.txt +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -34,38 +34,58 @@ endif() add_test(plain_homology ${CMAKE_CURRENT_BINARY_DIR}/plain_homology) add_test(persistence_from_simple_simplex_tree ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_simple_simplex_tree 1 0) -add_test(rips_distance_matrix ${CMAKE_CURRENT_BINARY_DIR}/rips_distance_matrix_persistence ${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv -r 1.0 -d 3 -p 3 -m 0) -add_test(rips_persistence_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3) -add_test(rips_persistence_step_by_step_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence_step_by_step ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3) -add_test(rips_persistence_via_boundary_matrix_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence_via_boundary_matrix ${CMAKE_SOURCE_DIR}/data/points/Kl.off -r 0.16 -d 3 -p 3 -m 100) -add_test(persistence_from_file_3_2_0 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file ${CMAKE_SOURCE_DIR}/data/filtered_simplicial_complex/bunny_5000_complex.fsc -p 2 -m 0) -add_test(persistence_from_file_3_3_100 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file ${CMAKE_SOURCE_DIR}/data/filtered_simplicial_complex/bunny_5000_complex.fsc -p 3 -m 100) +add_test(rips_distance_matrix ${CMAKE_CURRENT_BINARY_DIR}/rips_distance_matrix_persistence + ${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv -r 1.0 -d 3 -p 3 -m 0) +add_test(rips_persistence_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence + ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3) +add_test(rips_persistence_step_by_step_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence_step_by_step + ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3) +add_test(rips_persistence_via_boundary_matrix_3 ${CMAKE_CURRENT_BINARY_DIR}/rips_persistence_via_boundary_matrix + ${CMAKE_SOURCE_DIR}/data/points/Kl.off -r 0.16 -d 3 -p 3 -m 100) +add_test(persistence_from_file_3_2_0 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file + ${CMAKE_SOURCE_DIR}/data/filtered_simplicial_complex/bunny_5000_complex.fsc -p 2 -m 0) +add_test(persistence_from_file_3_3_100 ${CMAKE_CURRENT_BINARY_DIR}/persistence_from_file + ${CMAKE_SOURCE_DIR}/data/filtered_simplicial_complex/bunny_5000_complex.fsc -p 3 -m 100) if(GMP_FOUND) if(GMPXX_FOUND) 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}) + target_link_libraries(rips_multifield_persistence + ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY} ${GMPXX_LIBRARIES} ${GMP_LIBRARIES}) if (TBB_FOUND) target_link_libraries(rips_multifield_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/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 2 -q 71) + add_test(rips_multifield_persistence_2_71 ${CMAKE_CURRENT_BINARY_DIR}/rips_multifield_persistence + ${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 2 -q 71) endif(GMPXX_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} ${CGAL_LIBRARY}) + add_executable(exact_alpha_complex_3d_persistence exact_alpha_complex_3d_persistence.cpp) + target_link_libraries(exact_alpha_complex_3d_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) + add_executable(weighted_alpha_complex_3d_persistence weighted_alpha_complex_3d_persistence.cpp) + target_link_libraries(weighted_alpha_complex_3d_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY}) if (TBB_FOUND) target_link_libraries(alpha_complex_3d_persistence ${TBB_LIBRARIES}) + target_link_libraries(exact_alpha_complex_3d_persistence ${TBB_LIBRARIES}) + target_link_libraries(weighted_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) + 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) + add_test(exact_alpha_complex_3d_persistence_2_0_5 ${CMAKE_CURRENT_BINARY_DIR}/exact_alpha_complex_3d_persistence + ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off 2 0.45) + add_test(weighted_alpha_complex_3d_persistence_2_0_5 ${CMAKE_CURRENT_BINARY_DIR}/weighted_alpha_complex_3d_persistence + ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.weights 2 0.45) if (NOT CGAL_VERSION VERSION_LESS 4.7.0) if (EIGEN3_FOUND) add_executable (alpha_complex_persistence alpha_complex_persistence.cpp) - target_link_libraries(alpha_complex_persistence ${Boost_SYSTEM_LIBRARY} ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + 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} ${CGAL_LIBRARY}) @@ -78,8 +98,10 @@ if(CGAL_FOUND) 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 ${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt 2 0) + 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 ${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) endif(EIGEN3_FOUND) endif (NOT CGAL_VERSION VERSION_LESS 4.7.0) diff --git a/src/Persistent_cohomology/example/alpha_complex_3d_helper.h b/src/Persistent_cohomology/example/alpha_complex_3d_helper.h new file mode 100644 index 00000000..ceefd2ba --- /dev/null +++ b/src/Persistent_cohomology/example/alpha_complex_3d_helper.h @@ -0,0 +1,76 @@ +/* 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/>. + */ + +#ifndef EXAMPLE_PERSISTENT_COHOMOLOGY_ALPHA_COMPLEX_3D_HELPER_H_ +#define EXAMPLE_PERSISTENT_COHOMOLOGY_ALPHA_COMPLEX_3D_HELPER_H_ + +template<class Vertex_list, class Cell_handle> +Vertex_list from_cell(const Cell_handle& ch) { + Vertex_list the_list; + for (auto i = 0; i < 4; i++) { +#ifdef DEBUG_TRACES + std::cout << "from cell[" << i << "]=" << ch->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + the_list.push_back(ch->vertex(i)); + } + return the_list; +} + +template<class Vertex_list, class Facet> +Vertex_list from_facet(const Facet& fct) { + Vertex_list the_list; + for (auto i = 0; i < 4; i++) { + if (fct.second != i) { +#ifdef DEBUG_TRACES + std::cout << "from facet=[" << i << "]" << fct.first->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + the_list.push_back(fct.first->vertex(i)); + } + } + return the_list; +} + +template<class Vertex_list, class Edge_3> +Vertex_list from_edge(const Edge_3& edg) { + Vertex_list the_list; + for (auto i = 0; i < 4; i++) { + if ((edg.second == i) || (edg.third == i)) { +#ifdef DEBUG_TRACES + std::cout << "from edge[" << i << "]=" << edg.first->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + the_list.push_back(edg.first->vertex(i)); + } + } + return the_list; +} + +template<class Vertex_list, class Vertex_handle> +Vertex_list from_vertex(const Vertex_handle& vh) { + Vertex_list the_list; +#ifdef DEBUG_TRACES + std::cout << "from vertex=" << vh->point() << std::endl; +#endif // DEBUG_TRACES + the_list.push_back(vh); + return the_list; +} + +#endif // EXAMPLE_PERSISTENT_COHOMOLOGY_ALPHA_COMPLEX_3D_HELPER_H_ diff --git a/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp index 978dc942..ffb98283 100644 --- a/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp +++ b/src/Persistent_cohomology/example/alpha_complex_3d_persistence.cpp @@ -30,6 +30,8 @@ #include <CGAL/Alpha_shape_3.h> #include <CGAL/iterator.h> +#include "alpha_complex_3d_helper.h" + #include <fstream> #include <cmath> #include <string> @@ -40,84 +42,39 @@ #include <vector> // Alpha_shape_3 templates type definitions -typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; -typedef CGAL::Alpha_shape_vertex_base_3<Kernel> Vb; -typedef CGAL::Alpha_shape_cell_base_3<Kernel> Fb; -typedef CGAL::Triangulation_data_structure_3<Vb, Fb> Tds; -typedef CGAL::Delaunay_triangulation_3<Kernel, Tds> Triangulation_3; -typedef CGAL::Alpha_shape_3<Triangulation_3> Alpha_shape_3; +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel>; +using Fb = CGAL::Alpha_shape_cell_base_3<Kernel>; +using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>; +using Triangulation_3 = CGAL::Delaunay_triangulation_3<Kernel, Tds>; +using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3>; // From file type definition -typedef Kernel::Point_3 Point_3; +using Point_3 = Kernel::Point_3; // filtration with alpha values needed type definition -typedef Alpha_shape_3::FT Alpha_value_type; -typedef CGAL::Object Object; -typedef CGAL::Dispatch_output_iterator< -CGAL::cpp11::tuple<Object, Alpha_value_type>, -CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >, - std::back_insert_iterator< std::vector<Alpha_value_type> > > > Dispatch; -typedef Alpha_shape_3::Cell_handle Cell_handle; -typedef Alpha_shape_3::Facet Facet; -typedef Alpha_shape_3::Edge Edge_3; -typedef std::list<Alpha_shape_3::Vertex_handle> Vertex_list; +using Alpha_value_type = Alpha_shape_3::FT; +using Object = CGAL::Object; +using Dispatch = CGAL::Dispatch_output_iterator< + CGAL::cpp11::tuple<Object, Alpha_value_type>, + CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >, + std::back_insert_iterator< std::vector<Alpha_value_type> > > >; +using Cell_handle = Alpha_shape_3::Cell_handle; +using Facet = Alpha_shape_3::Facet; +using Edge_3 = Alpha_shape_3::Edge; +using Vertex_handle = Alpha_shape_3::Vertex_handle; +using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>; // gudhi type definition -typedef Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence> ST; -typedef ST::Filtration_value Filtration_value; -typedef ST::Vertex_handle Simplex_tree_vertex; -typedef std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex > Alpha_shape_simplex_tree_map; -typedef std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex> Alpha_shape_simplex_tree_pair; -typedef std::vector< Simplex_tree_vertex > Simplex_tree_vector_vertex; -typedef Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp > PCOH; - -Vertex_list from(const Cell_handle& ch) { - Vertex_list the_list; - for (auto i = 0; i < 4; i++) { -#ifdef DEBUG_TRACES - std::cout << "from cell[" << i << "]=" << ch->vertex(i)->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(ch->vertex(i)); - } - return the_list; -} - -Vertex_list from(const Facet& fct) { - Vertex_list the_list; - for (auto i = 0; i < 4; i++) { - if (fct.second != i) { -#ifdef DEBUG_TRACES - std::cout << "from facet=[" << i << "]" << fct.first->vertex(i)->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(fct.first->vertex(i)); - } - } - return the_list; -} - -Vertex_list from(const Edge_3& edg) { - Vertex_list the_list; - for (auto i = 0; i < 4; i++) { - if ((edg.second == i) || (edg.third == i)) { -#ifdef DEBUG_TRACES - std::cout << "from edge[" << i << "]=" << edg.first->vertex(i)->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(edg.first->vertex(i)); - } - } - return the_list; -} - -Vertex_list from(const Alpha_shape_3::Vertex_handle& vh) { - Vertex_list the_list; -#ifdef DEBUG_TRACES - std::cout << "from vertex=" << vh->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(vh); - return the_list; -} - -void usage(char * const progName) { +using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >; +using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; +using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >; +using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp >; + +void usage(const std::string& progName) { std::cerr << "Usage: " << progName << " path_to_file_graph coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n"; exit(-1); @@ -185,30 +142,29 @@ int main(int argc, char * const argv[]) { for (auto object_iterator : the_objects) { // Retrieve Alpha shape vertex list from object if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { - vertex_list = from(*cell); + vertex_list = from_cell<Vertex_list, Cell_handle>(*cell); count_cells++; if (dim_max < 3) { // Cell is of dim 3 dim_max = 3; } } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) { - vertex_list = from(*facet); + vertex_list = from_facet<Vertex_list, Facet>(*facet); count_facets++; if (dim_max < 2) { // Facet is of dim 2 dim_max = 2; } } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) { - vertex_list = from(*edge); + vertex_list = from_edge<Vertex_list, Edge_3>(*edge); count_edges++; if (dim_max < 1) { // Edge_3 is of dim 1 dim_max = 1; } - } else if (const Alpha_shape_3::Vertex_handle * vertex = - CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) { + } else if (const Vertex_handle * vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) { count_vertices++; - vertex_list = from(*vertex); + vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex); } // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex Simplex_tree_vector_vertex the_simplex_tree; diff --git a/src/Persistent_cohomology/example/exact_alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/exact_alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..c016a596 --- /dev/null +++ b/src/Persistent_cohomology/example/exact_alpha_complex_3d_persistence.cpp @@ -0,0 +1,244 @@ +/* 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 <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_3D_off_io.h> +#include <boost/variant.hpp> + +#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> +#include <CGAL/Delaunay_triangulation_3.h> +#include <CGAL/Alpha_shape_3.h> +#include <CGAL/iterator.h> + +#include "alpha_complex_3d_helper.h" + +#include <fstream> +#include <cmath> +#include <string> +#include <tuple> +#include <map> +#include <utility> +#include <list> +#include <vector> + +// Alpha_shape_3 templates type definitions +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using Exact_tag = CGAL::Tag_true; +using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel, CGAL::Default, Exact_tag>; +using Fb = CGAL::Alpha_shape_cell_base_3<Kernel, CGAL::Default, Exact_tag>; +using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>; +using Triangulation_3 = CGAL::Delaunay_triangulation_3<Kernel, Tds>; +using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3, Exact_tag>; + +// From file type definition +using Point_3 = Kernel::Point_3; + +// filtration with alpha values needed type definition +using Alpha_value_type = Alpha_shape_3::FT; +using Object = CGAL::Object; +using Dispatch = CGAL::Dispatch_output_iterator< + CGAL::cpp11::tuple<Object, Alpha_value_type>, + CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >, + std::back_insert_iterator< std::vector<Alpha_value_type> > > >; +using Cell_handle = Alpha_shape_3::Cell_handle; +using Facet = Alpha_shape_3::Facet; +using Edge_3 = Alpha_shape_3::Edge; +using Vertex_handle = Alpha_shape_3::Vertex_handle; +using Vertex_list = std::list<Vertex_handle>; + +// gudhi type definition +using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >; +using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; +using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >; +using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp >; + +void usage(char * const progName) { + std::cerr << "Usage: " << progName << + " path_to_file_graph 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) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; + usage(argv[0]); + } + + int coeff_field_characteristic = atoi(argv[2]); + + Filtration_value min_persistence = 0.0; + int returnedScanValue = sscanf(argv[3], "%f", &min_persistence); + if ((returnedScanValue == EOF) || (min_persistence < -1.0)) { + std::cerr << "Error: " << argv[3] << " is not correct\n"; + usage(argv[0]); + } + + // Read points from file + std::string offInputFile(argv[1]); + // Read the OFF file (input file name given as parameter) and triangulate points + Gudhi::Points_3D_off_reader<Point_3> off_reader(offInputFile); + // Check the read operation was correct + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << offInputFile << std::endl; + usage(argv[0]); + } + + // Retrieve the triangulation + std::vector<Point_3> lp = off_reader.get_point_cloud(); + + // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. + Alpha_shape_3 as(lp.begin(), lp.end(), 0, Alpha_shape_3::GENERAL); +#ifdef DEBUG_TRACES + std::cout << "Alpha shape computed in GENERAL mode" << std::endl; +#endif // DEBUG_TRACES + + // filtration with alpha values from alpha shape + std::vector<Object> the_objects; + std::vector<Alpha_value_type> the_alpha_values; + + Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects), + std::back_inserter(the_alpha_values)); + + as.filtration_with_alpha_values(disp); +#ifdef DEBUG_TRACES + std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl; +#endif // DEBUG_TRACES + + Alpha_shape_3::size_type count_vertices = 0; + Alpha_shape_3::size_type count_edges = 0; + Alpha_shape_3::size_type count_facets = 0; + Alpha_shape_3::size_type count_cells = 0; + + // Loop on objects vector + Vertex_list vertex_list; + ST simplex_tree; + Alpha_shape_simplex_tree_map map_cgal_simplex_tree; + std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin(); + int dim_max = 0; + Filtration_value filtration_max = 0.0; + for (auto object_iterator : the_objects) { + // Retrieve Alpha shape vertex list from object + if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { + vertex_list = from_cell<Vertex_list, Cell_handle>(*cell); + count_cells++; + if (dim_max < 3) { + // Cell is of dim 3 + dim_max = 3; + } + } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) { + vertex_list = from_facet<Vertex_list, Facet>(*facet); + count_facets++; + if (dim_max < 2) { + // Facet is of dim 2 + dim_max = 2; + } + } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) { + vertex_list = from_edge<Vertex_list, Edge_3>(*edge); + count_edges++; + if (dim_max < 1) { + // Edge_3 is of dim 1 + dim_max = 1; + } + } else if (const Vertex_handle * vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex); + } + // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex + Simplex_tree_vector_vertex the_simplex_tree; + for (auto the_alpha_shape_vertex : vertex_list) { + Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex); + if (the_map_iterator == map_cgal_simplex_tree.end()) { + // alpha shape not found + Simplex_tree_vertex vertex = map_cgal_simplex_tree.size(); +#ifdef DEBUG_TRACES + std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl; +#endif // DEBUG_TRACES + the_simplex_tree.push_back(vertex); + map_cgal_simplex_tree.insert(Alpha_shape_simplex_tree_pair(the_alpha_shape_vertex, vertex)); + } else { + // alpha shape found + Simplex_tree_vertex vertex = the_map_iterator->second; +#ifdef DEBUG_TRACES + std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl; +#endif // DEBUG_TRACES + the_simplex_tree.push_back(vertex); + } + } + // Construction of the simplex_tree + // you can also use the_alpha_value_iterator->exact() + Filtration_value filtr = /*std::sqrt*/CGAL::to_double(the_alpha_value_iterator->exact()); +#ifdef DEBUG_TRACES + std::cout << "filtration = " << filtr << std::endl; +#endif // DEBUG_TRACES + if (filtr > filtration_max) { + filtration_max = filtr; + } + simplex_tree.insert_simplex(the_simplex_tree, filtr); + if (the_alpha_value_iterator != the_alpha_values.end()) + ++the_alpha_value_iterator; + else + std::cout << "This shall not happen" << std::endl; + } + simplex_tree.set_filtration(filtration_max); + simplex_tree.set_dimension(dim_max); + +#ifdef DEBUG_TRACES + std::cout << "vertices \t\t" << count_vertices << std::endl; + std::cout << "edges \t\t" << count_edges << std::endl; + std::cout << "facets \t\t" << count_facets << std::endl; + std::cout << "cells \t\t" << count_cells << std::endl; + + + std::cout << "Information of the Simplex Tree: " << std::endl; + std::cout << " Number of vertices = " << simplex_tree.num_vertices() << " "; + std::cout << " Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl; + std::cout << " Dimension = " << simplex_tree.dimension() << " "; + std::cout << " filtration = " << simplex_tree.filtration() << std::endl << std::endl; +#endif // DEBUG_TRACES + +#ifdef DEBUG_TRACES + std::cout << "Iterator on vertices: " << std::endl; + for (auto vertex : simplex_tree.complex_vertex_range()) { + std::cout << vertex << " "; + } +#endif // DEBUG_TRACES + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl; + // Compute the persistence diagram of the complex + PCOH pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(coeff_field_characteristic); + + pcoh.compute_persistent_cohomology(min_persistence); + + pcoh.output_diagram(); + + 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 dbc42706..4f6ea40e 100644 --- a/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp +++ b/src/Persistent_cohomology/example/periodic_alpha_complex_3d_persistence.cpp @@ -31,6 +31,8 @@ #include <CGAL/Alpha_shape_3.h> #include <CGAL/iterator.h> +#include "alpha_complex_3d_helper.h" + #include <fstream> #include <cmath> #include <string> @@ -67,6 +69,7 @@ using Dispatch = CGAL::Dispatch_output_iterator< using Cell_handle = Alpha_shape_3::Cell_handle; using Facet = Alpha_shape_3::Facet; using Edge_3 = Alpha_shape_3::Edge; +using Vertex_handle = Alpha_shape_3::Vertex_handle; using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>; // gudhi type definition @@ -79,52 +82,6 @@ using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp >; -Vertex_list from(const Cell_handle& ch) { - Vertex_list the_list; - for (auto i = 0; i < 4; i++) { -#ifdef DEBUG_TRACES - std::cout << "from cell[" << i << "]=" << ch->vertex(i)->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(ch->vertex(i)); - } - return the_list; -} - -Vertex_list from(const Facet& fct) { - Vertex_list the_list; - for (auto i = 0; i < 4; i++) { - if (fct.second != i) { -#ifdef DEBUG_TRACES - std::cout << "from facet=[" << i << "]" << fct.first->vertex(i)->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(fct.first->vertex(i)); - } - } - return the_list; -} - -Vertex_list from(const Edge_3& edg) { - Vertex_list the_list; - for (auto i = 0; i < 4; i++) { - if ((edg.second == i) || (edg.third == i)) { -#ifdef DEBUG_TRACES - std::cout << "from edge[" << i << "]=" << edg.first->vertex(i)->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(edg.first->vertex(i)); - } - } - return the_list; -} - -Vertex_list from(const Alpha_shape_3::Vertex_handle& vh) { - Vertex_list the_list; -#ifdef DEBUG_TRACES - std::cout << "from vertex=" << vh->point() << std::endl; -#endif // DEBUG_TRACES - the_list.push_back(vh); - return the_list; -} - void usage(char * const progName) { std::cerr << "Usage: " << progName << " path_to_file_graph path_to_iso_cuboid_3_file coeff_field_characteristic[integer > 0] min_persistence[float >= -1.0]\n"; @@ -203,21 +160,21 @@ int main(int argc, char * const argv[]) { for (auto object_iterator : the_objects) { // Retrieve Alpha shape vertex list from object if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { - vertex_list = from(*cell); + vertex_list = from_cell<Vertex_list, Cell_handle>(*cell); count_cells++; if (dim_max < 3) { // Cell is of dim 3 dim_max = 3; } } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) { - vertex_list = from(*facet); + vertex_list = from_facet<Vertex_list, Facet>(*facet); count_facets++; if (dim_max < 2) { // Facet is of dim 2 dim_max = 2; } } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) { - vertex_list = from(*edge); + vertex_list = from_edge<Vertex_list, Edge_3>(*edge); count_edges++; if (dim_max < 1) { // Edge_3 is of dim 1 @@ -226,7 +183,7 @@ int main(int argc, char * const argv[]) { } else if (const Alpha_shape_3::Vertex_handle * vertex = CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) { count_vertices++; - vertex_list = from(*vertex); + vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex); } // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex Simplex_tree_vector_vertex the_simplex_tree; diff --git a/src/Persistent_cohomology/example/weighted_alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/weighted_alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..5203548a --- /dev/null +++ b/src/Persistent_cohomology/example/weighted_alpha_complex_3d_persistence.cpp @@ -0,0 +1,263 @@ +/* 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 + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_3D_off_io.h> +#include <boost/variant.hpp> + +#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> +#include <CGAL/Regular_triangulation_euclidean_traits_3.h> +#include <CGAL/Regular_triangulation_3.h> +#include <CGAL/Alpha_shape_3.h> +#include <CGAL/iterator.h> + +#include "alpha_complex_3d_helper.h" + +#include <fstream> +#include <cmath> +#include <string> +#include <tuple> +#include <map> +#include <utility> +#include <list> +#include <vector> +#include <cstdlib> + +// Traits +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using Gt = CGAL::Regular_triangulation_euclidean_traits_3<Kernel>; +using Vb = CGAL::Alpha_shape_vertex_base_3<Gt>; +using Fb = CGAL::Alpha_shape_cell_base_3<Gt>; +using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>; +using Triangulation_3 = CGAL::Regular_triangulation_3<Gt,Tds>; +using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3>; + +// From file type definition +using Point_3 = Gt::Bare_point; +using Weighted_point_3 = Gt::Weighted_point; + +// filtration with alpha values needed type definition +using Alpha_value_type = Alpha_shape_3::FT; +using Object = CGAL::Object; +using Dispatch = CGAL::Dispatch_output_iterator< + CGAL::cpp11::tuple<Object, Alpha_value_type>, + CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >, + std::back_insert_iterator< std::vector<Alpha_value_type> > > >; +using Cell_handle = Alpha_shape_3::Cell_handle; +using Facet = Alpha_shape_3::Facet; +using Edge_3 = Alpha_shape_3::Edge; +using Vertex_handle = Alpha_shape_3::Vertex_handle; +using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>; + +// gudhi type definition +using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >; +using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>; +using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology< + ST, Gudhi::persistent_cohomology::Field_Zp >; + +void usage(char * const progName) { + std::cerr << "Usage: " << progName << + " path_to_file_graph path_to_weight_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 != 5) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; + usage(argv[0]); + } + + int coeff_field_characteristic = atoi(argv[3]); + Filtration_value min_persistence = strtof(argv[4], nullptr); + + // Read points from file + std::string offInputFile(argv[1]); + // Read the OFF file (input file name given as parameter) and triangulate points + Gudhi::Points_3D_off_reader<Point_3> off_reader(offInputFile); + // Check the read operation was correct + if (!off_reader.is_valid()) { + std::cerr << "Unable to read file " << offInputFile << std::endl; + usage(argv[0]); + } + + // Retrieve the triangulation + std::vector<Point_3> lp = off_reader.get_point_cloud(); + + // Read weights information from file + std::ifstream weights_ifstr(argv[2]); + std::vector<Weighted_point_3> wp; + if (weights_ifstr.good()) { + double weight = 0.0; + std::size_t index = 0; + // Attempt read the weight in a double format, return false if it fails + while((weights_ifstr >> weight) && (index < lp.size())) + { + wp.push_back(Weighted_point_3(lp[index], weight)); + index++; + } + if (index != lp.size()) { + std::cerr << "Bad number of weights in file " << argv[2] << std::endl; + usage(argv[0]); + } + } else { + std::cerr << "Unable to read file " << argv[2] << std::endl; + usage(argv[0]); + } + + // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. + Alpha_shape_3 as(lp.begin(), lp.end(), 0, Alpha_shape_3::GENERAL); +#ifdef DEBUG_TRACES + std::cout << "Alpha shape computed in GENERAL mode" << std::endl; +#endif // DEBUG_TRACES + + // filtration with alpha values from alpha shape + std::vector<Object> the_objects; + std::vector<Alpha_value_type> the_alpha_values; + + Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects), + std::back_inserter(the_alpha_values)); + + as.filtration_with_alpha_values(disp); +#ifdef DEBUG_TRACES + std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl; +#endif // DEBUG_TRACES + + Alpha_shape_3::size_type count_vertices = 0; + Alpha_shape_3::size_type count_edges = 0; + Alpha_shape_3::size_type count_facets = 0; + Alpha_shape_3::size_type count_cells = 0; + + // Loop on objects vector + Vertex_list vertex_list; + ST simplex_tree; + Alpha_shape_simplex_tree_map map_cgal_simplex_tree; + std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin(); + int dim_max = 0; + Filtration_value filtration_max = 0.0; + for (auto object_iterator : the_objects) { + // Retrieve Alpha shape vertex list from object + if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { + vertex_list = from_cell<Vertex_list, Cell_handle>(*cell); + count_cells++; + if (dim_max < 3) { + // Cell is of dim 3 + dim_max = 3; + } + } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) { + vertex_list = from_facet<Vertex_list, Facet>(*facet); + count_facets++; + if (dim_max < 2) { + // Facet is of dim 2 + dim_max = 2; + } + } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) { + vertex_list = from_edge<Vertex_list, Edge_3>(*edge); + count_edges++; + if (dim_max < 1) { + // Edge_3 is of dim 1 + dim_max = 1; + } + } else if (const Alpha_shape_3::Vertex_handle * vertex = + CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex); + } + // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex + Simplex_tree_vector_vertex the_simplex_tree; + for (auto the_alpha_shape_vertex : vertex_list) { + Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex); + if (the_map_iterator == map_cgal_simplex_tree.end()) { + // alpha shape not found + Simplex_tree_vertex vertex = map_cgal_simplex_tree.size(); +#ifdef DEBUG_TRACES + std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl; +#endif // DEBUG_TRACES + the_simplex_tree.push_back(vertex); + map_cgal_simplex_tree.insert(Alpha_shape_simplex_tree_pair(the_alpha_shape_vertex, vertex)); + } else { + // alpha shape found + Simplex_tree_vertex vertex = the_map_iterator->second; +#ifdef DEBUG_TRACES + std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl; +#endif // DEBUG_TRACES + the_simplex_tree.push_back(vertex); + } + } + // Construction of the simplex_tree + Filtration_value filtr = /*std::sqrt*/(*the_alpha_value_iterator); +#ifdef DEBUG_TRACES + std::cout << "filtration = " << filtr << std::endl; +#endif // DEBUG_TRACES + if (filtr > filtration_max) { + filtration_max = filtr; + } + simplex_tree.insert_simplex(the_simplex_tree, filtr); + if (the_alpha_value_iterator != the_alpha_values.end()) + ++the_alpha_value_iterator; + else + std::cout << "This shall not happen" << std::endl; + } + simplex_tree.set_filtration(filtration_max); + simplex_tree.set_dimension(dim_max); + +#ifdef DEBUG_TRACES + std::cout << "vertices \t\t" << count_vertices << std::endl; + std::cout << "edges \t\t" << count_edges << std::endl; + std::cout << "facets \t\t" << count_facets << std::endl; + std::cout << "cells \t\t" << count_cells << std::endl; + + + std::cout << "Information of the Simplex Tree: " << std::endl; + std::cout << " Number of vertices = " << simplex_tree.num_vertices() << " "; + std::cout << " Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl; + std::cout << " Dimension = " << simplex_tree.dimension() << " "; + std::cout << " filtration = " << simplex_tree.filtration() << std::endl << std::endl; +#endif // DEBUG_TRACES + +#ifdef DEBUG_TRACES + std::cout << "Iterator on vertices: " << std::endl; + for (auto vertex : simplex_tree.complex_vertex_range()) { + std::cout << vertex << " "; + } +#endif // DEBUG_TRACES + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl; + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(simplex_tree, true); + // initializes the coefficient field for homology + pcoh.init_coefficients(coeff_field_characteristic); + + pcoh.compute_persistent_cohomology(min_persistence); + + pcoh.output_diagram(); + + return 0; +} |