From 73efba37dfcab3f9b8589cb3b0b80f6d7dd7bd0b Mon Sep 17 00:00:00 2001 From: cjamin Date: Wed, 4 Oct 2017 10:24:04 +0000 Subject: Document utilies in README files using Markdown + move/rename some utilities Only missing doc for now: Garland_heckbert git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/add_utils_in_gudhi_v2@2752 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 7ba0113762b6622130835dc9e372acfae29c2db8 --- src/Alpha_complex/utilities/CMakeLists.txt | 37 +++ src/Alpha_complex/utilities/README | 131 ++++++++++ .../utilities/alpha_complex_3d_helper.h | 76 ++++++ .../utilities/alpha_complex_3d_persistence.cpp | 242 +++++++++++++++++++ .../utilities/alpha_complex_persistence.cpp | 125 ++++++++++ .../periodic_alpha_complex_3d_persistence.cpp | 268 +++++++++++++++++++++ 6 files changed, 879 insertions(+) create mode 100644 src/Alpha_complex/utilities/CMakeLists.txt create mode 100644 src/Alpha_complex/utilities/README create mode 100644 src/Alpha_complex/utilities/alpha_complex_3d_helper.h create mode 100644 src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp create mode 100644 src/Alpha_complex/utilities/alpha_complex_persistence.cpp create mode 100644 src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp (limited to 'src/Alpha_complex') diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt new file mode 100644 index 00000000..69dbc28a --- /dev/null +++ b/src/Alpha_complex/utilities/CMakeLists.txt @@ -0,0 +1,37 @@ +cmake_minimum_required(VERSION 2.6) +project(Alpha_complex_utilities) + +if(CGAL_FOUND) + add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp) + target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY}) + + if (TBB_FOUND) + target_link_libraries(alpha_complex_3d_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) + add_test(NAME Alpha_complex_utilities_alpha_complex_3d_persistence COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45") + + install(TARGETS alpha_complex_3d_persistence DESTINATION bin) + + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) + add_executable (alpha_complex_persistence alpha_complex_persistence.cpp) + target_link_libraries(alpha_complex_persistence + ${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 ${CGAL_LIBRARY}) + + if (TBB_FOUND) + target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES}) + target_link_libraries(periodic_alpha_complex_3d_persistence ${TBB_LIBRARIES}) + endif(TBB_FOUND) + add_test(NAME Alpha_complex_utilities_alpha_complex_persistence COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45") + add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $ + "${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") + + install(TARGETS alpha_complex_persistence DESTINATION bin) + install(TARGETS periodic_alpha_complex_3d_persistence DESTINATION bin) + + endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) +endif(CGAL_FOUND) diff --git a/src/Alpha_complex/utilities/README b/src/Alpha_complex/utilities/README new file mode 100644 index 00000000..30e1b187 --- /dev/null +++ b/src/Alpha_complex/utilities/README @@ -0,0 +1,131 @@ +# Alpha_complex # + +## `alpha_complex_3d_persistence` ## +This program computes the persistent homology with coefficient field Z/pZ of the 3D alpha complex built from a 3D point cloud. The output diagram contains one bar per line, written with the convention: + +`p dim b d` + +where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p = p1*...*pr` is the product of prime numbers *pi* such that the homology feature exists in homology with *Z/piZ* coefficients). + +**Usage** +`alpha_complex_3d_persistence

` +where +`` is the path to the input point cloud in OFF format. +`

` is the characteristic p of the coefficient field *Z/pZ* for computing homology. It must be a stricly positive integer. +`` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. It must be a floating-point number >= -1. + +**Example** +`alpha_complex_3d_persistence ../../data/points/tore3D_300.off 2 0.45` + +outputs: +``` +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 +``` + +Here we retrieve expected Betti numbers on a tore 3D: +``` +Betti numbers[0] = 1 +Betti numbers[1] = 2 +Betti numbers[2] = 1 +``` + +N.B.: +* `alpha_complex_3d_persistence` only accepts OFF files in dimension 3. +* Filtration values are alpha square values. + + + + +## `periodic_alpha_complex_3d_persistence` ## +This program computes the persistent homology with coefficient field Z/pZ of the 3D periodic alpha complex built from a 3D point cloud. The output diagram contains one bar per line, written with the convention: + +`p dim b d` + +where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p = p1*...*pr` is the product of prime numbers *pi* such that the homology feature exists in homology with *Z/piZ* coefficients). + +**Usage** +`periodic_alpha_complex_3d_persistence

` +where +`` is the path to the input point cloud in OFF format. +`

` is the characteristic p of the coefficient field *Z/pZ* for computing homology. It must be a stricly positive integer. +`` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. It must be a floating-point number >= -1. + +**Example** +`periodic_alpha_complex_3d_persistence ../../data/points/grid_10_10_10_in_0_1.off ../../data/points/iso_cuboid_3_in_0_1.txt 3 1.0` + +outputs: +``` +Periodic Delaunay computed. +Simplex_tree dim: 3 +3 0 0 inf +3 1 0.0025 inf +3 1 0.0025 inf +3 1 0.0025 inf +3 2 0.005 inf +3 2 0.005 inf +3 2 0.005 inf +3 3 0.0075 inf +``` + +Here we retrieve expected Betti numbers on an 3D iso-oriented cuboids: +``` +Betti numbers[0] = 1 +Betti numbers[1] = 3 +Betti numbers[2] = 3 +Betti numbers[3] = 1 +``` + +N.B.: +* `periodic_alpha_complex_3d_persistence` only accepts OFF files in dimension 3. +* In this example, the periodic cube is hard coded to { x = [0,1]; y = [0,1]; z = [0,1] } +* Filtration values are alpha square values. + + + + + +## `alpha_complex_persistence` ## +This program computes the persistent homology with coefficient field Z/pZ of the dD alpha complex built from a dD point cloud. The output diagram contains one bar per line, written with the convention: + +`p dim b d` + +where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p = p1*...*pr` is the product of prime numbers pi such that the homology feature exists in homology with Z/piZ coefficients). + +**Usage** +`alpha_complex_persistence [options] ` + +**Allowed options** + +* `-h [ --help ]` Produce help message +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. By default, print in std::cout. +* `-r [ --max-alpha-square-value ]` (default = inf) Maximal alpha square value for the Alpha complex construction. +* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. + +**Example** +`alpha_complex_persistence -r 32 -p 2 -m 0.45 ../../data/points/tore3D_300.off` + +outputs: +``` +Alpha complex is of dimension 3 - 9273 simplices - 300 vertices. +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 +``` + +Here we retrieve expected Betti numbers on a tore 3D: +``` +Betti numbers[0] = 1 +Betti numbers[1] = 2 +Betti numbers[2] = 1 +``` + +N.B.: +* `alpha_complex_persistence` only accepts OFF files in dimension d. +* Filtration values are alpha square values. diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_helper.h b/src/Alpha_complex/utilities/alpha_complex_3d_helper.h new file mode 100644 index 00000000..7865e4ec --- /dev/null +++ b/src/Alpha_complex/utilities/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 . + */ + +#ifndef ALPHA_COMPLEX_3D_HELPER_H_ +#define ALPHA_COMPLEX_3D_HELPER_H_ + +template +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 +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 +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 +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 // ALPHA_COMPLEX_3D_HELPER_H_ diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..f63ff0f6 --- /dev/null +++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp @@ -0,0 +1,242 @@ +/* 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 . + */ + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "alpha_complex_3d_helper.h" + +// Alpha_shape_3 templates type definitions +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using Vb = CGAL::Alpha_shape_vertex_base_3; +using Fb = CGAL::Alpha_shape_cell_base_3; +using Tds = CGAL::Triangulation_data_structure_3; +using Triangulation_3 = CGAL::Delaunay_triangulation_3; +using Alpha_shape_3 = CGAL::Alpha_shape_3; + +// 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, + CGAL::cpp11::tuple >, + std::back_insert_iterator< std::vector > > >; +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; + +// gudhi type definition +using ST = Gudhi::Simplex_tree; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map; +using Alpha_shape_simplex_tree_pair = std::pair; +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:\n" << progName << " path_to_OFF_file coeff_field_characteristic[integer " << + "> 0] min_persistence[float >= -1.0]\n"; + std::cerr << " path_to_OFF_file is the path to your points cloud in OFF format.\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 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 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 the_objects; + std::vector the_alpha_values; + + Dispatch disp = CGAL::dispatch_output(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::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(&object_iterator)) { + vertex_list = from_cell(*cell); + count_cells++; + if (dim_max < 3) { + // Cell is of dim 3 + dim_max = 3; + } + } else if (const Facet * facet = CGAL::object_cast(&object_iterator)) { + vertex_list = from_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(&object_iterator)) { + vertex_list = from_edge(*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(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex(*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_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() << " "; +#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/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp new file mode 100644 index 00000000..9e84e91f --- /dev/null +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -0,0 +1,125 @@ +#include + +#include + +#include +#include +// to construct a simplex_tree from alpha complex +#include + +#include +#include +#include // for numeric_limits + +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & output_file_diag + , Filtration_value & alpha_square_max_value + , int & coeff_field_characteristic + , Filtration_value & min_persistence); + +int main(int argc, char **argv) { + std::string off_file_points; + std::string output_file_diag; + Filtration_value alpha_square_max_value; + int coeff_field_characteristic; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, output_file_diag, alpha_square_max_value, + coeff_field_characteristic, min_persistence); + + // ---------------------------------------------------------------------------- + // Init of an alpha complex from an OFF file + // ---------------------------------------------------------------------------- + using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >; + Gudhi::alpha_complex::Alpha_complex alpha_complex_from_file(off_file_points); + + Simplex_tree simplex; + if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) { + // ---------------------------------------------------------------------------- + // Display information about the alpha complex + // ---------------------------------------------------------------------------- + std::cout << "Simplicial complex is of dimension " << simplex.dimension() << + " - " << simplex.num_simplices() << " simplices - " << + simplex.num_vertices() << " vertices." << std::endl; + + // Sort the simplices in the order of the filtration + simplex.initialize_filtration(); + + std::cout << "Simplex_tree dim: " << simplex.dimension() << std::endl; + // Compute the persistence diagram of the complex + Gudhi::persistent_cohomology::Persistent_cohomology< Simplex_tree, + Gudhi::persistent_cohomology::Field_Zp > pcoh(simplex); + // initializes the coefficient field for homology + pcoh.init_coefficients(coeff_field_characteristic); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (output_file_diag.empty()) { + pcoh.output_diagram(); + } else { + std::cout << "Result in file: " << output_file_diag << std::endl; + std::ofstream out(output_file_diag); + pcoh.output_diagram(out); + out.close(); + } + } + + return 0; +} + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , std::string & output_file_diag + , Filtration_value & alpha_square_max_value + , int & coeff_field_characteristic + , 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), + "Name of file containing a point set. Format is one point per line: X1 ... Xd "); + + po::options_description visible("Allowed options", 100); + visible.add_options() + ("help,h", "produce help message") + ("output-file,o", po::value(&output_file_diag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout") + ("max-alpha-square-value,r", + po::value(&alpha_square_max_value)->default_value(std::numeric_limits::infinity()), + "Maximal alpha square value for the Alpha complex construction.") + ("field-charac,p", po::value(&coeff_field_characteristic)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.") + ("min-persistence,m", po::value(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map 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")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of an Alpha complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..8140a3c5 --- /dev/null +++ b/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp @@ -0,0 +1,268 @@ +/* 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 . + */ + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "alpha_complex_3d_helper.h" + +// Traits +using K = CGAL::Exact_predicates_inexact_constructions_kernel; +using PK = CGAL::Periodic_3_Delaunay_triangulation_traits_3; +// Vertex type +using DsVb = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>; +using Vb = CGAL::Triangulation_vertex_base_3; +using AsVb = CGAL::Alpha_shape_vertex_base_3; +// Cell type +using DsCb = CGAL::Periodic_3_triangulation_ds_cell_base_3<>; +using Cb = CGAL::Triangulation_cell_base_3; +using AsCb = CGAL::Alpha_shape_cell_base_3; +using Tds = CGAL::Triangulation_data_structure_3; +using P3DT3 = CGAL::Periodic_3_Delaunay_triangulation_3; +using Alpha_shape_3 = CGAL::Alpha_shape_3; +using Point_3 = PK::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, + CGAL::cpp11::tuple >, + std::back_insert_iterator< std::vector > > >; +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; + +// gudhi type definition +using ST = Gudhi::Simplex_tree; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map; +using Alpha_shape_simplex_tree_pair = std::pair; +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:\n" << progName << " path_to_OFF_file path_to_iso_cuboid_3_file coeff_field_characteristic[" << + "integer > 0] min_persistence[float >= -1.0]\n" << + " path_to_OFF_file is the path to your points cloud in OFF format.\n" << + " path_to_iso_cuboid_3_file is the path to the iso cuboid file with the following format :\n" << + " x_min y_min z_min x_max y_max z_max\n" << + " In this example, the periodic cube will be " << + "{ x = [x_min,x_max]; y = [y_min,y_max]; z = [z_min,z_max] }.\n" << + " For more information, please refer to\n" << + " https://doc.cgal.org/latest/Kernel_23/classCGAL_1_1Iso__cuboid__3.html\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 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]); + } + + // 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 lp = off_reader.get_point_cloud(); + + // Define the periodic cube + 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 + if (pdt.is_triangulation_in_1_sheet()) pdt.convert_to_1_sheeted_covering(); + std::cout << "Periodic Delaunay computed." << std::endl; + + // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode + // Maybe need to set it to GENERAL mode + Alpha_shape_3 as(pdt, 0, Alpha_shape_3::GENERAL); + + // filtration with alpha values from alpha shape + std::vector the_objects; + std::vector the_alpha_values; + + Dispatch disp = CGAL::dispatch_output(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::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(&object_iterator)) { + vertex_list = from_cell(*cell); + count_cells++; + if (dim_max < 3) { + // Cell is of dim 3 + dim_max = 3; + } + } else if (const Facet * facet = CGAL::object_cast(&object_iterator)) { + vertex_list = from_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(&object_iterator)) { + vertex_list = from_edge(*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(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex(*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_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() << " "; +#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; +} -- cgit v1.2.3 From 8428e315efa22a1b40eb996c98eff671f91a9b27 Mon Sep 17 00:00:00 2001 From: cjamin Date: Wed, 15 Nov 2017 14:20:30 +0000 Subject: Fixes according to Marc's review git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/add_utils_in_gudhi_v2@2881 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 84870193cb05cabd8c756a3bbc49083e4d739525 --- src/Alpha_complex/utilities/README | 10 +- src/Bottleneck_distance/example/CMakeLists.txt | 12 +- src/Bottleneck_distance/example/README | 19 +++ .../alpha_rips_persistence_bottleneck_distance.cpp | 190 +++++++++++++++++++++ .../example/bottleneck_read_file_example.cpp | 50 ------ src/Bottleneck_distance/utilities/CMakeLists.txt | 11 +- src/Bottleneck_distance/utilities/README | 19 +-- .../alpha_rips_persistence_bottleneck_distance.cpp | 190 --------------------- .../utilities/bottleneck_read_file_example.cpp | 50 ++++++ src/Rips_complex/utilities/README | 2 +- 10 files changed, 281 insertions(+), 272 deletions(-) create mode 100644 src/Bottleneck_distance/example/README create mode 100644 src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp delete mode 100644 src/Bottleneck_distance/example/bottleneck_read_file_example.cpp delete mode 100644 src/Bottleneck_distance/utilities/alpha_rips_persistence_bottleneck_distance.cpp create mode 100644 src/Bottleneck_distance/utilities/bottleneck_read_file_example.cpp (limited to 'src/Alpha_complex') diff --git a/src/Alpha_complex/utilities/README b/src/Alpha_complex/utilities/README index 30e1b187..abb17bf7 100644 --- a/src/Alpha_complex/utilities/README +++ b/src/Alpha_complex/utilities/README @@ -5,14 +5,14 @@ This program computes the persistent homology with coefficient field Z/pZ of the `p dim b d` -where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p = p1*...*pr` is the product of prime numbers *pi* such that the homology feature exists in homology with *Z/piZ* coefficients). +where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number). **Usage** `alpha_complex_3d_persistence

