diff options
author | Vincent Rouvreau <10407034+VincentRouvreau@users.noreply.github.com> | 2019-11-19 14:56:13 +0100 |
---|---|---|
committer | GitHub <noreply@github.com> | 2019-11-19 14:56:13 +0100 |
commit | 9471d51c322ccbff12482c9c9cf2df98b9bdb50c (patch) | |
tree | a48ff64e560c37af1bc504016873efbcc07e8226 /src | |
parent | 63616eba3031334d8e2c267e54f9aa4a801baeff (diff) | |
parent | c8be137d15f40f40c4150f0c0fb9776cd680a36c (diff) |
Merge pull request #127 from VincentRouvreau/alpha_complex_with_cgal_epeck_d
Alpha complex with cgal epeck d
Diffstat (limited to 'src')
-rw-r--r-- | src/Alpha_complex/doc/Intro_alpha_complex.h | 11 | ||||
-rw-r--r-- | src/Alpha_complex/example/Alpha_complex_from_off.cpp | 13 | ||||
-rw-r--r-- | src/Alpha_complex/example/Alpha_complex_from_points.cpp | 5 | ||||
-rw-r--r-- | src/Alpha_complex/example/CMakeLists.txt | 16 | ||||
-rw-r--r-- | src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp | 65 | ||||
-rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_complex.h | 55 | ||||
-rw-r--r-- | src/Alpha_complex/test/Alpha_complex_unit_test.cpp | 13 | ||||
-rw-r--r-- | src/Alpha_complex/utilities/CMakeLists.txt | 30 | ||||
-rw-r--r-- | src/Alpha_complex/utilities/alpha_complex_persistence.cpp | 107 | ||||
-rw-r--r-- | src/Alpha_complex/utilities/alphacomplex.md | 2 | ||||
-rw-r--r-- | src/Rips_complex/example/example_rips_complex_from_csv_distance_matrix_file.cpp | 8 | ||||
-rw-r--r-- | src/Rips_complex/example/example_rips_complex_from_off_file.cpp | 8 | ||||
-rw-r--r-- | src/common/doc/main_page.md | 1 |
13 files changed, 248 insertions, 86 deletions
diff --git a/src/Alpha_complex/doc/Intro_alpha_complex.h b/src/Alpha_complex/doc/Intro_alpha_complex.h index b075d1fc..adc1378f 100644 --- a/src/Alpha_complex/doc/Intro_alpha_complex.h +++ b/src/Alpha_complex/doc/Intro_alpha_complex.h @@ -51,6 +51,17 @@ namespace alpha_complex { * - For people only interested in the topology of the \ref alpha_complex (for instance persistence), * \ref alpha_complex is equivalent to the \ref cech_complex and much smaller if you do not bound the radii. * \ref cech_complex can still make sense in higher dimension precisely because you can bound the radii. + * - Using the default `CGAL::Epeck_d` makes the construction safe. If you pass exact=true to create_complex, the + * filtration values are the exact ones converted to the filtration value type of the simplicial complex. This can be + * very slow. If you pass exact=false (the default), the filtration values are only guaranteed to have a small + * multiplicative error compared to the exact value, see <code><a class="el" target="_blank" + * href="https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html"> + * CGAL::Lazy_exact_nt<NT>::set_relative_precision_of_to_double</a></code> for details. A drawback, when computing + * persistence, is that an empty exact interval [10^12,10^12] may become a non-empty approximate interval + * [10^12,10^12+10^6]. Using `CGAL::Epick_d` makes the computations slightly faster, and the combinatorics are still + * exact, but the computation of filtration values can exceptionally be arbitrarily bad. In all cases, we still + * guarantee that the output is a valid filtration (faces have a filtration value no larger than their cofaces). + * - For performances reasons, it is advised to use `Alpha_complex` with \ref cgal ≥ 5.0.0. * * \section pointsexample Example from points * diff --git a/src/Alpha_complex/example/Alpha_complex_from_off.cpp b/src/Alpha_complex/example/Alpha_complex_from_off.cpp index d411e90a..220a66de 100644 --- a/src/Alpha_complex/example/Alpha_complex_from_off.cpp +++ b/src/Alpha_complex/example/Alpha_complex_from_off.cpp @@ -2,8 +2,6 @@ // to construct a simplex_tree from alpha complex #include <gudhi/Simplex_tree.h> -#include <CGAL/Epick_d.h> - #include <iostream> #include <string> @@ -23,22 +21,21 @@ int main(int argc, char **argv) { // ---------------------------------------------------------------------------- // Init of an alpha complex from an OFF file // ---------------------------------------------------------------------------- - using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >; - Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_name); + Gudhi::alpha_complex::Alpha_complex<> alpha_complex_from_file(off_file_name); - std::streambuf* streambufffer; + std::streambuf* streambuffer; std::ofstream ouput_file_stream; if (argc == 4) { ouput_file_stream.open(std::string(argv[3])); - streambufffer = ouput_file_stream.rdbuf(); + streambuffer = ouput_file_stream.rdbuf(); } else { - streambufffer = std::cout.rdbuf(); + streambuffer = std::cout.rdbuf(); } Gudhi::Simplex_tree<> simplex; if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) { - std::ostream output_stream(streambufffer); + std::ostream output_stream(streambuffer); // ---------------------------------------------------------------------------- // Display information about the alpha complex diff --git a/src/Alpha_complex/example/Alpha_complex_from_points.cpp b/src/Alpha_complex/example/Alpha_complex_from_points.cpp index 981aa470..6526ca3a 100644 --- a/src/Alpha_complex/example/Alpha_complex_from_points.cpp +++ b/src/Alpha_complex/example/Alpha_complex_from_points.cpp @@ -2,12 +2,13 @@ // to construct a simplex_tree from alpha complex #include <gudhi/Simplex_tree.h> -#include <CGAL/Epick_d.h> +#include <CGAL/Epeck_d.h> #include <iostream> #include <vector> -using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<2> >; +// Explicit dimension 2 Epeck_d kernel +using Kernel = CGAL::Epeck_d< CGAL::Dimension_tag<2> >; using Point = Kernel::Point_d; using Vector_of_points = std::vector<Point>; diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt index b069b443..b0337934 100644 --- a/src/Alpha_complex/example/CMakeLists.txt +++ b/src/Alpha_complex/example/CMakeLists.txt @@ -5,9 +5,12 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) target_link_libraries(Alpha_complex_example_from_points ${CGAL_LIBRARY}) add_executable ( Alpha_complex_example_from_off Alpha_complex_from_off.cpp ) target_link_libraries(Alpha_complex_example_from_off ${CGAL_LIBRARY}) + add_executable ( Alpha_complex_example_fast_from_off Fast_alpha_complex_from_off.cpp ) + target_link_libraries(Alpha_complex_example_fast_from_off ${CGAL_LIBRARY}) if (TBB_FOUND) target_link_libraries(Alpha_complex_example_from_points ${TBB_LIBRARIES}) target_link_libraries(Alpha_complex_example_from_off ${TBB_LIBRARIES}) + target_link_libraries(Alpha_complex_example_fast_from_off ${TBB_LIBRARIES}) endif() add_test(NAME Alpha_complex_example_from_points COMMAND $<TARGET_FILE:Alpha_complex_example_from_points>) @@ -16,7 +19,13 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "60.0" "${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_60.txt") add_test(NAME Alpha_complex_example_from_off_32 COMMAND $<TARGET_FILE:Alpha_complex_example_from_off> "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "32.0" "${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_32.txt") - if (DIFF_PATH) + + add_test(NAME Alpha_complex_example_fast_from_off_60 COMMAND $<TARGET_FILE:Alpha_complex_example_fast_from_off> + "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "60.0" "${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_60.txt") + add_test(NAME Alpha_complex_example_fast_from_off_32 COMMAND $<TARGET_FILE:Alpha_complex_example_fast_from_off> + "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" "32.0" "${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_32.txt") + +if (DIFF_PATH) # Do not forget to copy test results files in current binary dir file(COPY "alphaoffreader_for_doc_32.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) file(COPY "alphaoffreader_for_doc_60.txt" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) @@ -25,6 +34,11 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_60.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_60.txt) add_test(Alpha_complex_example_from_off_32_diff_files ${DIFF_PATH} ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt) + + add_test(Alpha_complex_example_fast_from_off_60_diff_files ${DIFF_PATH} + ${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_60.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_60.txt) + add_test(Alpha_complex_example_fast_from_off_32_diff_files ${DIFF_PATH} + ${CMAKE_CURRENT_BINARY_DIR}/fastalphaoffreader_result_32.txt ${CMAKE_CURRENT_BINARY_DIR}/alphaoffreader_for_doc_32.txt) endif() add_executable ( Alpha_complex_example_weighted_3d_from_points Weighted_alpha_complex_3d_from_points.cpp ) diff --git a/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp b/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp new file mode 100644 index 00000000..f181005a --- /dev/null +++ b/src/Alpha_complex/example/Fast_alpha_complex_from_off.cpp @@ -0,0 +1,65 @@ +#include <gudhi/Alpha_complex.h> +// to construct a simplex_tree from alpha complex +#include <gudhi/Simplex_tree.h> + +#include <CGAL/Epick_d.h> + +#include <iostream> +#include <string> + +void usage(int nbArgs, char * const progName) { + std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + std::cerr << "Usage: " << progName << " filename.off alpha_square_max_value [ouput_file.txt]\n"; + std::cerr << " i.e.: " << progName << " ../../data/points/alphacomplexdoc.off 60.0\n"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 3) && (argc != 4)) usage(argc, (argv[0] - 1)); + + std::string off_file_name {argv[1]}; + double alpha_square_max_value {atof(argv[2])}; + + // WARNING : CGAL::Epick_d is fast but not safe (unlike CGAL::Epeck_d) + // (i.e. when the points are on a grid) + using Fast_kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; + // ---------------------------------------------------------------------------- + // Init of an alpha complex from an OFF file + // ---------------------------------------------------------------------------- + Gudhi::alpha_complex::Alpha_complex<Fast_kernel> alpha_complex_from_file(off_file_name); + + std::streambuf* streambuffer; + std::ofstream ouput_file_stream; + + if (argc == 4) { + ouput_file_stream.open(std::string(argv[3])); + streambuffer = ouput_file_stream.rdbuf(); + } else { + streambuffer = std::cout.rdbuf(); + } + + Gudhi::Simplex_tree<> simplex; + if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) { + std::ostream output_stream(streambuffer); + + // ---------------------------------------------------------------------------- + // Display information about the alpha complex + // ---------------------------------------------------------------------------- + output_stream << "Alpha complex is of dimension " << simplex.dimension() << + " - " << simplex.num_simplices() << " simplices - " << + simplex.num_vertices() << " vertices." << std::endl; + + output_stream << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << + std::endl; + for (auto f_simplex : simplex.filtration_simplex_range()) { + output_stream << " ( "; + for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { + output_stream << vertex << " "; + } + output_stream << ") -> " << "[" << simplex.filtration(f_simplex) << "] "; + output_stream << std::endl; + } + } + ouput_file_stream.close(); + return 0; +} diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 8919cdb9..0c2569c8 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -20,11 +20,12 @@ #include <math.h> // isnan, fmax #include <CGAL/Delaunay_triangulation.h> -#include <CGAL/Epick_d.h> +#include <CGAL/Epeck_d.h> // For EXACT or SAFE version +#include <CGAL/Epick_d.h> // For FAST version #include <CGAL/Spatial_sort_traits_adapter_d.h> #include <CGAL/property_map.h> // for CGAL::Identity_property_map -#include <CGAL/NT_converter.h> #include <CGAL/version.h> // for CGAL_VERSION_NR +#include <CGAL/NT_converter.h> #include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST @@ -39,17 +40,20 @@ // Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 #if CGAL_VERSION_NR < 1041101000 -# error Alpha_complex_3d is only available for CGAL >= 4.11 +# error Alpha_complex is only available for CGAL >= 4.11 #endif #if !EIGEN_VERSION_AT_LEAST(3,1,0) -# error Alpha_complex_3d is only available for Eigen3 >= 3.1.0 installed with CGAL +# error Alpha_complex is only available for Eigen3 >= 3.1.0 installed with CGAL #endif namespace Gudhi { namespace alpha_complex { +template<typename D> struct Is_Epeck_D { static const bool value = false; }; +template<typename D> struct Is_Epeck_D<CGAL::Epeck_d<D>> { static const bool value = true; }; + /** * \class Alpha_complex Alpha_complex.h gudhi/Alpha_complex.h * \brief Alpha complex data structure. @@ -63,17 +67,31 @@ namespace alpha_complex { * * Please refer to \ref alpha_complex for examples. * - * The complex is a template class requiring an Epick_d <a target="_blank" + * The complex is a template class requiring an <a target="_blank" + * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>, + * or an <a target="_blank" + * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epick__d.html">CGAL::Epick_d</a> <a target="_blank" * href="http://doc.cgal.org/latest/Kernel_d/index.html#Chapter_dD_Geometry_Kernel">dD Geometry Kernel</a> * \cite cgal:s-gkd-15b from CGAL as template, default value is <a target="_blank" - * href="http://doc.cgal.org/latest/Kernel_d/classCGAL_1_1Epick__d.html">CGAL::Epick_d</a> + * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a> * < <a target="_blank" href="http://doc.cgal.org/latest/Kernel_23/classCGAL_1_1Dynamic__dimension__tag.html"> * CGAL::Dynamic_dimension_tag </a> > * - * \remark When Alpha_complex is constructed with an infinite value of alpha, the complex is a Delaunay complex. - * + * \remark + * - When Alpha_complex is constructed with an infinite value of alpha, the complex is a Delaunay complex. + * - Using the default `CGAL::Epeck_d` makes the construction safe. If you pass exact=true to create_complex, the + * filtration values are the exact ones converted to the filtration value type of the simplicial complex. This can be + * very slow. If you pass exact=false (the default), the filtration values are only guaranteed to have a small + * multiplicative error compared to the exact value, see <code><a class="el" target="_blank" + * href="https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html"> + * CGAL::Lazy_exact_nt<NT>::set_relative_precision_of_to_double</a></code> for details. A drawback, when computing + * persistence, is that an empty exact interval [10^12,10^12] may become a non-empty approximate interval + * [10^12,10^12+10^6]. Using `CGAL::Epick_d` makes the computations slightly faster, and the combinatorics are still + * exact, but the computation of filtration values can exceptionally be arbitrarily bad. In all cases, we still + * guarantee that the output is a valid filtration (faces have a filtration value no larger than their cofaces). + * - For performances reasons, it is advised to use `Alpha_complex` with \ref cgal ≥ 5.0.0. */ -template<class Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>> +template<class Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>> class Alpha_complex { public: // Add an int in TDS to save point index in the structure @@ -184,6 +202,12 @@ class Alpha_complex { private: template<typename InputPointRange > void init_from_range(const InputPointRange& points) { + #if CGAL_VERSION_NR < 1050000000 + if (Is_Epeck_D<Kernel>::value) + std::cerr << "It is strongly advised to use a CGAL version >= 5.0 with Epeck_d Kernel for performance reasons." + << std::endl; + #endif + auto first = std::begin(points); auto last = std::end(points); @@ -237,6 +261,8 @@ class Alpha_complex { * @param[in] complex SimplicialComplexForAlpha to be created. * @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$, and there is very * little point using anything else since it does not save time. + * @param[in] exact Exact filtration values computation. Not exact if `Kernel` is not <a target="_blank" + * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>. * * @return true if creation succeeds, false otherwise. * @@ -248,7 +274,8 @@ class Alpha_complex { template <typename SimplicialComplexForAlpha, typename Filtration_value = typename SimplicialComplexForAlpha::Filtration_value> bool create_complex(SimplicialComplexForAlpha& complex, - Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity()) { + Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity(), + bool exact = false) { // From SimplicialComplexForAlpha type required to insert into a simplicial complex (with or without subfaces). typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle; typedef typename SimplicialComplexForAlpha::Simplex_handle Simplex_handle; @@ -324,9 +351,13 @@ class Alpha_complex { if (f_simplex_dim > 0) { // squared_radius function initialization Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object(); - CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv; - alpha_complex_filtration = cv(squared_radius(pointVector.begin(), pointVector.end())); + CGAL::NT_converter<typename Geom_traits::FT, Filtration_value> cv; + auto sqrad = squared_radius(pointVector.begin(), pointVector.end()); +#if CGAL_VERSION_NR >= 1050000000 + if(exact) CGAL::exact(sqrad); +#endif + alpha_complex_filtration = cv(sqrad); } complex.assign_filtration(f_simplex, alpha_complex_filtration); #ifdef DEBUG_TRACES diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp index 01e4cee3..40b3fe09 100644 --- a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp +++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp @@ -15,6 +15,7 @@ #include <CGAL/Delaunay_triangulation.h> #include <CGAL/Epick_d.h> +#include <CGAL/Epeck_d.h> #include <cmath> // float comparison #include <limits> @@ -28,12 +29,16 @@ #include <gudhi/Unitary_tests_utils.h> // Use dynamic_dimension_tag for the user to be able to set dimension -typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kernel_d; +typedef CGAL::Epeck_d< CGAL::Dynamic_dimension_tag > Exact_kernel_d; // Use static dimension_tag for the user not to be able to set dimension -typedef CGAL::Epick_d< CGAL::Dimension_tag<3> > Kernel_s; +typedef CGAL::Epeck_d< CGAL::Dimension_tag<3> > Exact_kernel_s; +// Use dynamic_dimension_tag for the user to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Inexact_kernel_d; +// Use static dimension_tag for the user not to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dimension_tag<3> > Inexact_kernel_s; // The triangulation uses the default instantiation of the TriangulationDataStructure template parameter -typedef boost::mpl::list<Kernel_d, Kernel_s> list_of_kernel_variants; +typedef boost::mpl::list<Exact_kernel_d, Exact_kernel_s, Inexact_kernel_d, Inexact_kernel_s> list_of_kernel_variants; BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of_kernel_variants) { // ---------------------------------------------------------------------------- @@ -86,7 +91,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of } // Use static dimension_tag for the user not to be able to set dimension -typedef CGAL::Epick_d< CGAL::Dimension_tag<4> > Kernel_4; +typedef CGAL::Epeck_d< CGAL::Dimension_tag<4> > Kernel_4; typedef Kernel_4::Point_d Point_4; typedef std::vector<Point_4> Vector_4_Points; diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt index 5295f3cd..57b92942 100644 --- a/src/Alpha_complex/utilities/CMakeLists.txt +++ b/src/Alpha_complex/utilities/CMakeLists.txt @@ -7,9 +7,19 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) if (TBB_FOUND) target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES}) endif(TBB_FOUND) - add_test(NAME Alpha_complex_utilities_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence> - "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45") - + add_test(NAME Alpha_complex_utilities_safe_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "safe.pers") + add_test(NAME Alpha_complex_utilities_fast_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "fast.pers" "-f") + add_test(NAME Alpha_complex_utilities_exact_alpha_complex_persistence COMMAND $<TARGET_FILE:alpha_complex_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45" "-o" "exact.pers" "-e") + if (DIFF_PATH) + add_test(Alpha_complex_utilities_diff_exact_alpha_complex ${DIFF_PATH} + "exact.pers" "safe.pers") + add_test(Alpha_complex_utilities_diff_fast_alpha_complex ${DIFF_PATH} + "fast.pers" "safe.pers") + endif() + install(TARGETS alpha_complex_persistence DESTINATION bin) add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp) @@ -20,21 +30,21 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) add_test(NAME Alpha_complex_utilities_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" - "-p" "2" "-m" "0.45" "-o" "safe.pers") + "-p" "2" "-m" "0.45" "-o" "safe_3d.pers") add_test(NAME Alpha_complex_utilities_exact_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" - "-p" "2" "-m" "0.45" "-o" "exact.pers" "-e") + "-p" "2" "-m" "0.45" "-o" "exact_3d.pers" "-e") add_test(NAME Alpha_complex_utilities_safe_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" - "-p" "2" "-m" "0.45" "-o" "fast.pers" "-f") + "-p" "2" "-m" "0.45" "-o" "fast_3d.pers" "-f") if (DIFF_PATH) - add_test(Alpha_complex_utilities_diff_alpha_complex_3d ${DIFF_PATH} - "exact.pers" "safe.pers") - add_test(Alpha_complex_utilities_diff_alpha_complex_3d ${DIFF_PATH} - "fast.pers" "safe.pers") + add_test(Alpha_complex_utilities_diff_exact_alpha_complex_3d ${DIFF_PATH} + "exact_3d.pers" "safe_3d.pers") + add_test(Alpha_complex_utilities_diff_fast_alpha_complex_3d ${DIFF_PATH} + "fast_3d.pers" "safe_3d.pers") endif() add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp index fab7bd30..486347cc 100644 --- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -24,63 +24,84 @@ 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, bool &exact, bool &fast, + 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; + bool exact_version = false; + bool fast_version = false; 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); + program_options(argc, argv, off_file_points, exact_version, fast_version, 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<Kernel> alpha_complex_from_file(off_file_points); + if ((exact_version) && (fast_version)) { + std::cerr << "You cannot set the exact and the fast version." << std::endl; + exit(-1); + } 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(); + if (fast_version) { + // WARNING : CGAL::Epick_d is fast but not safe (unlike CGAL::Epeck_d) + // (i.e. when the points are on a grid) + using Fast_kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; + + // Init of an alpha complex from an OFF file + Gudhi::alpha_complex::Alpha_complex<Fast_kernel> alpha_complex_from_file(off_file_points); + + if (!alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) { + std::cerr << "Fast Alpha complex simplicial complex creation failed." << std::endl; + exit(-1); } - } + } else { + using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>; + // Init of an alpha complex from an OFF file + Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_points); + + if (!alpha_complex_from_file.create_complex(simplex, alpha_square_max_value, exact_version)) { + std::cerr << "Alpha complex simplicial complex creation failed." << std::endl; + exit(-1); + } + } + // ---------------------------------------------------------------------------- + // 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) { +void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast, + 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<std::string>(&off_file_points), @@ -88,6 +109,10 @@ void program_options(int argc, char *argv[], std::string &off_file_points, std:: po::options_description visible("Allowed options", 100); visible.add_options()("help,h", "produce help message")( + "exact,e", po::bool_switch(&exact), + "To activate exact version of Alpha complex (default is false, not available if fast is set)")( + "fast,f", po::bool_switch(&fast), + "To activate fast version of Alpha complex (default is false, not available if exact is set)")( "output-file,o", po::value<std::string>(&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<Filtration_value>(&alpha_square_max_value) diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md index fcd16a3b..527598a9 100644 --- a/src/Alpha_complex/utilities/alphacomplex.md +++ b/src/Alpha_complex/utilities/alphacomplex.md @@ -46,6 +46,8 @@ for the Alpha complex construction. 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.
+* `-e [ --exact ]` for the exact computation version.
+* `-f [ --fast ]` for the fast computation version.
**Example**
diff --git a/src/Rips_complex/example/example_rips_complex_from_csv_distance_matrix_file.cpp b/src/Rips_complex/example/example_rips_complex_from_csv_distance_matrix_file.cpp index 9e182f1e..b7040453 100644 --- a/src/Rips_complex/example/example_rips_complex_from_csv_distance_matrix_file.cpp +++ b/src/Rips_complex/example/example_rips_complex_from_csv_distance_matrix_file.cpp @@ -35,19 +35,19 @@ int main(int argc, char **argv) { Distance_matrix distances = Gudhi::read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_file_name); Rips_complex rips_complex_from_file(distances, threshold); - std::streambuf* streambufffer; + std::streambuf* streambuffer; std::ofstream ouput_file_stream; if (argc == 5) { ouput_file_stream.open(std::string(argv[4])); - streambufffer = ouput_file_stream.rdbuf(); + streambuffer = ouput_file_stream.rdbuf(); } else { - streambufffer = std::cout.rdbuf(); + streambuffer = std::cout.rdbuf(); } Simplex_tree stree; rips_complex_from_file.create_complex(stree, dim_max); - std::ostream output_stream(streambufffer); + std::ostream output_stream(streambuffer); // ---------------------------------------------------------------------------- // Display information about the Rips complex diff --git a/src/Rips_complex/example/example_rips_complex_from_off_file.cpp b/src/Rips_complex/example/example_rips_complex_from_off_file.cpp index de2e4ea4..36b468a7 100644 --- a/src/Rips_complex/example/example_rips_complex_from_off_file.cpp +++ b/src/Rips_complex/example/example_rips_complex_from_off_file.cpp @@ -34,19 +34,19 @@ int main(int argc, char **argv) { Gudhi::Points_off_reader<Point> off_reader(off_file_name); Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); - std::streambuf* streambufffer; + std::streambuf* streambuffer; std::ofstream ouput_file_stream; if (argc == 5) { ouput_file_stream.open(std::string(argv[4])); - streambufffer = ouput_file_stream.rdbuf(); + streambuffer = ouput_file_stream.rdbuf(); } else { - streambufffer = std::cout.rdbuf(); + streambuffer = std::cout.rdbuf(); } Simplex_tree stree; rips_complex_from_file.create_complex(stree, dim_max); - std::ostream output_stream(streambufffer); + std::ostream output_stream(streambuffer); // ---------------------------------------------------------------------------- // Display information about the Rips complex diff --git a/src/common/doc/main_page.md b/src/common/doc/main_page.md index d8cbf97f..e8d11fdf 100644 --- a/src/common/doc/main_page.md +++ b/src/common/doc/main_page.md @@ -45,6 +45,7 @@ values of the codimension 1 cofaces that make it not Gabriel otherwise. All simplices that have a filtration value strictly greater than a given alpha squared value are not inserted into the complex.<br> + For performances reasons, it is advised to use \ref cgal ≥ 5.0.0. </td> <td width="15%"> <b>Author:</b> Vincent Rouvreau<br> |