diff options
Diffstat (limited to 'src/Cech_complex')
24 files changed, 2416 insertions, 0 deletions
diff --git a/src/Cech_complex/benchmark/CMakeLists.txt b/src/Cech_complex/benchmark/CMakeLists.txt new file mode 100644 index 00000000..b7697764 --- /dev/null +++ b/src/Cech_complex/benchmark/CMakeLists.txt @@ -0,0 +1,12 @@ +cmake_minimum_required(VERSION 2.6) +project(Cech_complex_benchmark) + +# Do not forget to copy test files in current binary dir +file(COPY "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + +add_executable(cech_complex_benchmark cech_complex_benchmark.cpp) +target_link_libraries(cech_complex_benchmark ${Boost_SYSTEM_LIBRARY} ${Boost_FILESYSTEM_LIBRARY}) + +if (TBB_FOUND) + target_link_libraries(cech_complex_benchmark ${TBB_LIBRARIES}) +endif() diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp new file mode 100644 index 00000000..d2d71dbf --- /dev/null +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -0,0 +1,132 @@ +/* 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 + */ + +#include <gudhi/Points_off_io.h> +#include <gudhi/distance_functions.h> +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/Clock.h> +#include <gudhi/Rips_complex.h> +#include <gudhi/Cech_complex.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Miniball.hpp> + +#include "boost/filesystem.hpp" // includes all needed Boost.Filesystem declarations + +#include <string> +#include <vector> + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Point = std::vector<Filtration_value>; +using Point_cloud = std::vector<Point>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; +using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>; + +class Minimal_enclosing_ball_radius { + public: + // boost::range_value is not SFINAE-friendly so we cannot use it in the return type + template <typename Point> + typename std::iterator_traits<typename boost::range_iterator<Point>::type>::value_type operator()( + const Point& p1, const Point& p2) const { + // Type def + using Point_cloud = std::vector<Point>; + using Point_iterator = typename Point_cloud::const_iterator; + using Coordinate_iterator = typename Point::const_iterator; + using Min_sphere = + typename Gudhi::Miniball::Miniball<Gudhi::Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>; + + Point_cloud point_cloud; + point_cloud.push_back(p1); + point_cloud.push_back(p2); + + GUDHI_CHECK((p1.end() - p1.begin()) != (p2.end() - p2.begin()), "inconsistent point dimensions"); + Min_sphere min_sphere(p1.end() - p1.begin(), point_cloud.begin(), point_cloud.end()); + + return std::sqrt(min_sphere.squared_radius()); + } +}; + +int main(int argc, char* argv[]) { + std::string off_file_points = "tore3D_1307.off"; + Filtration_value threshold = 1e20; + + // Extract the points from the file filepoints + Points_off_reader off_reader(off_file_points); + + Gudhi::Clock euclidean_clock("Gudhi::Euclidean_distance"); + // Compute the proximity graph of the points + Proximity_graph euclidean_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>( + off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + + std::cout << euclidean_clock << std::endl; + + Gudhi::Clock miniball_clock("Minimal_enclosing_ball_radius"); + // Compute the proximity graph of the points + Proximity_graph miniball_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>( + off_reader.get_point_cloud(), threshold, Minimal_enclosing_ball_radius()); + std::cout << miniball_clock << std::endl; + + Gudhi::Clock common_miniball_clock("Gudhi::Minimal_enclosing_ball_radius()"); + // Compute the proximity graph of the points + Proximity_graph common_miniball_prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>( + off_reader.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius()); + std::cout << common_miniball_clock << std::endl; + + boost::filesystem::path full_path(boost::filesystem::current_path()); + std::cout << "Current path is : " << full_path << std::endl; + + std::cout << "File name;Radius;Rips time;Cech time; Ratio Rips/Cech time;Rips nb simplices;Cech nb simplices;" + << std::endl; + boost::filesystem::directory_iterator end_itr; // default construction yields past-the-end + for (boost::filesystem::directory_iterator itr(boost::filesystem::current_path()); itr != end_itr; ++itr) { + if (!boost::filesystem::is_directory(itr->status())) { + if (itr->path().extension() == ".off") // see below + { + Points_off_reader off_reader(itr->path().string()); + Point p0 = off_reader.get_point_cloud()[0]; + + for (Filtration_value radius = 0.1; radius < 0.4; radius += 0.1) { + std::cout << itr->path().stem() << ";"; + std::cout << radius << ";"; + Gudhi::Clock rips_clock("Rips computation"); + Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), radius, + Gudhi::Minimal_enclosing_ball_radius()); + Simplex_tree rips_stree; + rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); + // ------------------------------------------ + // Display information about the Rips complex + // ------------------------------------------ + double rips_sec = rips_clock.num_seconds(); + std::cout << rips_sec << ";"; + + Gudhi::Clock cech_clock("Cech computation"); + Cech_complex cech_complex_from_points(off_reader.get_point_cloud(), radius); + Simplex_tree cech_stree; + cech_complex_from_points.create_complex(cech_stree, p0.size() - 1); + // ------------------------------------------ + // Display information about the Cech complex + // ------------------------------------------ + double cech_sec = cech_clock.num_seconds(); + std::cout << cech_sec << ";"; + std::cout << cech_sec / rips_sec << ";"; + + assert(rips_stree.num_simplices() >= cech_stree.num_simplices()); + std::cout << rips_stree.num_simplices() << ";"; + std::cout << cech_stree.num_simplices() << ";" << std::endl; + } + } + } + } + + return 0; +} diff --git a/src/Cech_complex/concept/SimplicialComplexForCech.h b/src/Cech_complex/concept/SimplicialComplexForCech.h new file mode 100644 index 00000000..00c7df3a --- /dev/null +++ b/src/Cech_complex/concept/SimplicialComplexForCech.h @@ -0,0 +1,54 @@ +/* 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_CECH_COMPLEX_SIMPLICIAL_COMPLEX_FOR_CECH_H_ +#define CONCEPT_CECH_COMPLEX_SIMPLICIAL_COMPLEX_FOR_CECH_H_ + +namespace Gudhi { + +namespace cech_complex { + +/** \brief The concept SimplicialComplexForCech describes the requirements for a type to implement a simplicial + * complex, that can be created from a `Cech_complex`. + */ +struct SimplicialComplexForCech { + /** 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; + + /** Assigns the 'simplex' with the given 'filtration' value. */ + int assign_filtration(Simplex_handle simplex, Filtration_value filtration); + + /** \brief Returns a range over vertices of a given + * simplex. */ + Simplex_vertex_range simplex_vertex_range(Simplex_handle const & simplex); + + /** \brief Inserts a given `Gudhi::ProximityGraph` in the simplicial complex. */ + template<class ProximityGraph> + void insert_graph(const ProximityGraph& proximity_graph); + + /** \brief Expands the simplicial complex containing only its one skeleton until a given maximal dimension. + * expansion can be blocked by the blocker oracle. */ + template< typename Blocker > + void expansion_with_blockers(int max_dim, Blocker block_simplex); + + /** Returns the number of vertices in the simplicial complex. */ + std::size_t num_vertices(); + +}; + +} // namespace alpha_complex + +} // namespace Gudhi + +#endif // CONCEPT_ALPHA_COMPLEX_SIMPLICIAL_COMPLEX_FOR_ALPHA_H_ diff --git a/src/Cech_complex/doc/COPYRIGHT b/src/Cech_complex/doc/COPYRIGHT new file mode 100644 index 00000000..61f17f6d --- /dev/null +++ b/src/Cech_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/Cech_complex/doc/Intro_cech_complex.h b/src/Cech_complex/doc/Intro_cech_complex.h new file mode 100644 index 00000000..90086de7 --- /dev/null +++ b/src/Cech_complex/doc/Intro_cech_complex.h @@ -0,0 +1,102 @@ +/* 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 DOC_CECH_COMPLEX_INTRO_CECH_COMPLEX_H_ +#define DOC_CECH_COMPLEX_INTRO_CECH_COMPLEX_H_ + +namespace Gudhi { + +namespace cech_complex { + +/** \defgroup cech_complex Čech complex + * + * \author Vincent Rouvreau + * + * @{ + * + * \section cechdefinition Čech complex definition + * + * Čech complex + * <a target="_blank" href="https://en.wikipedia.org/wiki/%C4%8Cech_cohomology">(Wikipedia)</a> is a + * <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a> constructed + * from a proximity graph. The set of all simplices is filtered by the radius of their minimal enclosing ball. + * + * The input shall be a point cloud in an Euclidean space. + * + * \remark For people only interested in the topology of the \ref cech_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. + * + * \subsection cechalgorithm Algorithm + * + * Cech_complex first builds a proximity graph from a point cloud. + * The filtration value of each edge of the `Gudhi::Proximity_graph` is computed from + * `Gudhi::Minimal_enclosing_ball_radius` function. + * + * All edges that have a filtration value strictly greater than a user given maximal radius value, \f$max\_radius\f$, + * are not inserted into the complex. + * + * Vertex name correspond to the index of the point in the given range (aka. the point cloud). + * + * \image html "cech_one_skeleton.png" "Čech complex proximity graph representation" + * + * When creating a simplicial complex from this proximity graph, Cech_complex inserts the proximity graph into the + * simplicial complex data structure, and then expands the simplicial complex when required. + * + * On this example, as edges \f$(x,y)\f$, \f$(y,z)\f$ and \f$(z,y)\f$ are in the complex, the minimal ball radius + * containing the points \f$(x,y,z)\f$ is computed. + * + * \f$(x,y,z)\f$ is inserted to the simplicial complex with the filtration value set with + * \f$mini\_ball\_radius(x,y,z))\f$ iff \f$mini\_ball\_radius(x,y,z)) \leq max\_radius\f$. + * + * And so on for higher dimensions. + * + * \image html "cech_complex_representation.png" "Čech complex expansion" + * + * The minimal ball radius computation is insured by + * <a target="_blank" href="https://people.inf.ethz.ch/gaertner/subdir/software/miniball.html"> + * the miniball software (V3.0)</a> - Smallest Enclosing Balls of Points - and distributed with GUDHI. + * Please refer to + * <a target="_blank" href="https://people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf"> + * the miniball software design description</a> for more information about this computation. + * + * This radius computation is the reason why the Cech_complex is taking much more time to be computed than the + * \ref rips_complex but it offers more topological guarantees. + * + * If the Cech_complex interfaces are not detailed enough for your need, please refer to + * <a href="_cech_complex_2cech_complex_step_by_step_8cpp-example.html"> + * cech_complex_step_by_step.cpp</a> example, where the graph construction over the Simplex_tree is more detailed. + * + * \subsection cechpointscloudexample Example from a point cloud + * + * This example builds the proximity graph from the given points, and maximal radius values. + * Then it creates a `Simplex_tree` with it. + * + * Then, it is asked to display information about the simplicial complex. + * + * \include Cech_complex/cech_complex_example_from_points.cpp + * + * When launching (maximal enclosing ball radius is 1., is expanded until dimension 2): + * + * \code $> ./Cech_complex_example_from_points + * \endcode + * + * the program output is: + * + * \include Cech_complex/cech_complex_example_from_points_for_doc.txt + * + */ +/** @} */ // end defgroup cech_complex + +} // namespace cech_complex + +} // namespace Gudhi + +#endif // DOC_CECH_COMPLEX_INTRO_CECH_COMPLEX_H_ diff --git a/src/Cech_complex/doc/cech_complex_representation.ipe b/src/Cech_complex/doc/cech_complex_representation.ipe new file mode 100644 index 00000000..377745a3 --- /dev/null +++ b/src/Cech_complex/doc/cech_complex_representation.ipe @@ -0,0 +1,330 @@ +<?xml version="1.0"?> +<!DOCTYPE ipe SYSTEM "ipe.dtd"> +<ipe version="70107" creator="Ipe 7.1.10"> +<info created="D:20150603143945" modified="D:20180305162524"/> +<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" stroke="black" fill="darkcyan"> +48 640 m +80 672 l +48 672 l +h +</path> +<text matrix="1 0 0 1 -222.178 174.178" transformations="translations" pos="380 530" stroke="seagreen" type="label" width="70.886" height="8.307" depth="2.32" valign="baseline" size="large">Cech complex</text> +<text matrix="1 0 0 1 -212.333 10.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 -314.178 58.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 -194.178 -13.8225" 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 -226.178 18.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 -218.178 -21.8225" 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 -89.478 -87.9762" 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 -302.178 -13.8225" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">6</text> +<use name="mark/circle(sx)" pos="80 544" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="48 576" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="112 576" size="normal" stroke="black"/> +<use name="mark/fdisk(sfx)" pos="48 672" size="normal" stroke="black" fill="white"/> +<use name="mark/circle(sx)" pos="48 640" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="48 672" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="80 672" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="144 672" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="144 608" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="200 640" size="normal" stroke="black"/> +<use matrix="1 0 0 1 -100 -96" name="mark/circle(sx)" pos="304 672" size="normal" stroke="darkgray"/> +<use matrix="1 0 0 1 -100 -96" name="mark/circle(sx)" pos="336 672" size="normal" stroke="darkgray"/> +<path matrix="1 0 0 1 -100 -96" stroke="darkgray"> +32 0 0 32 304 672 e +</path> +<path matrix="1 0 0 1 -100 -96" stroke="darkgray" pen="fat"> +304 672 m +336 672 l +</path> +<text matrix="1 0 0 1 -214.178 50.178" transformations="translations" pos="380 530" stroke="darkgray" type="label" width="80.052" height="8.302" depth="0" valign="baseline" size="large">Maximal radius</text> +<text matrix="1 0 0 1 -226.178 -13.8225" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">7</text> +<text matrix="1 0 0 1 -258.178 30.1775" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">8</text> +<text matrix="1 0 0 1 -334.178 -13.8225" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">9</text> +<path stroke="black"> +112 576 m +144 608 l +</path> +<path stroke="black"> +144 672 m +144 608 l +200 640 l +h +</path> +<path stroke="black" fill="darkcyan"> +48 576 m +112 576 l +80 544 l +h +</path> +<use name="mark/fdisk(sfx)" pos="112 728" size="normal" stroke="black"/> +<path stroke="black"> +80 672 m +144 672 l +112 728 l +h +</path> +<use name="mark/fdisk(sfx)" pos="112 728" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="80 672" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="144 672" size="normal" stroke="black" fill="white"/> +<path stroke="black" fill="darkcyan"> +48 576 m +48 640 l +32 608 l +h +</path> +<use name="mark/fdisk(sfx)" pos="200 640" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="144 608" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="112 576" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="80 544" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="48 576" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="48 640" size="normal" stroke="black" fill="white"/> +<path stroke="darkcyan"> +32 0 0 32 80 576 e +</path> +<path stroke="darkcyan"> +22.6274 0 0 22.6274 64 656 e +</path> +<path stroke="darkorange"> +37.1429 0 0 37.1429 112 690.857 e +</path> +<path stroke="darkorange"> +37.1429 0 0 37.1429 162.857 640 e +</path> +<use name="mark/fdisk(sfx)" pos="32 608" size="normal" stroke="black"/> +<text matrix="1 0 0 1 -334.178 94.1775" transformations="translations" pos="352.708 510.349" stroke="black" type="label" width="9.963" height="6.42" depth="0" valign="baseline">10</text> +<path stroke="darkcyan"> +32 0 0 32 48 608 e +</path> +<use name="mark/fdisk(sfx)" pos="204 576" size="normal" stroke="darkgray" fill="white"/> +<use name="mark/fdisk(sfx)" pos="236 576" size="normal" stroke="darkgray" fill="white"/> +</page> +</ipe> diff --git a/src/Cech_complex/doc/cech_complex_representation.png b/src/Cech_complex/doc/cech_complex_representation.png Binary files differnew file mode 100644 index 00000000..d0eb85a5 --- /dev/null +++ b/src/Cech_complex/doc/cech_complex_representation.png diff --git a/src/Cech_complex/doc/cech_one_skeleton.ipe b/src/Cech_complex/doc/cech_one_skeleton.ipe new file mode 100644 index 00000000..ed66e132 --- /dev/null +++ b/src/Cech_complex/doc/cech_one_skeleton.ipe @@ -0,0 +1,314 @@ +<?xml version="1.0"?> +<!DOCTYPE ipe SYSTEM "ipe.dtd"> +<ipe version="70107" creator="Ipe 7.1.10"> +<info created="D:20150603143945" modified="D:20180305162558"/> +<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"/> +<text layer="alpha" matrix="1 0 0 1 -222.178 174.178" transformations="translations" pos="380 530" stroke="seagreen" type="label" width="84.053" height="8.307" depth="2.32" valign="baseline" size="large">Proximity graph</text> +<text matrix="1 0 0 1 -212.333 10.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 -314.178 58.1775" transformations="translations" pos="352.708 510.349" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">1</text> +<path matrix="1 0 0 1 -100 -96" stroke="darkgray" pen="fat"> +304 672 m +336 672 l +</path> +<text matrix="1 0 0 1 -194.178 -13.8225" 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 -226.178 18.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 -218.178 -21.8225" 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 -89.478 -87.9762" 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 -302.178 -13.8225" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">6</text> +<use name="mark/circle(sx)" pos="80 544" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="48 576" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="112 576" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="48 640" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="48 672" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="80 672" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="144 672" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="144 608" size="normal" stroke="black"/> +<use name="mark/circle(sx)" pos="200 640" size="normal" stroke="black"/> +<use matrix="1 0 0 1 -100 -96" name="mark/circle(sx)" pos="336 672" size="normal" stroke="darkgray"/> +<path matrix="1 0 0 1 -100 -96" stroke="darkgray"> +32 0 0 32 304 672 e +</path> +<text matrix="1 0 0 1 -214.178 50.178" transformations="translations" pos="380 530" stroke="darkgray" type="label" width="80.052" height="8.302" depth="0" valign="baseline" size="large">Maximal radius</text> +<text matrix="1 0 0 1 -226.178 -13.8225" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">7</text> +<text matrix="1 0 0 1 -258.178 30.1775" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">8</text> +<text matrix="1 0 0 1 -334.178 -13.8225" transformations="translations" pos="375.332 689.453" stroke="black" type="label" width="4.981" height="6.42" depth="0" valign="baseline">9</text> +<path stroke="black"> +112 576 m +144 608 l +</path> +<path stroke="black"> +144 672 m +144 608 l +200 640 l +h +</path> +<path stroke="black"> +48 640 m +80 672 l +48 672 l +h +</path> +<path stroke="black"> +48 576 m +112 576 l +80 544 l +h +</path> +<use name="mark/fdisk(sfx)" pos="112 728" size="normal" stroke="black"/> +<path stroke="black"> +80 672 m +144 672 l +112 728 l +h +</path> +<use name="mark/fdisk(sfx)" pos="112 728" size="normal" stroke="black" fill="white"/> +<path stroke="black"> +48 576 m +48 640 l +32 608 l +h +</path> +<use name="mark/fdisk(sfx)" pos="80 672" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="144 672" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="200 640" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="144 608" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="112 576" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="80 544" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="48 576" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="48 640" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="48 672" size="normal" stroke="black" fill="white"/> +<use name="mark/fdisk(sfx)" pos="32 608" size="normal" stroke="black" fill="white"/> +<text matrix="1 0 0 1 -334.178 94.1775" transformations="translations" pos="352.708 510.349" stroke="black" type="label" width="9.963" height="6.42" depth="0" valign="baseline">10</text> +<use name="mark/fdisk(sfx)" pos="204 576" size="normal" stroke="darkgray" fill="white"/> +<use name="mark/fdisk(sfx)" pos="236 576" size="normal" stroke="darkgray" fill="white"/> +</page> +</ipe> diff --git a/src/Cech_complex/doc/cech_one_skeleton.png b/src/Cech_complex/doc/cech_one_skeleton.png Binary files differnew file mode 100644 index 00000000..cc636616 --- /dev/null +++ b/src/Cech_complex/doc/cech_one_skeleton.png diff --git a/src/Cech_complex/example/CMakeLists.txt b/src/Cech_complex/example/CMakeLists.txt new file mode 100644 index 00000000..ab391215 --- /dev/null +++ b/src/Cech_complex/example/CMakeLists.txt @@ -0,0 +1,16 @@ +cmake_minimum_required(VERSION 2.6) +project(Cech_complex_examples) + +add_executable ( Cech_complex_example_step_by_step cech_complex_step_by_step.cpp ) +target_link_libraries(Cech_complex_example_step_by_step ${Boost_PROGRAM_OPTIONS_LIBRARY}) +if (TBB_FOUND) + target_link_libraries(Cech_complex_example_step_by_step ${TBB_LIBRARIES}) +endif() +add_test(NAME Cech_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:Cech_complex_example_step_by_step> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-d" "3") + +add_executable ( Cech_complex_example_from_points cech_complex_example_from_points.cpp) +if (TBB_FOUND) + target_link_libraries(Cech_complex_example_from_points ${TBB_LIBRARIES}) +endif() +add_test(NAME Cech_complex_example_from_points COMMAND $<TARGET_FILE:Cech_complex_example_from_points>) diff --git a/src/Cech_complex/example/cech_complex_example_from_points.cpp b/src/Cech_complex/example/cech_complex_example_from_points.cpp new file mode 100644 index 00000000..3cc5a4df --- /dev/null +++ b/src/Cech_complex/example/cech_complex_example_from_points.cpp @@ -0,0 +1,54 @@ +#include <gudhi/Cech_complex.h> +#include <gudhi/Simplex_tree.h> + +#include <iostream> +#include <string> +#include <vector> +#include <array> + +int main() { + // Type definitions + using Point_cloud = std::vector<std::array<double, 2>>; + using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; + using Filtration_value = Simplex_tree::Filtration_value; + using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>; + + Point_cloud points; + points.push_back({1., 0.}); // 0 + points.push_back({0., 1.}); // 1 + points.push_back({2., 1.}); // 2 + points.push_back({3., 2.}); // 3 + points.push_back({0., 3.}); // 4 + points.push_back({3. + std::sqrt(3.), 3.}); // 5 + points.push_back({1., 4.}); // 6 + points.push_back({3., 4.}); // 7 + points.push_back({2., 4. + std::sqrt(3.)}); // 8 + points.push_back({0., 4.}); // 9 + points.push_back({-0.5, 2.}); // 10 + + // ---------------------------------------------------------------------------- + // Init of a Cech complex from points + // ---------------------------------------------------------------------------- + Filtration_value max_radius = 1.; + Cech_complex cech_complex_from_points(points, max_radius); + + Simplex_tree stree; + cech_complex_from_points.create_complex(stree, 2); + // ---------------------------------------------------------------------------- + // Display information about the one skeleton Cech complex + // ---------------------------------------------------------------------------- + std::cout << "Cech complex is of dimension " << stree.dimension() << " - " << stree.num_simplices() << " simplices - " + << stree.num_vertices() << " vertices." << std::endl; + + std::cout << "Iterator on Cech complex simplices in the filtration order, with [filtration value]:" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + std::cout << " ( "; + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " + << "[" << stree.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + return 0; +} diff --git a/src/Cech_complex/example/cech_complex_example_from_points_for_doc.txt b/src/Cech_complex/example/cech_complex_example_from_points_for_doc.txt new file mode 100644 index 00000000..be0afc76 --- /dev/null +++ b/src/Cech_complex/example/cech_complex_example_from_points_for_doc.txt @@ -0,0 +1,31 @@ +Iterator on Cech complex simplices in the filtration order, with [filtration value]: + ( 0 ) -> [0] + ( 1 ) -> [0] + ( 2 ) -> [0] + ( 3 ) -> [0] + ( 4 ) -> [0] + ( 5 ) -> [0] + ( 6 ) -> [0] + ( 7 ) -> [0] + ( 8 ) -> [0] + ( 9 ) -> [0] + ( 10 ) -> [0] + ( 9 4 ) -> [0.5] + ( 9 6 ) -> [0.5] + ( 10 1 ) -> [0.559017] + ( 10 4 ) -> [0.559017] + ( 1 0 ) -> [0.707107] + ( 2 0 ) -> [0.707107] + ( 3 2 ) -> [0.707107] + ( 6 4 ) -> [0.707107] + ( 9 6 4 ) -> [0.707107] + ( 2 1 ) -> [1] + ( 2 1 0 ) -> [1] + ( 4 1 ) -> [1] + ( 5 3 ) -> [1] + ( 7 3 ) -> [1] + ( 7 5 ) -> [1] + ( 7 6 ) -> [1] + ( 8 6 ) -> [1] + ( 8 7 ) -> [1] + ( 10 4 1 ) -> [1] diff --git a/src/Cech_complex/example/cech_complex_step_by_step.cpp b/src/Cech_complex/example/cech_complex_step_by_step.cpp new file mode 100644 index 00000000..b3d05697 --- /dev/null +++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp @@ -0,0 +1,154 @@ +/* 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 + */ + +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Points_off_io.h> + +#include <gudhi/Miniball.hpp> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity +#include <utility> // for pair +#include <map> + +// ---------------------------------------------------------------------------- +// rips_persistence_step_by_step is an example of each step that is required to +// build a Rips over a Simplex_tree. Please refer to rips_persistence to see +// how to do the same thing with the Rips_complex wrapper for less detailed +// steps. +// ---------------------------------------------------------------------------- + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<>; +using Simplex_handle = Simplex_tree::Simplex_handle; +using Filtration_value = Simplex_tree::Filtration_value; +using Point = std::vector<double>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; +using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>; + +class Cech_blocker { + private: + using Point_cloud = std::vector<Point>; + using Point_iterator = Point_cloud::const_iterator; + using Coordinate_iterator = Point::const_iterator; + using Min_sphere = Gudhi::Miniball::Miniball<Gudhi::Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>; + + public: + bool operator()(Simplex_handle sh) { + std::vector<Point> points; + for (auto vertex : simplex_tree_.simplex_vertex_range(sh)) { + points.push_back(point_cloud_[vertex]); +#ifdef DEBUG_TRACES + std::cout << "#(" << vertex << ")#"; +#endif // DEBUG_TRACES + } + Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); +#ifdef DEBUG_TRACES + std::cout << "radius = " << radius << " - " << (radius > max_radius_) << std::endl; +#endif // DEBUG_TRACES + simplex_tree_.assign_filtration(sh, radius); + return (radius > max_radius_); + } + Cech_blocker(Simplex_tree& simplex_tree, Filtration_value max_radius, const std::vector<Point>& point_cloud) + : simplex_tree_(simplex_tree), max_radius_(max_radius), point_cloud_(point_cloud) { + dimension_ = point_cloud_[0].size(); + } + + private: + Simplex_tree simplex_tree_; + Filtration_value max_radius_; + std::vector<Point> point_cloud_; + int dimension_; +}; + +void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& max_radius, int& dim_max); + +int main(int argc, char* argv[]) { + std::string off_file_points; + Filtration_value max_radius; + int dim_max; + + program_options(argc, argv, off_file_points, max_radius, dim_max); + + // Extract the points from the file filepoints + Points_off_reader off_reader(off_file_points); + + // Compute the proximity graph of the points + Proximity_graph prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(), max_radius, + Gudhi::Minimal_enclosing_ball_radius()); + + // Construct the Rips complex in a Simplex Tree + Simplex_tree st; + // insert the proximity graph in the simplex tree + st.insert_graph(prox_graph); + // expand the graph until dimension dim_max + st.expansion_with_blockers(dim_max, Cech_blocker(st, max_radius, off_reader.get_point_cloud())); + + std::cout << "The complex contains " << st.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << st.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + st.initialize_filtration(); + +#if DEBUG_TRACES + std::cout << "********************************************************************\n"; + std::cout << "* The complex contains " << st.num_simplices() << " simplices - dimension=" << st.dimension() << "\n"; + std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n"; + for (auto f_simplex : st.filtration_simplex_range()) { + std::cout << " " + << "[" << st.filtration(f_simplex) << "] "; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::cout << static_cast<int>(vertex) << " "; + } + std::cout << std::endl; + } +#endif // DEBUG_TRACES + + return 0; +} + +void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& max_radius, int& dim_max) { + 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 an OFF file containing a point set.\n"); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "max-radius,r", + po::value<Filtration_value>(&max_radius)->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal length of an edge for the Rips complex construction.")( + "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Rips complex we want to compute."); + + 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 << "Construct a Cech complex defined on a set of input points.\n \n"; + + std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl; + std::cout << visible << std::endl; + exit(-1); + } +} diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h new file mode 100644 index 00000000..b0871e10 --- /dev/null +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -0,0 +1,118 @@ +/* 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 CECH_COMPLEX_H_ +#define CECH_COMPLEX_H_ + +#include <gudhi/distance_functions.h> // for Gudhi::Minimal_enclosing_ball_radius +#include <gudhi/graph_simplicial_complex.h> // for Gudhi::Proximity_graph +#include <gudhi/Debug_utils.h> // for GUDHI_CHECK +#include <gudhi/Cech_complex_blocker.h> // for Gudhi::cech_complex::Cech_blocker + +#include <iostream> +#include <stdexcept> // for exception management +#include <vector> + +namespace Gudhi { + +namespace cech_complex { + +/** + * \class Cech_complex + * \brief Cech complex data structure. + * + * \ingroup cech_complex + * + * \details + * The data structure is a proximity graph, containing edges when the edge length is less or equal + * to a given max_radius. Edge length is computed from `Gudhi::Minimal_enclosing_ball_radius` distance function. + * + * \tparam SimplicialComplexForProximityGraph furnishes `Vertex_handle` and `Filtration_value` type definition required + * by `Gudhi::Proximity_graph`. + * + * \tparam ForwardPointRange must be a range for which `std::begin()` and `std::end()` methods return input + * iterators on a point. `std::begin()` and `std::end()` methods are also required for a point. + */ +template <typename SimplicialComplexForProximityGraph, typename ForwardPointRange> +class Cech_complex { + private: + // Required by compute_proximity_graph + using Vertex_handle = typename SimplicialComplexForProximityGraph::Vertex_handle; + using Filtration_value = typename SimplicialComplexForProximityGraph::Filtration_value; + using Proximity_graph = Gudhi::Proximity_graph<SimplicialComplexForProximityGraph>; + + // Retrieve Coordinate type from ForwardPointRange + using Point_from_range_iterator = typename boost::range_const_iterator<ForwardPointRange>::type; + using Point_from_range = typename std::iterator_traits<Point_from_range_iterator>::value_type; + using Coordinate_iterator = typename boost::range_const_iterator<Point_from_range>::type; + using Coordinate = typename std::iterator_traits<Coordinate_iterator>::value_type; + + public: + // Point and Point_cloud type definition + using Point = std::vector<Coordinate>; + using Point_cloud = std::vector<Point>; + + public: + /** \brief Cech_complex constructor from a list of points. + * + * @param[in] points Range of points. + * @param[in] max_radius Maximal radius value. + * + * \tparam ForwardPointRange must be a range of Point. Point must be a range of <b>copyable</b> Cartesian coordinates. + * + */ + Cech_complex(const ForwardPointRange& points, Filtration_value max_radius) : max_radius_(max_radius) { + // Point cloud deep copy + point_cloud_.reserve(boost::size(points)); + for (auto&& point : points) point_cloud_.emplace_back(std::begin(point), std::end(point)); + + cech_skeleton_graph_ = Gudhi::compute_proximity_graph<SimplicialComplexForProximityGraph>( + point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); + } + + /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal + * dimension, using the Cech blocker oracle. + * + * @param[in] complex SimplicialComplexForCech to be created. + * @param[in] dim_max graph expansion until this given maximal dimension. + * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0. + * + */ + template <typename SimplicialComplexForCechComplex> + void create_complex(SimplicialComplexForCechComplex& complex, int dim_max) { + GUDHI_CHECK(complex.num_vertices() == 0, + std::invalid_argument("Cech_complex::create_complex - simplicial complex is not empty")); + + // insert the proximity graph in the simplicial complex + complex.insert_graph(cech_skeleton_graph_); + // expand the graph until dimension dim_max + complex.expansion_with_blockers(dim_max, + Cech_blocker<SimplicialComplexForCechComplex, Cech_complex>(&complex, this)); + } + + /** @return max_radius value given at construction. */ + Filtration_value max_radius() const { return max_radius_; } + + /** @param[in] vertex Point position in the range. + * @return The point. + */ + const Point& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; } + + private: + Proximity_graph cech_skeleton_graph_; + Filtration_value max_radius_; + Point_cloud point_cloud_; +}; + +} // namespace cech_complex + +} // namespace Gudhi + +#endif // CECH_COMPLEX_H_ diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h new file mode 100644 index 00000000..068cdde3 --- /dev/null +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -0,0 +1,79 @@ +/* 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 CECH_COMPLEX_BLOCKER_H_ +#define CECH_COMPLEX_BLOCKER_H_ + +#include <gudhi/distance_functions.h> // for Gudhi::Minimal_enclosing_ball_radius + +#include <iostream> +#include <vector> +#include <cmath> // for std::sqrt + +namespace Gudhi { + +namespace cech_complex { + +/** \internal + * \class Cech_blocker + * \brief Čech complex blocker. + * + * \ingroup cech_complex + * + * \details + * Čech blocker is an oracle constructed from a Cech_complex and a simplicial complex. + * + * \tparam SimplicialComplexForProximityGraph furnishes `Simplex_handle` and `Filtration_value` type definition, + * `simplex_vertex_range(Simplex_handle sh)`and `assign_filtration(Simplex_handle sh, Filtration_value filt)` methods. + * + * \tparam Chech_complex is required by the blocker. + */ +template <typename SimplicialComplexForCech, typename Cech_complex> +class Cech_blocker { + private: + using Point_cloud = typename Cech_complex::Point_cloud; + + using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; + using Filtration_value = typename SimplicialComplexForCech::Filtration_value; + + public: + /** \internal \brief Čech complex blocker operator() - the oracle - assigns the filtration value from the simplex + * radius and returns if the simplex expansion must be blocked. + * \param[in] sh The Simplex_handle. + * \return true if the simplex radius is greater than the Cech_complex max_radius*/ + bool operator()(Simplex_handle sh) { + Point_cloud points; + for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { + points.push_back(cc_ptr_->get_point(vertex)); +#ifdef DEBUG_TRACES + std::cout << "#(" << vertex << ")#"; +#endif // DEBUG_TRACES + } + Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); +#ifdef DEBUG_TRACES + if (radius > cc_ptr_->max_radius()) std::cout << "radius > max_radius => expansion is blocked\n"; +#endif // DEBUG_TRACES + sc_ptr_->assign_filtration(sh, radius); + return (radius > cc_ptr_->max_radius()); + } + + /** \internal \brief Čech complex blocker constructor. */ + Cech_blocker(SimplicialComplexForCech* sc_ptr, Cech_complex* cc_ptr) : sc_ptr_(sc_ptr), cc_ptr_(cc_ptr) {} + + private: + SimplicialComplexForCech* sc_ptr_; + Cech_complex* cc_ptr_; +}; + +} // namespace cech_complex + +} // namespace Gudhi + +#endif // CECH_COMPLEX_BLOCKER_H_ diff --git a/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT b/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT new file mode 100644 index 00000000..dbe4c553 --- /dev/null +++ b/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT @@ -0,0 +1,4 @@ +The miniball software is available under the GNU General Public License (GPLv3 - https://www.gnu.org/copyleft/gpl.html). +If your intended use is not compliant with this license, please buy a commercial license (EUR 500 - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/license.html). +You need a license if the software that you develop using Miniball V3.0 is not open source. + diff --git a/src/Cech_complex/include/gudhi/Miniball.README b/src/Cech_complex/include/gudhi/Miniball.README new file mode 100644 index 00000000..033d8953 --- /dev/null +++ b/src/Cech_complex/include/gudhi/Miniball.README @@ -0,0 +1,26 @@ +https://people.inf.ethz.ch/gaertner/subdir/software/miniball.html + +Smallest Enclosing Balls of Points - Fast and Robust in C++. +(high-quality software for smallest enclosing balls of balls is available in the computational geometry algorithms library CGAL) + + +This is the miniball software (V3.0) for computing smallest enclosing balls of points in arbitrary dimensions. It consists of a C++ header file Miniball.hpp (around 500 lines of code) and two example programs miniball_example.cpp and miniball_example_containers.cpp that demonstrate the usage. The first example stores the coordinates of the input points in a two-dimensional array, the second example uses a list of vectors to show how generic containers can be used. + +Credits: Aditya Gupta and Alexandros Konstantinakis-Karmis have significantly contributed to this version of the software. + +Changes - https://people.inf.ethz.ch/gaertner/subdir/software/miniball/changes.txt - from previous versions. + +The theory - https://people.inf.ethz.ch/gaertner/subdir/texts/own_work/esa99_final.pdf - behind the miniball software (Proc. 7th Annual European Symposium on Algorithms (ESA), Lecture Notes in Computer Science 1643, Springer-Verlag, pp.325-338, 1999). + +Main Features: + + Very fast in low dimensions. 1 million points in 5-space are processed within 0.05 seconds on any recent machine. + + High numerical stability. Almost all input degeneracies (cospherical points, multiple points, points very close together) are routinely handled. + + Easily integrates into your code. You can freely choose the coordinate type of your points and the container to store the points. If you still need to adapt the code, the header is small and readable and contains documentation for all major methods. + + +Changes done for the GUDHI version of MiniBall: + - Add include guard + - Move Miniball namespace inside a new Gudhi namespace diff --git a/src/Cech_complex/include/gudhi/Miniball.hpp b/src/Cech_complex/include/gudhi/Miniball.hpp new file mode 100644 index 00000000..ce6cbb5b --- /dev/null +++ b/src/Cech_complex/include/gudhi/Miniball.hpp @@ -0,0 +1,523 @@ +// Copright (C) 1999-2013, Bernd Gaertner +// $Rev: 3581 $ +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. + +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. + +// You should have received a copy of the GNU General Public License +// along with this program. If not, see <http://www.gnu.org/licenses/>. +// +// Contact: +// -------- +// Bernd Gaertner +// Institute of Theoretical Computer Science +// ETH Zuerich +// CAB G31.1 +// CH-8092 Zuerich, Switzerland +// http://www.inf.ethz.ch/personal/gaertner + +#ifndef MINIBALL_HPP_ +#define MINIBALL_HPP_ + +#include <cassert> +#include <algorithm> +#include <list> +#include <ctime> +#include <limits> + +namespace Gudhi { + +namespace Miniball { + + // Global Functions + // ================ + template <typename NT> + inline NT mb_sqr (NT r) {return r*r;} + + // Functors + // ======== + + // functor to map a point iterator to the corresponding coordinate iterator; + // generic version for points whose coordinate containers have begin() + template < typename Pit_, typename Cit_ > + struct CoordAccessor { + typedef Pit_ Pit; + typedef Cit_ Cit; + inline Cit operator() (Pit it) const { return (*it).begin(); } + }; + + // partial specialization for points whose coordinate containers are arrays + template < typename Pit_, typename Cit_ > + struct CoordAccessor<Pit_, Cit_*> { + typedef Pit_ Pit; + typedef Cit_* Cit; + inline Cit operator() (Pit it) const { return *it; } + }; + + // Class Declaration + // ================= + + template <typename CoordAccessor> + class Miniball { + private: + // types + // The iterator type to go through the input points + typedef typename CoordAccessor::Pit Pit; + // The iterator type to go through the coordinates of a single point. + typedef typename CoordAccessor::Cit Cit; + // The coordinate type + typedef typename std::iterator_traits<Cit>::value_type NT; + // The iterator to go through the support points + typedef typename std::list<Pit>::iterator Sit; + + // data members... + const int d; // dimension + Pit points_begin; + Pit points_end; + CoordAccessor coord_accessor; + double time; + const NT nt0; // NT(0) + + //...for the algorithms + std::list<Pit> L; + Sit support_end; + int fsize; // number of forced points + int ssize; // number of support points + + // ...for the ball updates + NT* current_c; + NT current_sqr_r; + NT** c; + NT* sqr_r; + + // helper arrays + NT* q0; + NT* z; + NT* f; + NT** v; + NT** a; + + public: + // The iterator type to go through the support points + typedef typename std::list<Pit>::const_iterator SupportPointIterator; + + // PRE: [begin, end) is a nonempty range + // POST: computes the smallest enclosing ball of the points in the range + // [begin, end); the functor a maps a point iterator to an iterator + // through the d coordinates of the point + Miniball (int d_, Pit begin, Pit end, CoordAccessor ca = CoordAccessor()); + + // POST: returns a pointer to the first element of an array that holds + // the d coordinates of the center of the computed ball + const NT* center () const; + + // POST: returns the squared radius of the computed ball + NT squared_radius () const; + + // POST: returns the number of support points of the computed ball; + // the support points form a minimal set with the same smallest + // enclosing ball as the input set; in particular, the support + // points are on the boundary of the computed ball, and their + // number is at most d+1 + int nr_support_points () const; + + // POST: returns an iterator to the first support point + SupportPointIterator support_points_begin () const; + + // POST: returns a past-the-end iterator for the range of support points + SupportPointIterator support_points_end () const; + + // POST: returns the maximum excess of any input point w.r.t. the computed + // ball, divided by the squared radius of the computed ball. The + // excess of a point is the difference between its squared distance + // from the center and the squared radius; Ideally, the return value + // is 0. subopt is set to the absolute value of the most negative + // coefficient in the affine combination of the support points that + // yields the center. Ideally, this is a convex combination, and there + // is no negative coefficient in which case subopt is set to 0. + NT relative_error (NT& subopt) const; + + // POST: return true if the relative error is at most tol, and the + // suboptimality is 0; the default tolerance is 10 times the + // coordinate type's machine epsilon + bool is_valid (NT tol = NT(10) * std::numeric_limits<NT>::epsilon()) const; + + // POST: returns the time in seconds taken by the constructor call for + // computing the smallest enclosing ball + double get_time() const; + + // POST: deletes dynamically allocated arrays + ~Miniball(); + + private: + void mtf_mb (Sit n); + void mtf_move_to_front (Sit j); + void pivot_mb (Pit n); + void pivot_move_to_front (Pit j); + NT excess (Pit pit) const; + void pop (); + bool push (Pit pit); + NT suboptimality () const; + void create_arrays(); + void delete_arrays(); + }; + + // Class Definition + // ================ + template <typename CoordAccessor> + Miniball<CoordAccessor>::Miniball (int d_, Pit begin, Pit end, + CoordAccessor ca) + : d (d_), + points_begin (begin), + points_end (end), + coord_accessor (ca), + time (clock()), + nt0 (NT(0)), + L(), + support_end (L.begin()), + fsize(0), + ssize(0), + current_c (NULL), + current_sqr_r (NT(-1)), + c (NULL), + sqr_r (NULL), + q0 (NULL), + z (NULL), + f (NULL), + v (NULL), + a (NULL) + { + assert (points_begin != points_end); + create_arrays(); + + // set initial center + for (int j=0; j<d; ++j) c[0][j] = nt0; + current_c = c[0]; + + // compute miniball + pivot_mb (points_end); + + // update time + time = (clock() - time) / CLOCKS_PER_SEC; + } + + template <typename CoordAccessor> + Miniball<CoordAccessor>::~Miniball() + { + delete_arrays(); + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::create_arrays() + { + c = new NT*[d+1]; + v = new NT*[d+1]; + a = new NT*[d+1]; + for (int i=0; i<d+1; ++i) { + c[i] = new NT[d]; + v[i] = new NT[d]; + a[i] = new NT[d]; + } + sqr_r = new NT[d+1]; + q0 = new NT[d]; + z = new NT[d+1]; + f = new NT[d+1]; + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::delete_arrays() + { + delete[] f; + delete[] z; + delete[] q0; + delete[] sqr_r; + for (int i=0; i<d+1; ++i) { + delete[] a[i]; + delete[] v[i]; + delete[] c[i]; + } + delete[] a; + delete[] v; + delete[] c; + } + + template <typename CoordAccessor> + const typename Miniball<CoordAccessor>::NT* + Miniball<CoordAccessor>::center () const + { + return current_c; + } + + template <typename CoordAccessor> + typename Miniball<CoordAccessor>::NT + Miniball<CoordAccessor>::squared_radius () const + { + return current_sqr_r; + } + + template <typename CoordAccessor> + int Miniball<CoordAccessor>::nr_support_points () const + { + assert (ssize < d+2); + return ssize; + } + + template <typename CoordAccessor> + typename Miniball<CoordAccessor>::SupportPointIterator + Miniball<CoordAccessor>::support_points_begin () const + { + return L.begin(); + } + + template <typename CoordAccessor> + typename Miniball<CoordAccessor>::SupportPointIterator + Miniball<CoordAccessor>::support_points_end () const + { + return support_end; + } + + template <typename CoordAccessor> + typename Miniball<CoordAccessor>::NT + Miniball<CoordAccessor>::relative_error (NT& subopt) const + { + NT e, max_e = nt0; + // compute maximum absolute excess of support points + for (SupportPointIterator it = support_points_begin(); + it != support_points_end(); ++it) { + e = excess (*it); + if (e < nt0) e = -e; + if (e > max_e) { + max_e = e; + } + } + // compute maximum excess of any point + for (Pit i = points_begin; i != points_end; ++i) + if ((e = excess (i)) > max_e) + max_e = e; + + subopt = suboptimality(); + assert (current_sqr_r > nt0 || max_e == nt0); + return (current_sqr_r == nt0 ? nt0 : max_e / current_sqr_r); + } + + template <typename CoordAccessor> + bool Miniball<CoordAccessor>::is_valid (NT tol) const + { + NT suboptimality; + return ( (relative_error (suboptimality) <= tol) && (suboptimality == 0) ); + } + + template <typename CoordAccessor> + double Miniball<CoordAccessor>::get_time() const + { + return time; + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::mtf_mb (Sit n) + { + // Algorithm 1: mtf_mb (L_{n-1}, B), where L_{n-1} = [L.begin, n) + // B: the set of forced points, defining the current ball + // S: the superset of support points computed by the algorithm + // -------------------------------------------------------------- + // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, + // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf + + // PRE: B = S + assert (fsize == ssize); + + support_end = L.begin(); + if ((fsize) == d+1) return; + + // incremental construction + for (Sit i = L.begin(); i != n;) + { + // INV: (support_end - L.begin() == |S|-|B|) + assert (std::distance (L.begin(), support_end) == ssize - fsize); + + Sit j = i++; + if (excess(*j) > nt0) + if (push(*j)) { // B := B + p_i + mtf_mb (j); // mtf_mb (L_{i-1}, B + p_i) + pop(); // B := B - p_i + mtf_move_to_front(j); + } + } + // POST: the range [L.begin(), support_end) stores the set S\B + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::mtf_move_to_front (Sit j) + { + if (support_end == j) + support_end++; + L.splice (L.begin(), L, j); + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::pivot_mb (Pit n) + { + // Algorithm 2: pivot_mb (L_{n-1}), where L_{n-1} = [L.begin, n) + // -------------------------------------------------------------- + // from B. Gaertner, Fast and Robust Smallest Enclosing Balls, ESA 1999, + // http://www.inf.ethz.ch/personal/gaertner/texts/own_work/esa99_final.pdf + NT old_sqr_r; + const NT* c; + Pit pivot, k; + NT e, max_e, sqr_r; + Cit p; + do { + old_sqr_r = current_sqr_r; + sqr_r = current_sqr_r; + + pivot = points_begin; + max_e = nt0; + for (k = points_begin; k != n; ++k) { + p = coord_accessor(k); + e = -sqr_r; + c = current_c; + for (int j=0; j<d; ++j) + e += mb_sqr<NT>(*p++-*c++); + if (e > max_e) { + max_e = e; + pivot = k; + } + } + + if (max_e > nt0) { + // check if the pivot is already contained in the support set + if (std::find(L.begin(), support_end, pivot) == support_end) { + assert (fsize == 0); + if (push (pivot)) { + mtf_mb(support_end); + pop(); + pivot_move_to_front(pivot); + } + } + } + } while (old_sqr_r < current_sqr_r); + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::pivot_move_to_front (Pit j) + { + L.push_front(j); + if (std::distance(L.begin(), support_end) == d+2) + support_end--; + } + + template <typename CoordAccessor> + inline typename Miniball<CoordAccessor>::NT + Miniball<CoordAccessor>::excess (Pit pit) const + { + Cit p = coord_accessor(pit); + NT e = -current_sqr_r; + NT* c = current_c; + for (int k=0; k<d; ++k){ + e += mb_sqr<NT>(*p++-*c++); + } + return e; + } + + template <typename CoordAccessor> + void Miniball<CoordAccessor>::pop () + { + --fsize; + } + + template <typename CoordAccessor> + bool Miniball<CoordAccessor>::push (Pit pit) + { + int i, j; + NT eps = mb_sqr<NT>(std::numeric_limits<NT>::epsilon()); + + Cit cit = coord_accessor(pit); + Cit p = cit; + + if (fsize==0) { + for (i=0; i<d; ++i) + q0[i] = *p++; + for (i=0; i<d; ++i) + c[0][i] = q0[i]; + sqr_r[0] = nt0; + } + else { + // set v_fsize to Q_fsize + for (i=0; i<d; ++i) + //v[fsize][i] = p[i]-q0[i]; + v[fsize][i] = *p++-q0[i]; + + // compute the a_{fsize,i}, i< fsize + for (i=1; i<fsize; ++i) { + a[fsize][i] = nt0; + for (j=0; j<d; ++j) + a[fsize][i] += v[i][j] * v[fsize][j]; + a[fsize][i]*=(2/z[i]); + } + + // update v_fsize to Q_fsize-\bar{Q}_fsize + for (i=1; i<fsize; ++i) { + for (j=0; j<d; ++j) + v[fsize][j] -= a[fsize][i]*v[i][j]; + } + + // compute z_fsize + z[fsize]=nt0; + for (j=0; j<d; ++j) + z[fsize] += mb_sqr<NT>(v[fsize][j]); + z[fsize]*=2; + + // reject push if z_fsize too small + if (z[fsize]<eps*current_sqr_r) { + return false; + } + + // update c, sqr_r + p=cit; + NT e = -sqr_r[fsize-1]; + for (i=0; i<d; ++i) + e += mb_sqr<NT>(*p++-c[fsize-1][i]); + f[fsize]=e/z[fsize]; + + for (i=0; i<d; ++i) + c[fsize][i] = c[fsize-1][i]+f[fsize]*v[fsize][i]; + sqr_r[fsize] = sqr_r[fsize-1] + e*f[fsize]/2; + } + current_c = c[fsize]; + current_sqr_r = sqr_r[fsize]; + ssize = ++fsize; + return true; + } + + template <typename CoordAccessor> + typename Miniball<CoordAccessor>::NT + Miniball<CoordAccessor>::suboptimality () const + { + NT* l = new NT[d+1]; + NT min_l = nt0; + l[0] = NT(1); + for (int i=ssize-1; i>0; --i) { + l[i] = f[i]; + for (int k=ssize-1; k>i; --k) + l[i]-=a[k][i]*l[k]; + if (l[i] < min_l) min_l = l[i]; + l[0] -= l[i]; + } + if (l[0] < min_l) min_l = l[0]; + delete[] l; + if (min_l < nt0) + return -min_l; + return nt0; + } +} // namespace Miniball + +} // namespace Gudhi + +#endif // MINIBALL_HPP_ diff --git a/src/Cech_complex/test/CMakeLists.txt b/src/Cech_complex/test/CMakeLists.txt new file mode 100644 index 00000000..8db51173 --- /dev/null +++ b/src/Cech_complex/test/CMakeLists.txt @@ -0,0 +1,15 @@ +cmake_minimum_required(VERSION 2.6) +project(Cech_complex_tests) + +include(GUDHI_test_coverage) + +add_executable ( Cech_complex_test_unit test_cech_complex.cpp ) +target_link_libraries(Cech_complex_test_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +if (TBB_FOUND) + target_link_libraries(Cech_complex_test_unit ${TBB_LIBRARIES}) +endif() + +# 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}/) + +gudhi_add_coverage_test(Cech_complex_test_unit) diff --git a/src/Cech_complex/test/README b/src/Cech_complex/test/README new file mode 100644 index 00000000..adf704f7 --- /dev/null +++ b/src/Cech_complex/test/README @@ -0,0 +1,12 @@ +To compile: +*********** + +cmake . +make + +To launch with details: +*********************** + +./Cech_complex_test_unit --report_level=detailed --log_level=all + + ==> echo $? returns 0 in case of success (non-zero otherwise) diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp new file mode 100644 index 00000000..c6b15d7f --- /dev/null +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -0,0 +1,252 @@ +/* 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 + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "cech_complex" +#include <boost/test/unit_test.hpp> + +#include <cmath> // float comparison +#include <limits> +#include <string> +#include <vector> +#include <algorithm> // std::max + +#include <gudhi/Cech_complex.h> +// to construct Cech_complex from a OFF file of points +#include <gudhi/Points_off_io.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Unitary_tests_utils.h> +#include <gudhi/Miniball.hpp> + +// Type definitions +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Point = std::vector<Filtration_value>; +using Point_cloud = std::vector<Point>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; +using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>; + +using Point_iterator = Point_cloud::const_iterator; +using Coordinate_iterator = Point::const_iterator; +using Min_sphere = Gudhi::Miniball::Miniball<Gudhi::Miniball::CoordAccessor<Point_iterator, Coordinate_iterator>>; + +BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { + // ---------------------------------------------------------------------------- + // + // Init of a Cech complex from a point cloud + // + // ---------------------------------------------------------------------------- + Point_cloud points; + points.push_back({1., 0.}); // 0 + points.push_back({0., 1.}); // 1 + points.push_back({2., 1.}); // 2 + points.push_back({3., 2.}); // 3 + points.push_back({0., 3.}); // 4 + points.push_back({3. + std::sqrt(3.), 3.}); // 5 + points.push_back({1., 4.}); // 6 + points.push_back({3., 4.}); // 7 + points.push_back({2., 4. + std::sqrt(3.)}); // 8 + points.push_back({0., 4.}); // 9 + points.push_back({-0.5, 2.}); // 10 + + Filtration_value max_radius = 1.0; + std::cout << "========== NUMBER OF POINTS = " << points.size() << " - Cech max_radius = " << max_radius + << "==========" << std::endl; + + Cech_complex cech_complex_for_doc(points, max_radius); + + GUDHI_TEST_FLOAT_EQUALITY_CHECK(cech_complex_for_doc.max_radius(), max_radius); + std::size_t i = 0; + for (; i < points.size(); i++) { + BOOST_CHECK(points[i] == cech_complex_for_doc.get_point(i)); + } + + const int DIMENSION_1 = 1; + Simplex_tree st; + cech_complex_for_doc.create_complex(st, DIMENSION_1); + std::cout << "st.dimension()=" << st.dimension() << std::endl; + BOOST_CHECK(st.dimension() == DIMENSION_1); + + const int NUMBER_OF_VERTICES = 11; + std::cout << "st.num_vertices()=" << st.num_vertices() << std::endl; + BOOST_CHECK(st.num_vertices() == NUMBER_OF_VERTICES); + + std::cout << "st.num_simplices()=" << st.num_simplices() << std::endl; + BOOST_CHECK(st.num_simplices() == 27); + + // Check filtration values of vertices is 0.0 + for (auto f_simplex : st.skeleton_simplex_range(0)) { + BOOST_CHECK(st.filtration(f_simplex) == 0.0); + } + + // Check filtration values of edges + for (auto f_simplex : st.skeleton_simplex_range(DIMENSION_1)) { + if (DIMENSION_1 == st.dimension(f_simplex)) { + std::vector<Point> vp; + std::cout << "vertex = ("; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::cout << vertex << ","; + vp.push_back(points.at(vertex)); + } + std::cout << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1)) + << " - filtration =" << st.filtration(f_simplex) << std::endl; + BOOST_CHECK(vp.size() == 2); + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), + Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1))); + } + } + + const int DIMENSION_2 = 2; + +#ifdef GUDHI_DEBUG + BOOST_CHECK_THROW(cech_complex_for_doc.create_complex(st, DIMENSION_2), std::invalid_argument); +#endif + + Simplex_tree st2; + cech_complex_for_doc.create_complex(st2, DIMENSION_2); + std::cout << "st2.dimension()=" << st2.dimension() << std::endl; + BOOST_CHECK(st2.dimension() == DIMENSION_2); + + std::cout << "st2.num_vertices()=" << st2.num_vertices() << std::endl; + BOOST_CHECK(st2.num_vertices() == NUMBER_OF_VERTICES); + + std::cout << "st2.num_simplices()=" << st2.num_simplices() << std::endl; + BOOST_CHECK(st2.num_simplices() == 30); + + Point_cloud points012; + for (std::size_t vertex = 0; vertex <= 2; vertex++) { + points012.push_back(cech_complex_for_doc.get_point(vertex)); + } + std::size_t dimension = points[0].end() - points[0].begin(); + Min_sphere ms012(dimension, points012.begin(), points012.end()); + + Simplex_tree::Filtration_value f012 = st2.filtration(st2.find({0, 1, 2})); + std::cout << "f012= " << f012 << " | ms012_radius= " << std::sqrt(ms012.squared_radius()) << std::endl; + + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(ms012.squared_radius())); + + Point_cloud points1410; + points1410.push_back(cech_complex_for_doc.get_point(1)); + points1410.push_back(cech_complex_for_doc.get_point(4)); + points1410.push_back(cech_complex_for_doc.get_point(10)); + Min_sphere ms1410(dimension, points1410.begin(), points1410.end()); + + Simplex_tree::Filtration_value f1410 = st2.filtration(st2.find({1, 4, 10})); + std::cout << "f1410= " << f1410 << " | ms1410_radius= " << std::sqrt(ms1410.squared_radius()) << std::endl; + + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, std::sqrt(ms1410.squared_radius())); + + Point_cloud points469; + points469.push_back(cech_complex_for_doc.get_point(4)); + points469.push_back(cech_complex_for_doc.get_point(6)); + points469.push_back(cech_complex_for_doc.get_point(9)); + Min_sphere ms469(dimension, points469.begin(), points469.end()); + + Simplex_tree::Filtration_value f469 = st2.filtration(st2.find({4, 6, 9})); + std::cout << "f469= " << f469 << " | ms469_radius= " << std::sqrt(ms469.squared_radius()) << std::endl; + + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f469, std::sqrt(ms469.squared_radius())); + + BOOST_CHECK((st2.find({6, 7, 8}) == st2.null_simplex())); + BOOST_CHECK((st2.find({3, 5, 7}) == st2.null_simplex())); +} + +BOOST_AUTO_TEST_CASE(Cech_complex_from_points) { + // ---------------------------------------------------------------------------- + // Init of a list of points + // ---------------------------------------------------------------------------- + Point_cloud points; + std::vector<double> coords = {0.0, 0.0, 0.0, 1.0}; + points.push_back(Point(coords.begin(), coords.end())); + coords = {0.0, 0.0, 1.0, 0.0}; + points.push_back(Point(coords.begin(), coords.end())); + coords = {0.0, 1.0, 0.0, 0.0}; + points.push_back(Point(coords.begin(), coords.end())); + coords = {1.0, 0.0, 0.0, 0.0}; + points.push_back(Point(coords.begin(), coords.end())); + + // ---------------------------------------------------------------------------- + // Init of a Cech complex from the list of points + // ---------------------------------------------------------------------------- + Cech_complex cech_complex_from_points(points, 2.0); + + std::cout << "========== cech_complex_from_points ==========" << std::endl; + Simplex_tree st; + const int DIMENSION = 3; + cech_complex_from_points.create_complex(st, DIMENSION); + + // Another way to check num_simplices + std::cout << "Iterator on Cech complex simplices in the filtration order, with [filtration value]:" << std::endl; + int num_simplices = 0; + for (auto f_simplex : st.filtration_simplex_range()) { + num_simplices++; + std::cout << " ( "; + for (auto vertex : st.simplex_vertex_range(f_simplex)) { + std::cout << vertex << " "; + } + std::cout << ") -> " + << "[" << st.filtration(f_simplex) << "] "; + std::cout << std::endl; + } + BOOST_CHECK(num_simplices == 15); + std::cout << "st.num_simplices()=" << st.num_simplices() << std::endl; + BOOST_CHECK(st.num_simplices() == 15); + + std::cout << "st.dimension()=" << st.dimension() << std::endl; + BOOST_CHECK(st.dimension() == DIMENSION); + std::cout << "st.num_vertices()=" << st.num_vertices() << std::endl; + BOOST_CHECK(st.num_vertices() == 4); + + for (auto f_simplex : st.filtration_simplex_range()) { + std::cout << "dimension(" << st.dimension(f_simplex) << ") - f = " << st.filtration(f_simplex) << std::endl; + switch (st.dimension(f_simplex)) { + case 0: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.0); + break; + case 1: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.707107, .00001); + break; + case 2: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.816497, .00001); + break; + case 3: + GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.866025, .00001); + break; + default: + BOOST_CHECK(false); // Shall not happen + break; + } + } +} + +#ifdef GUDHI_DEBUG +BOOST_AUTO_TEST_CASE(Cech_create_complex_throw) { + // ---------------------------------------------------------------------------- + // + // Init of a Cech complex from a OFF file + // + // ---------------------------------------------------------------------------- + std::string off_file_name("alphacomplexdoc.off"); + double max_radius = 12.0; + std::cout << "========== OFF FILE NAME = " << off_file_name << " - Cech max_radius=" << max_radius + << "==========" << std::endl; + + Gudhi::Points_off_reader<Point> off_reader(off_file_name); + Cech_complex cech_complex_from_file(off_reader.get_point_cloud(), max_radius); + + Simplex_tree stree; + std::vector<int> simplex = {0, 1, 2}; + stree.insert_simplex_and_subfaces(simplex); + std::cout << "Check exception throw in debug mode" << std::endl; + // throw excpt because stree is not empty + BOOST_CHECK_THROW(cech_complex_from_file.create_complex(stree, 1), std::invalid_argument); +} +#endif diff --git a/src/Cech_complex/utilities/CMakeLists.txt b/src/Cech_complex/utilities/CMakeLists.txt new file mode 100644 index 00000000..30b99729 --- /dev/null +++ b/src/Cech_complex/utilities/CMakeLists.txt @@ -0,0 +1,14 @@ +cmake_minimum_required(VERSION 2.6) +project(Cech_complex_utilities) + +add_executable(cech_persistence cech_persistence.cpp) +target_link_libraries(cech_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY}) + +if (TBB_FOUND) + target_link_libraries(cech_persistence ${TBB_LIBRARIES}) +endif() + +add_test(NAME Cech_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:cech_persistence> + "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3") + +install(TARGETS cech_persistence DESTINATION bin) diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp new file mode 100644 index 00000000..8cfe018b --- /dev/null +++ b/src/Cech_complex/utilities/cech_persistence.cpp @@ -0,0 +1,124 @@ +/* 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 + */ + +#include <gudhi/Cech_complex.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Points_off_io.h> + +#include <boost/program_options.hpp> + +#include <string> +#include <vector> +#include <limits> // infinity + +// Types definition +using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>; +using Filtration_value = Simplex_tree::Filtration_value; +using Point = std::vector<double>; +using Point_cloud = std::vector<Point>; +using Points_off_reader = Gudhi::Points_off_reader<Point>; +using Cech_complex = Gudhi::cech_complex::Cech_complex<Simplex_tree, Point_cloud>; +using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; +using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>; + +void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, + Filtration_value& max_radius, int& dim_max, int& p, Filtration_value& min_persistence); + +int main(int argc, char* argv[]) { + std::string off_file_points; + std::string filediag; + Filtration_value max_radius; + int dim_max; + int p; + Filtration_value min_persistence; + + program_options(argc, argv, off_file_points, filediag, max_radius, dim_max, p, min_persistence); + + Points_off_reader off_reader(off_file_points); + Cech_complex cech_complex_from_file(off_reader.get_point_cloud(), max_radius); + + // Construct the Cech complex in a Simplex Tree + Simplex_tree simplex_tree; + + cech_complex_from_file.create_complex(simplex_tree, dim_max); + std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n"; + std::cout << " and has dimension " << simplex_tree.dimension() << " \n"; + + // Sort the simplices in the order of the filtration + simplex_tree.initialize_filtration(); + + // Compute the persistence diagram of the complex + Persistent_cohomology pcoh(simplex_tree); + // initializes the coefficient field for homology + pcoh.init_coefficients(p); + + pcoh.compute_persistent_cohomology(min_persistence); + + // Output the diagram in filediag + if (filediag.empty()) { + pcoh.output_diagram(); + } else { + std::ofstream out(filediag); + pcoh.output_diagram(out); + out.close(); + } + + return 0; +} + +void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, + Filtration_value& max_radius, int& dim_max, int& p, Filtration_value& min_persistence) { + namespace po = boost::program_options; + po::options_description hidden("Hidden options"); + hidden.add_options()("input-file", po::value<std::string>(&off_file_points), + "Name of an OFF file containing a point set.\n"); + + po::options_description visible("Allowed options", 100); + visible.add_options()("help,h", "produce help message")( + "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()), + "Name of file in which the persistence diagram is written. Default print in std::cout")( + "max-radius,r", + po::value<Filtration_value>(&max_radius)->default_value(std::numeric_limits<Filtration_value>::infinity()), + "Maximal length of an edge for the Cech complex construction.")( + "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1), + "Maximal dimension of the Cech complex we want to compute.")( + "field-charac,p", po::value<int>(&p)->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 a Cech 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/Cech_complex/utilities/cechcomplex.md b/src/Cech_complex/utilities/cechcomplex.md new file mode 100644 index 00000000..f7817dbb --- /dev/null +++ b/src/Cech_complex/utilities/cechcomplex.md @@ -0,0 +1,38 @@ + + +# Čech complex # + +## cech_persistence ## +This program computes the persistent homology with coefficient field *Z/pZ* of +a Čech complex defined on a set of input points, using Euclidean distance. 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** + +`cech_persistence [options] <OFF input file>` + +**Allowed options** + +* `-h [ --help ]` Produce help message +* `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. +* `-r [ --max-edge-length ]` (default = inf) Maximal length of an edge for the Čech complex construction. +* `-d [ --cpx-dimension ]` (default = 1) Maximal dimension of the Čech complex we want to compute. +* `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. +* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. + +Beware: this program may use a lot of RAM and take a lot of time if `max-edge-length` is set to a large value. + +**Example 1 with Z/2Z coefficients** + +`cech_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 2` + +**Example 2 with Z/3Z coefficients** + +`cech_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3` |
