diff options
Diffstat (limited to 'src/Alpha_complex')
34 files changed, 4767 insertions, 0 deletions
diff --git a/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp b/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp new file mode 100644 index 00000000..005a712a --- /dev/null +++ b/src/Alpha_complex/benchmark/Alpha_complex_3d_benchmark.cpp @@ -0,0 +1,278 @@ +#include <gudhi/Alpha_complex_3d.h> +#include <gudhi/Alpha_complex.h> +// to construct a simplex_tree from alpha complex +#include <gudhi/Simplex_tree.h> +#include <gudhi/random_point_generators.h> +#include <gudhi/Clock.h> + +#include <iostream> +#include <string> +#include <vector> +#include <limits> // for numeric limits +#include <fstream> + +#include <CGAL/Epick_d.h> +#include <CGAL/Epeck_d.h> +#include <CGAL/Random.h> + +std::ofstream results_csv("results.csv"); + +template <typename Kernel> +void benchmark_points_on_torus_dD(const std::string& msg) { + std::cout << "+ " << msg << std::endl; + + results_csv << "\"" << msg << "\";" << std::endl; + results_csv << "\"nb_points\";" + << "\"nb_simplices\";" + << "\"alpha_creation_time(sec.)\";" + << "\"complex_creation_time(sec.)\";" << std::endl; + + using K = CGAL::Epick_d<CGAL::Dimension_tag<3>>; + for (int nb_points = 1000; nb_points <= 125000; nb_points *= 5) { + std::cout << " Alpha complex dD on torus with " << nb_points << " points." << std::endl; + std::vector<K::Point_d> points_on_torus = Gudhi::generate_points_on_torus_3D<K>(nb_points, 1.0, 0.5); + std::vector<typename Kernel::Point_d> points; + + for (auto p : points_on_torus) { + points.push_back(typename Kernel::Point_d(p.begin(), p.end())); + } + + Gudhi::Clock ac_create_clock(" benchmark_points_on_torus_dD - Alpha complex 3d creation"); + ac_create_clock.begin(); + Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_points(points); + ac_create_clock.end(); + std::cout << ac_create_clock; + + Gudhi::Simplex_tree<> complex; + Gudhi::Clock st_create_clock(" benchmark_points_on_torus_dD - complex creation"); + st_create_clock.begin(); + alpha_complex_from_points.create_complex(complex); + st_create_clock.end(); + std::cout << st_create_clock; + + results_csv << nb_points << ";" << complex.num_simplices() << ";" << ac_create_clock.num_seconds() << ";" + << st_create_clock.num_seconds() << ";" << std::endl; + + std::cout << " benchmark_points_on_torus_dD - nb simplices = " << complex.num_simplices() << std::endl; + } +} + +template <typename Alpha_complex_3d> +void benchmark_points_on_torus_3D(const std::string& msg) { + using K = CGAL::Epick_d<CGAL::Dimension_tag<3>>; + std::cout << "+ " << msg << std::endl; + + results_csv << "\"" << msg << "\";" << std::endl; + results_csv << "\"nb_points\";" + << "\"nb_simplices\";" + << "\"alpha_creation_time(sec.)\";" + << "\"complex_creation_time(sec.)\";" << std::endl; + + for (int nb_points = 1000; nb_points <= 125000; nb_points *= 5) { + std::cout << " Alpha complex 3d on torus with " << nb_points << " points." << std::endl; + std::vector<K::Point_d> points_on_torus = Gudhi::generate_points_on_torus_3D<K>(nb_points, 1.0, 0.5); + std::vector<typename Alpha_complex_3d::Point_3> points; + + for (auto p : points_on_torus) { + points.push_back(typename Alpha_complex_3d::Point_3(p[0], p[1], p[2])); + } + + Gudhi::Clock ac_create_clock(" benchmark_points_on_torus_3D - Alpha complex 3d creation"); + ac_create_clock.begin(); + Alpha_complex_3d alpha_complex_from_points(points); + ac_create_clock.end(); + std::cout << ac_create_clock; + + Gudhi::Simplex_tree<> complex; + Gudhi::Clock st_create_clock(" benchmark_points_on_torus_3D - complex creation"); + st_create_clock.begin(); + alpha_complex_from_points.create_complex(complex); + st_create_clock.end(); + std::cout << st_create_clock; + + results_csv << nb_points << ";" << complex.num_simplices() << ";" << ac_create_clock.num_seconds() << ";" + << st_create_clock.num_seconds() << ";" << std::endl; + + std::cout << " benchmark_points_on_torus_3D - nb simplices = " << complex.num_simplices() << std::endl; + } +} + +template <typename Weighted_alpha_complex_3d> +void benchmark_weighted_points_on_torus_3D(const std::string& msg) { + using K = CGAL::Epick_d<CGAL::Dimension_tag<3>>; + + std::cout << "+ " << msg << std::endl; + + results_csv << "\"" << msg << "\";" << std::endl; + results_csv << "\"nb_points\";" + << "\"nb_simplices\";" + << "\"alpha_creation_time(sec.)\";" + << "\"complex_creation_time(sec.)\";" << std::endl; + + CGAL::Random random(8); + + for (int nb_points = 1000; nb_points <= 125000; nb_points *= 5) { + std::cout << " Alpha complex 3d on torus with " << nb_points << " points." << std::endl; + std::vector<K::Point_d> points_on_torus = Gudhi::generate_points_on_torus_3D<K>(nb_points, 1.0, 0.5); + + using Point = typename Weighted_alpha_complex_3d::Point_3; + using Weighted_point = typename Weighted_alpha_complex_3d::Weighted_point_3; + + std::vector<Weighted_point> points; + + for (auto p : points_on_torus) { + points.push_back(Weighted_point(Point(p[0], p[1], p[2]), 0.9 + random.get_double(0., 0.01))); + } + + Gudhi::Clock ac_create_clock(" benchmark_weighted_points_on_torus_3D - Alpha complex 3d creation"); + ac_create_clock.begin(); + Weighted_alpha_complex_3d alpha_complex_from_points(points); + ac_create_clock.end(); + std::cout << ac_create_clock; + + Gudhi::Simplex_tree<> complex; + Gudhi::Clock st_create_clock(" benchmark_weighted_points_on_torus_3D - complex creation"); + st_create_clock.begin(); + alpha_complex_from_points.create_complex(complex); + st_create_clock.end(); + std::cout << st_create_clock; + + results_csv << nb_points << ";" << complex.num_simplices() << ";" << ac_create_clock.num_seconds() << ";" + << st_create_clock.num_seconds() << ";" << std::endl; + + std::cout << " benchmark_weighted_points_on_torus_3D - nb simplices = " << complex.num_simplices() << std::endl; + } +} + +template <typename Periodic_alpha_complex_3d> +void benchmark_periodic_points(const std::string& msg) { + std::cout << "+ " << msg << std::endl; + + results_csv << "\"" << msg << "\";" << std::endl; + results_csv << "\"nb_points\";" + << "\"nb_simplices\";" + << "\"alpha_creation_time(sec.)\";" + << "\"complex_creation_time(sec.)\";" << std::endl; + + CGAL::Random random(8); + + for (double nb_points = 10.; nb_points <= 40.; nb_points += 10.) { + std::cout << " Periodic alpha complex 3d with " << nb_points * nb_points * nb_points << " points." << std::endl; + using Point = typename Periodic_alpha_complex_3d::Point_3; + std::vector<Point> points; + + for (double i = 0; i < nb_points; i++) { + for (double j = 0; j < nb_points; j++) { + for (double k = 0; k < nb_points; k++) { + points.push_back( + Point(i + random.get_double(0., 0.1), j + random.get_double(0., 0.1), k + random.get_double(0., 0.1))); + } + } + } + + Gudhi::Clock ac_create_clock(" benchmark_periodic_points - Alpha complex 3d creation"); + ac_create_clock.begin(); + Periodic_alpha_complex_3d alpha_complex_from_points(points, 0., 0., 0., nb_points, nb_points, nb_points); + ac_create_clock.end(); + std::cout << ac_create_clock; + + Gudhi::Simplex_tree<> complex; + Gudhi::Clock st_create_clock(" benchmark_periodic_points - complex creation"); + st_create_clock.begin(); + alpha_complex_from_points.create_complex(complex); + st_create_clock.end(); + std::cout << st_create_clock; + + results_csv << nb_points * nb_points * nb_points << ";" << complex.num_simplices() << ";" + << ac_create_clock.num_seconds() << ";" << st_create_clock.num_seconds() << ";" << std::endl; + + std::cout << " benchmark_periodic_points - nb simplices = " << complex.num_simplices() << std::endl; + } +} + +template <typename Weighted_periodic_alpha_complex_3d> +void benchmark_weighted_periodic_points(const std::string& msg) { + std::cout << "+ " << msg << std::endl; + + results_csv << "\"" << msg << "\";" << std::endl; + results_csv << "\"nb_points\";" + << "\"nb_simplices\";" + << "\"alpha_creation_time(sec.)\";" + << "\"complex_creation_time(sec.)\";" << std::endl; + + CGAL::Random random(8); + + for (double nb_points = 10.; nb_points <= 40.; nb_points += 10.) { + std::cout << " Weighted periodic alpha complex 3d with " << nb_points * nb_points * nb_points << " points." + << std::endl; + + using Point = typename Weighted_periodic_alpha_complex_3d::Point_3; + using Weighted_point = typename Weighted_periodic_alpha_complex_3d::Weighted_point_3; + std::vector<Weighted_point> points; + + for (double i = 0; i < nb_points; i++) { + for (double j = 0; j < nb_points; j++) { + for (double k = 0; k < nb_points; k++) { + points.push_back(Weighted_point( + Point(i + random.get_double(0., 0.1), j + random.get_double(0., 0.1), k + random.get_double(0., 0.1)), + random.get_double(0., (nb_points * nb_points) / 64.))); + } + } + } + + Gudhi::Clock ac_create_clock(" benchmark_weighted_periodic_points - Alpha complex 3d creation"); + ac_create_clock.begin(); + Weighted_periodic_alpha_complex_3d alpha_complex_from_points(points, 0., 0., 0., nb_points, nb_points, nb_points); + ac_create_clock.end(); + std::cout << ac_create_clock; + + Gudhi::Simplex_tree<> complex; + Gudhi::Clock st_create_clock(" benchmark_weighted_periodic_points - complex creation"); + st_create_clock.begin(); + alpha_complex_from_points.create_complex(complex); + st_create_clock.end(); + std::cout << st_create_clock; + + results_csv << nb_points * nb_points * nb_points << ";" << complex.num_simplices() << ";" + << ac_create_clock.num_seconds() << ";" << st_create_clock.num_seconds() << ";" << std::endl; + + std::cout << " benchmark_weighted_periodic_points - nb simplices = " << complex.num_simplices() << std::endl; + } +} + +int main(int argc, char** argv) { + benchmark_points_on_torus_dD<CGAL::Epick_d<CGAL::Dimension_tag<3>>>("Fast static dimension version"); + benchmark_points_on_torus_dD<CGAL::Epick_d<CGAL::Dynamic_dimension_tag>>("Fast dynamic dimension version"); + benchmark_points_on_torus_dD<CGAL::Epeck_d<CGAL::Dimension_tag<3>>>("Exact static dimension version"); + benchmark_points_on_torus_dD<CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>>("Exact dynamic dimension version"); + + benchmark_points_on_torus_3D< + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::FAST, false, false>>("Fast version"); + benchmark_points_on_torus_3D< + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, false, false>>("Safe version"); + benchmark_points_on_torus_3D< + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, false, false>>("Exact version"); + + benchmark_weighted_points_on_torus_3D< + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::FAST, true, false>>("Fast version"); + benchmark_weighted_points_on_torus_3D< + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, false>>("Safe version"); + benchmark_weighted_points_on_torus_3D< + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, true, false>>("Exact version"); + + benchmark_periodic_points< + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::FAST, false, true>>("Fast version"); + benchmark_periodic_points< + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, false, true>>("Safe version"); + benchmark_periodic_points< + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, false, true>>("Exact version"); + + benchmark_weighted_periodic_points< + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::FAST, true, true>>("Fast version"); + benchmark_weighted_periodic_points< + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, true>>("Safe version"); + benchmark_weighted_periodic_points< + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, true, true>>("Exact version"); + + return 0; +} diff --git a/src/Alpha_complex/benchmark/CMakeLists.txt b/src/Alpha_complex/benchmark/CMakeLists.txt new file mode 100644 index 00000000..622963dc --- /dev/null +++ b/src/Alpha_complex/benchmark/CMakeLists.txt @@ -0,0 +1,9 @@ +project(Alpha_complex_benchmark) + +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + add_executable(Alpha_complex_3d_benchmark Alpha_complex_3d_benchmark.cpp) + target_link_libraries(Alpha_complex_3d_benchmark ${CGAL_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Alpha_complex_3d_benchmark ${TBB_LIBRARIES}) + endif() +endif () diff --git a/src/Alpha_complex/concept/SimplicialComplexForAlpha.h b/src/Alpha_complex/concept/SimplicialComplexForAlpha.h new file mode 100644 index 00000000..1c6c3b0c --- /dev/null +++ b/src/Alpha_complex/concept/SimplicialComplexForAlpha.h @@ -0,0 +1,81 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONCEPT_ALPHA_COMPLEX_SIMPLICIAL_COMPLEX_FOR_ALPHA_H_ +#define CONCEPT_ALPHA_COMPLEX_SIMPLICIAL_COMPLEX_FOR_ALPHA_H_ + +namespace Gudhi { + +namespace alpha_complex { + +/** \brief The concept SimplicialComplexForAlpha describes the requirements for a type to implement a simplicial + * complex, that can be created from a `Alpha_complex`. + */ +struct SimplicialComplexForAlpha { + /** Handle to specify a simplex. */ + typedef unspecified Simplex_handle; + /** Handle to specify a vertex. Must be a non-negative integer. */ + typedef unspecified Vertex_handle; + /** Handle to specify the simplex filtration value. */ + typedef unspecified Filtration_value; + + /** Returns the number of vertices in the simplicial complex. */ + std::size_t num_vertices(); + + /** Gets the 'simplex' dimension. */ + int dimension(Simplex_handle simplex); + + /** Assigns the 'simplex' with the given 'filtration' value. */ + int assign_filtration(Simplex_handle simplex, Filtration_value filtration); + + /** \brief Inserts a simplex with vertices from a given simplex (represented by a vector of Vertex_handle) in the + * simplicial complex with the given 'filtration' value. */ + void insert_simplex_and_subfaces(std::vector<Vertex_handle> const & vertex_range, Filtration_value filtration); + + /** Browses the simplicial complex to make the filtration non-decreasing. */ + void make_filtration_non_decreasing(); + + /** Prune the simplicial complex above 'filtration' value given as parameter. */ + void prune_above_filtration(Filtration_value filtration); + + /** \brief Iterator over vertices of a simplex. + * + * 'value type' must be 'Vertex_handle'.*/ + typedef unspecified Simplex_vertex_range; + + /** \brief Returns a range over vertices of a given simplex. */ + Simplex_vertex_range simplex_vertex_range(Simplex_handle const & simplex); + + /** \brief Iterator over the boundaries of the complex, in an arbitrary order. + * + * 'value_type' must be 'Simplex_handle'.*/ + typedef unspecified Boundary_simplex_range; + + /** \brief Returns a range over boundaries of a given simplex. */ + Boundary_simplex_range boundary_simplex_range(Simplex_handle const & simplex); + + /** \brief Iterator over the simplices of the skeleton of the complex, for a given dimension. + * + * 'value_type' must be 'Simplex_handle'. */ + typedef unspecified Skeleton_simplex_range; + /** \brief Returns a range over the simplices of the skeleton of the simplicial complex, for a given + * dimension. */ + Skeleton_simplex_range skeleton_simplex_range; + + /** \brief Return type of an insertion of a simplex + */ + typedef unspecified Insertion_result_type; +}; + +} // namespace alpha_complex + +} // namespace Gudhi + +#endif // CONCEPT_ALPHA_COMPLEX_SIMPLICIAL_COMPLEX_FOR_ALPHA_H_ diff --git a/src/Alpha_complex/concept/SimplicialComplexForAlpha3d.h b/src/Alpha_complex/concept/SimplicialComplexForAlpha3d.h new file mode 100644 index 00000000..3a6830ff --- /dev/null +++ b/src/Alpha_complex/concept/SimplicialComplexForAlpha3d.h @@ -0,0 +1,45 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2018 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CONCEPT_ALPHA_COMPLEX_SIMPLICIAL_COMPLEX_FOR_ALPHA_3D_H_ +#define CONCEPT_ALPHA_COMPLEX_SIMPLICIAL_COMPLEX_FOR_ALPHA_3D_H_ + +namespace Gudhi { + +namespace alpha_complex { + +/** \brief The concept SimplicialComplexForAlpha3d describes the requirements for a type to implement a simplicial + * complex, that can be created from a `Alpha_complex_3d`. + */ +struct SimplicialComplexForAlpha3d { + /** Handle to specify a vertex. Must be a non-negative integer. */ + typedef unspecified Vertex_handle; + /** Handle to specify the simplex filtration value. */ + typedef unspecified Filtration_value; + + /** Returns the number of vertices in the simplicial complex. */ + std::size_t num_vertices(); + + /** \brief Inserts a simplex from a given simplex (represented by a vector of Vertex_handle) in the + * simplicial complex with the given 'filtration' value. */ + void insert_simplex(std::vector<Vertex_handle> const& vertex_range, Filtration_value filtration); + + /** Browses the simplicial complex to make the filtration non-decreasing. */ + void make_filtration_non_decreasing(); + + /** Prune the simplicial complex above 'filtration' value given as parameter. */ + void prune_above_filtration(Filtration_value filtration); +}; + +} // namespace alpha_complex + +} // namespace Gudhi + +#endif // CONCEPT_ALPHA_COMPLEX_SIMPLICIAL_COMPLEX_FOR_ALPHA_3D_H_ diff --git a/src/Alpha_complex/doc/COPYRIGHT b/src/Alpha_complex/doc/COPYRIGHT new file mode 100644 index 00000000..61f17f6d --- /dev/null +++ b/src/Alpha_complex/doc/COPYRIGHT @@ -0,0 +1,12 @@ +The files of this directory are part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. +See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + +Author(s): Vincent Rouvreau + +Copyright (C) 2015 Inria + +This gives everyone the freedoms to use openFrameworks in any context: +commercial or non-commercial, public or private, open or closed source. + +You should have received a copy of the MIT License along with this program. +If not, see https://opensource.org/licenses/MIT.
\ No newline at end of file diff --git a/src/Alpha_complex/doc/Intro_alpha_complex.h b/src/Alpha_complex/doc/Intro_alpha_complex.h new file mode 100644 index 00000000..b075d1fc --- /dev/null +++ b/src/Alpha_complex/doc/Intro_alpha_complex.h @@ -0,0 +1,190 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef DOC_ALPHA_COMPLEX_INTRO_ALPHA_COMPLEX_H_ +#define DOC_ALPHA_COMPLEX_INTRO_ALPHA_COMPLEX_H_ + +// needs namespace for Doxygen to link on classes +namespace Gudhi { +// needs namespace for Doxygen to link on classes +namespace alpha_complex { + +/** \defgroup alpha_complex Alpha complex + * + * \author Vincent Rouvreau + * + * @{ + * + * \section definition Definition + * + * Alpha_complex is a <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a> + * constructed from the finite cells of a Delaunay Triangulation. + * + * The filtration value of each simplex is computed as the square of the circumradius of the simplex if the + * circumsphere is empty (the simplex is then said to be Gabriel), and as the minimum of the filtration + * 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. + * + * \image html "alpha_complex_representation.png" "Alpha-complex representation" + * + * Alpha_complex is constructing a <a target="_blank" + * href="http://doc.cgal.org/latest/Triangulation/index.html#Chapter_Triangulations">Delaunay Triangulation</a> + * \cite cgal:hdj-t-15b from <a target="_blank" href="http://www.cgal.org/">CGAL</a> (the Computational Geometry + * Algorithms Library \cite cgal:eb-15b) and is able to create a `SimplicialComplexForAlpha`. + * + * The complex is a template class requiring an Epick_d <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 parameter. + * + * \remark + * - When the simplicial complex is constructed with an infinite value of alpha, the complex is a Delaunay + * 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. + * + * \section pointsexample Example from points + * + * This example builds the Delaunay triangulation from the given points in a 2D static kernel, and creates a + * `Simplex_tree` with it. + * + * Then, it is asked to display information about the simplicial complex. + * + * \include Alpha_complex/Alpha_complex_from_points.cpp + * + * When launching: + * + * \code $> ./Alpha_complex_example_from_points + * \endcode + * + * the program output is: + * + * \include Alpha_complex/alphaoffreader_for_doc_60.txt + * + * \section createcomplexalgorithm Create complex algorithm + * + * \subsection datastructure Data structure + * + * In order to create the simplicial complex, first, it is built from the cells of the Delaunay Triangulation. + * The filtration values are set to NaN, which stands for unknown value. + * + * In example, : + * \image html "alpha_complex_doc.png" "Simplicial complex structure construction example" + * + * \subsection filtrationcomputation Filtration value computation algorithm + * <br> + * \f$ + * \textbf{for } \text{i : dimension } \rightarrow 0 \textbf{ do}\\ + * \quad \textbf{for all } \sigma \text{ of dimension i}\\ + * \quad\quad \textbf{if } \text{filtration(} \sigma ) \text{ is NaN} \textbf{ then}\\ + * \quad\quad\quad \text{filtration(} \sigma ) = \alpha^2( \sigma )\\ + * \quad\quad \textbf{end if}\\ + * \quad\quad \textbf{for all } \tau \text{ face of } \sigma \textbf{ do}\quad\quad + * \textit{// propagate alpha filtration value}\\ + * \quad\quad\quad \textbf{if } \text{filtration(} \tau ) \text{ is not NaN} \textbf{ then}\\ + * \quad\quad\quad\quad \text{filtration(} \tau \text{) = min( filtration(} \tau \text{), filtration(} \sigma + * \text{) )}\\ + * \quad\quad\quad \textbf{else}\\ + * \quad\quad\quad\quad \textbf{if } \tau \text{ is not Gabriel for } \sigma \textbf{ then}\\ + * \quad\quad\quad\quad\quad \text{filtration(} \tau \text{) = filtration(} \sigma \text{)}\\ + * \quad\quad\quad\quad \textbf{end if}\\ + * \quad\quad\quad \textbf{end if}\\ + * \quad\quad \textbf{end for}\\ + * \quad \textbf{end for}\\ + * \textbf{end for}\\ + * \text{make_filtration_non_decreasing()}\\ + * \text{prune_above_filtration()}\\ + * \f$ + * + * \subsubsection dimension2 Dimension 2 + * + * From the example above, it means the algorithm looks into each triangle ([0,1,2], [0,2,4], [1,2,3], ...), + * computes the filtration value of the triangle, and then propagates the filtration value as described + * here : + * \image html "alpha_complex_doc_420.png" "Filtration value propagation example" + * + * \subsubsection dimension1 Dimension 1 + * + * Then, the algorithm looks into each edge ([0,1], [0,2], [1,2], ...), + * computes the filtration value of the edge (in this case, propagation will have no effect). + * + * \subsubsection dimension0 Dimension 0 + * + * Finally, the algorithm looks into each vertex ([0], [1], [2], [3], [4], [5] and [6]) and + * sets the filtration value (0 in case of a vertex - propagation will have no effect). + * + * \subsubsection nondecreasing Non decreasing filtration values + * + * As the squared radii computed by CGAL are an approximation, it might happen that these alpha squared values do not + * quite define a proper filtration (i.e. non-decreasing with respect to inclusion). + * We fix that up by calling `SimplicialComplexForAlpha::make_filtration_non_decreasing()`. + * + * \subsubsection pruneabove Prune above given filtration value + * + * The simplex tree is pruned from the given maximum alpha squared value (cf. + * `SimplicialComplexForAlpha::prune_above_filtration()`). + * In the following example, the value is given by the user as argument of the program. + * + * + * \section offexample Example from OFF file + * + * This example builds the Delaunay triangulation in a dynamic kernel, and initializes the alpha complex with it. + * + * + * Then, it is asked to display information about the alpha complex. + * + * \include Alpha_complex/Alpha_complex_from_off.cpp + * + * When launching: + * + * \code $> ./Alpha_complex_example_from_off ../../data/points/alphacomplexdoc.off 32.0 + * \endcode + * + * the program output is: + * + * \include Alpha_complex/alphaoffreader_for_doc_32.txt + * + * + * \section weighted3dexample 3d specific example + * + * A specific module for Alpha complex is available in 3d (cf. Alpha_complex_3d) and allows to construct standard, + * weighted, periodic or weighted and periodic versions of alpha complexes. Alpha values computation can be + * Gudhi::alpha_complex::complexity::FAST, Gudhi::alpha_complex::complexity::SAFE (default value) or + * Gudhi::alpha_complex::complexity::EXACT. + * + * This example builds the CGAL 3d weighted alpha shapes from a small molecule, and initializes the alpha complex with + * it. This example is taken from <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title13">CGAL 3d + * weighted alpha shapes</a>. + * + * Then, it is asked to display information about the alpha complex. + * + * \include Alpha_complex/Weighted_alpha_complex_3d_from_points.cpp + * + * When launching: + * + * \code $> ./Alpha_complex_example_weighted_3d_from_points + * \endcode + * + * the program output is: + * + * \include Alpha_complex/weightedalpha3dfrompoints_for_doc.txt + * + */ +/** @} */ // end defgroup alpha_complex + +} // namespace alpha_complex + +namespace alphacomplex = alpha_complex; + +} // namespace Gudhi + +#endif // DOC_ALPHA_COMPLEX_INTRO_ALPHA_COMPLEX_H_ diff --git a/src/Alpha_complex/doc/alpha_complex_doc.ipe b/src/Alpha_complex/doc/alpha_complex_doc.ipe new file mode 100644 index 00000000..71e5ce6c --- /dev/null +++ b/src/Alpha_complex/doc/alpha_complex_doc.ipe @@ -0,0 +1,296 @@ +<?xml version="1.0"?> +<!DOCTYPE ipe SYSTEM "ipe.dtd"> +<ipe version="70107" creator="Ipe 7.1.10"> +<info created="D:20150603143945" modified="D:20160921180211"/> +<ipestyle name="basic"> +<symbol name="arrow/arc(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/farc(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="mark/circle(sx)" transformations="translations"> +<path fill="sym-stroke"> +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e +</path> +</symbol> +<symbol name="mark/disk(sx)" transformations="translations"> +<path fill="sym-stroke"> +0.6 0 0 0.6 0 0 e +</path> +</symbol> +<symbol name="mark/fdisk(sfx)" transformations="translations"> +<group> +<path fill="sym-fill"> +0.5 0 0 0.5 0 0 e +</path> +<path fill="sym-stroke" fillrule="eofill"> +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e +</path> +</group> +</symbol> +<symbol name="mark/box(sx)" transformations="translations"> +<path fill="sym-stroke" fillrule="eofill"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h +</path> +</symbol> +<symbol name="mark/square(sx)" transformations="translations"> +<path fill="sym-stroke"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +</path> +</symbol> +<symbol name="mark/fsquare(sfx)" transformations="translations"> +<group> +<path fill="sym-fill"> +-0.5 -0.5 m +0.5 -0.5 l +0.5 0.5 l +-0.5 0.5 l +h +</path> +<path fill="sym-stroke" fillrule="eofill"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h +</path> +</group> +</symbol> +<symbol name="mark/cross(sx)" transformations="translations"> +<group> +<path fill="sym-stroke"> +-0.43 -0.57 m +0.57 0.43 l +0.43 0.57 l +-0.57 -0.43 l +h +</path> +<path fill="sym-stroke"> +-0.43 0.57 m +0.57 -0.43 l +0.43 -0.57 l +-0.57 0.43 l +h +</path> +</group> +</symbol> +<symbol name="arrow/fnormal(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/pointed(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/fpointed(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/linear(spx)"> +<path stroke="sym-stroke" pen="sym-pen"> +-1 0.333 m +0 0 l +-1 -0.333 l +</path> +</symbol> +<symbol name="arrow/fdouble(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/double(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h +</path> +</symbol> +<pen name="heavier" value="0.8"/> +<pen name="fat" value="1.2"/> +<pen name="ultrafat" value="2"/> +<symbolsize name="large" value="5"/> +<symbolsize name="small" value="2"/> +<symbolsize name="tiny" value="1.1"/> +<arrowsize name="large" value="10"/> +<arrowsize name="small" value="5"/> +<arrowsize name="tiny" value="3"/> +<color name="red" value="1 0 0"/> +<color name="green" value="0 1 0"/> +<color name="blue" value="0 0 1"/> +<color name="yellow" value="1 1 0"/> +<color name="orange" value="1 0.647 0"/> +<color name="gold" value="1 0.843 0"/> +<color name="purple" value="0.627 0.125 0.941"/> +<color name="gray" value="0.745"/> +<color name="brown" value="0.647 0.165 0.165"/> +<color name="navy" value="0 0 0.502"/> +<color name="pink" value="1 0.753 0.796"/> +<color name="seagreen" value="0.18 0.545 0.341"/> +<color name="turquoise" value="0.251 0.878 0.816"/> +<color name="violet" value="0.933 0.51 0.933"/> +<color name="darkblue" value="0 0 0.545"/> +<color name="darkcyan" value="0 0.545 0.545"/> +<color name="darkgray" value="0.663"/> +<color name="darkgreen" value="0 0.392 0"/> +<color name="darkmagenta" value="0.545 0 0.545"/> +<color name="darkorange" value="1 0.549 0"/> +<color name="darkred" value="0.545 0 0"/> +<color name="lightblue" value="0.678 0.847 0.902"/> +<color name="lightcyan" value="0.878 1 1"/> +<color name="lightgray" value="0.827"/> +<color name="lightgreen" value="0.565 0.933 0.565"/> +<color name="lightyellow" value="1 1 0.878"/> +<dashstyle name="dashed" value="[4] 0"/> +<dashstyle name="dotted" value="[1 3] 0"/> +<dashstyle name="dash dotted" value="[4 2 1 2] 0"/> +<dashstyle name="dash dot dotted" value="[4 2 1 2 1 2] 0"/> +<textsize name="large" value="\large"/> +<textsize name="small" value="\small"/> +<textsize name="tiny" value="\tiny"/> +<textsize name="Large" value="\Large"/> +<textsize name="LARGE" value="\LARGE"/> +<textsize name="huge" value="\huge"/> +<textsize name="Huge" value="\Huge"/> +<textsize name="footnote" value="\footnotesize"/> +<textstyle name="center" begin="\begin{center}" end="\end{center}"/> +<textstyle name="itemize" begin="\begin{itemize}" end="\end{itemize}"/> +<textstyle name="item" begin="\begin{itemize}\item{}" end="\end{itemize}"/> +<gridsize name="4 pts" value="4"/> +<gridsize name="8 pts (~3 mm)" value="8"/> +<gridsize name="16 pts (~6 mm)" value="16"/> +<gridsize name="32 pts (~12 mm)" value="32"/> +<gridsize name="10 pts (~3.5 mm)" value="10"/> +<gridsize name="20 pts (~7 mm)" value="20"/> +<gridsize name="14 pts (~5 mm)" value="14"/> +<gridsize name="28 pts (~10 mm)" value="28"/> +<gridsize name="56 pts (~20 mm)" value="56"/> +<anglesize name="90 deg" value="90"/> +<anglesize name="60 deg" value="60"/> +<anglesize name="45 deg" value="45"/> +<anglesize name="30 deg" value="30"/> +<anglesize name="22.5 deg" value="22.5"/> +<tiling name="falling" angle="-60" step="4" width="1"/> +<tiling name="rising" angle="30" step="4" width="1"/> +</ipestyle> +<page> +<layer name="alpha"/> +<view layers="alpha" active="alpha"/> +<path layer="alpha" matrix="1 0 0 1 -240 0" stroke="darkcyan"> +320 580 m +350 520 l +290 530 l +320 580 l +320 580 l +</path> +<path matrix="1 0 0 1 -240 0" stroke="darkcyan"> +320 580 m +280 660 l +290 530 l +320 580 l +320 580 l +</path> +<path matrix="1 0 0 1 -240 0" stroke="darkcyan"> +320 580 m +370 580 l +350 520 l +320 580 l +</path> +<text matrix="1 0 0 1 -260 0" transformations="translations" pos="380 530" stroke="darkcyan" type="label" width="118.196" height="8.307" depth="2.32" valign="baseline" size="large">Delaunay triangulation</text> +<text matrix="1 0 0 1 -242.155 -3.50128" transformations="translations" pos="282.952 524.893" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">0</text> +<text matrix="1 0 0 1 -240 0" transformations="translations" pos="352.708 510.349" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">1</text> +<text matrix="1 0 0 1 -240 0" transformations="translations" pos="310.693 578.759" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">2</text> +<text matrix="1 0 0 1 -240 0" transformations="translations" pos="375.332 578.49" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">3</text> +<text matrix="1 0 0 1 -240 0" transformations="translations" pos="272.179 660.635" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">4</text> +<text matrix="1 0 0 1 -239.3 -10.1537" transformations="translations" pos="296.419 724.197" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">5</text> +<text matrix="1 0 0 1 -240 0" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">6</text> +<path matrix="1 0 0 1 -240 0" stroke="darkcyan"> +280 660 m +300 710 l +370 690 l +280 660 l +</path> +<path matrix="1 0 0 1 -240 0" stroke="darkcyan"> +320 580 m +370 690 l +370 580 l +320 580 l +</path> +<path matrix="1 0 0 1 -240 0" stroke="darkcyan"> +280 660 m +370 690 l +320 580 l +280 660 l +</path> +<text matrix="1 0 0 1 76 36" transformations="translations" pos="180 620" stroke="black" type="label" width="153.148" height="6.926" depth="1.93" valign="baseline">Simplicial complex data structure :</text> +<use matrix="1 0 0 1 -239.3 -10.1537" name="mark/fdisk(sfx)" pos="300 720" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -240 0" name="mark/fdisk(sfx)" pos="370 690" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -240 0" name="mark/fdisk(sfx)" pos="280 660" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -240 0" name="mark/fdisk(sfx)" pos="320 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -240 0" name="mark/fdisk(sfx)" pos="370 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -240 0" name="mark/fdisk(sfx)" pos="350 520" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -240 0" name="mark/fdisk(sfx)" pos="290 530" size="normal" stroke="black" fill="white"/> +<text matrix="1 0 0 1 -20 -32" transformations="translations" pos="288 672" stroke="black" type="label" width="148.582" height="7.473" depth="2.49" valign="baseline">insert simplex and subfaces [0,1,2]</text> +<text matrix="1 0 0 1 -20 -56" transformations="translations" pos="288 672" stroke="black" type="label" width="148.582" height="7.473" depth="2.49" valign="baseline">insert simplex and subfaces [1,2,3]</text> +<text matrix="1 0 0 1 -20 -44" transformations="translations" pos="288 672" stroke="black" type="label" width="148.582" height="7.473" depth="2.49" valign="baseline">insert simplex and subfaces [0,2,4]</text> +<text matrix="1 0 0 1 -20 -68" transformations="translations" pos="288 672" stroke="black" type="label" width="148.582" height="7.473" depth="2.49" valign="baseline">insert simplex and subfaces [2,3,6]</text> +<text matrix="1 0 0 1 -20 -80" transformations="translations" pos="288 672" stroke="black" type="label" width="148.582" height="7.473" depth="2.49" valign="baseline">insert simplex and subfaces [2,4,6]</text> +<text matrix="1 0 0 1 -20 -92" transformations="translations" pos="288 672" stroke="black" type="label" width="148.582" height="7.473" depth="2.49" valign="baseline">insert simplex and subfaces [4,5,6]</text> +</page> +</ipe> diff --git a/src/Alpha_complex/doc/alpha_complex_doc.png b/src/Alpha_complex/doc/alpha_complex_doc.png Binary files differnew file mode 100644 index 00000000..170bae80 --- /dev/null +++ b/src/Alpha_complex/doc/alpha_complex_doc.png diff --git a/src/Alpha_complex/doc/alpha_complex_doc_420.ipe b/src/Alpha_complex/doc/alpha_complex_doc_420.ipe new file mode 100644 index 00000000..5d1d29d4 --- /dev/null +++ b/src/Alpha_complex/doc/alpha_complex_doc_420.ipe @@ -0,0 +1,514 @@ +<?xml version="1.0"?> +<!DOCTYPE ipe SYSTEM "ipe.dtd"> +<ipe version="70107" creator="Ipe 7.1.10"> +<info created="D:20150603143945" modified="D:20151130095019"/> +<ipestyle name="basic"> +<symbol name="arrow/arc(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/farc(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="mark/circle(sx)" transformations="translations"> +<path fill="sym-stroke"> +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e +</path> +</symbol> +<symbol name="mark/disk(sx)" transformations="translations"> +<path fill="sym-stroke"> +0.6 0 0 0.6 0 0 e +</path> +</symbol> +<symbol name="mark/fdisk(sfx)" transformations="translations"> +<group> +<path fill="sym-fill"> +0.5 0 0 0.5 0 0 e +</path> +<path fill="sym-stroke" fillrule="eofill"> +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e +</path> +</group> +</symbol> +<symbol name="mark/box(sx)" transformations="translations"> +<path fill="sym-stroke" fillrule="eofill"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h +</path> +</symbol> +<symbol name="mark/square(sx)" transformations="translations"> +<path fill="sym-stroke"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +</path> +</symbol> +<symbol name="mark/fsquare(sfx)" transformations="translations"> +<group> +<path fill="sym-fill"> +-0.5 -0.5 m +0.5 -0.5 l +0.5 0.5 l +-0.5 0.5 l +h +</path> +<path fill="sym-stroke" fillrule="eofill"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h +</path> +</group> +</symbol> +<symbol name="mark/cross(sx)" transformations="translations"> +<group> +<path fill="sym-stroke"> +-0.43 -0.57 m +0.57 0.43 l +0.43 0.57 l +-0.57 -0.43 l +h +</path> +<path fill="sym-stroke"> +-0.43 0.57 m +0.57 -0.43 l +0.43 -0.57 l +-0.57 0.43 l +h +</path> +</group> +</symbol> +<symbol name="arrow/fnormal(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/pointed(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/fpointed(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/linear(spx)"> +<path stroke="sym-stroke" pen="sym-pen"> +-1 0.333 m +0 0 l +-1 -0.333 l +</path> +</symbol> +<symbol name="arrow/fdouble(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/double(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h +</path> +</symbol> +<pen name="heavier" value="0.8"/> +<pen name="fat" value="1.2"/> +<pen name="ultrafat" value="2"/> +<symbolsize name="large" value="5"/> +<symbolsize name="small" value="2"/> +<symbolsize name="tiny" value="1.1"/> +<arrowsize name="large" value="10"/> +<arrowsize name="small" value="5"/> +<arrowsize name="tiny" value="3"/> +<color name="red" value="1 0 0"/> +<color name="green" value="0 1 0"/> +<color name="blue" value="0 0 1"/> +<color name="yellow" value="1 1 0"/> +<color name="orange" value="1 0.647 0"/> +<color name="gold" value="1 0.843 0"/> +<color name="purple" value="0.627 0.125 0.941"/> +<color name="gray" value="0.745"/> +<color name="brown" value="0.647 0.165 0.165"/> +<color name="navy" value="0 0 0.502"/> +<color name="pink" value="1 0.753 0.796"/> +<color name="seagreen" value="0.18 0.545 0.341"/> +<color name="turquoise" value="0.251 0.878 0.816"/> +<color name="violet" value="0.933 0.51 0.933"/> +<color name="darkblue" value="0 0 0.545"/> +<color name="darkcyan" value="0 0.545 0.545"/> +<color name="darkgray" value="0.663"/> +<color name="darkgreen" value="0 0.392 0"/> +<color name="darkmagenta" value="0.545 0 0.545"/> +<color name="darkorange" value="1 0.549 0"/> +<color name="darkred" value="0.545 0 0"/> +<color name="lightblue" value="0.678 0.847 0.902"/> +<color name="lightcyan" value="0.878 1 1"/> +<color name="lightgray" value="0.827"/> +<color name="lightgreen" value="0.565 0.933 0.565"/> +<color name="lightyellow" value="1 1 0.878"/> +<dashstyle name="dashed" value="[4] 0"/> +<dashstyle name="dotted" value="[1 3] 0"/> +<dashstyle name="dash dotted" value="[4 2 1 2] 0"/> +<dashstyle name="dash dot dotted" value="[4 2 1 2 1 2] 0"/> +<textsize name="large" value="\large"/> +<textsize name="small" value="\small"/> +<textsize name="tiny" value="\tiny"/> +<textsize name="Large" value="\Large"/> +<textsize name="LARGE" value="\LARGE"/> +<textsize name="huge" value="\huge"/> +<textsize name="Huge" value="\Huge"/> +<textsize name="footnote" value="\footnotesize"/> +<textstyle name="center" begin="\begin{center}" end="\end{center}"/> +<textstyle name="itemize" begin="\begin{itemize}" end="\end{itemize}"/> +<textstyle name="item" begin="\begin{itemize}\item{}" end="\end{itemize}"/> +<gridsize name="4 pts" value="4"/> +<gridsize name="8 pts (~3 mm)" value="8"/> +<gridsize name="16 pts (~6 mm)" value="16"/> +<gridsize name="32 pts (~12 mm)" value="32"/> +<gridsize name="10 pts (~3.5 mm)" value="10"/> +<gridsize name="20 pts (~7 mm)" value="20"/> +<gridsize name="14 pts (~5 mm)" value="14"/> +<gridsize name="28 pts (~10 mm)" value="28"/> +<gridsize name="56 pts (~20 mm)" value="56"/> +<anglesize name="90 deg" value="90"/> +<anglesize name="60 deg" value="60"/> +<anglesize name="45 deg" value="45"/> +<anglesize name="30 deg" value="30"/> +<anglesize name="22.5 deg" value="22.5"/> +<tiling name="falling" angle="-60" step="4" width="1"/> +<tiling name="rising" angle="30" step="4" width="1"/> +</ipestyle> +<page> +<layer name="alpha"/> +<view layers="alpha" active="alpha"/> +<path layer="alpha" matrix="1 0 0 1 0 80" stroke="lightgray"> +320 580 m +350 520 l +290 530 l +320 580 l +320 580 l +</path> +<path matrix="1 0 0 1 0 80" stroke="darkcyan" pen="heavier"> +320 580 m +280 660 l +290 530 l +320 580 l +320 580 l +</path> +<path matrix="1 0 0 1 0 80" stroke="lightgray"> +320 580 m +370 580 l +350 520 l +320 580 l +</path> +<text matrix="1 0 0 1 0 80" transformations="translations" pos="380 530" stroke="darkcyan" type="label" width="54.628" height="8.965" depth="2.99" valign="baseline" size="large">Cell [4,2,0]</text> +<text matrix="1 0 0 1 -2.15463 76.4987" transformations="translations" pos="282.952 524.893" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">0</text> +<text matrix="1 0 0 1 0 80" transformations="translations" pos="352.708 510.349" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">1</text> +<text matrix="1 0 0 1 0 80" transformations="translations" pos="310.693 578.759" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">2</text> +<text matrix="1 0 0 1 0 80" transformations="translations" pos="375.332 578.49" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">3</text> +<text matrix="1 0 0 1 0 80" transformations="translations" pos="272.179 660.635" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">4</text> +<text matrix="1 0 0 1 0.700256 69.8463" transformations="translations" pos="296.419 724.197" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">5</text> +<text matrix="1 0 0 1 0 80" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">6</text> +<path matrix="1 0 0 1 0 80" stroke="lightgray"> +280 660 m +300 710 l +370 690 l +280 660 l +</path> +<path matrix="1 0 0 1 0 80" stroke="lightgray"> +320 580 m +370 690 l +370 580 l +320 580 l +</path> +<path matrix="1 0 0 1 0 80" stroke="lightgray"> +280 660 m +370 690 l +320 580 l +280 660 l +</path> +<path matrix="1 0 0 1 0 80" stroke="darkcyan"> +77.2727 0 0 77.2727 243.636 591.818 e +</path> +<path matrix="1 0 0 1 0 80" stroke="darkcyan"> +243.428 591.569 m +186.061 643.28 l +</path> +<text matrix="1 0 0 1 0 80" transformations="translations" pos="212.724 627.389" stroke="darkcyan" type="label" width="18.785" height="4.294" depth="1.49" valign="baseline">$\alpha_{420}$</text> +<path matrix="1 0 0 1 -264 -162" stroke="lightgray"> +320 580 m +350 520 l +290 530 l +320 580 l +320 580 l +</path> +<path matrix="1 0 0 1 -264 -162" stroke="lightgray"> +320 580 m +280 660 l +290 530 l +320 580 l +320 580 l +</path> +<path matrix="1 0 0 1 -264 -162" stroke="lightgray"> +320 580 m +370 580 l +350 520 l +320 580 l +</path> +<text matrix="0.582962 0 0 1 -211.265 -209.555" transformations="translations" pos="380 530" stroke="darkcyan" type="label" width="231.798" height="8.965" depth="2.99" valign="baseline" size="large">[2,0] is Gabriel $\rightarrow$ $\alpha_{20}$ is not$\\$ +modified (NaN) +</text> +<text matrix="1 0 0 1 -266.155 -165.501" transformations="translations" pos="282.952 524.893" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">0</text> +<text matrix="1 0 0 1 -264 -162" transformations="translations" pos="310.693 578.759" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">2</text> +<text matrix="1 0 0 1 -264 -162" transformations="translations" pos="375.332 578.49" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">3</text> +<text matrix="1 0 0 1 -264 -172" transformations="translations" pos="272.179 660.635" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">4</text> +<text matrix="1 0 0 1 -263.3 -172.154" transformations="translations" pos="296.419 724.197" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">5</text> +<text matrix="1 0 0 1 -264 -162" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">6</text> +<path matrix="1 0 0 1 -264 -162" stroke="lightgray"> +280 660 m +300 710 l +370 690 l +280 660 l +</path> +<path matrix="1 0 0 1 -264 -162" stroke="lightgray"> +320 580 m +370 690 l +370 580 l +320 580 l +</path> +<path matrix="1 0 0 1 -264 -162" stroke="lightgray"> +280 660 m +370 690 l +320 580 l +280 660 l +</path> +<text matrix="1 0 0 1 -166.834 -240.52" transformations="translations" pos="212.724 627.389" stroke="darkcyan" type="label" width="14.814" height="4.294" depth="1.49" valign="baseline">$\alpha_{20}$</text> +<path matrix="1 0 0 1 -264 -162" stroke="darkcyan" pen="heavier"> +290 530 m +320 580 l +</path> +<path matrix="1 0 0 1 -264 -162" stroke="darkcyan"> +29.1548 0 0 29.1548 305 555 e +</path> +<path matrix="1 0 0 1 -264 -162" stroke="darkcyan"> +304.883 555.015 m +334.509 555.015 l +</path> +<path matrix="1 0 0 1 -37.2997 -163.65" stroke="lightgray"> +320 580 m +350 520 l +290 530 l +320 580 l +320 580 l +</path> +<path matrix="1 0 0 1 -38 -164" stroke="lightgray"> +320 580 m +280 660 l +290 530 l +320 580 l +320 580 l +</path> +<path matrix="1 0 0 1 -38 -164" stroke="lightgray"> +320 580 m +370 580 l +350 520 l +320 580 l +</path> +<text matrix="1 0 0 1 -199.21 -189.117" transformations="translations" pos="380 530" stroke="darkred" type="label" width="168.308" height="8.965" depth="2.99" valign="baseline" size="large">[0,4] is not Gabriel $\rightarrow$ $\alpha_{40} = \alpha_{420}$</text> +<text matrix="1 0 0 1 -40.1546 -167.501" transformations="translations" pos="282.952 524.893" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">0</text> +<text matrix="1 0 0 1 -38 -164" transformations="translations" pos="375.332 578.49" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">3</text> +<text matrix="1 0 0 1 -37.2997 -174.154" transformations="translations" pos="296.419 724.197" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">5</text> +<text matrix="1 0 0 1 -38 -164" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">6</text> +<path matrix="1 0 0 1 -38 -164" stroke="lightgray"> +280 660 m +300 710 l +370 690 l +280 660 l +</path> +<path matrix="1 0 0 1 -38 -164" stroke="lightgray"> +320 580 m +370 690 l +370 580 l +320 580 l +</path> +<path matrix="1 0 0 1 -38 -164" stroke="lightgray"> +280 660 m +370 690 l +320 580 l +280 660 l +</path> +<text matrix="1 0 0 1 52.4654 -193.97" transformations="translations" pos="212.724 627.389" stroke="darkcyan" type="label" width="14.814" height="4.294" depth="1.49" valign="baseline">$\alpha_{40}$</text> +<path matrix="1 0 0 1 -38 -164" stroke="darkcyan" pen="heavier"> +290 530 m +280 660 l +</path> +<path matrix="1 0 0 1 126 -162" stroke="lightgray"> +320 580 m +350 520 l +290 530 l +320 580 l +320 580 l +</path> +<path matrix="1 0 0 1 126 -162" stroke="lightgray"> +320 580 m +280 660 l +290 530 l +320 580 l +320 580 l +</path> +<path matrix="1 0 0 1 126 -162" stroke="lightgray"> +320 580 m +370 580 l +350 520 l +320 580 l +</path> +<text matrix="1 0 0 1 123.845 -165.501" transformations="translations" pos="282.952 524.893" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">0</text> +<text matrix="1 0 0 1 126 -162" transformations="translations" pos="352.708 510.349" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">1</text> +<text matrix="1 0 0 1 126 -162" transformations="translations" pos="310.693 578.759" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">2</text> +<text matrix="1 0 0 1 126 -162" transformations="translations" pos="375.332 578.49" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">3</text> +<text matrix="1 0 0 1 126.7 -172.154" transformations="translations" pos="296.419 724.197" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">5</text> +<text matrix="1 0 0 1 126 -162" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">6</text> +<path matrix="1 0 0 1 126 -162" stroke="lightgray"> +280 660 m +300 710 l +370 690 l +280 660 l +</path> +<path matrix="1 0 0 1 126 -162" stroke="lightgray"> +320 580 m +370 690 l +370 580 l +320 580 l +</path> +<path matrix="1 0 0 1 126 -162" stroke="lightgray"> +280 660 m +370 690 l +320 580 l +280 660 l +</path> +<text matrix="1 0 0 1 225.859 -165.729" transformations="translations" pos="212.724 627.389" stroke="darkcyan" type="label" width="14.814" height="4.294" depth="1.49" valign="baseline">$\alpha_{42}$</text> +<text matrix="1 0 0 1 122 -164" transformations="translations" pos="272.179 660.635" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">4</text> +<path stroke="darkcyan" pen="heavier"> +406.093 497.775 m +446.094 418.092 l +</path> +<path stroke="darkcyan"> +44.5799 0 0 44.5799 425.934 457.774 e +</path> +<path stroke="darkcyan"> +425.854 457.774 m +470.795 457.774 l +</path> +<text matrix="1 0 0 1 -48.9756 -209.799" transformations="translations" pos="380 530" stroke="darkcyan" type="label" width="231.798" height="8.965" depth="2.99" valign="baseline" size="large">[2,4] is Gabriel $\rightarrow$ $\alpha_{42}$ is not modified (NaN) +</text> +<path stroke="darkblue" arrow="normal/normal"> +205.028 596.091 m +110.946 544.02 l +</path> +<path stroke="darkblue" arrow="normal/normal"> +280.768 588.99 m +280.768 547.57 l +</path> +<path stroke="darkblue" arrow="normal/normal"> +341.123 594.316 m +413.904 554.079 l +</path> +<text matrix="1 0 0 1 39.645 -2.36686" transformations="translations" pos="199.703 569.464" stroke="darkblue" type="label" width="93.206" height="7.473" depth="2.49" valign="baseline">For all faces of [4,2,0]</text> +<text matrix="1 0 0 1 -93.391 2.68003" transformations="translations" pos="104.437 300.174" stroke="black" type="label" width="208.621" height="6.926" depth="1.93" valign="baseline">N.B. : is Gabriel on a single point has no sense.</text> +<text matrix="1 0 0 1 -36.9231 10" transformations="translations" pos="48 784" stroke="black" type="label" width="118.324" height="7.473" depth="2.49" valign="baseline">Dimension =2 - $\sigma$ = [4,2,0]</text> +<path stroke="darkcyan"> +247.333 430.892 m +311.764 430.892 l +</path> +<use matrix="1 0 0 1 0.700256 69.8463" name="mark/fdisk(sfx)" pos="300 720" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 0 80" name="mark/fdisk(sfx)" pos="370 690" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 0 80" name="mark/fdisk(sfx)" pos="280 660" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 0 80" name="mark/fdisk(sfx)" pos="243.636 591.818" size="normal" stroke="darkcyan" fill="white"/> +<use matrix="1 0 0 1 0 80" name="mark/fdisk(sfx)" pos="370 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 0 80" name="mark/fdisk(sfx)" pos="350 520" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 0 80" name="mark/fdisk(sfx)" pos="320 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 0 80" name="mark/fdisk(sfx)" pos="290 530" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -263.3 -172.154" name="mark/fdisk(sfx)" pos="300 720" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -264 -162" name="mark/fdisk(sfx)" pos="280 660" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -264 -162" name="mark/fdisk(sfx)" pos="370 690" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -264 -162" name="mark/fdisk(sfx)" pos="370 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -264 -162" name="mark/fdisk(sfx)" pos="350 520" size="normal" stroke="black" fill="white"/> +<text matrix="1 0 0 1 -264 -162" transformations="translations" pos="352.708 510.349" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">1</text> +<use matrix="1 0 0 1 -264 -162" name="mark/fdisk(sfx)" pos="305 555" size="normal" stroke="darkcyan" fill="white"/> +<use matrix="1 0 0 1 -264 -162" name="mark/fdisk(sfx)" pos="290 530" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -264 -162" name="mark/fdisk(sfx)" pos="320 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -37.2997 -174.154" name="mark/fdisk(sfx)" pos="300 720" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -38 -164" name="mark/fdisk(sfx)" pos="370 690" size="normal" stroke="black" fill="white"/> +<text matrix="1 0 0 1 -38 -164" transformations="translations" pos="272.179 660.635" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">4</text> +<use name="mark/fdisk(sfx)" pos="247 431" size="normal" stroke="darkcyan" fill="white"/> +<use matrix="1 0 0 1 -38 -164" name="mark/fdisk(sfx)" pos="350 520" size="normal" stroke="black" fill="white"/> +<text matrix="1 0 0 1 -38 -164" transformations="translations" pos="352.708 510.349" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">1</text> +<use matrix="1 0 0 1 -38 -164" name="mark/fdisk(sfx)" pos="370 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -38 -164" name="mark/fdisk(sfx)" pos="320 580" size="normal" stroke="darkred" fill="white"/> +<text matrix="1 0 0 1 -38 -164" transformations="translations" pos="310.693 578.759" stroke="darkred" type="label" width="4.981" height="6.42" depth="0" valign="baseline">2</text> +<path matrix="1 0 0 1 -38 -164" stroke="darkred" pen="heavier"> +65.192 0 0 65.192 285 595 e +</path> +<use matrix="1 0 0 1 -38 -164" name="mark/fdisk(sfx)" pos="290 530" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 126 -162" name="mark/fdisk(sfx)" pos="290 530" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 126.7 -172.154" name="mark/fdisk(sfx)" pos="300 720" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 126 -162" name="mark/fdisk(sfx)" pos="370 690" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 126 -162" name="mark/fdisk(sfx)" pos="280 660" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="425.934 457.774" size="normal" stroke="darkcyan" fill="white"/> +<use matrix="1 0 0 1 126 -162" name="mark/fdisk(sfx)" pos="320 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 126 -162" name="mark/fdisk(sfx)" pos="370 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 126 -162" name="mark/fdisk(sfx)" pos="350 520" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -38 -164" name="mark/fdisk(sfx)" pos="280 660" size="normal" stroke="black" fill="white"/> +</page> +</ipe> diff --git a/src/Alpha_complex/doc/alpha_complex_doc_420.png b/src/Alpha_complex/doc/alpha_complex_doc_420.png Binary files differnew file mode 100644 index 00000000..ef7187f7 --- /dev/null +++ b/src/Alpha_complex/doc/alpha_complex_doc_420.png diff --git a/src/Alpha_complex/doc/alpha_complex_representation.ipe b/src/Alpha_complex/doc/alpha_complex_representation.ipe new file mode 100644 index 00000000..e8096b93 --- /dev/null +++ b/src/Alpha_complex/doc/alpha_complex_representation.ipe @@ -0,0 +1,321 @@ +<?xml version="1.0"?> +<!DOCTYPE ipe SYSTEM "ipe.dtd"> +<ipe version="70107" creator="Ipe 7.1.10"> +<info created="D:20150603143945" modified="D:20160404172133"/> +<ipestyle name="basic"> +<symbol name="arrow/arc(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/farc(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="mark/circle(sx)" transformations="translations"> +<path fill="sym-stroke"> +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e +</path> +</symbol> +<symbol name="mark/disk(sx)" transformations="translations"> +<path fill="sym-stroke"> +0.6 0 0 0.6 0 0 e +</path> +</symbol> +<symbol name="mark/fdisk(sfx)" transformations="translations"> +<group> +<path fill="sym-fill"> +0.5 0 0 0.5 0 0 e +</path> +<path fill="sym-stroke" fillrule="eofill"> +0.6 0 0 0.6 0 0 e +0.4 0 0 0.4 0 0 e +</path> +</group> +</symbol> +<symbol name="mark/box(sx)" transformations="translations"> +<path fill="sym-stroke" fillrule="eofill"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h +</path> +</symbol> +<symbol name="mark/square(sx)" transformations="translations"> +<path fill="sym-stroke"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +</path> +</symbol> +<symbol name="mark/fsquare(sfx)" transformations="translations"> +<group> +<path fill="sym-fill"> +-0.5 -0.5 m +0.5 -0.5 l +0.5 0.5 l +-0.5 0.5 l +h +</path> +<path fill="sym-stroke" fillrule="eofill"> +-0.6 -0.6 m +0.6 -0.6 l +0.6 0.6 l +-0.6 0.6 l +h +-0.4 -0.4 m +0.4 -0.4 l +0.4 0.4 l +-0.4 0.4 l +h +</path> +</group> +</symbol> +<symbol name="mark/cross(sx)" transformations="translations"> +<group> +<path fill="sym-stroke"> +-0.43 -0.57 m +0.57 0.43 l +0.43 0.57 l +-0.57 -0.43 l +h +</path> +<path fill="sym-stroke"> +-0.43 0.57 m +0.57 -0.43 l +0.43 -0.57 l +-0.57 0.43 l +h +</path> +</group> +</symbol> +<symbol name="arrow/fnormal(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/pointed(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/fpointed(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-0.8 0 l +-1 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/linear(spx)"> +<path stroke="sym-stroke" pen="sym-pen"> +-1 0.333 m +0 0 l +-1 -0.333 l +</path> +</symbol> +<symbol name="arrow/fdouble(spx)"> +<path stroke="sym-stroke" fill="white" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h +</path> +</symbol> +<symbol name="arrow/double(spx)"> +<path stroke="sym-stroke" fill="sym-stroke" pen="sym-pen"> +0 0 m +-1 0.333 l +-1 -0.333 l +h +-1 0 m +-2 0.333 l +-2 -0.333 l +h +</path> +</symbol> +<pen name="heavier" value="0.8"/> +<pen name="fat" value="1.2"/> +<pen name="ultrafat" value="2"/> +<symbolsize name="large" value="5"/> +<symbolsize name="small" value="2"/> +<symbolsize name="tiny" value="1.1"/> +<arrowsize name="large" value="10"/> +<arrowsize name="small" value="5"/> +<arrowsize name="tiny" value="3"/> +<color name="red" value="1 0 0"/> +<color name="green" value="0 1 0"/> +<color name="blue" value="0 0 1"/> +<color name="yellow" value="1 1 0"/> +<color name="orange" value="1 0.647 0"/> +<color name="gold" value="1 0.843 0"/> +<color name="purple" value="0.627 0.125 0.941"/> +<color name="gray" value="0.745"/> +<color name="brown" value="0.647 0.165 0.165"/> +<color name="navy" value="0 0 0.502"/> +<color name="pink" value="1 0.753 0.796"/> +<color name="seagreen" value="0.18 0.545 0.341"/> +<color name="turquoise" value="0.251 0.878 0.816"/> +<color name="violet" value="0.933 0.51 0.933"/> +<color name="darkblue" value="0 0 0.545"/> +<color name="darkcyan" value="0 0.545 0.545"/> +<color name="darkgray" value="0.663"/> +<color name="darkgreen" value="0 0.392 0"/> +<color name="darkmagenta" value="0.545 0 0.545"/> +<color name="darkorange" value="1 0.549 0"/> +<color name="darkred" value="0.545 0 0"/> +<color name="lightblue" value="0.678 0.847 0.902"/> +<color name="lightcyan" value="0.878 1 1"/> +<color name="lightgray" value="0.827"/> +<color name="lightgreen" value="0.565 0.933 0.565"/> +<color name="lightyellow" value="1 1 0.878"/> +<dashstyle name="dashed" value="[4] 0"/> +<dashstyle name="dotted" value="[1 3] 0"/> +<dashstyle name="dash dotted" value="[4 2 1 2] 0"/> +<dashstyle name="dash dot dotted" value="[4 2 1 2 1 2] 0"/> +<textsize name="large" value="\large"/> +<textsize name="small" value="\small"/> +<textsize name="tiny" value="\tiny"/> +<textsize name="Large" value="\Large"/> +<textsize name="LARGE" value="\LARGE"/> +<textsize name="huge" value="\huge"/> +<textsize name="Huge" value="\Huge"/> +<textsize name="footnote" value="\footnotesize"/> +<textstyle name="center" begin="\begin{center}" end="\end{center}"/> +<textstyle name="itemize" begin="\begin{itemize}" end="\end{itemize}"/> +<textstyle name="item" begin="\begin{itemize}\item{}" end="\end{itemize}"/> +<gridsize name="4 pts" value="4"/> +<gridsize name="8 pts (~3 mm)" value="8"/> +<gridsize name="16 pts (~6 mm)" value="16"/> +<gridsize name="32 pts (~12 mm)" value="32"/> +<gridsize name="10 pts (~3.5 mm)" value="10"/> +<gridsize name="20 pts (~7 mm)" value="20"/> +<gridsize name="14 pts (~5 mm)" value="14"/> +<gridsize name="28 pts (~10 mm)" value="28"/> +<gridsize name="56 pts (~20 mm)" value="56"/> +<anglesize name="90 deg" value="90"/> +<anglesize name="60 deg" value="60"/> +<anglesize name="45 deg" value="45"/> +<anglesize name="30 deg" value="30"/> +<anglesize name="22.5 deg" value="22.5"/> +<tiling name="falling" angle="-60" step="4" width="1"/> +<tiling name="rising" angle="30" step="4" width="1"/> +</ipestyle> +<page> +<layer name="alpha"/> +<view layers="alpha" active="alpha"/> +<path layer="alpha" fill="lightblue"> +109.771 601.912 m +159.595 601.797 l +140.058 541.915 l +h +</path> +<path fill="lightblue"> +79.8776 552.169 m +109.756 601.699 l +139.812 542.209 l +h +</path> +<path fill="lightblue"> +69.8453 682.419 m +159.925 712.208 l +90.12 732.039 l +h +</path> +<text matrix="1 0 0 1 -230.178 22.1775" transformations="translations" pos="380 530" stroke="seagreen" type="label" width="76.735" height="8.307" depth="2.32" valign="baseline" size="large">Alpha complex</text> +<text matrix="1 0 0 1 -212.333 18.6762" transformations="translations" pos="282.952 524.893" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">0</text> +<text matrix="1 0 0 1 -210.178 22.1775" transformations="translations" pos="352.708 510.349" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">1</text> +<text matrix="1 0 0 1 -210.178 22.1775" transformations="translations" pos="310.693 578.759" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">2</text> +<text matrix="1 0 0 1 -210.178 22.1775" transformations="translations" pos="375.332 578.49" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">3</text> +<text matrix="1 0 0 1 -210.178 22.1775" transformations="translations" pos="272.179 660.635" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">4</text> +<text matrix="1 0 0 1 -209.478 12.0238" transformations="translations" pos="296.419 724.197" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">5</text> +<text matrix="1 0 0 1 -210.178 22.1775" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">6</text> +<path matrix="1 0 0 1 31.9779 -58.7483" stroke="darkgray"> +58.1341 0 0 58.1341 218.925 692.601 e +</path> +<path matrix="1 0 0 1 29.8225 22.1775" stroke="black" pen="heavier"> +60 710 m +40 660 l +</path> +<path matrix="1 0 0 1 29.8225 22.1775" stroke="black" pen="heavier"> +40 660 m +130 690 l +</path> +<path matrix="1 0 0 1 29.8225 22.1775" stroke="black" pen="heavier"> +130 690 m +60 710 l +</path> +<path matrix="1 0 0 1 29.8225 22.1775" stroke="black" pen="heavier"> +40 660 m +80 580 l +</path> +<path matrix="1 0 0 1 29.8225 22.1775" stroke="black" pen="heavier"> +80 580 m +130 580 l +130 580 l +</path> +<path matrix="1 0 0 1 29.8225 22.1775" stroke="black" pen="heavier"> +130 580 m +110 520 l +</path> +<path matrix="1 0 0 1 29.8225 22.1775" stroke="black" pen="heavier"> +110 520 m +50 530 l +50 530 l +50 530 l +</path> +<path matrix="1 0 0 1 29.8225 22.1775" stroke="black" pen="heavier"> +50 530 m +80 580 l +</path> +<path matrix="1 0 0 1 29.8225 22.1775" stroke="black" pen="heavier"> +130 580 m +130 690 l +</path> +<use matrix="1 0 0 1 142.618 -109.867" name="mark/fdisk(sfx)" pos="108.285 743.72" size="normal" stroke="darkgray" fill="white"/> +<path matrix="1 0 0 1 142.618 -109.867" stroke="darkgray"> +108.275 743.531 m +166.45 743.531 l +</path> +<text matrix="1 0 0 1 142.618 -109.867" transformations="translations" pos="127.397 746.763" stroke="darkgray" type="label" width="6.41" height="4.289" depth="0" valign="baseline">$\alpha$</text> +<use matrix="1 0 0 1 -209.478 12.0238" name="mark/fdisk(sfx)" pos="300 720" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 22.1775" name="mark/fdisk(sfx)" pos="280 660" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 22.1775" name="mark/fdisk(sfx)" pos="370 690" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 22.1775" name="mark/fdisk(sfx)" pos="370 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 22.1775" name="mark/fdisk(sfx)" pos="290 530" size="normal" stroke="black" fill="white"/> +<path matrix="1 0 0 1 -40 -8" stroke="black" pen="heavier"> +150.038 609.9 m +179.929 549.727 l +</path> +<use matrix="1 0 0 1 -210.178 22.1775" name="mark/fdisk(sfx)" pos="320 580" size="normal" stroke="black" fill="white"/> +<use matrix="1 0 0 1 -210.178 22.1775" name="mark/fdisk(sfx)" pos="350 520" size="normal" stroke="black" fill="white"/> +</page> +</ipe> diff --git a/src/Alpha_complex/doc/alpha_complex_representation.png b/src/Alpha_complex/doc/alpha_complex_representation.png Binary files differnew file mode 100644 index 00000000..7b81cd69 --- /dev/null +++ b/src/Alpha_complex/doc/alpha_complex_representation.png diff --git a/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp b/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp new file mode 100644 index 00000000..0e359a27 --- /dev/null +++ b/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp @@ -0,0 +1,56 @@ +#include <gudhi/Alpha_complex_3d.h> +// to construct a simplex_tree from alpha complex +#include <gudhi/Simplex_tree.h> + +#include <iostream> +#include <string> +#include <vector> +#include <limits> // for numeric limits + +using Alpha_complex_3d = Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, false, false>; +using Point = Alpha_complex_3d::Point_3; +using Vector_of_points = std::vector<Point>; + +int main(int argc, char **argv) { + if (argc != 1) { + std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n"; + std::cerr << "Usage: " << (argv[0] - 1) << " \n"; + exit(-1); // ----- >> + } + + // ---------------------------------------------------------------------------- + // Init of a list of points from a small molecule + // ---------------------------------------------------------------------------- + Vector_of_points points; + points.push_back(Point(1, -1, -1)); + points.push_back(Point(-1, 1, -1)); + points.push_back(Point(-1, -1, 1)); + points.push_back(Point(1, 1, 1)); + points.push_back(Point(2, 2, 2)); + + // ---------------------------------------------------------------------------- + // Init of an alpha complex from the list of points + // ---------------------------------------------------------------------------- + Alpha_complex_3d alpha_complex_from_points(points); + + Gudhi::Simplex_tree<> simplex; + if (alpha_complex_from_points.create_complex(simplex)) { + // ---------------------------------------------------------------------------- + // Display information about the alpha complex + // ---------------------------------------------------------------------------- + std::cout << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() + << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + for (auto f_simplex : simplex.filtration_simplex_range()) { + std::cout << " ( "; + for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " + << "[" << simplex.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + } + return 0; +} diff --git a/src/Alpha_complex/example/Alpha_complex_from_off.cpp b/src/Alpha_complex/example/Alpha_complex_from_off.cpp new file mode 100644 index 00000000..d411e90a --- /dev/null +++ b/src/Alpha_complex/example/Alpha_complex_from_off.cpp @@ -0,0 +1,63 @@ +#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])}; + + // ---------------------------------------------------------------------------- + // 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); + + std::streambuf* streambufffer; + std::ofstream ouput_file_stream; + + if (argc == 4) { + ouput_file_stream.open(std::string(argv[3])); + streambufffer = ouput_file_stream.rdbuf(); + } else { + streambufffer = std::cout.rdbuf(); + } + + Gudhi::Simplex_tree<> simplex; + if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) { + std::ostream output_stream(streambufffer); + + // ---------------------------------------------------------------------------- + // 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/example/Alpha_complex_from_points.cpp b/src/Alpha_complex/example/Alpha_complex_from_points.cpp new file mode 100644 index 00000000..981aa470 --- /dev/null +++ b/src/Alpha_complex/example/Alpha_complex_from_points.cpp @@ -0,0 +1,52 @@ +#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 <vector> + +using Kernel = CGAL::Epick_d< CGAL::Dimension_tag<2> >; +using Point = Kernel::Point_d; +using Vector_of_points = std::vector<Point>; + +int main() { + // ---------------------------------------------------------------------------- + // Init of a list of points + // ---------------------------------------------------------------------------- + Vector_of_points points; + points.push_back(Point(1.0, 1.0)); + points.push_back(Point(7.0, 0.0)); + points.push_back(Point(4.0, 6.0)); + points.push_back(Point(9.0, 6.0)); + points.push_back(Point(0.0, 14.0)); + points.push_back(Point(2.0, 19.0)); + points.push_back(Point(9.0, 17.0)); + + // ---------------------------------------------------------------------------- + // Init of an alpha complex from the list of points + // ---------------------------------------------------------------------------- + Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_points(points); + + Gudhi::Simplex_tree<> simplex; + if (alpha_complex_from_points.create_complex(simplex)) { + // ---------------------------------------------------------------------------- + // Display information about the alpha complex + // ---------------------------------------------------------------------------- + std::cout << "Alpha complex is of dimension " << simplex.dimension() << + " - " << simplex.num_simplices() << " simplices - " << + simplex.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + for (auto f_simplex : simplex.filtration_simplex_range()) { + std::cout << " ( "; + for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " << "[" << simplex.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + } + return 0; +} diff --git a/src/Alpha_complex/example/CMakeLists.txt b/src/Alpha_complex/example/CMakeLists.txt new file mode 100644 index 00000000..b069b443 --- /dev/null +++ b/src/Alpha_complex/example/CMakeLists.txt @@ -0,0 +1,46 @@ +project(Alpha_complex_examples) + +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + add_executable ( Alpha_complex_example_from_points Alpha_complex_from_points.cpp ) + 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}) + if (TBB_FOUND) + target_link_libraries(Alpha_complex_example_from_points ${TBB_LIBRARIES}) + target_link_libraries(Alpha_complex_example_from_off ${TBB_LIBRARIES}) + endif() + + add_test(NAME Alpha_complex_example_from_points COMMAND $<TARGET_FILE:Alpha_complex_example_from_points>) + + add_test(NAME Alpha_complex_example_from_off_60 COMMAND $<TARGET_FILE:Alpha_complex_example_from_off> + "${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) + # 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}/) + + add_test(Alpha_complex_example_from_off_60_diff_files ${DIFF_PATH} + ${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) + endif() + + add_executable ( Alpha_complex_example_weighted_3d_from_points Weighted_alpha_complex_3d_from_points.cpp ) + target_link_libraries(Alpha_complex_example_weighted_3d_from_points ${CGAL_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Alpha_complex_example_weighted_3d_from_points ${TBB_LIBRARIES}) + endif() + add_test(NAME Alpha_complex_example_weighted_3d_from_points + COMMAND $<TARGET_FILE:Alpha_complex_example_weighted_3d_from_points>) + + add_executable ( Alpha_complex_example_3d_from_points Alpha_complex_3d_from_points.cpp ) + target_link_libraries(Alpha_complex_example_3d_from_points ${CGAL_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Alpha_complex_example_3d_from_points ${TBB_LIBRARIES}) + endif() + add_test(NAME Alpha_complex_example_3d_from_points + COMMAND $<TARGET_FILE:Alpha_complex_example_3d_from_points>) + +endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) diff --git a/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp b/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp new file mode 100644 index 00000000..ac11b68c --- /dev/null +++ b/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp @@ -0,0 +1,52 @@ +#include <gudhi/Alpha_complex_3d.h> +// to construct a simplex_tree from alpha complex +#include <gudhi/Simplex_tree.h> + +#include <iostream> +#include <string> +#include <vector> +#include <limits> // for numeric limits + +// Complexity = FAST, weighted = true, periodic = false +using Weighted_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, false>; +using Point = Weighted_alpha_complex_3d::Point_3; +using Weighted_point = Weighted_alpha_complex_3d::Weighted_point_3; + +int main(int argc, char **argv) { + // ---------------------------------------------------------------------------- + // Init of a list of points and weights from a small molecule + // ---------------------------------------------------------------------------- + std::vector<Weighted_point> weighted_points; + weighted_points.push_back(Weighted_point(Point(1, -1, -1), 4.)); + weighted_points.push_back(Weighted_point(Point(-1, 1, -1), 4.)); + weighted_points.push_back(Weighted_point(Point(-1, -1, 1), 4.)); + weighted_points.push_back(Weighted_point(Point(1, 1, 1), 4.)); + weighted_points.push_back(Weighted_point(Point(2, 2, 2), 1.)); + + // ---------------------------------------------------------------------------- + // Init of an alpha complex from the list of points + // ---------------------------------------------------------------------------- + Weighted_alpha_complex_3d alpha_complex_from_points(weighted_points); + + Gudhi::Simplex_tree<> simplex; + if (alpha_complex_from_points.create_complex(simplex)) { + // ---------------------------------------------------------------------------- + // Display information about the alpha complex + // ---------------------------------------------------------------------------- + std::cout << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() + << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + for (auto f_simplex : simplex.filtration_simplex_range()) { + std::cout << " ( "; + for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " + << "[" << simplex.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + } + return 0; +} diff --git a/src/Alpha_complex/example/alphaoffreader_for_doc_32.txt b/src/Alpha_complex/example/alphaoffreader_for_doc_32.txt new file mode 100644 index 00000000..13183e86 --- /dev/null +++ b/src/Alpha_complex/example/alphaoffreader_for_doc_32.txt @@ -0,0 +1,22 @@ +Alpha complex is of dimension 2 - 20 simplices - 7 vertices. +Iterator on alpha complex simplices in the filtration order, with [filtration value]: + ( 0 ) -> [0] + ( 1 ) -> [0] + ( 2 ) -> [0] + ( 3 ) -> [0] + ( 4 ) -> [0] + ( 5 ) -> [0] + ( 6 ) -> [0] + ( 3 2 ) -> [6.25] + ( 5 4 ) -> [7.25] + ( 2 0 ) -> [8.5] + ( 1 0 ) -> [9.25] + ( 3 1 ) -> [10] + ( 2 1 ) -> [11.25] + ( 3 2 1 ) -> [12.5] + ( 2 1 0 ) -> [12.9959] + ( 6 5 ) -> [13.25] + ( 4 2 ) -> [20] + ( 6 4 ) -> [22.7367] + ( 6 5 4 ) -> [22.7367] + ( 6 3 ) -> [30.25] diff --git a/src/Alpha_complex/example/alphaoffreader_for_doc_60.txt b/src/Alpha_complex/example/alphaoffreader_for_doc_60.txt new file mode 100644 index 00000000..71f29a00 --- /dev/null +++ b/src/Alpha_complex/example/alphaoffreader_for_doc_60.txt @@ -0,0 +1,27 @@ +Alpha complex is of dimension 2 - 25 simplices - 7 vertices. +Iterator on alpha complex simplices in the filtration order, with [filtration value]: + ( 0 ) -> [0] + ( 1 ) -> [0] + ( 2 ) -> [0] + ( 3 ) -> [0] + ( 4 ) -> [0] + ( 5 ) -> [0] + ( 6 ) -> [0] + ( 3 2 ) -> [6.25] + ( 5 4 ) -> [7.25] + ( 2 0 ) -> [8.5] + ( 1 0 ) -> [9.25] + ( 3 1 ) -> [10] + ( 2 1 ) -> [11.25] + ( 3 2 1 ) -> [12.5] + ( 2 1 0 ) -> [12.9959] + ( 6 5 ) -> [13.25] + ( 4 2 ) -> [20] + ( 6 4 ) -> [22.7367] + ( 6 5 4 ) -> [22.7367] + ( 6 3 ) -> [30.25] + ( 6 2 ) -> [36.5] + ( 6 3 2 ) -> [36.5] + ( 6 4 2 ) -> [37.2449] + ( 4 0 ) -> [59.7107] + ( 4 2 0 ) -> [59.7107] diff --git a/src/Alpha_complex/example/weightedalpha3dfrompoints_for_doc.txt b/src/Alpha_complex/example/weightedalpha3dfrompoints_for_doc.txt new file mode 100644 index 00000000..7a09998d --- /dev/null +++ b/src/Alpha_complex/example/weightedalpha3dfrompoints_for_doc.txt @@ -0,0 +1,31 @@ +Alpha complex is of dimension 3 - 29 simplices - 5 vertices. +Iterator on alpha complex simplices in the filtration order, with [filtration value]: + ( 0 ) -> [-4] + ( 1 ) -> [-4] + ( 2 ) -> [-4] + ( 3 ) -> [-4] + ( 1 0 ) -> [-2] + ( 2 0 ) -> [-2] + ( 2 1 ) -> [-2] + ( 3 0 ) -> [-2] + ( 3 1 ) -> [-2] + ( 3 2 ) -> [-2] + ( 2 1 0 ) -> [-1.33333] + ( 3 1 0 ) -> [-1.33333] + ( 3 2 0 ) -> [-1.33333] + ( 3 2 1 ) -> [-1.33333] + ( 3 2 1 0 ) -> [-1] + ( 4 ) -> [-1] + ( 4 2 ) -> [-1] + ( 4 0 ) -> [23] + ( 4 1 ) -> [23] + ( 4 2 0 ) -> [23] + ( 4 2 1 ) -> [23] + ( 4 3 ) -> [23] + ( 4 3 2 ) -> [23] + ( 4 1 0 ) -> [95] + ( 4 2 1 0 ) -> [95] + ( 4 3 0 ) -> [95] + ( 4 3 1 ) -> [95] + ( 4 3 2 0 ) -> [95] + ( 4 3 2 1 ) -> [95] diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h new file mode 100644 index 00000000..8919cdb9 --- /dev/null +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -0,0 +1,432 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL and Eigen3 + * - YYYY/MM Author: Description of the modification + */ + +#ifndef ALPHA_COMPLEX_H_ +#define ALPHA_COMPLEX_H_ + +#include <gudhi/Debug_utils.h> +// to construct Alpha_complex from a OFF file of points +#include <gudhi/Points_off_io.h> + +#include <stdlib.h> +#include <math.h> // isnan, fmax + +#include <CGAL/Delaunay_triangulation.h> +#include <CGAL/Epick_d.h> +#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 <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST + +#include <iostream> +#include <vector> +#include <string> +#include <limits> // NaN +#include <map> +#include <utility> // std::pair +#include <stdexcept> +#include <numeric> // for std::iota + +// 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 +#endif + +#if !EIGEN_VERSION_AT_LEAST(3,1,0) +# error Alpha_complex_3d is only available for Eigen3 >= 3.1.0 installed with CGAL +#endif + +namespace Gudhi { + +namespace alpha_complex { + +/** + * \class Alpha_complex Alpha_complex.h gudhi/Alpha_complex.h + * \brief Alpha complex data structure. + * + * \ingroup alpha_complex + * + * \details + * The data structure is constructing a CGAL Delaunay triangulation (for more informations on CGAL Delaunay + * triangulation, please refer to the corresponding chapter in page http://doc.cgal.org/latest/Triangulation/) from a + * range of points or from an OFF file (cf. Points_off_reader). + * + * Please refer to \ref alpha_complex for examples. + * + * The complex is a template class requiring an Epick_d <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> + * < <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. + * + */ +template<class Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>> +class Alpha_complex { + public: + // Add an int in TDS to save point index in the structure + typedef CGAL::Triangulation_data_structure<typename Kernel::Dimension, + CGAL::Triangulation_vertex<Kernel, std::ptrdiff_t>, + CGAL::Triangulation_full_cell<Kernel> > TDS; + /** \brief A Delaunay triangulation of a set of points in \f$ \mathbb{R}^D\f$.*/ + typedef CGAL::Delaunay_triangulation<Kernel, TDS> Delaunay_triangulation; + + /** \brief A point in Euclidean space.*/ + typedef typename Kernel::Point_d Point_d; + /** \brief Geometric traits class that provides the geometric types and predicates needed by Delaunay + * triangulations.*/ + typedef Kernel Geom_traits; + + private: + typedef typename Kernel::Compute_squared_radius_d Squared_Radius; + typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel; + typedef typename Kernel::Point_dimension_d Point_Dimension; + + // Type required to compute squared radius, or side of bounded sphere on a vector of points. + typedef typename std::vector<Point_d> Vector_of_CGAL_points; + + // Vertex_iterator type from CGAL. + typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator; + + // size_type type from CGAL. + typedef typename Delaunay_triangulation::size_type size_type; + + // Map type to switch from simplex tree vertex handle to CGAL vertex iterator. + typedef typename std::map< std::size_t, CGAL_vertex_iterator > Vector_vertex_iterator; + + private: + /** \brief Vertex iterator vector to switch from simplex tree vertex handle to CGAL vertex iterator. + * Vertex handles are inserted sequentially, starting at 0.*/ + Vector_vertex_iterator vertex_handle_to_iterator_; + /** \brief Pointer on the CGAL Delaunay triangulation.*/ + Delaunay_triangulation* triangulation_; + /** \brief Kernel for triangulation_ functions access.*/ + Kernel kernel_; + + public: + /** \brief Alpha_complex constructor from an OFF file name. + * + * Uses the Points_off_reader to construct the Delaunay triangulation required to initialize + * the Alpha_complex. + * + * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous. + * + * @param[in] off_file_name OFF file [path and] name. + */ + Alpha_complex(const std::string& off_file_name) + : triangulation_(nullptr) { + Gudhi::Points_off_reader<Point_d> off_reader(off_file_name); + if (!off_reader.is_valid()) { + std::cerr << "Alpha_complex - Unable to read file " << off_file_name << "\n"; + exit(-1); // ----- >> + } + + init_from_range(off_reader.get_point_cloud()); + } + + /** \brief Alpha_complex constructor from a list of points. + * + * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous. + * + * @param[in] points Range of points to triangulate. Points must be in Kernel::Point_d + * + * The type InputPointRange must be a range for which std::begin and + * std::end return input iterators on a Kernel::Point_d. + */ + template<typename InputPointRange > + Alpha_complex(const InputPointRange& points) + : triangulation_(nullptr) { + init_from_range(points); + } + + /** \brief Alpha_complex destructor deletes the Delaunay triangulation. + */ + ~Alpha_complex() { + delete triangulation_; + } + + // Forbid copy/move constructor/assignment operator + Alpha_complex(const Alpha_complex& other) = delete; + Alpha_complex& operator= (const Alpha_complex& other) = delete; + Alpha_complex (Alpha_complex&& other) = delete; + Alpha_complex& operator= (Alpha_complex&& other) = delete; + + /** \brief get_point returns the point corresponding to the vertex given as parameter. + * + * @param[in] vertex Vertex handle of the point to retrieve. + * @return The point found. + * @exception std::out_of_range In case vertex is not found (cf. std::vector::at). + */ + const Point_d& get_point(std::size_t vertex) const { + return vertex_handle_to_iterator_.at(vertex)->point(); + } + + /** \brief number_of_vertices returns the number of vertices (same as the number of points). + * + * @return The number of vertices. + */ + std::size_t number_of_vertices() const { + return vertex_handle_to_iterator_.size(); + } + + private: + template<typename InputPointRange > + void init_from_range(const InputPointRange& points) { + auto first = std::begin(points); + auto last = std::end(points); + + if (first != last) { + // point_dimension function initialization + Point_Dimension point_dimension = kernel_.point_dimension_d_object(); + + // Delaunay triangulation is point dimension. + triangulation_ = new Delaunay_triangulation(point_dimension(*first)); + + std::vector<Point_d> point_cloud(first, last); + + // Creates a vector {0, 1, ..., N-1} + std::vector<std::ptrdiff_t> indices(boost::counting_iterator<std::ptrdiff_t>(0), + boost::counting_iterator<std::ptrdiff_t>(point_cloud.size())); + + typedef boost::iterator_property_map<typename std::vector<Point_d>::iterator, + CGAL::Identity_property_map<std::ptrdiff_t>> Point_property_map; + typedef CGAL::Spatial_sort_traits_adapter_d<Kernel, Point_property_map> Search_traits_d; + + CGAL::spatial_sort(indices.begin(), indices.end(), Search_traits_d(std::begin(point_cloud))); + + typename Delaunay_triangulation::Full_cell_handle hint; + for (auto index : indices) { + typename Delaunay_triangulation::Vertex_handle pos = triangulation_->insert(point_cloud[index], hint); + // Save index value as data to retrieve it after insertion + pos->data() = index; + hint = pos->full_cell(); + } + // -------------------------------------------------------------------------------------------- + // double map to retrieve simplex tree vertex handles from CGAL vertex iterator and vice versa + // Loop on triangulation vertices list + for (CGAL_vertex_iterator vit = triangulation_->vertices_begin(); vit != triangulation_->vertices_end(); ++vit) { + if (!triangulation_->is_infinite(*vit)) { +#ifdef DEBUG_TRACES + std::cout << "Vertex insertion - " << vit->data() << " -> " << vit->point() << std::endl; +#endif // DEBUG_TRACES + vertex_handle_to_iterator_.emplace(vit->data(), vit); + } + } + // -------------------------------------------------------------------------------------------- + } + } + + public: + /** \brief Inserts all Delaunay triangulation into the simplicial complex. + * It also computes the filtration values accordingly to the \ref createcomplexalgorithm + * + * \tparam SimplicialComplexForAlpha must meet `SimplicialComplexForAlpha` concept. + * + * @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. + * + * @return true if creation succeeds, false otherwise. + * + * @pre Delaunay triangulation must be already constructed with dimension strictly greater than 0. + * @pre The simplicial complex must be empty (no vertices) + * + * Initialization can be launched once. + */ + 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()) { + // 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; + typedef std::vector<Vertex_handle> Vector_vertex; + + if (triangulation_ == nullptr) { + std::cerr << "Alpha_complex cannot create_complex from a NULL triangulation\n"; + return false; // ----- >> + } + if (triangulation_->maximal_dimension() < 1) { + std::cerr << "Alpha_complex cannot create_complex from a zero-dimension triangulation\n"; + return false; // ----- >> + } + if (complex.num_vertices() > 0) { + std::cerr << "Alpha_complex create_complex - complex is not empty\n"; + return false; // ----- >> + } + + // -------------------------------------------------------------------------------------------- + // Simplex_tree construction from loop on triangulation finite full cells list + if (triangulation_->number_of_vertices() > 0) { + for (auto cit = triangulation_->finite_full_cells_begin(); + cit != triangulation_->finite_full_cells_end(); + ++cit) { + Vector_vertex vertexVector; +#ifdef DEBUG_TRACES + std::cout << "Simplex_tree insertion "; +#endif // DEBUG_TRACES + for (auto vit = cit->vertices_begin(); vit != cit->vertices_end(); ++vit) { + if (*vit != nullptr) { +#ifdef DEBUG_TRACES + std::cout << " " << (*vit)->data(); +#endif // DEBUG_TRACES + // Vector of vertex construction for simplex_tree structure + vertexVector.push_back((*vit)->data()); + } + } +#ifdef DEBUG_TRACES + std::cout << std::endl; +#endif // DEBUG_TRACES + // Insert each simplex and its subfaces in the simplex tree - filtration is NaN + complex.insert_simplex_and_subfaces(vertexVector, std::numeric_limits<double>::quiet_NaN()); + } + } + // -------------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------------- + // Will be re-used many times + Vector_of_CGAL_points pointVector; + // ### For i : d -> 0 + for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) { + // ### Foreach Sigma of dim i + for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) { + int f_simplex_dim = complex.dimension(f_simplex); + if (decr_dim == f_simplex_dim) { + pointVector.clear(); +#ifdef DEBUG_TRACES + std::cout << "Sigma of dim " << decr_dim << " is"; +#endif // DEBUG_TRACES + for (auto vertex : complex.simplex_vertex_range(f_simplex)) { + pointVector.push_back(get_point(vertex)); +#ifdef DEBUG_TRACES + std::cout << " " << vertex; +#endif // DEBUG_TRACES + } +#ifdef DEBUG_TRACES + std::cout << std::endl; +#endif // DEBUG_TRACES + // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma) + if (std::isnan(complex.filtration(f_simplex))) { + Filtration_value alpha_complex_filtration = 0.0; + // No need to compute squared_radius on a single point - alpha is 0.0 + 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())); + } + complex.assign_filtration(f_simplex, alpha_complex_filtration); +#ifdef DEBUG_TRACES + std::cout << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl; +#endif // DEBUG_TRACES + } + // No need to propagate further, unweighted points all have value 0 + if (decr_dim > 1) + propagate_alpha_filtration(complex, f_simplex); + } + } + } + // -------------------------------------------------------------------------------------------- + + // -------------------------------------------------------------------------------------------- + // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension + complex.make_filtration_non_decreasing(); + // Remove all simplices that have a filtration value greater than max_alpha_square + complex.prune_above_filtration(max_alpha_square); + // -------------------------------------------------------------------------------------------- + return true; + } + + private: + template <typename SimplicialComplexForAlpha, typename Simplex_handle> + void propagate_alpha_filtration(SimplicialComplexForAlpha& complex, Simplex_handle f_simplex) { + // From SimplicialComplexForAlpha type required to assign filtration values. + typedef typename SimplicialComplexForAlpha::Filtration_value Filtration_value; +#ifdef DEBUG_TRACES + typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle; +#endif // DEBUG_TRACES + + // ### Foreach Tau face of Sigma + for (auto f_boundary : complex.boundary_simplex_range(f_simplex)) { +#ifdef DEBUG_TRACES + std::cout << " | --------------------------------------------------\n"; + std::cout << " | Tau "; + for (auto vertex : complex.simplex_vertex_range(f_boundary)) { + std::cout << vertex << " "; + } + std::cout << "is a face of Sigma\n"; + std::cout << " | isnan(complex.filtration(Tau)=" << std::isnan(complex.filtration(f_boundary)) << std::endl; +#endif // DEBUG_TRACES + // ### If filt(Tau) is not NaN + if (!std::isnan(complex.filtration(f_boundary))) { + // ### filt(Tau) = fmin(filt(Tau), filt(Sigma)) + Filtration_value alpha_complex_filtration = fmin(complex.filtration(f_boundary), + complex.filtration(f_simplex)); + complex.assign_filtration(f_boundary, alpha_complex_filtration); +#ifdef DEBUG_TRACES + std::cout << " | filt(Tau) = fmin(filt(Tau), filt(Sigma)) = " << complex.filtration(f_boundary) << std::endl; +#endif // DEBUG_TRACES + // ### Else + } else { + // insert the Tau points in a vector for is_gabriel function + Vector_of_CGAL_points pointVector; +#ifdef DEBUG_TRACES + Vertex_handle vertexForGabriel = Vertex_handle(); +#endif // DEBUG_TRACES + for (auto vertex : complex.simplex_vertex_range(f_boundary)) { + pointVector.push_back(get_point(vertex)); + } + // Retrieve the Sigma point that is not part of Tau - parameter for is_gabriel function + Point_d point_for_gabriel; + for (auto vertex : complex.simplex_vertex_range(f_simplex)) { + point_for_gabriel = get_point(vertex); + if (std::find(pointVector.begin(), pointVector.end(), point_for_gabriel) == pointVector.end()) { +#ifdef DEBUG_TRACES + // vertex is not found in Tau + vertexForGabriel = vertex; +#endif // DEBUG_TRACES + // No need to continue loop + break; + } + } + // is_gabriel function initialization + Is_Gabriel is_gabriel = kernel_.side_of_bounded_sphere_d_object(); + bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), point_for_gabriel) + != CGAL::ON_BOUNDED_SIDE; +#ifdef DEBUG_TRACES + std::cout << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << vertexForGabriel << std::endl; +#endif // DEBUG_TRACES + // ### If Tau is not Gabriel of Sigma + if (false == is_gab) { + // ### filt(Tau) = filt(Sigma) + Filtration_value alpha_complex_filtration = complex.filtration(f_simplex); + complex.assign_filtration(f_boundary, alpha_complex_filtration); +#ifdef DEBUG_TRACES + std::cout << " | filt(Tau) = filt(Sigma) = " << complex.filtration(f_boundary) << std::endl; +#endif // DEBUG_TRACES + } + } + } + } +}; + +} // namespace alpha_complex + +namespace alphacomplex = alpha_complex; + +} // namespace Gudhi + +#endif // ALPHA_COMPLEX_H_ diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h new file mode 100644 index 00000000..13ebb9c1 --- /dev/null +++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h @@ -0,0 +1,579 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2018 Inria + * + * Modification(s): + * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL and Eigen3 + * - YYYY/MM Author: Description of the modification + */ + +#ifndef ALPHA_COMPLEX_3D_H_ +#define ALPHA_COMPLEX_3D_H_ + +#include <boost/version.hpp> +#include <boost/variant.hpp> + +#include <gudhi/Debug_utils.h> +#include <gudhi/Alpha_complex_options.h> + +#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> +#include <CGAL/Exact_predicates_exact_constructions_kernel.h> +#include <CGAL/Delaunay_triangulation_3.h> +#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h> +#include <CGAL/Periodic_3_Delaunay_triangulation_3.h> +#include <CGAL/Periodic_3_regular_triangulation_traits_3.h> +#include <CGAL/Periodic_3_regular_triangulation_3.h> +#include <CGAL/Regular_triangulation_3.h> +#include <CGAL/Alpha_shape_3.h> +#include <CGAL/Alpha_shape_cell_base_3.h> +#include <CGAL/Alpha_shape_vertex_base_3.h> + +#include <CGAL/Object.h> +#include <CGAL/tuple.h> +#include <CGAL/iterator.h> +#include <CGAL/version.h> // for CGAL_VERSION_NR + +#include <Eigen/src/Core/util/Macros.h> // for EIGEN_VERSION_AT_LEAST + +#include <boost/container/static_vector.hpp> + +#include <iostream> +#include <vector> +#include <unordered_map> +#include <stdexcept> +#include <cstddef> +#include <memory> // for std::unique_ptr +#include <type_traits> // for std::conditional and std::enable_if +#include <limits> // for numeric_limits<> + +// 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 +#endif + +#if !EIGEN_VERSION_AT_LEAST(3,1,0) +# error Alpha_complex_3d is only available for Eigen3 >= 3.1.0 installed with CGAL +#endif + +namespace Gudhi { + +namespace alpha_complex { + +#ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL +thread_local +#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL + double RELATIVE_PRECISION_OF_TO_DOUBLE = 0.00001; + +// Value_from_iterator returns the filtration value from an iterator on alpha shapes values +// +// FAST SAFE EXACT +// CGAL::to_double(*iterator) CGAL::to_double(*iterator) CGAL::to_double(iterator->exact()) + +template <complexity Complexity> +struct Value_from_iterator { + template <typename Iterator> + static double perform(Iterator it) { + // Default behaviour + return CGAL::to_double(*it); + } +}; + +template <> +struct Value_from_iterator<complexity::EXACT> { + template <typename Iterator> + static double perform(Iterator it) { + return CGAL::to_double(it->exact()); + } +}; + +/** + * \class Alpha_complex_3d + * \brief Alpha complex data structure for 3d specific case. + * + * \ingroup alpha_complex + * + * \details + * The data structure is constructing a <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html">CGAL 3D Alpha + * Shapes</a> from a range of points (can be read from an OFF file, cf. Points_off_reader). + * Duplicate points are inserted once in the Alpha_complex. This is the reason why the vertices may be not contiguous. + * + * \tparam Complexity shall be `Gudhi::alpha_complex::complexity` type. Default value is + * `Gudhi::alpha_complex::complexity::SAFE`. + * + * \tparam Weighted Boolean used to set/unset the weighted version of Alpha_complex_3d. Default value is false. + * + * \tparam Periodic Boolean used to set/unset the periodic version of Alpha_complex_3d. Default value is false. + * + * For the weighted version, weights values are explained on CGAL + * <a href="https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0">Alpha shapes 3d</a> and + * <a href="https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation">Regular + * triangulation</a> documentation. + * + * For the periodic version, refer to the + * <a href="https://doc.cgal.org/latest/Periodic_3_triangulation_3/index.html">CGAL’s 3D Periodic Triangulations User + * Manual </a> for more details. + * The periodicity is defined by an iso-oriented cuboid with diagonal opposite vertices (x_min, y_min, z_min) and + * (x_max, y_max, z_max). + * + * Please refer to \ref alpha_complex for examples. + * + * \remark When Alpha_complex_3d is constructed with an infinite value of alpha (default value), the complex is a + * 3d Delaunay complex. + * + */ +template <complexity Complexity = complexity::SAFE, bool Weighted = false, bool Periodic = false> +class Alpha_complex_3d { + // Epick = Exact_predicates_inexact_constructions_kernel + // Epeck = Exact_predicates_exact_constructions_kernel + // Exact_alpha_comparison_tag = exact version of CGAL Alpha_shape_3 and of its objects (Alpha_shape_vertex_base_3 and + // Alpha_shape_cell_base_3). Not available if weighted or periodic. + // Can be CGAL::Tag_false or CGAL::Tag_true. Default is False. + // cf. https://doc.cgal.org/latest/Alpha_shapes_3/classCGAL_1_1Alpha__shape__3.html + // + // We could use Epick + CGAL::Tag_true for not weighted nor periodic, but during benchmark, we found a bug + // https://github.com/CGAL/cgal/issues/3460 + // This is the reason we only use Epick + CGAL::Tag_false, or Epeck + // + // FAST SAFE EXACT + // Epick + CGAL::Tag_false Epeck Epeck + using Predicates = typename std::conditional<(Complexity == complexity::FAST), + CGAL::Exact_predicates_inexact_constructions_kernel, + CGAL::Exact_predicates_exact_constructions_kernel>::type; + + // The other way to do a conditional type. Here there are 3 possibilities + template <typename Predicates, bool Weighted_version, bool Periodic_version> + struct Kernel_3 {}; + + template <typename Predicates, bool Is_periodic> + struct Kernel_3<Predicates, Is_periodic, false> { + using Kernel = Predicates; + }; + + template <typename Predicates> + struct Kernel_3<Predicates, false, true> { + using Kernel = CGAL::Periodic_3_Delaunay_triangulation_traits_3<Predicates>; + }; + template <typename Predicates> + struct Kernel_3<Predicates, true, true> { + using Kernel = CGAL::Periodic_3_regular_triangulation_traits_3<Predicates>; + }; + + using Kernel = typename Kernel_3<Predicates, Weighted, Periodic>::Kernel; + + using TdsVb = typename std::conditional<Periodic, CGAL::Periodic_3_triangulation_ds_vertex_base_3<>, + CGAL::Triangulation_ds_vertex_base_3<>>::type; + + using Tvb = typename std::conditional<Weighted, CGAL::Regular_triangulation_vertex_base_3<Kernel, TdsVb>, + CGAL::Triangulation_vertex_base_3<Kernel, TdsVb>>::type; + + using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel, Tvb>; + + using TdsCb = typename std::conditional<Periodic, CGAL::Periodic_3_triangulation_ds_cell_base_3<>, + CGAL::Triangulation_ds_cell_base_3<>>::type; + + using Tcb = typename std::conditional<Weighted, CGAL::Regular_triangulation_cell_base_3<Kernel, TdsCb>, + CGAL::Triangulation_cell_base_3<Kernel, TdsCb>>::type; + + using Cb = CGAL::Alpha_shape_cell_base_3<Kernel, Tcb>; + using Tds = CGAL::Triangulation_data_structure_3<Vb, Cb>; + + // The other way to do a conditional type. Here there 4 possibilities, cannot use std::conditional + template <typename Kernel, typename Tds, bool Weighted_version, bool Periodic_version> + struct Triangulation_3 {}; + + template <typename Kernel, typename Tds> + struct Triangulation_3<Kernel, Tds, false, false> { + using Dt = CGAL::Delaunay_triangulation_3<Kernel, Tds>; + using Weighted_point_3 = void; + }; + template <typename Kernel, typename Tds> + struct Triangulation_3<Kernel, Tds, true, false> { + using Dt = CGAL::Regular_triangulation_3<Kernel, Tds>; + using Weighted_point_3 = typename Dt::Weighted_point; + }; + template <typename Kernel, typename Tds> + struct Triangulation_3<Kernel, Tds, false, true> { + using Dt = CGAL::Periodic_3_Delaunay_triangulation_3<Kernel, Tds>; + using Weighted_point_3 = void; + }; + template <typename Kernel, typename Tds> + struct Triangulation_3<Kernel, Tds, true, true> { + using Dt = CGAL::Periodic_3_regular_triangulation_3<Kernel, Tds>; + using Weighted_point_3 = typename Dt::Weighted_point; + }; + + /** \brief Is either Delaunay_triangulation_3 (Weighted = false and Periodic = false), + * Regular_triangulation_3 (Weighted = true and Periodic = false), + * Periodic_3_Delaunay_triangulation_3 (Weighted = false and Periodic = true) + * or Periodic_3_regular_triangulation_3 (Weighted = true and Periodic = true). + * + * This type is required by `Gudhi::alpha_complex::Alpha_complex_3d::Alpha_shape_3`. + * */ + using Dt = typename Triangulation_3<Kernel, Tds, Weighted, Periodic>::Dt; + + public: + /** \brief The <a href="https://doc.cgal.org/latest/Alpha_shapes_3/classCGAL_1_1Alpha__shape__3.html">CGAL 3D Alpha + * Shapes</a> type. + * + * The `Gudhi::alpha_complex::Alpha_complex_3d` is a wrapper on top of this class to ease the standard, weighted + * and/or periodic build of the Alpha complex 3d.*/ + using Alpha_shape_3 = CGAL::Alpha_shape_3<Dt>; + + /** \brief The alpha values type. + * Must be compatible with double. */ + using FT = typename Alpha_shape_3::FT; + + /** \brief Gives public access to the Point_3 type. Here is a Point_3 constructor example: +\code{.cpp} +using Alpha_complex_3d = Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, false, false>; + +// x0 = 1., y0 = -1.1, z0 = -1.. +Alpha_complex_3d::Point_3 p0(1., -1.1, -1.); +\endcode + * */ + using Point_3 = typename Kernel::Point_3; + + /** \brief Gives public access to the Weighted_point_3 type. A Weighted point can be constructed as follows: +\code{.cpp} +using Weighted_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, false>; + +// x0 = 1., y0 = -1.1, z0 = -1., weight = 4. +Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Point_3(1., -1.1, -1.), 4.); +\endcode + * + * Note: This type is defined to void if Alpha complex is not weighted. + * + * */ + using Weighted_point_3 = typename Triangulation_3<Kernel, Tds, Weighted, Periodic>::Weighted_point_3; + + private: + using Dispatch = + CGAL::Dispatch_output_iterator<CGAL::cpp11::tuple<CGAL::Object, FT>, + CGAL::cpp11::tuple<std::back_insert_iterator<std::vector<CGAL::Object>>, + std::back_insert_iterator<std::vector<FT>>>>; + + using Cell_handle = typename Alpha_shape_3::Cell_handle; + using Facet = typename Alpha_shape_3::Facet; + using Edge = typename Alpha_shape_3::Edge; + using Alpha_vertex_handle = typename Alpha_shape_3::Vertex_handle; + using Vertex_list = boost::container::static_vector<Alpha_vertex_handle, 4>; + + public: + /** \brief Alpha_complex constructor from a list of points. + * + * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3` or + * `Alpha_complex_3d::Weighted_point_3`. + * + * @pre Available if Alpha_complex_3d is not Periodic. + * + * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a + * `Alpha_complex_3d::Point_3` or a `Alpha_complex_3d::Weighted_point_3`. + */ + template <typename InputPointRange> + Alpha_complex_3d(const InputPointRange& points) { + static_assert(!Periodic, "This constructor is not available for periodic versions of Alpha_complex_3d"); + + alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>( + new Alpha_shape_3(std::begin(points), std::end(points), 0, Alpha_shape_3::GENERAL)); + } + + /** \brief Alpha_complex constructor from a list of points and associated weights. + * + * @exception std::invalid_argument In debug mode, if points and weights do not have the same size. + * + * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3`. + * @param[in] weights Range of weights on points. Weights shall be in double. + * + * @pre Available if Alpha_complex_3d is Weighted and not Periodic. + * + * The type InputPointRange must be a range for which std::begin and + * std::end return input iterators on a `Alpha_complex_3d::Point_3`. + * The type WeightRange must be a range for which std::begin and + * std::end return an input iterator on a double. + */ + template <typename InputPointRange, typename WeightRange> + Alpha_complex_3d(const InputPointRange& points, WeightRange weights) { + static_assert(Weighted, "This constructor is not available for non-weighted versions of Alpha_complex_3d"); + static_assert(!Periodic, "This constructor is not available for periodic versions of Alpha_complex_3d"); + GUDHI_CHECK((weights.size() == points.size()), + std::invalid_argument("Points number in range different from weights range number")); + + std::vector<Weighted_point_3> weighted_points_3; + + std::size_t index = 0; + weighted_points_3.reserve(points.size()); + while ((index < weights.size()) && (index < points.size())) { + weighted_points_3.push_back(Weighted_point_3(points[index], weights[index])); + index++; + } + + alpha_shape_3_ptr_ = std::unique_ptr<Alpha_shape_3>( + new Alpha_shape_3(std::begin(weighted_points_3), std::end(weighted_points_3), 0, Alpha_shape_3::GENERAL)); + } + + /** \brief Alpha_complex constructor from a list of points and an iso-cuboid coordinates. + * + * @exception std::invalid_argument In debug mode, if the size of the cuboid in every directions is not the same. + * + * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3` or + * `Alpha_complex_3d::Weighted_point_3`. + * @param[in] x_min Iso-oriented cuboid x_min. + * @param[in] y_min Iso-oriented cuboid y_min. + * @param[in] z_min Iso-oriented cuboid z_min. + * @param[in] x_max Iso-oriented cuboid x_max. + * @param[in] y_max Iso-oriented cuboid y_max. + * @param[in] z_max Iso-oriented cuboid z_max. + * + * @pre Available if Alpha_complex_3d is Periodic. + * + * The type InputPointRange must be a range for which std::begin and std::end return input iterators on a + * `Alpha_complex_3d::Point_3` or a `Alpha_complex_3d::Weighted_point_3`. + * + * @note In weighted version, please check weights are greater than zero, and lower than 1/64*cuboid length + * squared. + */ + template <typename InputPointRange> + Alpha_complex_3d(const InputPointRange& points, FT x_min, FT y_min, + FT z_min, FT x_max, FT y_max, FT z_max) { + static_assert(Periodic, "This constructor is not available for non-periodic versions of Alpha_complex_3d"); + // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it. + GUDHI_CHECK( + (x_max - x_min == y_max - y_min) && (x_max - x_min == z_max - z_min) && (z_max - z_min == y_max - y_min), + std::invalid_argument("The size of the cuboid in every directions is not the same.")); + + // Define the periodic cube + Dt pdt(typename Kernel::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(std::begin(points), std::end(points), true); + // As pdt won't be modified anymore switch to 1-sheeted cover if possible + if (!pdt.is_triangulation_in_1_sheet()) { + throw std::invalid_argument("Unable to construct a triangulation within a single periodic domain."); + } + pdt.convert_to_1_sheeted_covering(); + + // 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_ptr_ = std::unique_ptr<Alpha_shape_3>(new Alpha_shape_3(pdt, 0, Alpha_shape_3::GENERAL)); + } + + /** \brief Alpha_complex constructor from a list of points, associated weights and an iso-cuboid coordinates. + * + * @exception std::invalid_argument In debug mode, if points and weights do not have the same size. + * @exception std::invalid_argument In debug mode, if the size of the cuboid in every directions is not the same. + * @exception std::invalid_argument In debug mode, if a weight is negative, zero, or greater than 1/64*cuboid length + * squared. + * + * @param[in] points Range of points to triangulate. Points must be in `Alpha_complex_3d::Point_3`. + * @param[in] weights Range of weights on points. Weights shall be in double. + * @param[in] x_min Iso-oriented cuboid x_min. + * @param[in] y_min Iso-oriented cuboid y_min. + * @param[in] z_min Iso-oriented cuboid z_min. + * @param[in] x_max Iso-oriented cuboid x_max. + * @param[in] y_max Iso-oriented cuboid y_max. + * @param[in] z_max Iso-oriented cuboid z_max. + * + * @pre Available if Alpha_complex_3d is Weighted and Periodic. + * + * The type InputPointRange must be a range for which std::begin and + * std::end return input iterators on a `Alpha_complex_3d::Point_3`. + * The type WeightRange must be a range for which std::begin and + * std::end return an input iterator on a double. + * The type of x_min, y_min, z_min, x_max, y_max and z_max must be a double. + */ + template <typename InputPointRange, typename WeightRange> + Alpha_complex_3d(const InputPointRange& points, WeightRange weights, FT x_min, FT y_min, + FT z_min, FT x_max, FT y_max, FT z_max) { + static_assert(Weighted, "This constructor is not available for non-weighted versions of Alpha_complex_3d"); + static_assert(Periodic, "This constructor is not available for non-periodic versions of Alpha_complex_3d"); + GUDHI_CHECK((weights.size() == points.size()), + std::invalid_argument("Points number in range different from weights range number")); + // Checking if the cuboid is the same in x,y and z direction. If not, CGAL will not process it. + GUDHI_CHECK( + (x_max - x_min == y_max - y_min) && (x_max - x_min == z_max - z_min) && (z_max - z_min == y_max - y_min), + std::invalid_argument("The size of the cuboid in every directions is not the same.")); + + std::vector<Weighted_point_3> weighted_points_3; + + std::size_t index = 0; + weighted_points_3.reserve(points.size()); + +#ifdef GUDHI_DEBUG + // Defined in GUDHI_DEBUG to avoid unused variable warning for GUDHI_CHECK + FT maximal_possible_weight = 0.015625 * (x_max - x_min) * (x_max - x_min); +#endif + + while ((index < weights.size()) && (index < points.size())) { + GUDHI_CHECK((weights[index] < maximal_possible_weight) && (weights[index] >= 0), + std::invalid_argument("Invalid weight at index " + std::to_string(index + 1) + + ". Must be positive and less than maximal possible weight = 1/64*cuboid length " + "squared, which is not an acceptable input.")); + weighted_points_3.push_back(Weighted_point_3(points[index], weights[index])); + index++; + } + + // Define the periodic cube + Dt pdt(typename Kernel::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(std::begin(weighted_points_3), std::end(weighted_points_3), true); + // As pdt won't be modified anymore switch to 1-sheeted cover if possible + if (!pdt.is_triangulation_in_1_sheet()) { + throw std::invalid_argument("Unable to construct a triangulation within a single periodic domain."); + } + pdt.convert_to_1_sheeted_covering(); + + // 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_ptr_ = std::unique_ptr<Alpha_shape_3>(new Alpha_shape_3(pdt, 0, Alpha_shape_3::GENERAL)); + } + + /** \brief Inserts all Delaunay triangulation into the simplicial complex. + * It also computes the filtration values accordingly to the \ref createcomplexalgorithm + * + * \tparam SimplicialComplexForAlpha3d must meet `SimplicialComplexForAlpha3d` concept. + * + * @param[in] complex SimplicialComplexForAlpha3d 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. + * + * @return true if creation succeeds, false otherwise. + * + * @pre The simplicial complex must be empty (no vertices). + * + */ + template <typename SimplicialComplexForAlpha3d, + typename Filtration_value = typename SimplicialComplexForAlpha3d::Filtration_value> + bool create_complex(SimplicialComplexForAlpha3d& complex, + Filtration_value max_alpha_square = std::numeric_limits<Filtration_value>::infinity()) { + if (complex.num_vertices() > 0) { + std::cerr << "Alpha_complex_3d create_complex - complex is not empty\n"; + return false; // ----- >> + } + + // using Filtration_value = typename SimplicialComplexForAlpha3d::Filtration_value; + using Complex_vertex_handle = typename SimplicialComplexForAlpha3d::Vertex_handle; + using Alpha_shape_simplex_tree_map = std::unordered_map<Alpha_vertex_handle, Complex_vertex_handle>; + using Simplex_tree_vector_vertex = std::vector<Complex_vertex_handle>; + +#ifdef DEBUG_TRACES + std::size_t count_vertices = 0; + std::size_t count_edges = 0; + std::size_t count_facets = 0; + std::size_t count_cells = 0; +#endif // DEBUG_TRACES + std::vector<CGAL::Object> objects; + std::vector<FT> alpha_values; + + Dispatch dispatcher = CGAL::dispatch_output<CGAL::Object, FT>(std::back_inserter(objects), + std::back_inserter(alpha_values)); + + alpha_shape_3_ptr_->filtration_with_alpha_values(dispatcher); +#ifdef DEBUG_TRACES + std::cout << "filtration_with_alpha_values returns : " << objects.size() << " objects" << std::endl; +#endif // DEBUG_TRACES + + Alpha_shape_simplex_tree_map map_cgal_simplex_tree; + using Alpha_value_iterator = typename std::vector<FT>::const_iterator; + Alpha_value_iterator alpha_value_iterator = alpha_values.begin(); + for (auto object_iterator : objects) { + Vertex_list vertex_list; + + // Retrieve Alpha shape vertex list from object + if (const Cell_handle* cell = CGAL::object_cast<Cell_handle>(&object_iterator)) { + for (auto i = 0; i < 4; i++) { +#ifdef DEBUG_TRACES + std::cout << "from cell[" << i << "]=" << (*cell)->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + vertex_list.push_back((*cell)->vertex(i)); + } +#ifdef DEBUG_TRACES + count_cells++; +#endif // DEBUG_TRACES + } else if (const Facet* facet = CGAL::object_cast<Facet>(&object_iterator)) { + for (auto i = 0; i < 4; i++) { + if ((*facet).second != i) { +#ifdef DEBUG_TRACES + std::cout << "from facet=[" << i << "]" << (*facet).first->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + vertex_list.push_back((*facet).first->vertex(i)); + } + } +#ifdef DEBUG_TRACES + count_facets++; +#endif // DEBUG_TRACES + } else if (const Edge* edge = CGAL::object_cast<Edge>(&object_iterator)) { + for (auto i : {(*edge).second, (*edge).third}) { +#ifdef DEBUG_TRACES + std::cout << "from edge[" << i << "]=" << (*edge).first->vertex(i)->point() << std::endl; +#endif // DEBUG_TRACES + vertex_list.push_back((*edge).first->vertex(i)); + } +#ifdef DEBUG_TRACES + count_edges++; +#endif // DEBUG_TRACES + } else if (const Alpha_vertex_handle* vertex = CGAL::object_cast<Alpha_vertex_handle>(&object_iterator)) { +#ifdef DEBUG_TRACES + count_vertices++; + std::cout << "from vertex=" << (*vertex)->point() << std::endl; +#endif // DEBUG_TRACES + vertex_list.push_back((*vertex)); + } + // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex + Simplex_tree_vector_vertex the_simplex; + for (auto the_alpha_shape_vertex : vertex_list) { + auto 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 + Complex_vertex_handle 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.push_back(vertex); + map_cgal_simplex_tree.emplace(the_alpha_shape_vertex, vertex); + } else { + // alpha shape found + Complex_vertex_handle 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.push_back(vertex); + } + } + // Construction of the simplex_tree + Filtration_value filtr = Value_from_iterator<Complexity>::perform(alpha_value_iterator); + +#ifdef DEBUG_TRACES + std::cout << "filtration = " << filtr << std::endl; +#endif // DEBUG_TRACES + complex.insert_simplex(the_simplex, static_cast<Filtration_value>(filtr)); + GUDHI_CHECK(alpha_value_iterator != alpha_values.end(), "CGAL provided more simplices than values"); + ++alpha_value_iterator; + } + +#ifdef DEBUG_TRACES + std::cout << "vertices \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; +#endif // DEBUG_TRACES + // -------------------------------------------------------------------------------------------- + // As Alpha value is an approximation, we have to make filtration non decreasing while increasing the dimension + complex.make_filtration_non_decreasing(); + // Remove all simplices that have a filtration value greater than max_alpha_square + complex.prune_above_filtration(max_alpha_square); + // -------------------------------------------------------------------------------------------- + return true; + } + + private: + // use of a unique_ptr on cgal Alpha_shape_3, as copy and default constructor is not available - no need to be freed + std::unique_ptr<Alpha_shape_3> alpha_shape_3_ptr_; +}; + +} // namespace alpha_complex + +} // namespace Gudhi + +#endif // ALPHA_COMPLEX_3D_H_ diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_options.h b/src/Alpha_complex/include/gudhi/Alpha_complex_options.h new file mode 100644 index 00000000..85c83672 --- /dev/null +++ b/src/Alpha_complex/include/gudhi/Alpha_complex_options.h @@ -0,0 +1,33 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2018 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef ALPHA_COMPLEX_OPTIONS_H_ +#define ALPHA_COMPLEX_OPTIONS_H_ + +namespace Gudhi { + +namespace alpha_complex { + +/** + * \brief Alpha complex complexity template parameter possible values. + * + * \ingroup alpha_complex + */ +enum class complexity : char { + FAST = 'f', ///< Fast version. + SAFE = 's', ///< Safe version. + EXACT = 'e', ///< Exact version. +}; + +} // namespace alpha_complex + +} // namespace Gudhi + +#endif // ALPHA_COMPLEX_OPTIONS_H_ diff --git a/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp new file mode 100644 index 00000000..1102838a --- /dev/null +++ b/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp @@ -0,0 +1,152 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "alpha_complex_3d" +#include <boost/test/unit_test.hpp> + +#include <cmath> // float comparison +#include <limits> +#include <string> +#include <vector> +#include <random> +#include <cstddef> // for std::size_t + +#include <gudhi/Alpha_complex_3d.h> +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Unitary_tests_utils.h> +// to construct Alpha_complex from a OFF file of points +#include <gudhi/Points_3D_off_io.h> + +#include <CGAL/Random.h> +#include <CGAL/point_generators_3.h> + +using Fast_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::FAST, false, false>; +using Safe_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, false, false>; +using Exact_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, false, false>; + +template <typename Point> +std::vector<Point> get_points() { + std::vector<Point> points; + points.push_back(Point(0.0, 0.0, 0.0)); + points.push_back(Point(0.0, 0.0, 0.2)); + points.push_back(Point(0.2, 0.0, 0.2)); + points.push_back(Point(0.6, 0.6, 0.0)); + points.push_back(Point(0.8, 0.8, 0.2)); + points.push_back(Point(0.2, 0.8, 0.6)); + + return points; +} + + +BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) { + // ----------------- + // Fast version + // ----------------- + std::cout << "Fast alpha complex 3d" << std::endl; + + Fast_alpha_complex_3d alpha_complex(get_points<Fast_alpha_complex_3d::Point_3>()); + + Gudhi::Simplex_tree<> stree; + alpha_complex.create_complex(stree); + + // ----------------- + // Exact version + // ----------------- + std::cout << "Exact alpha complex 3d" << std::endl; + + Exact_alpha_complex_3d exact_alpha_complex(get_points<Exact_alpha_complex_3d::Point_3>()); + + Gudhi::Simplex_tree<> exact_stree; + exact_alpha_complex.create_complex(exact_stree); + + // --------------------- + // Compare both versions + // --------------------- + std::cout << "Exact Alpha complex 3d is of dimension " << exact_stree.dimension() << " - Fast is " + << stree.dimension() << std::endl; + BOOST_CHECK(exact_stree.dimension() == stree.dimension()); + std::cout << "Exact Alpha complex 3d num_simplices " << exact_stree.num_simplices() << " - Fast is " + << stree.num_simplices() << std::endl; + BOOST_CHECK(exact_stree.num_simplices() == stree.num_simplices()); + std::cout << "Exact Alpha complex 3d num_vertices " << exact_stree.num_vertices() << " - Fast is " + << stree.num_vertices() << std::endl; + BOOST_CHECK(exact_stree.num_vertices() == stree.num_vertices()); + + auto sh = stree.filtration_simplex_range().begin(); + while (sh != stree.filtration_simplex_range().end()) { + std::vector<int> simplex; + std::vector<int> exact_simplex; + std::cout << "Fast ( "; + for (auto vertex : stree.simplex_vertex_range(*sh)) { + simplex.push_back(vertex); + std::cout << vertex << " "; + } + std::cout << ") -> [" << stree.filtration(*sh) << "] "; + + // Find it in the exact structure + auto sh_exact = exact_stree.find(simplex); + BOOST_CHECK(sh_exact != exact_stree.null_simplex()); + + std::cout << " versus [" << exact_stree.filtration(sh_exact) << "] " << std::endl; + // Exact and non-exact version is not exactly the same due to float comparison + GUDHI_TEST_FLOAT_EQUALITY_CHECK(exact_stree.filtration(sh_exact), stree.filtration(*sh)); + + ++sh; + } + // ----------------- + // Safe version + // ----------------- + std::cout << "Safe alpha complex 3d" << std::endl; + + Safe_alpha_complex_3d safe_alpha_complex(get_points<Safe_alpha_complex_3d::Point_3>()); + + Gudhi::Simplex_tree<> safe_stree; + safe_alpha_complex.create_complex(safe_stree); + + // --------------------- + // Compare both versions + // --------------------- + std::cout << "Safe Alpha complex 3d is of dimension " << safe_stree.dimension() << " - Fast is " + << stree.dimension() << std::endl; + BOOST_CHECK(safe_stree.dimension() == stree.dimension()); + std::cout << "Safe Alpha complex 3d num_simplices " << safe_stree.num_simplices() << " - Fast is " + << stree.num_simplices() << std::endl; + BOOST_CHECK(safe_stree.num_simplices() == stree.num_simplices()); + std::cout << "Safe Alpha complex 3d num_vertices " << safe_stree.num_vertices() << " - Fast is " + << stree.num_vertices() << std::endl; + BOOST_CHECK(safe_stree.num_vertices() == stree.num_vertices()); + + auto safe_sh = stree.filtration_simplex_range().begin(); + while (safe_sh != stree.filtration_simplex_range().end()) { + std::vector<int> simplex; + std::vector<int> exact_simplex; + std::cout << "Fast ( "; + for (auto vertex : stree.simplex_vertex_range(*safe_sh)) { + simplex.push_back(vertex); + std::cout << vertex << " "; + } + std::cout << ") -> [" << stree.filtration(*safe_sh) << "] "; + + // Find it in the exact structure + auto sh_exact = safe_stree.find(simplex); + BOOST_CHECK(sh_exact != safe_stree.null_simplex()); + + std::cout << " versus [" << safe_stree.filtration(sh_exact) << "] " << std::endl; + // Exact and non-exact version is not exactly the same due to float comparison + GUDHI_TEST_FLOAT_EQUALITY_CHECK(safe_stree.filtration(sh_exact), stree.filtration(*safe_sh), 1e-15); + + ++safe_sh; + } +} diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp new file mode 100644 index 00000000..01e4cee3 --- /dev/null +++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp @@ -0,0 +1,271 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "alpha_complex" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <CGAL/Delaunay_triangulation.h> +#include <CGAL/Epick_d.h> + +#include <cmath> // float comparison +#include <limits> +#include <string> +#include <vector> + +#include <gudhi/Alpha_complex.h> +// to construct a simplex_tree from Delaunay_triangulation +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Simplex_tree.h> +#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; +// Use static dimension_tag for the user not to be able to set dimension +typedef CGAL::Epick_d< CGAL::Dimension_tag<3> > 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; + +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of_kernel_variants) { + // ---------------------------------------------------------------------------- + // + // Init of an alpha-complex from a OFF file + // + // ---------------------------------------------------------------------------- + std::string off_file_name("alphacomplexdoc.off"); + double max_alpha_square_value = 60.0; + std::cout << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << + max_alpha_square_value << "==========" << std::endl; + + Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_file(off_file_name); + + std::cout << "alpha_complex_from_points.number_of_vertices()=" << alpha_complex_from_file.number_of_vertices() + << std::endl; + BOOST_CHECK(alpha_complex_from_file.number_of_vertices() == 7); + + Gudhi::Simplex_tree<> simplex_tree_60; + BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_60, max_alpha_square_value)); + + std::cout << "simplex_tree_60.dimension()=" << simplex_tree_60.dimension() << std::endl; + BOOST_CHECK(simplex_tree_60.dimension() == 2); + + std::cout << "alpha_complex_from_points.number_of_vertices()=" << alpha_complex_from_file.number_of_vertices() + << std::endl; + BOOST_CHECK(alpha_complex_from_file.number_of_vertices() == 7); + + std::cout << "simplex_tree_60.num_vertices()=" << simplex_tree_60.num_vertices() << std::endl; + BOOST_CHECK(simplex_tree_60.num_vertices() == 7); + + std::cout << "simplex_tree_60.num_simplices()=" << simplex_tree_60.num_simplices() << std::endl; + BOOST_CHECK(simplex_tree_60.num_simplices() == 25); + + max_alpha_square_value = 59.0; + std::cout << "========== OFF FILE NAME = " << off_file_name << " - alpha²=" << + max_alpha_square_value << "==========" << std::endl; + + Gudhi::Simplex_tree<> simplex_tree_59; + BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_59, max_alpha_square_value)); + + std::cout << "simplex_tree_59.dimension()=" << simplex_tree_59.dimension() << std::endl; + BOOST_CHECK(simplex_tree_59.dimension() == 2); + + std::cout << "simplex_tree_59.num_vertices()=" << simplex_tree_59.num_vertices() << std::endl; + BOOST_CHECK(simplex_tree_59.num_vertices() == 7); + + std::cout << "simplex_tree_59.num_simplices()=" << simplex_tree_59.num_simplices() << std::endl; + BOOST_CHECK(simplex_tree_59.num_simplices() == 23); +} + +// 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 Kernel_4::Point_d Point_4; +typedef std::vector<Point_4> Vector_4_Points; + +bool is_point_in_list(Vector_4_Points points_list, Point_4 point) { + for (auto& point_in_list : points_list) { + if (point_in_list == point) { + return true; // point found + } + } + return false; // point not found +} + +BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { + // ---------------------------------------------------------------------------- + // Init of a list of points + // ---------------------------------------------------------------------------- + Vector_4_Points points; + std::vector<double> coords = { 0.0, 0.0, 0.0, 1.0 }; + points.push_back(Point_4(coords.begin(), coords.end())); + coords = { 0.0, 0.0, 1.0, 0.0 }; + points.push_back(Point_4(coords.begin(), coords.end())); + coords = { 0.0, 1.0, 0.0, 0.0 }; + points.push_back(Point_4(coords.begin(), coords.end())); + coords = { 1.0, 0.0, 0.0, 0.0 }; + points.push_back(Point_4(coords.begin(), coords.end())); + + // ---------------------------------------------------------------------------- + // Init of an alpha complex from the list of points + // ---------------------------------------------------------------------------- + Gudhi::alpha_complex::Alpha_complex<Kernel_4> alpha_complex_from_points(points); + + std::cout << "========== Alpha_complex_from_points ==========" << std::endl; + + Gudhi::Simplex_tree<> simplex_tree; + BOOST_CHECK(alpha_complex_from_points.create_complex(simplex_tree)); + + std::cout << "alpha_complex_from_points.number_of_vertices()=" << alpha_complex_from_points.number_of_vertices() + << std::endl; + BOOST_CHECK(alpha_complex_from_points.number_of_vertices() == points.size()); + + // Another way to check num_simplices + std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + int num_simplices = 0; + for (auto f_simplex : simplex_tree.filtration_simplex_range()) { + num_simplices++; + std::cout << " ( "; + for (auto vertex : simplex_tree.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " << "[" << simplex_tree.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + BOOST_CHECK(num_simplices == 15); + std::cout << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl; + BOOST_CHECK(simplex_tree.num_simplices() == 15); + + std::cout << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl; + BOOST_CHECK(simplex_tree.dimension() == 3); + std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; + BOOST_CHECK(simplex_tree.num_vertices() == 4); + + for (auto f_simplex : simplex_tree.filtration_simplex_range()) { + switch (simplex_tree.dimension(f_simplex)) { + case 0: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(f_simplex), 0.0); + break; + case 1: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(f_simplex), 1.0/2.0); + break; + case 2: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(f_simplex), 2.0/3.0); + break; + case 3: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(f_simplex), 3.0/4.0); + break; + default: + BOOST_CHECK(false); // Shall not happen + break; + } + } + + Point_4 p0 = alpha_complex_from_points.get_point(0); + std::cout << "alpha_complex_from_points.get_point(0)=" << p0 << std::endl; + BOOST_CHECK(4 == p0.dimension()); + BOOST_CHECK(is_point_in_list(points, p0)); + + Point_4 p1 = alpha_complex_from_points.get_point(1); + std::cout << "alpha_complex_from_points.get_point(1)=" << p1 << std::endl; + BOOST_CHECK(4 == p1.dimension()); + BOOST_CHECK(is_point_in_list(points, p1)); + + Point_4 p2 = alpha_complex_from_points.get_point(2); + std::cout << "alpha_complex_from_points.get_point(2)=" << p2 << std::endl; + BOOST_CHECK(4 == p2.dimension()); + BOOST_CHECK(is_point_in_list(points, p2)); + + Point_4 p3 = alpha_complex_from_points.get_point(3); + std::cout << "alpha_complex_from_points.get_point(3)=" << p3 << std::endl; + BOOST_CHECK(4 == p3.dimension()); + BOOST_CHECK(is_point_in_list(points, p3)); + + // Test to the limit + BOOST_CHECK_THROW (alpha_complex_from_points.get_point(4), std::out_of_range); + BOOST_CHECK_THROW (alpha_complex_from_points.get_point(-1), std::out_of_range); + BOOST_CHECK_THROW (alpha_complex_from_points.get_point(1234), std::out_of_range); + + // Test after prune_above_filtration + bool modified = simplex_tree.prune_above_filtration(0.6); + if (modified) { + simplex_tree.initialize_filtration(); + } + BOOST_CHECK(modified); + + // Another way to check num_simplices + std::cout << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + num_simplices = 0; + for (auto f_simplex : simplex_tree.filtration_simplex_range()) { + num_simplices++; + std::cout << " ( "; + for (auto vertex : simplex_tree.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " << "[" << simplex_tree.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + BOOST_CHECK(num_simplices == 10); + std::cout << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl; + BOOST_CHECK(simplex_tree.num_simplices() == 10); + + std::cout << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl; + BOOST_CHECK(simplex_tree.dimension() == 1); + std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; + BOOST_CHECK(simplex_tree.num_vertices() == 4); + + for (auto f_simplex : simplex_tree.filtration_simplex_range()) { + switch (simplex_tree.dimension(f_simplex)) { + case 0: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(f_simplex), 0.0); + break; + case 1: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(simplex_tree.filtration(f_simplex), 1.0/2.0); + break; + default: + BOOST_CHECK(false); // Shall not happen + break; + } + } + +} + +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_empty_points, TestedKernel, list_of_kernel_variants) { + std::cout << "========== Alpha_complex_from_empty_points ==========" << std::endl; + + // ---------------------------------------------------------------------------- + // Init of an empty list of points + // ---------------------------------------------------------------------------- + std::vector<typename TestedKernel::Point_d> points; + + // ---------------------------------------------------------------------------- + // Init of an alpha complex from the list of points + // ---------------------------------------------------------------------------- + Gudhi::alpha_complex::Alpha_complex<TestedKernel> alpha_complex_from_points(points); + + // Test to the limit + BOOST_CHECK_THROW (alpha_complex_from_points.get_point(0), std::out_of_range); + + Gudhi::Simplex_tree<> simplex_tree; + BOOST_CHECK(!alpha_complex_from_points.create_complex(simplex_tree)); + + std::cout << "alpha_complex_from_points.number_of_vertices()=" << alpha_complex_from_points.number_of_vertices() + << std::endl; + BOOST_CHECK(alpha_complex_from_points.number_of_vertices() == points.size()); + + std::cout << "simplex_tree.num_simplices()=" << simplex_tree.num_simplices() << std::endl; + BOOST_CHECK(simplex_tree.num_simplices() == 0); + + std::cout << "simplex_tree.dimension()=" << simplex_tree.dimension() << std::endl; + BOOST_CHECK(simplex_tree.dimension() == -1); + + std::cout << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; + BOOST_CHECK(simplex_tree.num_vertices() == 0); +} diff --git a/src/Alpha_complex/test/CMakeLists.txt b/src/Alpha_complex/test/CMakeLists.txt new file mode 100644 index 00000000..ad5b6314 --- /dev/null +++ b/src/Alpha_complex/test/CMakeLists.txt @@ -0,0 +1,37 @@ +project(Alpha_complex_tests) + +include(GUDHI_test_coverage) +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + + # Do not forget to copy test files in current binary dir + file(COPY "${CMAKE_SOURCE_DIR}/data/points/alphacomplexdoc.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + + add_executable ( Alpha_complex_test_unit Alpha_complex_unit_test.cpp ) + target_link_libraries(Alpha_complex_test_unit ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Alpha_complex_test_unit ${TBB_LIBRARIES}) + endif() + + gudhi_add_coverage_test(Alpha_complex_test_unit) + + add_executable ( Alpha_complex_3d_test_unit Alpha_complex_3d_unit_test.cpp ) + target_link_libraries(Alpha_complex_3d_test_unit ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + add_executable ( Weighted_alpha_complex_3d_test_unit Weighted_alpha_complex_3d_unit_test.cpp ) + target_link_libraries(Weighted_alpha_complex_3d_test_unit ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + add_executable ( Periodic_alpha_complex_3d_test_unit Periodic_alpha_complex_3d_unit_test.cpp ) + target_link_libraries(Periodic_alpha_complex_3d_test_unit ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + add_executable ( Weighted_periodic_alpha_complex_3d_test_unit Weighted_periodic_alpha_complex_3d_unit_test.cpp ) + target_link_libraries(Weighted_periodic_alpha_complex_3d_test_unit ${CGAL_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + if (TBB_FOUND) + target_link_libraries(Alpha_complex_3d_test_unit ${TBB_LIBRARIES}) + target_link_libraries(Weighted_alpha_complex_3d_test_unit ${TBB_LIBRARIES}) + target_link_libraries(Periodic_alpha_complex_3d_test_unit ${TBB_LIBRARIES}) + target_link_libraries(Weighted_periodic_alpha_complex_3d_test_unit ${TBB_LIBRARIES}) + endif() + + gudhi_add_coverage_test(Alpha_complex_3d_test_unit) + gudhi_add_coverage_test(Weighted_alpha_complex_3d_test_unit) + gudhi_add_coverage_test(Periodic_alpha_complex_3d_test_unit) + gudhi_add_coverage_test(Weighted_periodic_alpha_complex_3d_test_unit) + +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) diff --git a/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp new file mode 100644 index 00000000..ac3791a4 --- /dev/null +++ b/src/Alpha_complex/test/Periodic_alpha_complex_3d_unit_test.cpp @@ -0,0 +1,178 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "alpha_complex_3d" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <cmath> // float comparison +#include <limits> +#include <string> +#include <vector> +#include <random> +#include <cstddef> // for std::size_t + +#include <gudhi/Alpha_complex_3d.h> +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Unitary_tests_utils.h> +// to construct Alpha_complex from a OFF file of points +#include <gudhi/Points_3D_off_io.h> + +#include <CGAL/Random.h> +#include <CGAL/point_generators_3.h> + +using Fast_periodic_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::FAST, false, true>; +using Safe_periodic_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, false, true>; +using Exact_periodic_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, false, true>; + +#ifdef GUDHI_DEBUG +typedef boost::mpl::list<Fast_periodic_alpha_complex_3d, Safe_periodic_alpha_complex_3d, + Exact_periodic_alpha_complex_3d> + periodic_variants_type_list; + +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_periodic_throw, Periodic_alpha_complex_3d, periodic_variants_type_list) { + std::cout << "Periodic alpha complex 3d exception throw" << std::endl; + using Point_3 = typename Periodic_alpha_complex_3d::Point_3; + std::vector<Point_3> p_points; + + // Not important, this is not what we want to check + p_points.push_back(Point_3(0.0, 0.0, 0.0)); + + std::cout << "Check exception throw in debug mode" << std::endl; + // Check it throws an exception when the cuboid is not iso + BOOST_CHECK_THROW(Periodic_alpha_complex_3d periodic_alpha_complex(p_points, 0., 0., 0., 0.9, 1., 1.), + std::invalid_argument); + BOOST_CHECK_THROW(Periodic_alpha_complex_3d periodic_alpha_complex(p_points, 0., 0., 0., 1., 0.9, 1.), + std::invalid_argument); + BOOST_CHECK_THROW(Periodic_alpha_complex_3d periodic_alpha_complex(p_points, 0., 0., 0., 1., 1., 0.9), + std::invalid_argument); + BOOST_CHECK_THROW(Periodic_alpha_complex_3d periodic_alpha_complex(p_points, 0., 0., 0., 1.1, 1., 1.), + std::invalid_argument); + BOOST_CHECK_THROW(Periodic_alpha_complex_3d periodic_alpha_complex(p_points, 0., 0., 0., 1., 1.1, 1.), + std::invalid_argument); + BOOST_CHECK_THROW(Periodic_alpha_complex_3d periodic_alpha_complex(p_points, 0., 0., 0., 1., 1., 1.1), + std::invalid_argument); +} +#endif + +BOOST_AUTO_TEST_CASE(Alpha_complex_periodic) { + // --------------------- + // Fast periodic version + // --------------------- + std::cout << "Fast periodic alpha complex 3d" << std::endl; + + using Creator = CGAL::Creator_uniform_3<double, Fast_periodic_alpha_complex_3d::Point_3>; + CGAL::Random random(7); + CGAL::Random_points_in_cube_3<Fast_periodic_alpha_complex_3d::Point_3, Creator> in_cube(1, random); + std::vector<Fast_periodic_alpha_complex_3d::Point_3> p_points; + + for (int i = 0; i < 50; i++) { + Fast_periodic_alpha_complex_3d::Point_3 p = *in_cube++; + p_points.push_back(p); + } + + Fast_periodic_alpha_complex_3d periodic_alpha_complex(p_points, -1., -1., -1., 1., 1., 1.); + + Gudhi::Simplex_tree<> stree; + periodic_alpha_complex.create_complex(stree); + + // ---------------------- + // Exact periodic version + // ---------------------- + std::cout << "Exact periodic alpha complex 3d" << std::endl; + + std::vector<Exact_periodic_alpha_complex_3d::Point_3> e_p_points; + + for (auto p : p_points) { + e_p_points.push_back(Exact_periodic_alpha_complex_3d::Point_3(p[0], p[1], p[2])); + } + + Exact_periodic_alpha_complex_3d exact_alpha_complex(e_p_points, -1., -1., -1., 1., 1., 1.); + + Gudhi::Simplex_tree<> exact_stree; + exact_alpha_complex.create_complex(exact_stree); + + // --------------------- + // Compare both versions + // --------------------- + std::cout << "Exact periodic alpha complex 3d is of dimension " << exact_stree.dimension() << " - Non exact is " + << stree.dimension() << std::endl; + BOOST_CHECK(exact_stree.dimension() == stree.dimension()); + std::cout << "Exact periodic alpha complex 3d num_simplices " << exact_stree.num_simplices() << " - Non exact is " + << stree.num_simplices() << std::endl; + BOOST_CHECK(exact_stree.num_simplices() == stree.num_simplices()); + std::cout << "Exact periodic alpha complex 3d num_vertices " << exact_stree.num_vertices() << " - Non exact is " + << stree.num_vertices() << std::endl; + BOOST_CHECK(exact_stree.num_vertices() == stree.num_vertices()); + + // We cannot compare as objects from dispatcher on the alpha shape is not deterministic. + // cf. https://github.com/CGAL/cgal/issues/3346 + auto sh = stree.filtration_simplex_range().begin(); + auto sh_exact = exact_stree.filtration_simplex_range().begin(); + + while (sh != stree.filtration_simplex_range().end() || sh_exact != exact_stree.filtration_simplex_range().end()) { + GUDHI_TEST_FLOAT_EQUALITY_CHECK(stree.filtration(*sh), exact_stree.filtration(*sh_exact), 1e-14); + + std::vector<int> vh(stree.simplex_vertex_range(*sh).begin(), stree.simplex_vertex_range(*sh).end()); + std::vector<int> exact_vh(exact_stree.simplex_vertex_range(*sh_exact).begin(), + exact_stree.simplex_vertex_range(*sh_exact).end()); + + BOOST_CHECK(vh.size() == exact_vh.size()); + ++sh; + ++sh_exact; + } + + BOOST_CHECK(sh == stree.filtration_simplex_range().end()); + BOOST_CHECK(sh_exact == exact_stree.filtration_simplex_range().end()); + + // ---------------------- + // Safe periodic version + // ---------------------- + std::cout << "Safe periodic alpha complex 3d" << std::endl; + + std::vector<Safe_periodic_alpha_complex_3d::Point_3> s_p_points; + + for (auto p : p_points) { + s_p_points.push_back(Safe_periodic_alpha_complex_3d::Point_3(p[0], p[1], p[2])); + } + + Safe_periodic_alpha_complex_3d safe_alpha_complex(s_p_points, -1., -1., -1., 1., 1., 1.); + + Gudhi::Simplex_tree<> safe_stree; + safe_alpha_complex.create_complex(safe_stree); + + // --------------------- + // Compare both versions + // --------------------- + // We cannot compare as objects from dispatcher on the alpha shape is not deterministic. + // cf. https://github.com/CGAL/cgal/issues/3346 + sh = stree.filtration_simplex_range().begin(); + auto sh_safe = safe_stree.filtration_simplex_range().begin(); + + while (sh != stree.filtration_simplex_range().end() || sh_safe != safe_stree.filtration_simplex_range().end()) { + GUDHI_TEST_FLOAT_EQUALITY_CHECK(stree.filtration(*sh), safe_stree.filtration(*sh_safe), 1e-14); + + std::vector<int> vh(stree.simplex_vertex_range(*sh).begin(), stree.simplex_vertex_range(*sh).end()); + std::vector<int> safe_vh(safe_stree.simplex_vertex_range(*sh_safe).begin(), + safe_stree.simplex_vertex_range(*sh_safe).end()); + + BOOST_CHECK(vh.size() == safe_vh.size()); + ++sh; + ++sh_safe; + } + + BOOST_CHECK(sh == stree.filtration_simplex_range().end()); + BOOST_CHECK(sh_safe == safe_stree.filtration_simplex_range().end()); +} diff --git a/src/Alpha_complex/test/README b/src/Alpha_complex/test/README new file mode 100644 index 00000000..0e5b9eb1 --- /dev/null +++ b/src/Alpha_complex/test/README @@ -0,0 +1,12 @@ +To compile: +*********** + +cmake . +make + +To launch with details: +*********************** + +./Alpha_complex_unit_test --report_level=detailed --log_level=all + + ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp new file mode 100644 index 00000000..44deb930 --- /dev/null +++ b/src/Alpha_complex/test/Weighted_alpha_complex_3d_unit_test.cpp @@ -0,0 +1,135 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "alpha_complex_3d" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <cmath> // float comparison +#include <limits> +#include <string> +#include <vector> +#include <random> +#include <cstddef> // for std::size_t + +#include <gudhi/Alpha_complex_3d.h> +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Unitary_tests_utils.h> +// to construct Alpha_complex from a OFF file of points +#include <gudhi/Points_3D_off_io.h> + +#include <CGAL/Random.h> +#include <CGAL/point_generators_3.h> + +using Fast_weighted_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::FAST, true, false>; +using Safe_weighted_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, false>; +using Exact_weighted_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, true, false>; + +typedef boost::mpl::list<Fast_weighted_alpha_complex_3d, Safe_weighted_alpha_complex_3d, + Exact_weighted_alpha_complex_3d> + weighted_variants_type_list; + +#ifdef GUDHI_DEBUG +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_throw, Weighted_alpha_complex_3d, weighted_variants_type_list) { + using Point_3 = typename Weighted_alpha_complex_3d::Point_3; + std::vector<Point_3> w_points; + w_points.push_back(Point_3(0.0, 0.0, 0.0)); + w_points.push_back(Point_3(0.0, 0.0, 0.2)); + w_points.push_back(Point_3(0.2, 0.0, 0.2)); + // w_points.push_back(Point_3(0.6, 0.6, 0.0)); + // w_points.push_back(Point_3(0.8, 0.8, 0.2)); + // w_points.push_back(Point_3(0.2, 0.8, 0.6)); + + // weights size is different from w_points size to make weighted Alpha_complex_3d throw in debug mode + std::vector<double> weights = {0.01, 0.005, 0.006, 0.01, 0.009, 0.001}; + + std::cout << "Check exception throw in debug mode" << std::endl; + BOOST_CHECK_THROW(Weighted_alpha_complex_3d wac(w_points, weights), std::invalid_argument); +} +#endif + +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted, Weighted_alpha_complex_3d, weighted_variants_type_list) { + std::cout << "Weighted alpha complex 3d from points and weights" << std::endl; + using Point_3 = typename Weighted_alpha_complex_3d::Point_3; + std::vector<Point_3> w_points; + w_points.push_back(Point_3(0.0, 0.0, 0.0)); + w_points.push_back(Point_3(0.0, 0.0, 0.2)); + w_points.push_back(Point_3(0.2, 0.0, 0.2)); + w_points.push_back(Point_3(0.6, 0.6, 0.0)); + w_points.push_back(Point_3(0.8, 0.8, 0.2)); + w_points.push_back(Point_3(0.2, 0.8, 0.6)); + + // weights size is different from w_points size to make weighted Alpha_complex_3d throw in debug mode + std::vector<double> weights = {0.01, 0.005, 0.006, 0.01, 0.009, 0.001}; + + Weighted_alpha_complex_3d alpha_complex_p_a_w(w_points, weights); + Gudhi::Simplex_tree<> stree; + alpha_complex_p_a_w.create_complex(stree); + + std::cout << "Weighted alpha complex 3d from weighted points" << std::endl; + using Weighted_point_3 = typename Weighted_alpha_complex_3d::Weighted_point_3; + + std::vector<Weighted_point_3> weighted_points; + + for (std::size_t i = 0; i < w_points.size(); i++) { + weighted_points.push_back(Weighted_point_3(w_points[i], weights[i])); + } + Weighted_alpha_complex_3d alpha_complex_w_p(weighted_points); + + Gudhi::Simplex_tree<> stree_bis; + alpha_complex_w_p.create_complex(stree_bis); + + // --------------------- + // Compare both versions + // --------------------- + std::cout << "Weighted alpha complex 3d is of dimension " << stree_bis.dimension() << " - versus " + << stree.dimension() << std::endl; + BOOST_CHECK(stree_bis.dimension() == stree.dimension()); + std::cout << "Weighted alpha complex 3d num_simplices " << stree_bis.num_simplices() << " - versus " + << stree.num_simplices() << std::endl; + BOOST_CHECK(stree_bis.num_simplices() == stree.num_simplices()); + std::cout << "Weighted alpha complex 3d num_vertices " << stree_bis.num_vertices() << " - versus " + << stree.num_vertices() << std::endl; + BOOST_CHECK(stree_bis.num_vertices() == stree.num_vertices()); + + auto sh = stree.filtration_simplex_range().begin(); + while (sh != stree.filtration_simplex_range().end()) { + std::vector<int> simplex; + std::vector<int> exact_simplex; +#ifdef DEBUG_TRACES + std::cout << " ( "; +#endif + for (auto vertex : stree.simplex_vertex_range(*sh)) { + simplex.push_back(vertex); +#ifdef DEBUG_TRACES + std::cout << vertex << " "; +#endif + } +#ifdef DEBUG_TRACES + std::cout << ") -> " + << "[" << stree.filtration(*sh) << "] "; + std::cout << std::endl; +#endif + + // Find it in the exact structure + auto sh_exact = stree_bis.find(simplex); + BOOST_CHECK(sh_exact != stree_bis.null_simplex()); + + // Exact and non-exact version is not exactly the same due to float comparison + GUDHI_TEST_FLOAT_EQUALITY_CHECK(stree_bis.filtration(sh_exact), stree.filtration(*sh)); + + ++sh; + } +} diff --git a/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp new file mode 100644 index 00000000..670c7799 --- /dev/null +++ b/src/Alpha_complex/test/Weighted_periodic_alpha_complex_3d_unit_test.cpp @@ -0,0 +1,227 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "alpha_complex_3d" +#include <boost/test/unit_test.hpp> +#include <boost/mpl/list.hpp> + +#include <cmath> // float comparison +#include <limits> +#include <string> +#include <vector> +#include <random> +#include <cstddef> // for std::size_t + +#include <gudhi/Alpha_complex_3d.h> +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Unitary_tests_utils.h> +// to construct Alpha_complex from a OFF file of points +#include <gudhi/Points_3D_off_io.h> + +#include <CGAL/Random.h> +#include <CGAL/point_generators_3.h> + +using Fast_weighted_periodic_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::FAST, true, true>; +using Safe_weighted_periodic_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, true>; +using Exact_weighted_periodic_alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, true, true>; + + +typedef boost::mpl::list<Fast_weighted_periodic_alpha_complex_3d, Exact_weighted_periodic_alpha_complex_3d, + Safe_weighted_periodic_alpha_complex_3d> + wp_variants_type_list; + +#ifdef GUDHI_DEBUG +BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_weighted_periodic_throw, Weighted_periodic_alpha_complex_3d, + wp_variants_type_list) { + std::cout << "Weighted periodic alpha complex 3d exception throw" << std::endl; + + using Creator = CGAL::Creator_uniform_3<double, typename Weighted_periodic_alpha_complex_3d::Point_3>; + CGAL::Random random(7); + CGAL::Random_points_in_cube_3<typename Weighted_periodic_alpha_complex_3d::Point_3, Creator> in_cube(1, random); + std::vector<typename Weighted_periodic_alpha_complex_3d::Point_3> wp_points; + + for (int i = 0; i < 50; i++) { + typename Weighted_periodic_alpha_complex_3d::Point_3 p = *in_cube++; + wp_points.push_back(p); + } + std::vector<double> p_weights; + // Weights must be in range ]0, 1/64 = 0.015625[ + for (std::size_t i = 0; i < wp_points.size(); ++i) { + p_weights.push_back(random.get_double(0., 0.01)); + } + + std::cout << "Cuboid is not iso exception" << std::endl; + // Check it throws an exception when the cuboid is not iso + BOOST_CHECK_THROW( + Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, -1., -1., -1., 0.9, 1., 1.), + std::invalid_argument); + BOOST_CHECK_THROW( + Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, -1., -1., -1., 1., 0.9, 1.), + std::invalid_argument); + BOOST_CHECK_THROW( + Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, -1., -1., -1., 1., 1., 0.9), + std::invalid_argument); + BOOST_CHECK_THROW( + Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, -1., -1., -1., 1.1, 1., 1.), + std::invalid_argument); + BOOST_CHECK_THROW( + Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, -1., -1., -1., 1., 1.1, 1.), + std::invalid_argument); + BOOST_CHECK_THROW( + Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, -1., -1., -1., 1., 1., 1.1), + std::invalid_argument); + + std::cout << "0 <= point.weight() < 1/64 * domain_size * domain_size exception" << std::endl; + // Weights must be in range ]0, 1/64 = 0.015625[ + double temp = p_weights[25]; + p_weights[25] = 1.0; + BOOST_CHECK_THROW(Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, 0., 0., 0., 1., 1., 1.), + std::invalid_argument); + // Weights must be in range ]0, 1/64 = 0.015625[ + p_weights[25] = temp; + temp = p_weights[14]; + p_weights[14] = -1e-10; + BOOST_CHECK_THROW(Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, 0., 0., 0., 1., 1., 1.), + std::invalid_argument); + p_weights[14] = temp; + + std::cout << "wp_points and p_weights size exception" << std::endl; + // Weights and points must have the same size + // + 1 + p_weights.push_back(1e-10); + BOOST_CHECK_THROW(Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, 0., 0., 0., 1., 1., 1.), + std::invalid_argument); + // - 1 + p_weights.pop_back(); + p_weights.pop_back(); + BOOST_CHECK_THROW(Weighted_periodic_alpha_complex_3d wp_alpha_complex(wp_points, p_weights, 0., 0., 0., 1., 1., 1.), + std::invalid_argument); +} +#endif + +BOOST_AUTO_TEST_CASE(Alpha_complex_weighted_periodic) { + // --------------------- + // Fast weighted periodic version + // --------------------- + std::cout << "Fast weighted periodic alpha complex 3d" << std::endl; + + using Creator = CGAL::Creator_uniform_3<double, Fast_weighted_periodic_alpha_complex_3d::Point_3>; + CGAL::Random random(7); + CGAL::Random_points_in_cube_3<Fast_weighted_periodic_alpha_complex_3d::Point_3, Creator> in_cube(1, random); + std::vector<Fast_weighted_periodic_alpha_complex_3d::Point_3> p_points; + + for (int i = 0; i < 50; i++) { + Fast_weighted_periodic_alpha_complex_3d::Point_3 p = *in_cube++; + p_points.push_back(p); + } + std::vector<double> p_weights; + // Weights must be in range ]0, 1/64 = 0.015625[ + for (std::size_t i = 0; i < p_points.size(); ++i) { + p_weights.push_back(random.get_double(0., 0.01)); + } + + Fast_weighted_periodic_alpha_complex_3d periodic_alpha_complex(p_points, p_weights, -1., -1., -1., 1., 1., 1.); + + Gudhi::Simplex_tree<> stree; + periodic_alpha_complex.create_complex(stree); + + // ---------------------- + // Exact weighted periodic version + // ---------------------- + std::cout << "Exact weighted periodic alpha complex 3d" << std::endl; + + std::vector<Exact_weighted_periodic_alpha_complex_3d::Point_3> e_p_points; + + for (auto p : p_points) { + e_p_points.push_back(Exact_weighted_periodic_alpha_complex_3d::Point_3(p[0], p[1], p[2])); + } + + Exact_weighted_periodic_alpha_complex_3d exact_alpha_complex(e_p_points, p_weights, -1., -1., -1., 1., 1., 1.); + + Gudhi::Simplex_tree<> exact_stree; + exact_alpha_complex.create_complex(exact_stree); + + // --------------------- + // Compare both versions + // --------------------- + std::cout << "Exact weighted periodic alpha complex 3d is of dimension " << exact_stree.dimension() + << " - Non exact is " << stree.dimension() << std::endl; + BOOST_CHECK(exact_stree.dimension() == stree.dimension()); + std::cout << "Exact weighted periodic alpha complex 3d num_simplices " << exact_stree.num_simplices() + << " - Non exact is " << stree.num_simplices() << std::endl; + BOOST_CHECK(exact_stree.num_simplices() == stree.num_simplices()); + std::cout << "Exact weighted periodic alpha complex 3d num_vertices " << exact_stree.num_vertices() + << " - Non exact is " << stree.num_vertices() << std::endl; + BOOST_CHECK(exact_stree.num_vertices() == stree.num_vertices()); + + // We cannot compare as objects from dispatcher on the alpha shape is not deterministic. + // cf. https://github.com/CGAL/cgal/issues/3346 + auto sh = stree.filtration_simplex_range().begin(); + auto sh_exact = exact_stree.filtration_simplex_range().begin(); + + while (sh != stree.filtration_simplex_range().end() || sh_exact != exact_stree.filtration_simplex_range().end()) { + GUDHI_TEST_FLOAT_EQUALITY_CHECK(stree.filtration(*sh), exact_stree.filtration(*sh_exact), 1e-14); + + std::vector<int> vh(stree.simplex_vertex_range(*sh).begin(), stree.simplex_vertex_range(*sh).end()); + std::vector<int> exact_vh(exact_stree.simplex_vertex_range(*sh_exact).begin(), + exact_stree.simplex_vertex_range(*sh_exact).end()); + + BOOST_CHECK(vh.size() == exact_vh.size()); + ++sh; + ++sh_exact; + } + + BOOST_CHECK(sh == stree.filtration_simplex_range().end()); + BOOST_CHECK(sh_exact == exact_stree.filtration_simplex_range().end()); + + // ---------------------- + // Safe weighted periodic version + // ---------------------- + std::cout << "Safe weighted periodic alpha complex 3d" << std::endl; + + std::vector<Safe_weighted_periodic_alpha_complex_3d::Point_3> s_p_points; + + for (auto p : p_points) { + s_p_points.push_back(Safe_weighted_periodic_alpha_complex_3d::Point_3(p[0], p[1], p[2])); + } + + Safe_weighted_periodic_alpha_complex_3d safe_alpha_complex(s_p_points, p_weights, -1., -1., -1., 1., 1., 1.); + + Gudhi::Simplex_tree<> safe_stree; + safe_alpha_complex.create_complex(safe_stree); + + // --------------------- + // Compare both versions + // --------------------- + // We cannot compare as objects from dispatcher on the alpha shape is not deterministic. + // cf. https://github.com/CGAL/cgal/issues/3346 + sh = stree.filtration_simplex_range().begin(); + auto sh_safe = safe_stree.filtration_simplex_range().begin(); + + while (sh != stree.filtration_simplex_range().end() || sh_safe != safe_stree.filtration_simplex_range().end()) { + GUDHI_TEST_FLOAT_EQUALITY_CHECK(stree.filtration(*sh), safe_stree.filtration(*sh_safe), 1e-14); + + std::vector<int> vh(stree.simplex_vertex_range(*sh).begin(), stree.simplex_vertex_range(*sh).end()); + std::vector<int> safe_vh(safe_stree.simplex_vertex_range(*sh_safe).begin(), + safe_stree.simplex_vertex_range(*sh_safe).end()); + + BOOST_CHECK(vh.size() == safe_vh.size()); + ++sh; + ++sh_safe; + } + + BOOST_CHECK(sh == stree.filtration_simplex_range().end()); + BOOST_CHECK(sh_safe == safe_stree.filtration_simplex_range().end()); +} diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt new file mode 100644 index 00000000..5295f3cd --- /dev/null +++ b/src/Alpha_complex/utilities/CMakeLists.txt @@ -0,0 +1,58 @@ +project(Alpha_complex_utilities) + +if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + add_executable (alpha_complex_persistence alpha_complex_persistence.cpp) + target_link_libraries(alpha_complex_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY}) + + 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") + + install(TARGETS alpha_complex_persistence DESTINATION bin) + + add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp) + target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_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 COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" + "-p" "2" "-m" "0.45" "-o" "safe.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") + + 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") + + 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") + endif() + + add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" + "-c" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" + "-p" "2" "-m" "0") + + add_test(NAME Alpha_complex_utilities_weighted_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" + "-w" "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.weights" + "-p" "2" "-m" "0") + + add_test(NAME Alpha_complex_utilities_weighted_periodic_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence> + "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" + "-w" "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.weights" + "-c" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" + "-p" "2" "-m" "0" "-e") + + install(TARGETS alpha_complex_3d_persistence DESTINATION bin) + +endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) 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..2272576e --- /dev/null +++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp @@ -0,0 +1,305 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2014 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <boost/program_options.hpp> +#include <boost/variant.hpp> + +#include <gudhi/Alpha_complex_3d.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_3D_off_io.h> + +#include <fstream> +#include <string> +#include <vector> +#include <limits> // for numeric_limits<> + +// gudhi type definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Persistent_cohomology = + Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp>; + +void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &safe, + std::string &weight_file, std::string &cuboid_file, std::string &output_file_diag, + Filtration_value &alpha_square_max_value, int &coeff_field_characteristic, + Filtration_value &min_persistence); + +bool read_weight_file(const std::string &weight_file, std::vector<double> &weights) { + // Read weights information from file + std::ifstream weights_ifstr(weight_file); + if (weights_ifstr.good()) { + double weight = 0.0; + // Attempt read the weight in a double format, return false if it fails + while (weights_ifstr >> weight) { + weights.push_back(weight); + } + } else { + return false; + } + return true; +} + +bool read_cuboid_file(const std::string &cuboid_file, double &x_min, double &y_min, double &z_min, double &x_max, + double &y_max, double &z_max) { + // Read weights information from file + std::ifstream iso_cuboid_str(cuboid_file); + if (iso_cuboid_str.is_open()) { + if (!(iso_cuboid_str >> x_min >> y_min >> z_min >> x_max >> y_max >> z_max)) { + return false; + } + } else { + return false; + } + return true; +} + +template <typename AlphaComplex3d> +std::vector<typename AlphaComplex3d::Point_3> read_off(const std::string &off_file_points) { + // Read the OFF file (input file name given as parameter) and triangulate points + Gudhi::Points_3D_off_reader<typename AlphaComplex3d::Point_3> off_reader(off_file_points); + // Check the read operation was correct + if (!off_reader.is_valid()) { + std::cerr << "Unable to read OFF file " << off_file_points << std::endl; + exit(-1); + } + return off_reader.get_point_cloud(); +} + +int main(int argc, char **argv) { + std::string off_file_points; + std::string weight_file; + std::string cuboid_file; + std::string output_file_diag; + Filtration_value alpha_square_max_value = 0.; + int coeff_field_characteristic = 0; + Filtration_value min_persistence = 0.; + bool exact_version = false; + bool fast_version = false; + bool weighted_version = false; + bool periodic_version = false; + + program_options(argc, argv, off_file_points, exact_version, fast_version, weight_file, cuboid_file, output_file_diag, + alpha_square_max_value, coeff_field_characteristic, min_persistence); + + std::vector<double> weights; + if (weight_file != std::string()) { + if (!read_weight_file(weight_file, weights)) { + std::cerr << "Unable to read weights file " << weight_file << std::endl; + exit(-1); + } + weighted_version = true; + } + + double x_min = 0., y_min = 0., z_min = 0., x_max = 0., y_max = 0., z_max = 0.; + std::ifstream iso_cuboid_str(argv[3]); + if (cuboid_file != std::string()) { + if (!read_cuboid_file(cuboid_file, x_min, y_min, z_min, x_max, y_max, z_max)) { + std::cerr << "Unable to read cuboid file " << cuboid_file << std::endl; + exit(-1); + } + periodic_version = true; + } + + Gudhi::alpha_complex::complexity complexity = Gudhi::alpha_complex::complexity::SAFE; + if (exact_version) { + if (fast_version) { + std::cerr << "You cannot set the exact and the fast version." << std::endl; + exit(-1); + } + complexity = Gudhi::alpha_complex::complexity::EXACT; + } + if (fast_version) { + complexity = Gudhi::alpha_complex::complexity::FAST; + } + + Simplex_tree simplex_tree; + + switch (complexity) { + case Gudhi::alpha_complex::complexity::FAST: + if (weighted_version) { + if (periodic_version) { + using Alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::FAST, true, true>; + auto points = read_off<Alpha_complex_3d>(off_file_points); + Alpha_complex_3d alpha_complex(points, weights, x_min, y_min, z_min, x_max, y_max, z_max); + alpha_complex.create_complex(simplex_tree, alpha_square_max_value); + } else { + using Alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::FAST, true, false>; + auto points = read_off<Alpha_complex_3d>(off_file_points); + Alpha_complex_3d alpha_complex(points, weights); + alpha_complex.create_complex(simplex_tree, alpha_square_max_value); + } + } else { + if (periodic_version) { + using Alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::FAST, false, true>; + auto points = read_off<Alpha_complex_3d>(off_file_points); + Alpha_complex_3d alpha_complex(points, x_min, y_min, z_min, x_max, y_max, z_max); + alpha_complex.create_complex(simplex_tree, alpha_square_max_value); + } else { + using Alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::FAST, false, false>; + auto points = read_off<Alpha_complex_3d>(off_file_points); + Alpha_complex_3d alpha_complex(points); + alpha_complex.create_complex(simplex_tree, alpha_square_max_value); + } + } + break; + case Gudhi::alpha_complex::complexity::EXACT: + if (weighted_version) { + if (periodic_version) { + using Alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, true, true>; + auto points = read_off<Alpha_complex_3d>(off_file_points); + Alpha_complex_3d alpha_complex(points, weights, x_min, y_min, z_min, x_max, y_max, z_max); + alpha_complex.create_complex(simplex_tree, alpha_square_max_value); + } else { + using Alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, true, false>; + auto points = read_off<Alpha_complex_3d>(off_file_points); + Alpha_complex_3d alpha_complex(points, weights); + alpha_complex.create_complex(simplex_tree, alpha_square_max_value); + } + } else { + if (periodic_version) { + using Alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, false, true>; + auto points = read_off<Alpha_complex_3d>(off_file_points); + Alpha_complex_3d alpha_complex(points, x_min, y_min, z_min, x_max, y_max, z_max); + alpha_complex.create_complex(simplex_tree, alpha_square_max_value); + } else { + using Alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::EXACT, false, false>; + auto points = read_off<Alpha_complex_3d>(off_file_points); + Alpha_complex_3d alpha_complex(points); + alpha_complex.create_complex(simplex_tree, alpha_square_max_value); + } + } + break; + case Gudhi::alpha_complex::complexity::SAFE: + if (weighted_version) { + if (periodic_version) { + using Alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, true>; + auto points = read_off<Alpha_complex_3d>(off_file_points); + Alpha_complex_3d alpha_complex(points, weights, x_min, y_min, z_min, x_max, y_max, z_max); + alpha_complex.create_complex(simplex_tree, alpha_square_max_value); + } else { + using Alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, true, false>; + auto points = read_off<Alpha_complex_3d>(off_file_points); + Alpha_complex_3d alpha_complex(points, weights); + alpha_complex.create_complex(simplex_tree, alpha_square_max_value); + } + } else { + if (periodic_version) { + using Alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, false, true>; + auto points = read_off<Alpha_complex_3d>(off_file_points); + Alpha_complex_3d alpha_complex(points, x_min, y_min, z_min, x_max, y_max, z_max); + alpha_complex.create_complex(simplex_tree, alpha_square_max_value); + } else { + using Alpha_complex_3d = + Gudhi::alpha_complex::Alpha_complex_3d<Gudhi::alpha_complex::complexity::SAFE, false, false>; + auto points = read_off<Alpha_complex_3d>(off_file_points); + Alpha_complex_3d alpha_complex(points); + alpha_complex.create_complex(simplex_tree, alpha_square_max_value); + } + } + break; + default: + std::cerr << "Unknown complexity value " << std::endl; + exit(-1); + break; + } + + // 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); + + // 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, bool &exact, bool &fast, + std::string &weight_file, std::string &cuboid_file, 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), + "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")( + "exact,e", po::bool_switch(&exact), + "To activate exact version of Alpha complex 3d (default is false, not available if fast is set)")( + "fast,f", po::bool_switch(&fast), + "To activate fast version of Alpha complex 3d (default is false, not available if exact is set)")( + "weight-file,w", po::value<std::string>(&weight_file)->default_value(std::string()), + "Name of file containing a point weights. Format is one weight per line:\n W1\n ...\n Wn ")( + "cuboid-file,c", po::value<std::string>(&cuboid_file), + "Name of file describing the periodic domain. Format is:\n min_hx min_hy min_hz\n max_hx max_hy max_hz")( + "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) + ->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal alpha square value for the Alpha complex construction.")( + "field-charac,p", po::value<int>(&coeff_field_characteristic)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::value<Filtration_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") || !vm.count("weight-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"; + std::cout << "3D Alpha complex can be safe (by default) exact or fast, weighted and/or periodic\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.\n\n"; + + std::cout << "Usage: " << argv[0] << " [options] input-file weight-file\n\n"; + std::cout << visible << std::endl; + exit(-1); + } +} 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..fab7bd30 --- /dev/null +++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp @@ -0,0 +1,126 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Vincent Rouvreau + * + * Copyright (C) 2016 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#include <boost/program_options.hpp> + +#include <CGAL/Epick_d.h> + +#include <gudhi/Alpha_complex.h> +#include <gudhi/Persistent_cohomology.h> +// to construct a simplex_tree from alpha complex +#include <gudhi/Simplex_tree.h> + +#include <iostream> +#include <string> +#include <limits> // 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<Kernel> 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<std::string>(&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<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) + ->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal alpha square value for the Alpha complex construction.")( + "field-charac,p", po::value<int>(&coeff_field_characteristic)->default_value(11), + "Characteristic p of the coefficient field Z/pZ for computing homology.")( + "min-persistence,m", po::value<Filtration_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; + exit(-1); + } +} diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md new file mode 100644 index 00000000..fcd16a3b --- /dev/null +++ b/src/Alpha_complex/utilities/alphacomplex.md @@ -0,0 +1,127 @@ +---
+layout: page
+title: "Alpha complex"
+meta_title: "Alpha complex"
+teaser: ""
+permalink: /alphacomplex/
+---
+{::comment}
+Leave the lines above as it is required by the web site generator 'Jekyll'
+{:/comment}
+
+
+## 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 birth death
+```
+
+where `dim` is the dimension of the homological feature, `birth` and `death`
+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] <input OFF file>
+```
+
+where
+`<input OFF file>` is the path to the input point cloud in
+[nOFF ASCII format]({{ site.officialurl }}/doc/latest/fileformats.html#FileFormatsOFF).
+
+**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-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
+```
+
+N.B.:
+
+* Filtration values are alpha square values.
+
+
+## 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.
+One can use exact computation. It is slower, but it is necessary when points
+are on a grid for instance.
+Alpha complex 3d can be weighted and/or periodic (refer to the
+[CGAL's 3D Periodic Triangulations User Manual](
+https://doc.cgal.org/latest/Periodic_3_triangulation_3/index.html)
+for more details).
+
+The output diagram contains
+one bar per line, written with the convention:
+
+```
+p dim birth death
+```
+
+where `dim` is the dimension of the homological feature, `birth` and `death`
+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 [options] <input OFF file>
+```
+
+where `<input OFF file>` is the path to the input point cloud in
+[nOFF ASCII format]({{ site.officialurl }}/doc/latest/fileformats.html#FileFormatsOFF).
+
+**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-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.
+* `-c [ --cuboid-file ]` is the path to the file describing the periodic domain.
+It must be in the format described
+[here]({{ site.officialurl }}/doc/latest/fileformats.html#FileFormatsIsoCuboid).
+Default version is not periodic.
+* `-w [ --weight-file ]` is the path to the file containing the weights of the
+points (one value per line).
+Default version is not weighted.
+* `-e [ --exact ]` for the exact computation version (not compatible with
+weight and periodic version).
+* `-f [ --fast ]` for the fast computation version.
+
+**Example**
+
+```
+alpha_complex_3d_persistence ../../data/points/tore3D_300.off -p 2 -m 0.45
+```
+
+N.B.:
+
+* `alpha_complex_3d_persistence` only accepts OFF files in dimension 3.
+* Filtration values are alpha square values.
+* 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.
|