` where `` is the path to the input point cloud in OFF format. `

` is the characteristic p of the coefficient field *Z/pZ* for computing homology. It must be a stricly positive integer. -`` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. It must be a floating-point number >= -1. +`` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. **Example** `alpha_complex_3d_persistence ../../data/points/tore3D_300.off 2 0.45` @@ -45,14 +45,14 @@ This program computes the persistent homology with coefficient field Z/pZ of the `p dim b d` -where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p = p1*...*pr` is the product of prime numbers *pi* such that the homology feature exists in homology with *Z/piZ* coefficients). +where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number). **Usage** `periodic_alpha_complex_3d_persistence

` where `` is the path to the input point cloud in OFF format. `

` is the characteristic p of the coefficient field *Z/pZ* for computing homology. It must be a stricly positive integer. -`` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. It must be a floating-point number >= -1. +`` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. **Example** `periodic_alpha_complex_3d_persistence ../../data/points/grid_10_10_10_in_0_1.off ../../data/points/iso_cuboid_3_in_0_1.txt 3 1.0` @@ -93,7 +93,7 @@ This program computes the persistent homology with coefficient field Z/pZ of the `p dim b d` -where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p = p1*...*pr` is the product of prime numbers pi such that the homology feature exists in homology with Z/piZ coefficients). +where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number). **Usage** `alpha_complex_persistence [options] ` diff --git a/src/Bottleneck_distance/example/CMakeLists.txt b/src/Bottleneck_distance/example/CMakeLists.txt index 6a602dbb..9677f5c5 100644 --- a/src/Bottleneck_distance/example/CMakeLists.txt +++ b/src/Bottleneck_distance/example/CMakeLists.txt @@ -3,15 +3,15 @@ project(Bottleneck_distance_examples) if (NOT CGAL_VERSION VERSION_LESS 4.8.1) add_executable (bottleneck_basic_example bottleneck_basic_example.cpp) - add_executable (bottleneck_read_file_example bottleneck_read_file_example.cpp) + add_executable (alpha_rips_persistence_bottleneck_distance alpha_rips_persistence_bottleneck_distance.cpp) + target_link_libraries(alpha_rips_persistence_bottleneck_distance ${Boost_PROGRAM_OPTIONS_LIBRARY}) add_test(NAME Bottleneck_distance_example_basic COMMAND $) + add_test(NAME Bottleneck_distance_example_alpha_rips_persistence_bottleneck + COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.15" "-m" "0.12" "-d" "3" "-p" "3") - add_test(NAME Bottleneck_read_file_example - COMMAND $ - "${CMAKE_SOURCE_DIR}/data/persistence_diagram/first.pers" "${CMAKE_SOURCE_DIR}/data/persistence_diagram/second.pers") - - install(TARGETS bottleneck_read_file_example DESTINATION bin) install(TARGETS bottleneck_basic_example DESTINATION bin) + install(TARGETS alpha_rips_persistence_bottleneck_distance DESTINATION bin) endif (NOT CGAL_VERSION VERSION_LESS 4.8.1) diff --git a/src/Bottleneck_distance/example/README b/src/Bottleneck_distance/example/README new file mode 100644 index 00000000..0e314608 --- /dev/null +++ b/src/Bottleneck_distance/example/README @@ -0,0 +1,19 @@ +# Bottleneck_distance # + +## `alpha_rips_persistence_bottleneck_distance` ## +This program computes the persistent homology with coefficient field Z/pZ of a Rips complex defined on a set of input points. The output diagram contains one bar per line, written with the convention: + +`p dim b d` + +where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients. + +Usage: +`alpha_rips_persistence_bottleneck_distance [options] ` + +Allowed options: + +* `-h [ --help ]` Produce help message +* `-r [ --max-edge-length ]` (default = inf) Maximal length of an edge for the Rips complex construction.` +* `-d [ --cpx-dimension ]` (default = 1) Maximal dimension of the Rips complex we want to compute.` +* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. diff --git a/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp b/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp new file mode 100644 index 00000000..fd164b22 --- /dev/null +++ b/src/Bottleneck_distance/example/alpha_rips_persistence_bottleneck_distance.cpp @@ -0,0 +1,190 @@ +/* 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) 2017 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 + +#include +#include +#include // infinity +#include // for pair +#include // for transform + + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; +using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >; +using Point_d = Kernel::Point_d; +using Points_off_reader = Gudhi::Points_off_reader; + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , Filtration_value & threshold + , int & dim_max + , int & p + , Filtration_value & min_persistence); + +static inline std::pair compute_root_square(std::pair input) { + return std::make_pair(std::sqrt(input.first), std::sqrt(input.second)); +} + +int main(int argc, char * argv[]) { + std::string off_file_points; + Filtration_value threshold; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, threshold, dim_max, p, min_persistence); + + Points_off_reader off_reader(off_file_points); + + // -------------------------------------------- + // Rips persistence + // -------------------------------------------- + Rips_complex rips_complex(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree rips_stree; + + rips_complex.create_complex(rips_stree, dim_max); + std::cout << "The Rips complex contains " << rips_stree.num_simplices() << " simplices and has dimension " + << rips_stree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + rips_stree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology rips_pcoh(rips_stree); + // initializes the coefficient field for homology + rips_pcoh.init_coefficients(p); + rips_pcoh.compute_persistent_cohomology(min_persistence); + + // rips_pcoh.output_diagram(); + + // -------------------------------------------- + // Alpha persistence + // -------------------------------------------- + Gudhi::alpha_complex::Alpha_complex alpha_complex(off_reader.get_point_cloud()); + + Simplex_tree alpha_stree; + alpha_complex.create_complex(alpha_stree, threshold * threshold); + std::cout << "The Alpha complex contains " << alpha_stree.num_simplices() << " simplices and has dimension " + << alpha_stree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + alpha_stree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology alpha_pcoh(alpha_stree); + // initializes the coefficient field for homology + alpha_pcoh.init_coefficients(p); + alpha_pcoh.compute_persistent_cohomology(min_persistence * min_persistence); + + // alpha_pcoh.output_diagram(); + + // -------------------------------------------- + // Bottleneck distance between both persistence + // -------------------------------------------- + double max_b_distance {}; + for (int dim = 0; dim < dim_max; dim ++) { + std::vector< std::pair< Filtration_value , Filtration_value > > rips_intervals; + std::vector< std::pair< Filtration_value , Filtration_value > > alpha_intervals; + rips_intervals = rips_pcoh.intervals_in_dimension(dim); + alpha_intervals = alpha_pcoh.intervals_in_dimension(dim); + std::transform(alpha_intervals.begin(), alpha_intervals.end(), alpha_intervals.begin(), compute_root_square); + + double bottleneck_distance = Gudhi::persistence_diagram::bottleneck_distance(rips_intervals, alpha_intervals); + std::cout << "In dimension " << dim << ", bottleneck distance = " << bottleneck_distance << std::endl; + if (bottleneck_distance > max_b_distance) + max_b_distance = bottleneck_distance; + } + std::cout << "================================================================================" << std::endl; + std::cout << "Bottleneck distance is " << max_b_distance << std::endl; + + return 0; +} + +void program_options(int argc, char * argv[] + , std::string & off_file_points + , Filtration_value & threshold + , 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), + "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-edge-length,r", + po::value(&threshold)->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.") + ("field-charac,p", po::value(&p)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.") + ("min-persistence,m", po::value(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map 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")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a Rips complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp b/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp deleted file mode 100644 index 24d73c57..00000000 --- a/src/Bottleneck_distance/example/bottleneck_read_file_example.cpp +++ /dev/null @@ -1,50 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Authors: Francois Godi, small modifications by Pawel Dlotko - * - * 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 . - */ - -#include -#include -#include -#include -#include // for pair -#include -#include // for numeric_limits - -int main(int argc, char** argv) { - if (argc < 3) { - std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file" << - " should contain a birth-death pair per line. Third, optional parameter is an error bound on a bottleneck" << - " distance (set by default to the smallest positive double value). If you set the error bound to 0, be" << - " aware this version is exact but expensive. The program will now terminate \n"; - return -1; - } - std::vector> diag1 = Gudhi::read_persistence_intervals_in_dimension(argv[1]); - std::vector> diag2 = Gudhi::read_persistence_intervals_in_dimension(argv[2]); - - double tolerance = std::numeric_limits::min(); - if (argc == 4) { - tolerance = atof(argv[3]); - } - double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance); - std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl; - - return 0; -} diff --git a/src/Bottleneck_distance/utilities/CMakeLists.txt b/src/Bottleneck_distance/utilities/CMakeLists.txt index 634b2479..decdc789 100644 --- a/src/Bottleneck_distance/utilities/CMakeLists.txt +++ b/src/Bottleneck_distance/utilities/CMakeLists.txt @@ -2,16 +2,15 @@ cmake_minimum_required(VERSION 2.6) project(Bottleneck_distance_utilities) if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) - add_executable (alpha_rips_persistence_bottleneck_distance alpha_rips_persistence_bottleneck_distance.cpp) - target_link_libraries(alpha_rips_persistence_bottleneck_distance ${Boost_PROGRAM_OPTIONS_LIBRARY}) + add_executable (bottleneck_read_file_example bottleneck_read_file_example.cpp) if (TBB_FOUND) target_link_libraries(alpha_rips_persistence_bottleneck_distance ${TBB_LIBRARIES}) endif(TBB_FOUND) - add_test(NAME Bottleneck_distance_example_alpha_rips_persistence_bottleneck - COMMAND $ - "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.15" "-m" "0.12" "-d" "3" "-p" "3") + add_test(NAME Bottleneck_read_file_example + COMMAND $ + "${CMAKE_SOURCE_DIR}/data/persistence_diagram/first.pers" "${CMAKE_SOURCE_DIR}/data/persistence_diagram/second.pers") - install(TARGETS alpha_rips_persistence_bottleneck_distance DESTINATION bin) + install(TARGETS bottleneck_read_file_example DESTINATION bin) endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) diff --git a/src/Bottleneck_distance/utilities/README b/src/Bottleneck_distance/utilities/README index 4c9b0735..cfdccf0b 100644 --- a/src/Bottleneck_distance/utilities/README +++ b/src/Bottleneck_distance/utilities/README @@ -1,19 +1,10 @@ # Bottleneck_distance # -## `alpha_rips_persistence_bottleneck_distance` ## -This program computes the persistent homology with coefficient field Z/pZ of a Rips complex defined on a set of input points. The output diagram contains one bar per line, written with the convention: - -`p dim b d` - -where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients. +## `bottleneck_read_file_example` ## +This program computes the Bottleneck distance between two persistence diagram files. Usage: -`alpha_rips_persistence_bottleneck_distance [options] ` - -Allowed options: +`bottleneck_read_file_example []` -* `-h [ --help ]` Produce help message -* `-r [ --max-edge-length ]` (default = inf) Maximal length of an edge for the Rips complex construction.` -* `-d [ --cpx-dimension ]` (default = 1) Maximal dimension of the Rips complex we want to compute.` -* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. -* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. \ No newline at end of file + and must be in the format described [here](http://gudhi.gforge.inria.fr/doc/latest/fileformats.html#FileFormatsPers). + is an error bound on the bottleneck distance (set by default to the smallest positive double value). diff --git a/src/Bottleneck_distance/utilities/alpha_rips_persistence_bottleneck_distance.cpp b/src/Bottleneck_distance/utilities/alpha_rips_persistence_bottleneck_distance.cpp deleted file mode 100644 index fd164b22..00000000 --- a/src/Bottleneck_distance/utilities/alpha_rips_persistence_bottleneck_distance.cpp +++ /dev/null @@ -1,190 +0,0 @@ -/* 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) 2017 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 - -#include -#include -#include // infinity -#include // for pair -#include // for transform - - -// Types definition -using Simplex_tree = Gudhi::Simplex_tree; -using Filtration_value = Simplex_tree::Filtration_value; -using Rips_complex = Gudhi::rips_complex::Rips_complex; -using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; -using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; -using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >; -using Point_d = Kernel::Point_d; -using Points_off_reader = Gudhi::Points_off_reader; - -void program_options(int argc, char * argv[] - , std::string & off_file_points - , Filtration_value & threshold - , int & dim_max - , int & p - , Filtration_value & min_persistence); - -static inline std::pair compute_root_square(std::pair input) { - return std::make_pair(std::sqrt(input.first), std::sqrt(input.second)); -} - -int main(int argc, char * argv[]) { - std::string off_file_points; - Filtration_value threshold; - int dim_max; - int p; - Filtration_value min_persistence; - - program_options(argc, argv, off_file_points, threshold, dim_max, p, min_persistence); - - Points_off_reader off_reader(off_file_points); - - // -------------------------------------------- - // Rips persistence - // -------------------------------------------- - Rips_complex rips_complex(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); - - // Construct the Rips complex in a Simplex Tree - Simplex_tree rips_stree; - - rips_complex.create_complex(rips_stree, dim_max); - std::cout << "The Rips complex contains " << rips_stree.num_simplices() << " simplices and has dimension " - << rips_stree.dimension() << " \n"; - - // Sort the simplices in the order of the filtration - rips_stree.initialize_filtration(); - - // Compute the persistence diagram of the complex - Persistent_cohomology rips_pcoh(rips_stree); - // initializes the coefficient field for homology - rips_pcoh.init_coefficients(p); - rips_pcoh.compute_persistent_cohomology(min_persistence); - - // rips_pcoh.output_diagram(); - - // -------------------------------------------- - // Alpha persistence - // -------------------------------------------- - Gudhi::alpha_complex::Alpha_complex alpha_complex(off_reader.get_point_cloud()); - - Simplex_tree alpha_stree; - alpha_complex.create_complex(alpha_stree, threshold * threshold); - std::cout << "The Alpha complex contains " << alpha_stree.num_simplices() << " simplices and has dimension " - << alpha_stree.dimension() << " \n"; - - // Sort the simplices in the order of the filtration - alpha_stree.initialize_filtration(); - - // Compute the persistence diagram of the complex - Persistent_cohomology alpha_pcoh(alpha_stree); - // initializes the coefficient field for homology - alpha_pcoh.init_coefficients(p); - alpha_pcoh.compute_persistent_cohomology(min_persistence * min_persistence); - - // alpha_pcoh.output_diagram(); - - // -------------------------------------------- - // Bottleneck distance between both persistence - // -------------------------------------------- - double max_b_distance {}; - for (int dim = 0; dim < dim_max; dim ++) { - std::vector< std::pair< Filtration_value , Filtration_value > > rips_intervals; - std::vector< std::pair< Filtration_value , Filtration_value > > alpha_intervals; - rips_intervals = rips_pcoh.intervals_in_dimension(dim); - alpha_intervals = alpha_pcoh.intervals_in_dimension(dim); - std::transform(alpha_intervals.begin(), alpha_intervals.end(), alpha_intervals.begin(), compute_root_square); - - double bottleneck_distance = Gudhi::persistence_diagram::bottleneck_distance(rips_intervals, alpha_intervals); - std::cout << "In dimension " << dim << ", bottleneck distance = " << bottleneck_distance << std::endl; - if (bottleneck_distance > max_b_distance) - max_b_distance = bottleneck_distance; - } - std::cout << "================================================================================" << std::endl; - std::cout << "Bottleneck distance is " << max_b_distance << std::endl; - - return 0; -} - -void program_options(int argc, char * argv[] - , std::string & off_file_points - , Filtration_value & threshold - , 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), - "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-edge-length,r", - po::value(&threshold)->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.") - ("field-charac,p", po::value(&p)->default_value(11), - "Characteristic p of the coefficient field Z/pZ for computing homology.") - ("min-persistence,m", po::value(&min_persistence), - "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); - - po::positional_options_description pos; - pos.add("input-file", 1); - - po::options_description all; - all.add(visible).add(hidden); - - po::variables_map 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")) { - std::cout << std::endl; - std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; - std::cout << "of a Rips complex defined on a set of input points.\n \n"; - std::cout << "The output diagram contains one bar per line, written with the convention: \n"; - std::cout << " p dim b d \n"; - std::cout << "where dim is the dimension of the homological feature,\n"; - std::cout << "b and d are respectively the birth and death of the feature and \n"; - std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; - - std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; - std::cout << visible << std::endl; - std::abort(); - } -} diff --git a/src/Bottleneck_distance/utilities/bottleneck_read_file_example.cpp b/src/Bottleneck_distance/utilities/bottleneck_read_file_example.cpp new file mode 100644 index 00000000..9dd52b31 --- /dev/null +++ b/src/Bottleneck_distance/utilities/bottleneck_read_file_example.cpp @@ -0,0 +1,50 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Authors: Francois Godi, small modifications by Pawel Dlotko + * + * 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 . + */ + +#include +#include +#include +#include +#include // for pair +#include +#include // for numeric_limits + +int main(int argc, char** argv) { + if (argc < 3) { + std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file" << + " should contain a birth-death pair per line. Third, optional parameter is an error bound on the bottleneck" << + " distance (set by default to the smallest positive double value). If you set the error bound to 0, be" << + " aware this version is exact but expensive. The program will now terminate \n"; + return -1; + } + std::vector> diag1 = Gudhi::read_persistence_intervals_in_dimension(argv[1]); + std::vector> diag2 = Gudhi::read_persistence_intervals_in_dimension(argv[2]); + + double tolerance = std::numeric_limits::min(); + if (argc == 4) { + tolerance = atof(argv[3]); + } + double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance); + std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl; + + return 0; +} diff --git a/src/Rips_complex/utilities/README b/src/Rips_complex/utilities/README index ddb7860f..93ec3658 100644 --- a/src/Rips_complex/utilities/README +++ b/src/Rips_complex/utilities/README @@ -5,7 +5,7 @@ This program computes the persistent homology with coefficient field *Z/pZ* of a `p dim b d` -where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p = p1*...*pr` is the product of prime numbers *pi* such that the homology feature exists in homology with *Z/piZ* coefficients). +where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number). **Usage** `rips_persistence [options] ` -- cgit v1.2.3 From ae41883215cd89e2ca71fd7e2a46b66d5075dc2d Mon Sep 17 00:00:00 2001 From: cjamin Date: Wed, 15 Nov 2017 14:39:41 +0000 Subject: Add iso-cuboid file format git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/add_utils_in_gudhi_v2@2882 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: ad82867ef7523c4a03fdb89054e5bd7b8d40f7da --- src/Alpha_complex/utilities/README | 3 ++- src/Bottleneck_distance/utilities/README | 6 +++--- src/common/doc/file_formats.h | 15 +++++++++++++++ 3 files changed, 20 insertions(+), 4 deletions(-) (limited to 'src/Alpha_complex') diff --git a/src/Alpha_complex/utilities/README b/src/Alpha_complex/utilities/README index abb17bf7..67ea1875 100644 --- a/src/Alpha_complex/utilities/README +++ b/src/Alpha_complex/utilities/README @@ -48,9 +48,10 @@ This program computes the persistent homology with coefficient field Z/pZ of the where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number). **Usage** -`periodic_alpha_complex_3d_persistence

` +`periodic_alpha_complex_3d_persistence

` where `` is the path to the input point cloud in OFF format. +`` is the path to the file describing the periodic domain. It must be in the format described [here](http://gudhi.gforge.inria.fr/doc/latest/fileformats.html#FileFormatsIsoCuboid). `

` is the characteristic p of the coefficient field *Z/pZ* for computing homology. It must be a stricly positive integer. `` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. diff --git a/src/Bottleneck_distance/utilities/README b/src/Bottleneck_distance/utilities/README index cfdccf0b..d9fdd252 100644 --- a/src/Bottleneck_distance/utilities/README +++ b/src/Bottleneck_distance/utilities/README @@ -5,6 +5,6 @@ This program computes the Bottleneck distance between two persistence diagram fi Usage: `bottleneck_read_file_example []` - - and must be in the format described [here](http://gudhi.gforge.inria.fr/doc/latest/fileformats.html#FileFormatsPers). - is an error bound on the bottleneck distance (set by default to the smallest positive double value). +where +`` and `` must be in the format described [here](http://gudhi.gforge.inria.fr/doc/latest/fileformats.html#FileFormatsPers). +`` is an error bound on the bottleneck distance (set by default to the smallest positive double value). diff --git a/src/common/doc/file_formats.h b/src/common/doc/file_formats.h index d715aa4d..3b3d9248 100644 --- a/src/common/doc/file_formats.h +++ b/src/common/doc/file_formats.h @@ -53,6 +53,21 @@ namespace Gudhi { Such files can be generated with `Gudhi::persistent_cohomology::Persistent_cohomology::output_diagram()` and read with `Gudhi::read_persistence_intervals_and_dimension()`, `Gudhi::read_persistence_intervals_grouped_by_dimension()` or `Gudhi::read_persistence_intervals_in_dimension()`. + + + \section FileFormatsIsoCuboid Iso-cuboid + + Such a file describes an iso-oriented cuboid with diagonal opposite vertices (min_hx, min_hy, min_hz,...) and (max_hx, max_hy, max_hz, ...). The format is:
+ \verbatim + min_hx min_hy [min_hz ...] + max_hx max_hy [max_hz ...] + \endverbatim + + Here is a simple sample file in the 3D case: + \verbatim + -1. -1. -1. + 1. 1. 1. + \endverbatim */ } // namespace Gudhi -- cgit v1.2.3 From 048108b77120bf7b607a5093977b1da294bc472b Mon Sep 17 00:00:00 2001 From: cjamin Date: Wed, 15 Nov 2017 14:50:29 +0000 Subject: Move 2 more examples into utilities git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/add_utils_in_gudhi_v2@2883 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: b9fd71a4833c84640be2a0f792f1eaea5109874d --- src/Alpha_complex/utilities/CMakeLists.txt | 13 + .../exact_alpha_complex_3d_persistence.cpp | 244 +++++++++++++++++++ .../weighted_alpha_complex_3d_persistence.cpp | 267 +++++++++++++++++++++ src/Persistent_cohomology/example/CMakeLists.txt | 17 -- .../example/exact_alpha_complex_3d_persistence.cpp | 244 ------------------- .../weighted_alpha_complex_3d_persistence.cpp | 267 --------------------- 6 files changed, 524 insertions(+), 528 deletions(-) create mode 100644 src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp create mode 100644 src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp delete mode 100644 src/Persistent_cohomology/example/exact_alpha_complex_3d_persistence.cpp delete mode 100644 src/Persistent_cohomology/example/weighted_alpha_complex_3d_persistence.cpp (limited to 'src/Alpha_complex') diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt index 69dbc28a..120ca416 100644 --- a/src/Alpha_complex/utilities/CMakeLists.txt +++ b/src/Alpha_complex/utilities/CMakeLists.txt @@ -4,14 +4,27 @@ project(Alpha_complex_utilities) if(CGAL_FOUND) add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp) target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY}) + add_executable(exact_alpha_complex_3d_persistence exact_alpha_complex_3d_persistence.cpp) + target_link_libraries(exact_alpha_complex_3d_persistence ${CGAL_LIBRARY}) + add_executable(weighted_alpha_complex_3d_persistence weighted_alpha_complex_3d_persistence.cpp) + target_link_libraries(weighted_alpha_complex_3d_persistence ${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(NAME Alpha_complex_utilities_alpha_complex_3d_persistence COMMAND $ "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45") + add_test(NAME Persistent_cohomology_example_exact_alpha_complex_3d COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45") + add_test(NAME Persistent_cohomology_example_weighted_alpha_complex_3d COMMAND $ + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.weights" "2" "0.45") install(TARGETS alpha_complex_3d_persistence DESTINATION bin) + install(TARGETS exact_alpha_complex_3d_persistence DESTINATION bin) + install(TARGETS weighted_alpha_complex_3d_persistence DESTINATION bin) if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) add_executable (alpha_complex_persistence alpha_complex_persistence.cpp) diff --git a/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..8ce68406 --- /dev/null +++ b/src/Alpha_complex/utilities/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 . + */ + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../utilities/alpha_complex_3d_helper.h" + +// 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; +using Fb = CGAL::Alpha_shape_cell_base_3; +using Tds = CGAL::Triangulation_data_structure_3; +using Triangulation_3 = CGAL::Delaunay_triangulation_3; +using Alpha_shape_3 = CGAL::Alpha_shape_3; + +// 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, + CGAL::cpp11::tuple >, + std::back_insert_iterator< std::vector > > >; +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; + +// gudhi type definition +using ST = Gudhi::Simplex_tree; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map; +using Alpha_shape_simplex_tree_pair = std::pair; +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:\n" << progName << " path_to_OFF_file coeff_field_characteristic[integer " << + "> 0] min_persistence[float >= -1.0]\n"; + std::cerr << " path_to_OFF_file is the path to your points cloud in OFF format.\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 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 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 the_objects; + std::vector the_alpha_values; + + Dispatch disp = CGAL::dispatch_output(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::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(&object_iterator)) { + vertex_list = from_cell(*cell); + count_cells++; + if (dim_max < 3) { + // Cell is of dim 3 + dim_max = 3; + } + } else if (const Facet * facet = CGAL::object_cast(&object_iterator)) { + vertex_list = from_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(&object_iterator)) { + vertex_list = from_edge(*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(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex(*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_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() << " "; +#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/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp new file mode 100644 index 00000000..a4ecf9da --- /dev/null +++ b/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp @@ -0,0 +1,267 @@ +/* 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 . + */ + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../utilities/alpha_complex_3d_helper.h" + +// Traits +using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; +using Gt = CGAL::Regular_triangulation_euclidean_traits_3; +using Vb = CGAL::Alpha_shape_vertex_base_3; +using Fb = CGAL::Alpha_shape_cell_base_3; +using Tds = CGAL::Triangulation_data_structure_3; +using Triangulation_3 = CGAL::Regular_triangulation_3; +using Alpha_shape_3 = CGAL::Alpha_shape_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, + CGAL::cpp11::tuple >, + std::back_insert_iterator< std::vector > > >; +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; + +// gudhi type definition +using ST = Gudhi::Simplex_tree; +using Filtration_value = ST::Filtration_value; +using Simplex_tree_vertex = ST::Vertex_handle; +using Alpha_shape_simplex_tree_map = std::map; +using Alpha_shape_simplex_tree_pair = std::pair; +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:\n" << progName << " path_to_OFF_file path_to_weight_file coeff_field_characteristic[integer " << + "> 0] min_persistence[float >= -1.0]\n"; + std::cerr << " path_to_OFF_file is the path to your points cloud in OFF format.\n"; + std::cerr << " path_to_weight_file is the path to the weights of your points cloud (one value per line.)\n"; + std::cerr << " Weights values are explained on CGAL documentation:\n"; + std::cerr << " https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0\n"; + std::cerr << " https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation\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 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 lp = off_reader.get_point_cloud(); + + // Read weights information from file + std::ifstream weights_ifstr(argv[2]); + std::vector wp; + if (weights_ifstr.good()) { + double weight = 0.0; + std::size_t index = 0; + wp.reserve(lp.size()); + // 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(wp.begin(), wp.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 the_objects; + std::vector the_alpha_values; + + Dispatch disp = CGAL::dispatch_output(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::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(&object_iterator)) { + vertex_list = from_cell(*cell); + count_cells++; + if (dim_max < 3) { + // Cell is of dim 3 + dim_max = 3; + } + } else if (const Facet * facet = CGAL::object_cast(&object_iterator)) { + vertex_list = from_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(&object_iterator)) { + vertex_list = from_edge(*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(&object_iterator)) { + count_vertices++; + vertex_list = from_vertex(*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_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() << " "; +#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; +} diff --git a/src/Persistent_cohomology/example/CMakeLists.txt b/src/Persistent_cohomology/example/CMakeLists.txt index 98ed753d..18e2913b 100644 --- a/src/Persistent_cohomology/example/CMakeLists.txt +++ b/src/Persistent_cohomology/example/CMakeLists.txt @@ -55,23 +55,6 @@ if(GMP_FOUND) endif(GMP_FOUND) if(CGAL_FOUND) - add_executable(exact_alpha_complex_3d_persistence exact_alpha_complex_3d_persistence.cpp) - target_link_libraries(exact_alpha_complex_3d_persistence ${CGAL_LIBRARY}) - add_executable(weighted_alpha_complex_3d_persistence weighted_alpha_complex_3d_persistence.cpp) - target_link_libraries(weighted_alpha_complex_3d_persistence ${CGAL_LIBRARY}) - - if (TBB_FOUND) - 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(NAME Persistent_cohomology_example_exact_alpha_complex_3d COMMAND $ - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45") - add_test(NAME Persistent_cohomology_example_weighted_alpha_complex_3d COMMAND $ - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.weights" "2" "0.45") - - install(TARGETS exact_alpha_complex_3d_persistence DESTINATION bin) - install(TARGETS weighted_alpha_complex_3d_persistence DESTINATION bin) - if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) add_executable(custom_persistence_sort custom_persistence_sort.cpp) target_link_libraries(custom_persistence_sort ${CGAL_LIBRARY}) diff --git a/src/Persistent_cohomology/example/exact_alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/exact_alpha_complex_3d_persistence.cpp deleted file mode 100644 index 8ce68406..00000000 --- a/src/Persistent_cohomology/example/exact_alpha_complex_3d_persistence.cpp +++ /dev/null @@ -1,244 +0,0 @@ -/* 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 . - */ - -#include - -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#include "../utilities/alpha_complex_3d_helper.h" - -// 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; -using Fb = CGAL::Alpha_shape_cell_base_3; -using Tds = CGAL::Triangulation_data_structure_3; -using Triangulation_3 = CGAL::Delaunay_triangulation_3; -using Alpha_shape_3 = CGAL::Alpha_shape_3; - -// 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, - CGAL::cpp11::tuple >, - std::back_insert_iterator< std::vector > > >; -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; - -// gudhi type definition -using ST = Gudhi::Simplex_tree; -using Filtration_value = ST::Filtration_value; -using Simplex_tree_vertex = ST::Vertex_handle; -using Alpha_shape_simplex_tree_map = std::map; -using Alpha_shape_simplex_tree_pair = std::pair; -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:\n" << progName << " path_to_OFF_file coeff_field_characteristic[integer " << - "> 0] min_persistence[float >= -1.0]\n"; - std::cerr << " path_to_OFF_file is the path to your points cloud in OFF format.\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 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 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 the_objects; - std::vector the_alpha_values; - - Dispatch disp = CGAL::dispatch_output(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::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(&object_iterator)) { - vertex_list = from_cell(*cell); - count_cells++; - if (dim_max < 3) { - // Cell is of dim 3 - dim_max = 3; - } - } else if (const Facet * facet = CGAL::object_cast(&object_iterator)) { - vertex_list = from_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(&object_iterator)) { - vertex_list = from_edge(*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(&object_iterator)) { - count_vertices++; - vertex_list = from_vertex(*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_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() << " "; -#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/weighted_alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/weighted_alpha_complex_3d_persistence.cpp deleted file mode 100644 index a4ecf9da..00000000 --- a/src/Persistent_cohomology/example/weighted_alpha_complex_3d_persistence.cpp +++ /dev/null @@ -1,267 +0,0 @@ -/* 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 . - */ - -#include - -#include -#include -#include - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "../utilities/alpha_complex_3d_helper.h" - -// Traits -using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; -using Gt = CGAL::Regular_triangulation_euclidean_traits_3; -using Vb = CGAL::Alpha_shape_vertex_base_3; -using Fb = CGAL::Alpha_shape_cell_base_3; -using Tds = CGAL::Triangulation_data_structure_3; -using Triangulation_3 = CGAL::Regular_triangulation_3; -using Alpha_shape_3 = CGAL::Alpha_shape_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, - CGAL::cpp11::tuple >, - std::back_insert_iterator< std::vector > > >; -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; - -// gudhi type definition -using ST = Gudhi::Simplex_tree; -using Filtration_value = ST::Filtration_value; -using Simplex_tree_vertex = ST::Vertex_handle; -using Alpha_shape_simplex_tree_map = std::map; -using Alpha_shape_simplex_tree_pair = std::pair; -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:\n" << progName << " path_to_OFF_file path_to_weight_file coeff_field_characteristic[integer " << - "> 0] min_persistence[float >= -1.0]\n"; - std::cerr << " path_to_OFF_file is the path to your points cloud in OFF format.\n"; - std::cerr << " path_to_weight_file is the path to the weights of your points cloud (one value per line.)\n"; - std::cerr << " Weights values are explained on CGAL documentation:\n"; - std::cerr << " https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0\n"; - std::cerr << " https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation\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 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 lp = off_reader.get_point_cloud(); - - // Read weights information from file - std::ifstream weights_ifstr(argv[2]); - std::vector wp; - if (weights_ifstr.good()) { - double weight = 0.0; - std::size_t index = 0; - wp.reserve(lp.size()); - // 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(wp.begin(), wp.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 the_objects; - std::vector the_alpha_values; - - Dispatch disp = CGAL::dispatch_output(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::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(&object_iterator)) { - vertex_list = from_cell(*cell); - count_cells++; - if (dim_max < 3) { - // Cell is of dim 3 - dim_max = 3; - } - } else if (const Facet * facet = CGAL::object_cast(&object_iterator)) { - vertex_list = from_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(&object_iterator)) { - vertex_list = from_edge(*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(&object_iterator)) { - count_vertices++; - vertex_list = from_vertex(*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_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() << " "; -#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; -} -- cgit v1.2.3 From 6dbf0ab8f84fd57d4fb80cd3c1bcf0aca6f0bb05 Mon Sep 17 00:00:00 2001 From: cjamin Date: Wed, 15 Nov 2017 15:57:48 +0000 Subject: Doc git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/add_utils_in_gudhi_v2@2884 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: c5e0a65a484a71d8a9fc1fc2b695dd20a367fd0c --- src/Alpha_complex/utilities/CMakeLists.txt | 4 +- src/Alpha_complex/utilities/README | 17 ++++++++ src/Bottleneck_distance/example/CMakeLists.txt | 4 ++ src/Bottleneck_distance/utilities/CMakeLists.txt | 11 ++--- .../utilities/bottleneck_read_file.cpp | 50 ++++++++++++++++++++++ .../utilities/bottleneck_read_file_example.cpp | 50 ---------------------- 6 files changed, 77 insertions(+), 59 deletions(-) create mode 100644 src/Bottleneck_distance/utilities/bottleneck_read_file.cpp delete mode 100644 src/Bottleneck_distance/utilities/bottleneck_read_file_example.cpp (limited to 'src/Alpha_complex') diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt index 120ca416..18699650 100644 --- a/src/Alpha_complex/utilities/CMakeLists.txt +++ b/src/Alpha_complex/utilities/CMakeLists.txt @@ -17,9 +17,9 @@ if(CGAL_FOUND) add_test(NAME Alpha_complex_utilities_alpha_complex_3d_persistence COMMAND $ "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45") - add_test(NAME Persistent_cohomology_example_exact_alpha_complex_3d COMMAND $ + add_test(NAME Alpha_complex_utilities_exact_alpha_complex_3d COMMAND $ "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45") - add_test(NAME Persistent_cohomology_example_weighted_alpha_complex_3d COMMAND $ + add_test(NAME Alpha_complex_utilities_weighted_alpha_complex_3d COMMAND $ "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.weights" "2" "0.45") install(TARGETS alpha_complex_3d_persistence DESTINATION bin) diff --git a/src/Alpha_complex/utilities/README b/src/Alpha_complex/utilities/README index 67ea1875..6e7d16b7 100644 --- a/src/Alpha_complex/utilities/README +++ b/src/Alpha_complex/utilities/README @@ -39,6 +39,23 @@ N.B.: +## `exact_alpha_complex_3d_persistence` ## +Same as `alpha_complex_3d_persistence`, but using exact computation. It is slower, but it is necessary when points are on a grid for instance. + + + +## `weighted_alpha_complex_3d_persistence` ## +Same as `alpha_complex_3d_persistence`, but using weighted points. + +**Usage** +`weighted_alpha_complex_3d_persistence

` +where +`` is the path to the input point cloud in OFF format. +`` is the path to the file containing the weights of the points (one value per line). +`

` is the characteristic p of the coefficient field *Z/pZ* for computing homology. It must be a stricly positive integer. +`` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. + + ## `periodic_alpha_complex_3d_persistence` ## This program computes the persistent homology with coefficient field Z/pZ of the 3D periodic alpha complex built from a 3D point cloud. The output diagram contains one bar per line, written with the convention: diff --git a/src/Bottleneck_distance/example/CMakeLists.txt b/src/Bottleneck_distance/example/CMakeLists.txt index 9677f5c5..2ae45dc5 100644 --- a/src/Bottleneck_distance/example/CMakeLists.txt +++ b/src/Bottleneck_distance/example/CMakeLists.txt @@ -6,6 +6,10 @@ if (NOT CGAL_VERSION VERSION_LESS 4.8.1) add_executable (alpha_rips_persistence_bottleneck_distance alpha_rips_persistence_bottleneck_distance.cpp) target_link_libraries(alpha_rips_persistence_bottleneck_distance ${Boost_PROGRAM_OPTIONS_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(alpha_rips_persistence_bottleneck_distance ${TBB_LIBRARIES}) + endif(TBB_FOUND) + add_test(NAME Bottleneck_distance_example_basic COMMAND $) add_test(NAME Bottleneck_distance_example_alpha_rips_persistence_bottleneck COMMAND $ diff --git a/src/Bottleneck_distance/utilities/CMakeLists.txt b/src/Bottleneck_distance/utilities/CMakeLists.txt index decdc789..cdf76c85 100644 --- a/src/Bottleneck_distance/utilities/CMakeLists.txt +++ b/src/Bottleneck_distance/utilities/CMakeLists.txt @@ -2,15 +2,12 @@ cmake_minimum_required(VERSION 2.6) project(Bottleneck_distance_utilities) if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) - add_executable (bottleneck_read_file_example bottleneck_read_file_example.cpp) - if (TBB_FOUND) - target_link_libraries(alpha_rips_persistence_bottleneck_distance ${TBB_LIBRARIES}) - endif(TBB_FOUND) + add_executable (bottleneck_read_file bottleneck_read_file.cpp) - add_test(NAME Bottleneck_read_file_example - COMMAND $ + add_test(NAME Bottleneck_distance_utilities_Bottleneck_read_file + COMMAND $ "${CMAKE_SOURCE_DIR}/data/persistence_diagram/first.pers" "${CMAKE_SOURCE_DIR}/data/persistence_diagram/second.pers") - install(TARGETS bottleneck_read_file_example DESTINATION bin) + install(TARGETS bottleneck_read_file DESTINATION bin) endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) diff --git a/src/Bottleneck_distance/utilities/bottleneck_read_file.cpp b/src/Bottleneck_distance/utilities/bottleneck_read_file.cpp new file mode 100644 index 00000000..9dd52b31 --- /dev/null +++ b/src/Bottleneck_distance/utilities/bottleneck_read_file.cpp @@ -0,0 +1,50 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Authors: Francois Godi, small modifications by Pawel Dlotko + * + * 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 . + */ + +#include +#include +#include +#include +#include // for pair +#include +#include // for numeric_limits + +int main(int argc, char** argv) { + if (argc < 3) { + std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file" << + " should contain a birth-death pair per line. Third, optional parameter is an error bound on the bottleneck" << + " distance (set by default to the smallest positive double value). If you set the error bound to 0, be" << + " aware this version is exact but expensive. The program will now terminate \n"; + return -1; + } + std::vector> diag1 = Gudhi::read_persistence_intervals_in_dimension(argv[1]); + std::vector> diag2 = Gudhi::read_persistence_intervals_in_dimension(argv[2]); + + double tolerance = std::numeric_limits::min(); + if (argc == 4) { + tolerance = atof(argv[3]); + } + double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance); + std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl; + + return 0; +} diff --git a/src/Bottleneck_distance/utilities/bottleneck_read_file_example.cpp b/src/Bottleneck_distance/utilities/bottleneck_read_file_example.cpp deleted file mode 100644 index 9dd52b31..00000000 --- a/src/Bottleneck_distance/utilities/bottleneck_read_file_example.cpp +++ /dev/null @@ -1,50 +0,0 @@ -/* This file is part of the Gudhi Library. The Gudhi library - * (Geometric Understanding in Higher Dimensions) is a generic C++ - * library for computational topology. - * - * Authors: Francois Godi, small modifications by Pawel Dlotko - * - * 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 . - */ - -#include -#include -#include -#include -#include // for pair -#include -#include // for numeric_limits - -int main(int argc, char** argv) { - if (argc < 3) { - std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file" << - " should contain a birth-death pair per line. Third, optional parameter is an error bound on the bottleneck" << - " distance (set by default to the smallest positive double value). If you set the error bound to 0, be" << - " aware this version is exact but expensive. The program will now terminate \n"; - return -1; - } - std::vector> diag1 = Gudhi::read_persistence_intervals_in_dimension(argv[1]); - std::vector> diag2 = Gudhi::read_persistence_intervals_in_dimension(argv[2]); - - double tolerance = std::numeric_limits::min(); - if (argc == 4) { - tolerance = atof(argv[3]); - } - double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance); - std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl; - - return 0; -} -- cgit v1.2.3 From aafe601c247c78a474c6af110ca5b34540a34f44 Mon Sep 17 00:00:00 2001 From: vrouvrea Date: Thu, 16 Nov 2017 08:31:26 +0000 Subject: Code review : Homogenize alpha shape examples command line git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/add_utils_in_gudhi_v2@2887 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 8c692434eb48c525047bc34a6fab7d7c1af64d56 --- src/Alpha_complex/utilities/CMakeLists.txt | 19 ++- src/Alpha_complex/utilities/README | 89 +++++++++------ .../utilities/alpha_complex_3d_helper.h | 8 +- .../utilities/alpha_complex_3d_persistence.cpp | 115 ++++++++++++------- .../utilities/alpha_complex_persistence.cpp | 63 +++++----- .../exact_alpha_complex_3d_persistence.cpp | 107 ++++++++++------- .../periodic_alpha_complex_3d_persistence.cpp | 124 ++++++++++++-------- .../weighted_alpha_complex_3d_persistence.cpp | 127 +++++++++++++-------- src/Bottleneck_distance/utilities/CMakeLists.txt | 3 + src/Rips_complex/utilities/README | 3 +- 10 files changed, 406 insertions(+), 252 deletions(-) (limited to 'src/Alpha_complex') diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt index 18699650..cefc1ad6 100644 --- a/src/Alpha_complex/utilities/CMakeLists.txt +++ b/src/Alpha_complex/utilities/CMakeLists.txt @@ -3,11 +3,11 @@ project(Alpha_complex_utilities) if(CGAL_FOUND) add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp) - target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY}) + target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) add_executable(exact_alpha_complex_3d_persistence exact_alpha_complex_3d_persistence.cpp) - target_link_libraries(exact_alpha_complex_3d_persistence ${CGAL_LIBRARY}) + target_link_libraries(exact_alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) add_executable(weighted_alpha_complex_3d_persistence weighted_alpha_complex_3d_persistence.cpp) - target_link_libraries(weighted_alpha_complex_3d_persistence ${CGAL_LIBRARY}) + target_link_libraries(weighted_alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) if (TBB_FOUND) target_link_libraries(alpha_complex_3d_persistence ${TBB_LIBRARIES}) @@ -16,11 +16,11 @@ if(CGAL_FOUND) endif(TBB_FOUND) add_test(NAME Alpha_complex_utilities_alpha_complex_3d_persistence COMMAND $ - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45") + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45") add_test(NAME Alpha_complex_utilities_exact_alpha_complex_3d COMMAND $ - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45") + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45") add_test(NAME Alpha_complex_utilities_weighted_alpha_complex_3d COMMAND $ - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.weights" "2" "0.45") + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.weights" "-p" "2" "-m" "0.45") install(TARGETS alpha_complex_3d_persistence DESTINATION bin) install(TARGETS exact_alpha_complex_3d_persistence DESTINATION bin) @@ -28,11 +28,10 @@ if(CGAL_FOUND) if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0) add_executable (alpha_complex_persistence alpha_complex_persistence.cpp) - target_link_libraries(alpha_complex_persistence - ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + target_link_libraries(alpha_complex_persistence ${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 ${CGAL_LIBRARY}) + target_link_libraries(periodic_alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) if (TBB_FOUND) target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES}) @@ -41,7 +40,7 @@ if(CGAL_FOUND) add_test(NAME Alpha_complex_utilities_alpha_complex_persistence COMMAND $ "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45") add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $ - "${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") + "${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" "-p" "2" "-m" "0") install(TARGETS alpha_complex_persistence DESTINATION bin) install(TARGETS periodic_alpha_complex_3d_persistence DESTINATION bin) diff --git a/src/Alpha_complex/utilities/README b/src/Alpha_complex/utilities/README index 6e7d16b7..c3dd170b 100644 --- a/src/Alpha_complex/utilities/README +++ b/src/Alpha_complex/utilities/README @@ -8,22 +8,24 @@ This program computes the persistent homology with coefficient field Z/pZ of the where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number). **Usage** -`alpha_complex_3d_persistence

` -where -`` is the path to the input point cloud in OFF format. -`

` is the characteristic p of the coefficient field *Z/pZ* for computing homology. It must be a stricly positive integer. -`` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. +`alpha_complex_3d_persistence [options] ` + +**Allowed options** +* `-h [ --help ]` Produce help message +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. +* `-p [ --field-charac ]` (default=11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. **Example** -`alpha_complex_3d_persistence ../../data/points/tore3D_300.off 2 0.45` +`alpha_complex_3d_persistence ../../data/points/tore3D_300.off -p 2 -m 0.45` outputs: ``` 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 +2 0 0 inf +2 1 0.0682162 1.0001 +2 1 0.0934117 1.00003 +2 2 0.56444 1.03938 ``` Here we retrieve expected Betti numbers on a tore 3D: @@ -48,21 +50,34 @@ Same as `alpha_complex_3d_persistence`, but using exact computation. It is slowe Same as `alpha_complex_3d_persistence`, but using weighted points. **Usage** -`weighted_alpha_complex_3d_persistence

` -where -`` is the path to the input point cloud in OFF format. -`` is the path to the file containing the weights of the points (one value per line). -`

` is the characteristic p of the coefficient field *Z/pZ* for computing homology. It must be a stricly positive integer. -`` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. +`weighted_alpha_complex_3d_persistence [options] ` +**Allowed options** +* `-h [ --help ]` Produce help message +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. +* `-p [ --field-charac ]` (default=11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. +**Example** +`weighted_alpha_complex_3d_persistence ../../data/points/tore3D_300.off ../../data/points/tore3D_300.weights -p 2 -m 0.45` -## `periodic_alpha_complex_3d_persistence` ## -This program computes the persistent homology with coefficient field Z/pZ of the 3D periodic alpha complex built from a 3D point cloud. The output diagram contains one bar per line, written with the convention: +outputs: +``` +Simplex_tree dim: 3 +2 0 -1 inf +2 1 -0.931784 0.000103311 +2 1 -0.906588 2.60165e-05 +2 2 -0.43556 0.0393798 +``` -`p dim b d` +N.B.: +* Weights values are explained on CGAL [Alpha shape](https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0) +and [Regular triangulation](https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation) documentation. +* Filtration values are alpha square values. -where `dim` is the dimension of the homological feature, `b` and `d` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number). + +## `periodic_alpha_complex_3d_persistence` ## +Same as `alpha_complex_3d_persistence`, but using periodic alpha shape 3d. **Usage** `periodic_alpha_complex_3d_persistence

` @@ -72,21 +87,31 @@ where `

` is the characteristic p of the coefficient field *Z/pZ* for computing homology. It must be a stricly positive integer. `` is the minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. +**Usage** +`./periodic_alpha_complex_3d_persistence [options] input-file cuboid-file` + +**Allowed options** +* `-h [ --help ]` Produce help message +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. +* `-p [ --field-charac ]` (default=11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals + + **Example** -`periodic_alpha_complex_3d_persistence ../../data/points/grid_10_10_10_in_0_1.off ../../data/points/iso_cuboid_3_in_0_1.txt 3 1.0` +`periodic_alpha_complex_3d_persistence ../../data/points/grid_10_10_10_in_0_1.off ../../data/points/iso_cuboid_3_in_0_1.txt -p 3 -m 1.0` outputs: ``` Periodic Delaunay computed. Simplex_tree dim: 3 -3 0 0 inf -3 1 0.0025 inf -3 1 0.0025 inf -3 1 0.0025 inf -3 2 0.005 inf -3 2 0.005 inf -3 2 0.005 inf -3 3 0.0075 inf +3 0 0 inf +3 1 0.0025 inf +3 1 0.0025 inf +3 1 0.0025 inf +3 2 0.005 inf +3 2 0.005 inf +3 2 0.005 inf +3 3 0.0075 inf ``` Here we retrieve expected Betti numbers on an 3D iso-oriented cuboids: @@ -98,14 +123,12 @@ Betti numbers[3] = 1 ``` N.B.: -* `periodic_alpha_complex_3d_persistence` only accepts OFF files in dimension 3. -* In this example, the periodic cube is hard coded to { x = [0,1]; y = [0,1]; z = [0,1] } +* Cuboid file must be in the format described [here](http://gudhi.gforge.inria.fr/doc/latest/fileformats.html#FileFormatsIsoCuboid). * Filtration values are alpha square values. - ## `alpha_complex_persistence` ## This program computes the persistent homology with coefficient field Z/pZ of the dD alpha complex built from a dD point cloud. The output diagram contains one bar per line, written with the convention: @@ -119,7 +142,7 @@ where `dim` is the dimension of the homological feature, `b` and `d` are respect **Allowed options** * `-h [ --help ]` Produce help message -* `-o [ --output-file ]` Name of file in which the persistence diagram is written. By default, print in std::cout. +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. * `-r [ --max-alpha-square-value ]` (default = inf) Maximal alpha square value for the Alpha complex construction. * `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. * `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_helper.h b/src/Alpha_complex/utilities/alpha_complex_3d_helper.h index 7865e4ec..6b3b7d5d 100644 --- a/src/Alpha_complex/utilities/alpha_complex_3d_helper.h +++ b/src/Alpha_complex/utilities/alpha_complex_3d_helper.h @@ -23,7 +23,7 @@ #ifndef ALPHA_COMPLEX_3D_HELPER_H_ #define ALPHA_COMPLEX_3D_HELPER_H_ -template +template Vertex_list from_cell(const Cell_handle& ch) { Vertex_list the_list; for (auto i = 0; i < 4; i++) { @@ -35,7 +35,7 @@ Vertex_list from_cell(const Cell_handle& ch) { return the_list; } -template +template Vertex_list from_facet(const Facet& fct) { Vertex_list the_list; for (auto i = 0; i < 4; i++) { @@ -49,7 +49,7 @@ Vertex_list from_facet(const Facet& fct) { return the_list; } -template +template Vertex_list from_edge(const Edge_3& edg) { Vertex_list the_list; for (auto i = 0; i < 4; i++) { @@ -63,7 +63,7 @@ Vertex_list from_edge(const Edge_3& edg) { return the_list; } -template +template Vertex_list from_vertex(const Vertex_handle& vh) { Vertex_list the_list; #ifdef DEBUG_TRACES diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp index f63ff0f6..f2085ab2 100644 --- a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp @@ -20,6 +20,7 @@ * along with this program. If not, see . */ +#include #include #include @@ -56,10 +57,10 @@ 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, - CGAL::cpp11::tuple >, - std::back_insert_iterator< std::vector > > >; +using Dispatch = + CGAL::Dispatch_output_iterator, + CGAL::cpp11::tuple >, + std::back_insert_iterator > > >; using Cell_handle = Alpha_shape_3::Cell_handle; using Facet = Alpha_shape_3::Facet; using Edge_3 = Alpha_shape_3::Edge; @@ -70,42 +71,28 @@ using Vertex_list = std::list; using ST = Gudhi::Simplex_tree; using Filtration_value = ST::Filtration_value; using Simplex_tree_vertex = ST::Vertex_handle; -using Alpha_shape_simplex_tree_map = std::map; +using Alpha_shape_simplex_tree_map = std::map; using Alpha_shape_simplex_tree_pair = std::pair; -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:\n" << progName << " path_to_OFF_file coeff_field_characteristic[integer " << - "> 0] min_persistence[float >= -1.0]\n"; - std::cerr << " path_to_OFF_file is the path to your points cloud in OFF format.\n"; - exit(-1); -} +using Simplex_tree_vector_vertex = std::vector; +using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology; -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]); - } +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag, + int &coeff_field_characteristic, Filtration_value &min_persistence); - int coeff_field_characteristic = atoi(argv[2]); +int main(int argc, char **argv) { + std::string off_file_points; + std::string output_file_diag; + int coeff_field_characteristic; + Filtration_value min_persistence; - 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]); - } + program_options(argc, argv, off_file_points, output_file_diag, coeff_field_characteristic, min_persistence); - // 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 off_reader(offInputFile); + Gudhi::Points_3D_off_reader off_reader(off_file_points); // Check the read operation was correct if (!off_reader.is_valid()) { - std::cerr << "Unable to read file " << offInputFile << std::endl; - usage(argv[0]); + std::cerr << "Unable to read file " << off_file_points << std::endl; + exit(-1); } // Retrieve the triangulation @@ -143,28 +130,28 @@ int main(int argc, char * const argv[]) { 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(&object_iterator)) { + if (const Cell_handle *cell = CGAL::object_cast(&object_iterator)) { vertex_list = from_cell(*cell); count_cells++; if (dim_max < 3) { // Cell is of dim 3 dim_max = 3; } - } else if (const Facet * facet = CGAL::object_cast(&object_iterator)) { + } else if (const Facet *facet = CGAL::object_cast(&object_iterator)) { vertex_list = from_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(&object_iterator)) { + } else if (const Edge_3 *edge = CGAL::object_cast(&object_iterator)) { vertex_list = from_edge(*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(&object_iterator)) { + } else if (const Vertex_handle *vertex = CGAL::object_cast(&object_iterator)) { count_vertices++; vertex_list = from_vertex(*vertex); } @@ -190,7 +177,7 @@ int main(int argc, char * const argv[]) { } } // Construction of the simplex_tree - Filtration_value filtr = /*std::sqrt*/(*the_alpha_value_iterator); + Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator); #ifdef DEBUG_TRACES std::cout << "filtration = " << filtr << std::endl; #endif // DEBUG_TRACES @@ -211,7 +198,6 @@ int main(int argc, char * const argv[]) { 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; @@ -236,7 +222,58 @@ int main(int argc, char * const argv[]) { pcoh.compute_persistent_cohomology(min_persistence); - pcoh.output_diagram(); + // Output the diagram in filediag + if (output_file_diag.empty()) { + pcoh.output_diagram(); + } else { + std::cout << "Result in file: " << output_file_diag << std::endl; + std::ofstream out(output_file_diag); + pcoh.output_diagram(out); + out.close(); + } return 0; } + +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag, + int &coeff_field_characteristic, 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), + "Name of file containing a point set. Format is one point per line: X1 ... Xd "); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value(&output_file_diag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "field-charac,p", po::value(&coeff_field_characteristic)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::value(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length " + "intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map 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")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a 3D Alpha complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp index 9e84e91f..2105220a 100644 --- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -14,12 +14,9 @@ using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; -void program_options(int argc, char * argv[] - , std::string & off_file_points - , std::string & output_file_diag - , Filtration_value & alpha_square_max_value - , int & coeff_field_characteristic - , Filtration_value & min_persistence); +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag, + Filtration_value &alpha_square_max_value, int &coeff_field_characteristic, + Filtration_value &min_persistence); int main(int argc, char **argv) { std::string off_file_points; @@ -28,13 +25,13 @@ int main(int argc, char **argv) { int coeff_field_characteristic; Filtration_value min_persistence; - program_options(argc, argv, off_file_points, output_file_diag, alpha_square_max_value, - coeff_field_characteristic, min_persistence); + program_options(argc, argv, off_file_points, output_file_diag, alpha_square_max_value, coeff_field_characteristic, + min_persistence); // ---------------------------------------------------------------------------- // Init of an alpha complex from an OFF file // ---------------------------------------------------------------------------- - using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >; + using Kernel = CGAL::Epick_d; Gudhi::alpha_complex::Alpha_complex alpha_complex_from_file(off_file_points); Simplex_tree simplex; @@ -42,17 +39,16 @@ int main(int argc, char **argv) { // ---------------------------------------------------------------------------- // Display information about the alpha complex // ---------------------------------------------------------------------------- - std::cout << "Simplicial complex is of dimension " << simplex.dimension() << - " - " << simplex.num_simplices() << " simplices - " << - simplex.num_vertices() << " vertices." << std::endl; + std::cout << "Simplicial complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() + << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; // Sort the simplices in the order of the filtration simplex.initialize_filtration(); std::cout << "Simplex_tree dim: " << simplex.dimension() << std::endl; // Compute the persistence diagram of the complex - Gudhi::persistent_cohomology::Persistent_cohomology< Simplex_tree, - Gudhi::persistent_cohomology::Field_Zp > pcoh(simplex); + Gudhi::persistent_cohomology::Persistent_cohomology pcoh( + simplex); // initializes the coefficient field for homology pcoh.init_coefficients(coeff_field_characteristic); @@ -72,30 +68,26 @@ int main(int argc, char **argv) { return 0; } -void program_options(int argc, char * argv[] - , std::string & off_file_points - , std::string & output_file_diag - , Filtration_value & alpha_square_max_value - , int & coeff_field_characteristic - , Filtration_value & min_persistence) { +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag, + Filtration_value &alpha_square_max_value, int &coeff_field_characteristic, + 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), - "Name of file containing a point set. Format is one point per line: X1 ... Xd "); + hidden.add_options()("input-file", po::value(&off_file_points), + "Name of file containing a point set. Format is one point per line: X1 ... Xd "); po::options_description visible("Allowed options", 100); - visible.add_options() - ("help,h", "produce help message") - ("output-file,o", po::value(&output_file_diag)->default_value(std::string()), - "Name of file in which the persistence diagram is written. Default print in std::cout") - ("max-alpha-square-value,r", - po::value(&alpha_square_max_value)->default_value(std::numeric_limits::infinity()), - "Maximal alpha square value for the Alpha complex construction.") - ("field-charac,p", po::value(&coeff_field_characteristic)->default_value(11), - "Characteristic p of the coefficient field Z/pZ for computing homology.") - ("min-persistence,m", po::value(&min_persistence), - "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals"); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value(&output_file_diag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "max-alpha-square-value,r", po::value(&alpha_square_max_value) + ->default_value(std::numeric_limits::infinity()), + "Maximal alpha square value for the Alpha complex construction.")( + "field-charac,p", po::value(&coeff_field_characteristic)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::value(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length " + "intervals"); po::positional_options_description pos; pos.add("input-file", 1); @@ -104,8 +96,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/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp index 8ce68406..7919c7fd 100644 --- a/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp +++ b/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp @@ -20,6 +20,7 @@ * along with this program. If not, see . */ +#include #include #include @@ -40,7 +41,7 @@ #include #include -#include "../utilities/alpha_complex_3d_helper.h" +#include "alpha_complex_3d_helper.h" // Alpha_shape_3 templates type definitions using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; @@ -57,10 +58,10 @@ 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, - CGAL::cpp11::tuple >, - std::back_insert_iterator< std::vector > > >; +using Dispatch = + CGAL::Dispatch_output_iterator, + CGAL::cpp11::tuple >, + std::back_insert_iterator > > >; using Cell_handle = Alpha_shape_3::Cell_handle; using Facet = Alpha_shape_3::Facet; using Edge_3 = Alpha_shape_3::Edge; @@ -71,42 +72,28 @@ using Vertex_list = std::list; using ST = Gudhi::Simplex_tree; using Filtration_value = ST::Filtration_value; using Simplex_tree_vertex = ST::Vertex_handle; -using Alpha_shape_simplex_tree_map = std::map; +using Alpha_shape_simplex_tree_map = std::map; using Alpha_shape_simplex_tree_pair = std::pair; -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:\n" << progName << " path_to_OFF_file coeff_field_characteristic[integer " << - "> 0] min_persistence[float >= -1.0]\n"; - std::cerr << " path_to_OFF_file is the path to your points cloud in OFF format.\n"; - exit(-1); -} +using Simplex_tree_vector_vertex = std::vector; +using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology; -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]); - } +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag, + int &coeff_field_characteristic, Filtration_value &min_persistence); - int coeff_field_characteristic = atoi(argv[2]); +int main(int argc, char **argv) { + std::string off_file_points; + std::string output_file_diag; + int coeff_field_characteristic; + Filtration_value min_persistence; - 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]); - } + program_options(argc, argv, off_file_points, output_file_diag, coeff_field_characteristic, min_persistence); - // 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 off_reader(offInputFile); + Gudhi::Points_3D_off_reader off_reader(off_file_points); // Check the read operation was correct if (!off_reader.is_valid()) { - std::cerr << "Unable to read file " << offInputFile << std::endl; - usage(argv[0]); + std::cerr << "Unable to read file " << off_file_points << std::endl; + exit(-1); } // Retrieve the triangulation @@ -144,28 +131,28 @@ int main(int argc, char * const argv[]) { 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(&object_iterator)) { + if (const Cell_handle *cell = CGAL::object_cast(&object_iterator)) { vertex_list = from_cell(*cell); count_cells++; if (dim_max < 3) { // Cell is of dim 3 dim_max = 3; } - } else if (const Facet * facet = CGAL::object_cast(&object_iterator)) { + } else if (const Facet *facet = CGAL::object_cast(&object_iterator)) { vertex_list = from_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(&object_iterator)) { + } else if (const Edge_3 *edge = CGAL::object_cast(&object_iterator)) { vertex_list = from_edge(*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(&object_iterator)) { + } else if (const Vertex_handle *vertex = CGAL::object_cast(&object_iterator)) { count_vertices++; vertex_list = from_vertex(*vertex); } @@ -192,7 +179,7 @@ int main(int argc, char * const argv[]) { } // 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()); + 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 @@ -213,7 +200,6 @@ int main(int argc, char * const argv[]) { 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; @@ -242,3 +228,46 @@ int main(int argc, char * const argv[]) { return 0; } + +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &output_file_diag, + int &coeff_field_characteristic, 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), + "Name of file containing a point set. Format is one point per line: X1 ... Xd "); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value(&output_file_diag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "field-charac,p", po::value(&coeff_field_characteristic)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::value(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length " + "intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map 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")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a 3D Alpha complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp index 8140a3c5..10b62f75 100644 --- a/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp +++ b/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp @@ -20,6 +20,7 @@ * along with this program. If not, see . */ +#include #include #include @@ -63,10 +64,10 @@ using Point_3 = PK::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, - CGAL::cpp11::tuple >, - std::back_insert_iterator< std::vector > > >; +using Dispatch = + CGAL::Dispatch_output_iterator, + CGAL::cpp11::tuple >, + std::back_insert_iterator > > >; using Cell_handle = Alpha_shape_3::Cell_handle; using Facet = Alpha_shape_3::Facet; using Edge_3 = Alpha_shape_3::Edge; @@ -77,54 +78,41 @@ using Vertex_list = std::list; using ST = Gudhi::Simplex_tree; using Filtration_value = ST::Filtration_value; using Simplex_tree_vertex = ST::Vertex_handle; -using Alpha_shape_simplex_tree_map = std::map; +using Alpha_shape_simplex_tree_map = std::map; using Alpha_shape_simplex_tree_pair = std::pair; -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:\n" << progName << " path_to_OFF_file path_to_iso_cuboid_3_file coeff_field_characteristic[" << - "integer > 0] min_persistence[float >= -1.0]\n" << - " path_to_OFF_file is the path to your points cloud in OFF format.\n" << - " path_to_iso_cuboid_3_file is the path to the iso cuboid file with the following format :\n" << - " x_min y_min z_min x_max y_max z_max\n" << - " In this example, the periodic cube will be " << - "{ x = [x_min,x_max]; y = [y_min,y_max]; z = [z_min,z_max] }.\n" << - " For more information, please refer to\n" << - " https://doc.cgal.org/latest/Kernel_23/classCGAL_1_1Iso__cuboid__3.html\n"; - - exit(-1); -} +using Simplex_tree_vector_vertex = std::vector; +using Persistent_cohomology = + Gudhi::persistent_cohomology::Persistent_cohomology; -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]); - } +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &cuboid_file, + std::string &output_file_diag, int &coeff_field_characteristic, Filtration_value &min_persistence); - int coeff_field_characteristic = atoi(argv[3]); - Filtration_value min_persistence = strtof(argv[4], nullptr); +int main(int argc, char **argv) { + std::string off_file_points; + std::string cuboid_file; + std::string output_file_diag; + int coeff_field_characteristic; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, cuboid_file, output_file_diag, coeff_field_characteristic, + min_persistence); - // 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 off_reader(offInputFile); + Gudhi::Points_3D_off_reader off_reader(off_file_points); // Check the read operation was correct if (!off_reader.is_valid()) { - std::cerr << "Unable to read file " << offInputFile << std::endl; - usage(argv[0]); + std::cerr << "Unable to read OFF file " << off_file_points << std::endl; + exit(-1); } // Read iso_cuboid_3 information from file - std::ifstream iso_cuboid_str(argv[2]); + std::ifstream iso_cuboid_str(cuboid_file); 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]); + std::cerr << "Unable to read file " << cuboid_file << std::endl; + exit(-1); } // Retrieve the triangulation @@ -168,29 +156,29 @@ int main(int argc, char * const argv[]) { 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(&object_iterator)) { + if (const Cell_handle *cell = CGAL::object_cast(&object_iterator)) { vertex_list = from_cell(*cell); count_cells++; if (dim_max < 3) { // Cell is of dim 3 dim_max = 3; } - } else if (const Facet * facet = CGAL::object_cast(&object_iterator)) { + } else if (const Facet *facet = CGAL::object_cast(&object_iterator)) { vertex_list = from_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(&object_iterator)) { + } else if (const Edge_3 *edge = CGAL::object_cast(&object_iterator)) { vertex_list = from_edge(*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(&object_iterator)) { + } else if (const Alpha_shape_3::Vertex_handle *vertex = + CGAL::object_cast(&object_iterator)) { count_vertices++; vertex_list = from_vertex(*vertex); } @@ -216,7 +204,7 @@ int main(int argc, char * const argv[]) { } } // Construction of the simplex_tree - Filtration_value filtr = /*std::sqrt*/(*the_alpha_value_iterator); + Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator); #ifdef DEBUG_TRACES std::cout << "filtration = " << filtr << std::endl; #endif // DEBUG_TRACES @@ -237,7 +225,6 @@ int main(int argc, char * const argv[]) { 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; @@ -266,3 +253,50 @@ int main(int argc, char * const argv[]) { return 0; } + +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &cuboid_file, + std::string &output_file_diag, int &coeff_field_characteristic, + 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), + "Name of file containing a point set. Format is one point per line: X1 ... Xd ")( + "cuboid-file", po::value(&cuboid_file), + "Name of file describing the periodic domain. Format is: min_hx min_hy min_hz\nmax_hx max_hy max_hz"); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value(&output_file_diag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "field-charac,p", po::value(&coeff_field_characteristic)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::value(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length " + "intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + pos.add("cuboid-file", 2); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map 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") || !vm.count("cuboid-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a periodic 3D Alpha complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file cuboid-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp index a4ecf9da..293170f7 100644 --- a/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp +++ b/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp @@ -20,6 +20,7 @@ * along with this program. If not, see . */ +#include #include #include @@ -42,7 +43,7 @@ #include #include -#include "../utilities/alpha_complex_3d_helper.h" +#include "alpha_complex_3d_helper.h" // Traits using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel; @@ -60,10 +61,10 @@ 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, - CGAL::cpp11::tuple >, - std::back_insert_iterator< std::vector > > >; +using Dispatch = + CGAL::Dispatch_output_iterator, + CGAL::cpp11::tuple >, + std::back_insert_iterator > > >; using Cell_handle = Alpha_shape_3::Cell_handle; using Facet = Alpha_shape_3::Facet; using Edge_3 = Alpha_shape_3::Edge; @@ -74,48 +75,38 @@ using Vertex_list = std::list; using ST = Gudhi::Simplex_tree; using Filtration_value = ST::Filtration_value; using Simplex_tree_vertex = ST::Vertex_handle; -using Alpha_shape_simplex_tree_map = std::map; +using Alpha_shape_simplex_tree_map = std::map; using Alpha_shape_simplex_tree_pair = std::pair; -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:\n" << progName << " path_to_OFF_file path_to_weight_file coeff_field_characteristic[integer " << - "> 0] min_persistence[float >= -1.0]\n"; - std::cerr << " path_to_OFF_file is the path to your points cloud in OFF format.\n"; - std::cerr << " path_to_weight_file is the path to the weights of your points cloud (one value per line.)\n"; - std::cerr << " Weights values are explained on CGAL documentation:\n"; - std::cerr << " https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0\n"; - std::cerr << " https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation\n"; - exit(-1); -} +using Simplex_tree_vector_vertex = std::vector; +using Persistent_cohomology = + Gudhi::persistent_cohomology::Persistent_cohomology; -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]); - } +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &weight_file, + std::string &output_file_diag, int &coeff_field_characteristic, Filtration_value &min_persistence); - int coeff_field_characteristic = atoi(argv[3]); - Filtration_value min_persistence = strtof(argv[4], nullptr); +int main(int argc, char **argv) { + std::string off_file_points; + std::string weight_file; + std::string output_file_diag; + int coeff_field_characteristic; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, weight_file, output_file_diag, coeff_field_characteristic, + min_persistence); - // 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 off_reader(offInputFile); + Gudhi::Points_3D_off_reader off_reader(off_file_points); // Check the read operation was correct if (!off_reader.is_valid()) { - std::cerr << "Unable to read file " << offInputFile << std::endl; - usage(argv[0]); + std::cerr << "Unable to read OFF file " << off_file_points << std::endl; + exit(-1); } // Retrieve the triangulation std::vector lp = off_reader.get_point_cloud(); // Read weights information from file - std::ifstream weights_ifstr(argv[2]); + std::ifstream weights_ifstr(weight_file); std::vector wp; if (weights_ifstr.good()) { double weight = 0.0; @@ -127,12 +118,12 @@ int main(int argc, char * const argv[]) { index++; } if (index != lp.size()) { - std::cerr << "Bad number of weights in file " << argv[2] << std::endl; - usage(argv[0]); + std::cerr << "Bad number of weights in file " << weight_file << std::endl; + exit(-1); } } else { - std::cerr << "Unable to read file " << argv[2] << std::endl; - usage(argv[0]); + std::cerr << "Unable to read weights file " << weight_file << std::endl; + exit(-1); } // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. @@ -167,29 +158,29 @@ int main(int argc, char * const argv[]) { 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(&object_iterator)) { + if (const Cell_handle *cell = CGAL::object_cast(&object_iterator)) { vertex_list = from_cell(*cell); count_cells++; if (dim_max < 3) { // Cell is of dim 3 dim_max = 3; } - } else if (const Facet * facet = CGAL::object_cast(&object_iterator)) { + } else if (const Facet *facet = CGAL::object_cast(&object_iterator)) { vertex_list = from_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(&object_iterator)) { + } else if (const Edge_3 *edge = CGAL::object_cast(&object_iterator)) { vertex_list = from_edge(*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(&object_iterator)) { + } else if (const Alpha_shape_3::Vertex_handle *vertex = + CGAL::object_cast(&object_iterator)) { count_vertices++; vertex_list = from_vertex(*vertex); } @@ -215,7 +206,7 @@ int main(int argc, char * const argv[]) { } } // Construction of the simplex_tree - Filtration_value filtr = /*std::sqrt*/(*the_alpha_value_iterator); + Filtration_value filtr = /*std::sqrt*/ (*the_alpha_value_iterator); #ifdef DEBUG_TRACES std::cout << "filtration = " << filtr << std::endl; #endif // DEBUG_TRACES @@ -236,7 +227,6 @@ int main(int argc, char * const argv[]) { 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; @@ -265,3 +255,50 @@ int main(int argc, char * const argv[]) { return 0; } + +void program_options(int argc, char *argv[], std::string &off_file_points, std::string &weight_file, + std::string &output_file_diag, int &coeff_field_characteristic, + 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), + "Name of file containing a point set. Format is one point per line: X1 ... Xd ")( + "weight-file", po::value(&weight_file), + "Name of file containing a point weights. Format is one weigt per line: W1\n...\nWn "); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value(&output_file_diag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "field-charac,p", po::value(&coeff_field_characteristic)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::value(&min_persistence), + "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length " + "intervals"); + + po::positional_options_description pos; + pos.add("input-file", 1); + pos.add("weight-file", 2); + + po::options_description all; + all.add(visible).add(hidden); + + po::variables_map 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") || !vm.count("weight-file")) { + std::cout << std::endl; + std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; + std::cout << "of a weighted 3D Alpha complex defined on a set of input points.\n \n"; + std::cout << "The output diagram contains one bar per line, written with the convention: \n"; + std::cout << " p dim b d \n"; + std::cout << "where dim is the dimension of the homological feature,\n"; + std::cout << "b and d are respectively the birth and death of the feature and \n"; + std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl; + + std::cout << "Usage: " << argv[0] << " [options] input-file weight-file" << std::endl << std::endl; + std::cout << visible << std::endl; + std::abort(); + } +} diff --git a/src/Bottleneck_distance/utilities/CMakeLists.txt b/src/Bottleneck_distance/utilities/CMakeLists.txt index cdf76c85..063b6ae3 100644 --- a/src/Bottleneck_distance/utilities/CMakeLists.txt +++ b/src/Bottleneck_distance/utilities/CMakeLists.txt @@ -3,6 +3,9 @@ project(Bottleneck_distance_utilities) if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1) add_executable (bottleneck_read_file bottleneck_read_file.cpp) + if (TBB_FOUND) + target_link_libraries(bottleneck_read_file ${TBB_LIBRARIES}) + endif(TBB_FOUND) add_test(NAME Bottleneck_distance_utilities_Bottleneck_read_file COMMAND $ diff --git a/src/Rips_complex/utilities/README b/src/Rips_complex/utilities/README index 93ec3658..ca10424d 100644 --- a/src/Rips_complex/utilities/README +++ b/src/Rips_complex/utilities/README @@ -13,6 +13,7 @@ where `dim` is the dimension of the homological feature, `b` and `d` are respect **Allowed options** * `-h [ --help ]` Produce help message +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. * `-r [ --max-edge-length ]` (default = inf) Maximal length of an edge for the Rips complex construction. * `-d [ --cpx-dimension ]` (default = 1) Maximal dimension of the Rips complex we want to compute. * `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. @@ -45,7 +46,7 @@ outputs: ## `rips_distance_matrix_persistence` ## -Same as `rips_persistence` but taking an distance matrix as input. +Same as `rips_persistence` but taking a distance matrix as input. **Example** `rips_distance_matrix_persistence data/distance_matrix/full_square_distance_matrix.csv -r 15 -d 3 -p 3 -m 0` -- cgit v1.2.3 From 1b895dd0d76cbe13b92b68b7198f58631678200f Mon Sep 17 00:00:00 2001 From: cjamin Date: Thu, 16 Nov 2017 11:54:36 +0000 Subject: Add missing includes git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/add_utils_in_gudhi_v2@2890 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 9ddda864abff7b9509984dcbc657225ea1c80b1a --- src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp | 2 ++ src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp | 2 ++ src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp | 2 ++ src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp | 2 ++ .../example/weighted_periodic_alpha_complex_3d_persistence.cpp | 2 ++ .../example/example_alpha_shapes_3_simplex_tree_from_off_file.cpp | 2 ++ 6 files changed, 12 insertions(+) (limited to 'src/Alpha_complex') diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp index f2085ab2..28daca5e 100644 --- a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp @@ -30,6 +30,8 @@ #include #include #include +#include +#include #include #include diff --git a/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp index 7919c7fd..c69482ab 100644 --- a/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp +++ b/src/Alpha_complex/utilities/exact_alpha_complex_3d_persistence.cpp @@ -30,6 +30,8 @@ #include #include #include +#include +#include #include #include diff --git a/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp index 10b62f75..b480a9ec 100644 --- a/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp +++ b/src/Alpha_complex/utilities/periodic_alpha_complex_3d_persistence.cpp @@ -31,6 +31,8 @@ #include #include #include +#include +#include #include #include diff --git a/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp index 293170f7..91747a8c 100644 --- a/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp +++ b/src/Alpha_complex/utilities/weighted_alpha_complex_3d_persistence.cpp @@ -31,6 +31,8 @@ #include #include #include +#include +#include #include #include diff --git a/src/Persistent_cohomology/example/weighted_periodic_alpha_complex_3d_persistence.cpp b/src/Persistent_cohomology/example/weighted_periodic_alpha_complex_3d_persistence.cpp index 13634ff7..67d830dc 100644 --- a/src/Persistent_cohomology/example/weighted_periodic_alpha_complex_3d_persistence.cpp +++ b/src/Persistent_cohomology/example/weighted_periodic_alpha_complex_3d_persistence.cpp @@ -30,6 +30,8 @@ #include #include #include +#include +#include #include #include diff --git a/src/Simplex_tree/example/example_alpha_shapes_3_simplex_tree_from_off_file.cpp b/src/Simplex_tree/example/example_alpha_shapes_3_simplex_tree_from_off_file.cpp index ff2eebcb..d8289ba9 100644 --- a/src/Simplex_tree/example/example_alpha_shapes_3_simplex_tree_from_off_file.cpp +++ b/src/Simplex_tree/example/example_alpha_shapes_3_simplex_tree_from_off_file.cpp @@ -28,6 +28,8 @@ #include #include #include +#include +#include #include #include -- cgit v1.2.3