diff options
author | Vincent Rouvreau <10407034+VincentRouvreau@users.noreply.github.com> | 2020-07-02 00:32:34 -0700 |
---|---|---|
committer | GitHub <noreply@github.com> | 2020-07-02 00:32:34 -0700 |
commit | 5131abd569ae5c4f11e753f1a6dc1ee232bcb96f (patch) | |
tree | 786e04a0b9ead4620d5fcdf7cb67befeaf4646f2 /src | |
parent | 8d316e831c6af51efb9c362a5b203528a9fd3b15 (diff) | |
parent | f39473b7c1de0fe42b3f4ebf5cb37bca0a84a247 (diff) |
Merge pull request #355 from VincentRouvreau/alpha_complex_3d_python
Alpha complex is using 3d version in dim 3
Diffstat (limited to 'src')
-rw-r--r-- | src/Alpha_complex/include/gudhi/Alpha_complex_3d.h | 56 | ||||
-rw-r--r-- | src/python/CMakeLists.txt | 28 | ||||
-rw-r--r-- | src/python/doc/alpha_complex_user.rst | 51 | ||||
-rw-r--r-- | src/python/gudhi/alpha_complex.pyx | 29 | ||||
-rw-r--r-- | src/python/include/Alpha_complex_factory.h | 139 | ||||
-rw-r--r-- | src/python/include/Alpha_complex_interface.h | 92 | ||||
-rwxr-xr-x | src/python/test/test_alpha_complex.py | 83 |
7 files changed, 312 insertions, 166 deletions
diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h index c19ebb79..f56e12d0 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h @@ -14,6 +14,9 @@ #include <boost/version.hpp> #include <boost/variant.hpp> +#include <boost/range/size.hpp> +#include <boost/range/combine.hpp> +#include <boost/range/adaptor/transformed.hpp> #include <gudhi/Debug_utils.h> #include <gudhi/Alpha_complex_options.h> @@ -277,8 +280,8 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ 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)); + alpha_shape_3_ptr_ = std::make_unique<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. @@ -299,20 +302,15 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ 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()), + // FIXME: this test is only valid if we have a forward range + GUDHI_CHECK(boost::size(weights) == boost::size(points), std::invalid_argument("Points number in range different from weights range number")); - std::vector<Weighted_point_3> weighted_points_3; + auto weighted_points_3 = boost::range::combine(points, weights) + | boost::adaptors::transformed([](auto const&t){return Weighted_point_3(boost::get<0>(t), boost::get<1>(t));}); - 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)); + alpha_shape_3_ptr_ = std::make_unique<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. @@ -356,7 +354,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ // 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)); + alpha_shape_3_ptr_ = std::make_unique<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. @@ -388,31 +386,27 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ 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()), + // FIXME: this test is only valid if we have a forward range + GUDHI_CHECK(boost::size(weights) == boost::size(points), 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++; - } + auto weighted_points_3 = boost::range::combine(points, weights) + | boost::adaptors::transformed([=](auto const&t){ + auto w = boost::get<1>(t); + GUDHI_CHECK((w < maximal_possible_weight) && (w >= 0), + std::invalid_argument("Invalid weight " + std::to_string(w) + + ". Must be non-negative and less than maximal possible weight = 1/64*cuboid length squared.")); + return Weighted_point_3(boost::get<0>(t), w); + }); // Define the periodic cube Dt pdt(typename Kernel::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max)); @@ -426,7 +420,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ // 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)); + alpha_shape_3_ptr_ = std::make_unique<Alpha_shape_3>(pdt, 0, Alpha_shape_3::GENERAL); } /** \brief Inserts all Delaunay triangulation into the simplicial complex. @@ -471,6 +465,10 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ #ifdef DEBUG_TRACES std::clog << "filtration_with_alpha_values returns : " << objects.size() << " objects" << std::endl; #endif // DEBUG_TRACES + if (objects.size() == 0) { + std::cerr << "Alpha_complex_3d create_complex - no triangulation as points are on a 2d plane\n"; + return false; // ----- >> + } using Alpha_value_iterator = typename std::vector<FT>::const_iterator; Alpha_value_iterator alpha_value_iterator = alpha_values.begin(); diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 04f64ce3..4f26481e 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -14,6 +14,15 @@ function( add_GUDHI_PYTHON_lib THE_LIB ) endif(EXISTS ${THE_LIB}) endfunction( add_GUDHI_PYTHON_lib ) +function( add_GUDHI_PYTHON_lib_dir THE_LIB_DIR ) + # deals when it is not set - error on windows + if(EXISTS ${THE_LIB_DIR}) + set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${THE_LIB_DIR}', " PARENT_SCOPE) + else() + message("add_GUDHI_PYTHON_lib_dir - '${THE_LIB_DIR}' does not exist") + endif() +endfunction( add_GUDHI_PYTHON_lib_dir ) + # THE_TEST is the python test file name (without .py extension) containing tests functions function( add_gudhi_py_test THE_TEST ) if(PYTEST_FOUND) @@ -154,7 +163,7 @@ if(PYTHONINTERP_FOUND) else(CGAL_HEADER_ONLY) add_gudhi_debug_info("CGAL version ${CGAL_VERSION}") add_GUDHI_PYTHON_lib("${CGAL_LIBRARY}") - set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${CGAL_LIBRARIES_DIR}', ") + add_GUDHI_PYTHON_lib_dir("${CGAL_LIBRARIES_DIR}") message("** Add CGAL ${CGAL_LIBRARIES_DIR}") # If CGAL is not header only, CGAL library may link with boost system, if(CMAKE_BUILD_TYPE MATCHES Debug) @@ -162,7 +171,7 @@ if(PYTHONINTERP_FOUND) else() add_GUDHI_PYTHON_lib("${Boost_SYSTEM_LIBRARY_RELEASE}") endif() - set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${Boost_LIBRARY_DIRS}', ") + add_GUDHI_PYTHON_lib_dir("${Boost_LIBRARY_DIRS}") message("** Add Boost ${Boost_LIBRARY_DIRS}") endif(CGAL_HEADER_ONLY) # GMP and GMPXX are not required, but if present, CGAL will link with them. @@ -170,13 +179,17 @@ if(PYTHONINTERP_FOUND) add_gudhi_debug_info("GMP_LIBRARIES = ${GMP_LIBRARIES}") set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_USE_GMP', ") add_GUDHI_PYTHON_lib("${GMP_LIBRARIES}") - set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${GMP_LIBRARIES_DIR}', ") + if(NOT GMP_LIBRARIES_DIR) + get_filename_component(GMP_LIBRARIES_DIR ${GMP_LIBRARIES} PATH) + message("GMP_LIBRARIES_DIR from GMP_LIBRARIES set to ${GMP_LIBRARIES_DIR}") + endif(NOT GMP_LIBRARIES_DIR) + add_GUDHI_PYTHON_lib_dir("${GMP_LIBRARIES_DIR}") message("** Add gmp ${GMP_LIBRARIES_DIR}") if(GMPXX_FOUND) add_gudhi_debug_info("GMPXX_LIBRARIES = ${GMPXX_LIBRARIES}") set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_USE_GMPXX', ") add_GUDHI_PYTHON_lib("${GMPXX_LIBRARIES}") - set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${GMPXX_LIBRARIES_DIR}', ") + add_GUDHI_PYTHON_lib_dir("${GMPXX_LIBRARIES_DIR}") message("** Add gmpxx ${GMPXX_LIBRARIES_DIR}") endif(GMPXX_FOUND) endif(GMP_FOUND) @@ -187,9 +200,10 @@ if(PYTHONINTERP_FOUND) # In case CGAL is not header only, all MPFR variables are set except MPFR_LIBRARIES_DIR - Just set it if(NOT MPFR_LIBRARIES_DIR) get_filename_component(MPFR_LIBRARIES_DIR ${MPFR_LIBRARIES} PATH) + message("MPFR_LIBRARIES_DIR from MPFR_LIBRARIES set to ${MPFR_LIBRARIES_DIR}") endif(NOT MPFR_LIBRARIES_DIR) - set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${MPFR_LIBRARIES_DIR}', ") - message("** Add mpfr ${MPFR_LIBRARIES}") + add_GUDHI_PYTHON_lib_dir("${MPFR_LIBRARIES_DIR}") + message("** Add mpfr ${MPFR_LIBRARIES_DIR}") endif(MPFR_FOUND) endif(CGAL_FOUND) @@ -216,7 +230,7 @@ if(PYTHONINTERP_FOUND) add_GUDHI_PYTHON_lib("${TBB_RELEASE_LIBRARY}") add_GUDHI_PYTHON_lib("${TBB_MALLOC_RELEASE_LIBRARY}") endif() - set(GUDHI_PYTHON_LIBRARY_DIRS "${GUDHI_PYTHON_LIBRARY_DIRS}'${TBB_LIBRARY_DIRS}', ") + add_GUDHI_PYTHON_lib_dir("${TBB_LIBRARY_DIRS}") message("** Add tbb ${TBB_LIBRARY_DIRS}") set(GUDHI_PYTHON_INCLUDE_DIRS "${GUDHI_PYTHON_INCLUDE_DIRS}'${TBB_INCLUDE_DIRS}', ") endif() diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst index 097470c1..fffcb3db 100644 --- a/src/python/doc/alpha_complex_user.rst +++ b/src/python/doc/alpha_complex_user.rst @@ -34,8 +34,8 @@ Remarks the computation of filtration values can exceptionally be arbitrarily bad. In all cases, we still guarantee that the output is a valid filtration (faces have a filtration value no larger than their cofaces). * For performances reasons, it is advised to use Alpha_complex with `CGAL <installation.html#cgal>`_ :math:`\geq` 5.0.0. - -For performances reasons, it is advised to use CGAL :math:`\geq` 5.0.0. +* The vertices in the output simplex tree are not guaranteed to match the order of the input points. One can use + :func:`~gudhi.AlphaComplex.get_point` to get the initial point back. Example from points ------------------- @@ -178,49 +178,22 @@ In the following example, a threshold of :math:`\alpha^2 = 32.0` is used. Example from OFF file ^^^^^^^^^^^^^^^^^^^^^ -This example builds the Delaunay triangulation from the points given by an OFF file, and initializes the alpha complex -with it. - +This example builds the alpha complex from 300 random points on a 2-torus. -Then, it is asked to display information about the alpha complex: +Then, it computes the persistence diagram and displays it: -.. testcode:: +.. plot:: + :include-source: + import matplotlib.pyplot as plt import gudhi alpha_complex = gudhi.AlphaComplex(off_file=gudhi.__root_source_dir__ + \ - '/data/points/alphacomplexdoc.off') - simplex_tree = alpha_complex.create_simplex_tree(max_alpha_square=32.0) + '/data/points/tore3D_300.off') + simplex_tree = alpha_complex.create_simplex_tree() result_str = 'Alpha complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \ repr(simplex_tree.num_simplices()) + ' simplices - ' + \ repr(simplex_tree.num_vertices()) + ' vertices.' print(result_str) - fmt = '%s -> %.2f' - for filtered_value in simplex_tree.get_filtration(): - print(fmt % tuple(filtered_value)) - -the program output is: - -.. testoutput:: - - Alpha complex is of dimension 2 - 20 simplices - 7 vertices. - [0] -> 0.00 - [1] -> 0.00 - [2] -> 0.00 - [3] -> 0.00 - [4] -> 0.00 - [5] -> 0.00 - [6] -> 0.00 - [2, 3] -> 6.25 - [4, 5] -> 7.25 - [0, 2] -> 8.50 - [0, 1] -> 9.25 - [1, 3] -> 10.00 - [1, 2] -> 11.25 - [1, 2, 3] -> 12.50 - [0, 1, 2] -> 13.00 - [5, 6] -> 13.25 - [2, 4] -> 20.00 - [4, 6] -> 22.74 - [4, 5, 6] -> 22.74 - [3, 6] -> 30.25 - + diag = simplex_tree.persistence() + gudhi.plot_persistence_diagram(diag) + plt.show() diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx index a356384d..ea128743 100644 --- a/src/python/gudhi/alpha_complex.pyx +++ b/src/python/gudhi/alpha_complex.pyx @@ -20,6 +20,7 @@ import os from gudhi.simplex_tree cimport * from gudhi.simplex_tree import SimplexTree +from gudhi import read_points_from_off_file __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" @@ -27,11 +28,9 @@ __license__ = "GPL v3" cdef extern from "Alpha_complex_interface.h" namespace "Gudhi": cdef cppclass Alpha_complex_interface "Gudhi::alpha_complex::Alpha_complex_interface": - Alpha_complex_interface(vector[vector[double]] points, bool fast_version) nogil except + - # bool from_file is a workaround for cython to find the correct signature - Alpha_complex_interface(string off_file, bool fast_version, bool from_file) nogil except + + Alpha_complex_interface(vector[vector[double]] points, bool fast_version, bool exact_version) nogil except + vector[double] get_point(int vertex) nogil except + - void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool exact_version, bool default_filtration_value) nogil except + + void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool default_filtration_value) nogil except + # AlphaComplex python interface cdef class AlphaComplex: @@ -54,7 +53,6 @@ cdef class AlphaComplex: """ cdef Alpha_complex_interface * this_ptr - cdef bool exact # Fake constructor that does nothing but documenting the constructor def __init__(self, points=None, off_file='', precision='safe'): @@ -76,21 +74,20 @@ cdef class AlphaComplex: def __cinit__(self, points = None, off_file = '', precision = 'safe'): assert precision in ['fast', 'safe', 'exact'], "Alpha complex precision can only be 'fast', 'safe' or 'exact'" cdef bool fast = precision == 'fast' - self.exact = precision == 'exact' + cdef bool exact = precision == 'exact' cdef vector[vector[double]] pts if off_file: if os.path.isfile(off_file): - self.this_ptr = new Alpha_complex_interface(off_file.encode('utf-8'), fast, True) + points = read_points_from_off_file(off_file = off_file) else: print("file " + off_file + " not found.") - else: - if points is None: - # Empty Alpha construction - points=[] - pts = points - with nogil: - self.this_ptr = new Alpha_complex_interface(pts, fast) + if points is None: + # Empty Alpha construction + points=[] + pts = points + with nogil: + self.this_ptr = new Alpha_complex_interface(pts, fast, exact) def __dealloc__(self): if self.this_ptr != NULL: @@ -102,7 +99,7 @@ cdef class AlphaComplex: return self.this_ptr != NULL def get_point(self, vertex): - """This function returns the point corresponding to a given vertex. + """This function returns the point corresponding to a given vertex from the :class:`~gudhi.SimplexTree`. :param vertex: The vertex. :type vertex: int @@ -128,5 +125,5 @@ cdef class AlphaComplex: cdef bool compute_filtration = default_filtration_value == True with nogil: self.this_ptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr, - mas, self.exact, compute_filtration) + mas, compute_filtration) return stree diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h new file mode 100644 index 00000000..d699ad9b --- /dev/null +++ b/src/python/include/Alpha_complex_factory.h @@ -0,0 +1,139 @@ +/* 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) 2020 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef INCLUDE_ALPHA_COMPLEX_FACTORY_H_ +#define INCLUDE_ALPHA_COMPLEX_FACTORY_H_ + +#include <gudhi/Simplex_tree.h> +#include <gudhi/Alpha_complex.h> +#include <gudhi/Alpha_complex_3d.h> +#include <gudhi/Alpha_complex_options.h> +#include <CGAL/Epeck_d.h> +#include <CGAL/Epick_d.h> + +#include <boost/range/adaptor/transformed.hpp> + +#include "Simplex_tree_interface.h" + +#include <iostream> +#include <vector> +#include <string> +#include <memory> // for std::unique_ptr + +namespace Gudhi { + +namespace alpha_complex { + +template <typename CgalPointType> +std::vector<double> pt_cgal_to_cython(CgalPointType const& point) { + std::vector<double> vd; + vd.reserve(point.dimension()); + for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++) + vd.push_back(CGAL::to_double(*coord)); + return vd; +} + +template <typename CgalPointType> +static CgalPointType pt_cython_to_cgal(std::vector<double> const& vec) { + return CgalPointType(vec.size(), vec.begin(), vec.end()); +} + +class Abstract_alpha_complex { + public: + virtual std::vector<double> get_point(int vh) = 0; + virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, + bool default_filtration_value) = 0; +}; + +class Exact_Alphacomplex_dD : public Abstract_alpha_complex { + private: + using Kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>; + using Point = typename Kernel::Point_d; + + public: + Exact_Alphacomplex_dD(const std::vector<std::vector<double>>& points, bool exact_version) + : exact_version_(exact_version), + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Point>)) { + } + + virtual std::vector<double> get_point(int vh) override { + Point const& point = alpha_complex_.get_point(vh); + return pt_cgal_to_cython(point); + } + + virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, + bool default_filtration_value) override { + return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value); + } + + private: + bool exact_version_; + Alpha_complex<Kernel> alpha_complex_; +}; + +class Inexact_Alphacomplex_dD : public Abstract_alpha_complex { + private: + using Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; + using Point = typename Kernel::Point_d; + + public: + Inexact_Alphacomplex_dD(const std::vector<std::vector<double>>& points, bool exact_version) + : exact_version_(exact_version), + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal<Point>)) { + } + + virtual std::vector<double> get_point(int vh) override { + Point const& point = alpha_complex_.get_point(vh); + return pt_cgal_to_cython(point); + } + virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, + bool default_filtration_value) override { + return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value); + } + + private: + bool exact_version_; + Alpha_complex<Kernel> alpha_complex_; +}; + +template <complexity Complexity> +class Alphacomplex_3D : public Abstract_alpha_complex { + private: + using Point = typename Alpha_complex_3d<Complexity, false, false>::Bare_point_3; + + static Point pt_cython_to_cgal_3(std::vector<double> const& vec) { + return Point(vec[0], vec[1], vec[2]); + } + + public: + Alphacomplex_3D(const std::vector<std::vector<double>>& points) + : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3)) { + } + + virtual std::vector<double> get_point(int vh) override { + Point const& point = alpha_complex_.get_point(vh); + return pt_cgal_to_cython(point); + } + + virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, + bool default_filtration_value) override { + return alpha_complex_.create_complex(*simplex_tree, max_alpha_square); + } + + private: + Alpha_complex_3d<Complexity, false, false> alpha_complex_; +}; + + +} // namespace alpha_complex + +} // namespace Gudhi + +#endif // INCLUDE_ALPHA_COMPLEX_FACTORY_H_ diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index 3ac5db1f..23be194d 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -11,12 +11,8 @@ #ifndef INCLUDE_ALPHA_COMPLEX_INTERFACE_H_ #define INCLUDE_ALPHA_COMPLEX_INTERFACE_H_ -#include <gudhi/Simplex_tree.h> -#include <gudhi/Alpha_complex.h> -#include <CGAL/Epeck_d.h> -#include <CGAL/Epick_d.h> - -#include <boost/range/adaptor/transformed.hpp> +#include "Alpha_complex_factory.h" +#include <gudhi/Alpha_complex_options.h> #include "Simplex_tree_interface.h" @@ -30,67 +26,51 @@ namespace Gudhi { namespace alpha_complex { class Alpha_complex_interface { - private: - using Exact_kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>; - using Inexact_kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>; - using Point_exact_kernel = typename Exact_kernel::Point_d; - using Point_inexact_kernel = typename Inexact_kernel::Point_d; - - template <typename CgalPointType> - std::vector<double> pt_cgal_to_cython(CgalPointType& point) { - std::vector<double> vd; - for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++) - vd.push_back(CGAL::to_double(*coord)); - return vd; - } - - template <typename CgalPointType> - static CgalPointType pt_cython_to_cgal(std::vector<double> const& vec) { - return CgalPointType(vec.size(), vec.begin(), vec.end()); - } - public: - Alpha_complex_interface(const std::vector<std::vector<double>>& points, bool fast_version) - : fast_version_(fast_version) { - if (fast_version_) { - ac_inexact_ptr_ = std::make_unique<Alpha_complex<Inexact_kernel>>( - boost::adaptors::transform(points, pt_cython_to_cgal<Point_inexact_kernel>)); - } else { - ac_exact_ptr_ = std::make_unique<Alpha_complex<Exact_kernel>>( - boost::adaptors::transform(points, pt_cython_to_cgal<Point_exact_kernel>)); - } - } - - Alpha_complex_interface(const std::string& off_file_name, bool fast_version, bool from_file = true) - : fast_version_(fast_version) { - if (fast_version_) - ac_inexact_ptr_ = std::make_unique<Alpha_complex<Inexact_kernel>>(off_file_name); - else - ac_exact_ptr_ = std::make_unique<Alpha_complex<Exact_kernel>>(off_file_name); + Alpha_complex_interface(const std::vector<std::vector<double>>& points, bool fast_version, bool exact_version) + : points_(points), + fast_version_(fast_version), + exact_version_(exact_version) { } std::vector<double> get_point(int vh) { - if (fast_version_) { - Point_inexact_kernel const& point = ac_inexact_ptr_->get_point(vh); - return pt_cgal_to_cython(point); - } else { - Point_exact_kernel const& point = ac_exact_ptr_->get_point(vh); - return pt_cgal_to_cython(point); - } + return alpha_ptr_->get_point(vh); } - void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool exact_version, + void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) { - if (fast_version_) - ac_inexact_ptr_->create_complex(*simplex_tree, max_alpha_square, exact_version, default_filtration_value); - else - ac_exact_ptr_->create_complex(*simplex_tree, max_alpha_square, exact_version, default_filtration_value); + if (points_.size() > 0) { + std::size_t dimension = points_[0].size(); + if (dimension == 3 && !default_filtration_value) { + if (fast_version_) + alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::FAST>>(points_); + else if (exact_version_) + alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::EXACT>>(points_); + else + alpha_ptr_ = std::make_unique<Alphacomplex_3D<Gudhi::alpha_complex::complexity::SAFE>>(points_); + if (!alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value)) { + // create_simplex_tree will fail if all points are on a plane - Retry with dD by setting dimension to 2 + dimension--; + alpha_ptr_.reset(); + } + } + // Not ** else ** because we have to take into account if 3d fails + if (dimension != 3 || default_filtration_value) { + if (fast_version_) { + alpha_ptr_ = std::make_unique<Inexact_Alphacomplex_dD>(points_, exact_version_); + } else { + alpha_ptr_ = std::make_unique<Exact_Alphacomplex_dD>(points_, exact_version_); + } + alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value); + } + } } private: + std::unique_ptr<Abstract_alpha_complex> alpha_ptr_; + std::vector<std::vector<double>> points_; bool fast_version_; - std::unique_ptr<Alpha_complex<Exact_kernel>> ac_exact_ptr_; - std::unique_ptr<Alpha_complex<Inexact_kernel>> ac_inexact_ptr_; + bool exact_version_; }; } // namespace alpha_complex diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py index 943ad2c4..a4ee260b 100755 --- a/src/python/test/test_alpha_complex.py +++ b/src/python/test/test_alpha_complex.py @@ -8,7 +8,7 @@ - YYYY/MM Author: Description of the modification """ -from gudhi import AlphaComplex, SimplexTree +import gudhi as gd import math import numpy as np import pytest @@ -25,17 +25,16 @@ __license__ = "MIT" def _empty_alpha(precision): - alpha_complex = AlphaComplex(points=[[0, 0]], precision = precision) + alpha_complex = gd.AlphaComplex(points=[[0, 0]], precision = precision) assert alpha_complex.__is_defined() == True def test_empty_alpha(): - _empty_alpha('fast') - _empty_alpha('safe') - _empty_alpha('exact') + for precision in ['fast', 'safe', 'exact']: + _empty_alpha(precision) def _infinite_alpha(precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - alpha_complex = AlphaComplex(points=point_list, precision = precision) + alpha_complex = gd.AlphaComplex(points=point_list, precision = precision) assert alpha_complex.__is_defined() == True simplex_tree = alpha_complex.create_simplex_tree() @@ -84,13 +83,12 @@ def _infinite_alpha(precision): assert False def test_infinite_alpha(): - _infinite_alpha('fast') - _infinite_alpha('safe') - _infinite_alpha('exact') + for precision in ['fast', 'safe', 'exact']: + _infinite_alpha(precision) def _filtered_alpha(precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - filtered_alpha = AlphaComplex(points=point_list, precision = precision) + filtered_alpha = gd.AlphaComplex(points=point_list, precision = precision) simplex_tree = filtered_alpha.create_simplex_tree(max_alpha_square=0.25) @@ -128,9 +126,8 @@ def _filtered_alpha(precision): assert simplex_tree.get_cofaces([0], 1) == [([0, 1], 0.25), ([0, 2], 0.25)] def test_filtered_alpha(): - _filtered_alpha('fast') - _filtered_alpha('safe') - _filtered_alpha('exact') + for precision in ['fast', 'safe', 'exact']: + _filtered_alpha(precision) def _safe_alpha_persistence_comparison(precision): #generate periodic signal @@ -144,10 +141,10 @@ def _safe_alpha_persistence_comparison(precision): embedding2 = [[signal[i], delayed[i]] for i in range(len(time))] #build alpha complex and simplex tree - alpha_complex1 = AlphaComplex(points=embedding1, precision = precision) + alpha_complex1 = gd.AlphaComplex(points=embedding1, precision = precision) simplex_tree1 = alpha_complex1.create_simplex_tree() - alpha_complex2 = AlphaComplex(points=embedding2, precision = precision) + alpha_complex2 = gd.AlphaComplex(points=embedding2, precision = precision) simplex_tree2 = alpha_complex2.create_simplex_tree() diag1 = simplex_tree1.persistence() @@ -165,7 +162,7 @@ def test_safe_alpha_persistence_comparison(): def _delaunay_complex(precision): point_list = [[0, 0], [1, 0], [0, 1], [1, 1]] - filtered_alpha = AlphaComplex(points=point_list, precision = precision) + filtered_alpha = gd.AlphaComplex(points=point_list, precision = precision) simplex_tree = filtered_alpha.create_simplex_tree(default_filtration_value = True) @@ -197,6 +194,54 @@ def _delaunay_complex(precision): assert math.isnan(filtered_value[1]) def test_delaunay_complex(): - _delaunay_complex('fast') - _delaunay_complex('safe') - _delaunay_complex('exact') + for precision in ['fast', 'safe', 'exact']: + _delaunay_complex(precision) + +def _3d_points_on_a_plane(precision, default_filtration_value): + alpha = gd.AlphaComplex(off_file=gd.__root_source_dir__ + '/data/points/alphacomplexdoc.off', + precision = precision) + + simplex_tree = alpha.create_simplex_tree(default_filtration_value = default_filtration_value) + assert simplex_tree.dimension() == 2 + assert simplex_tree.num_vertices() == 7 + assert simplex_tree.num_simplices() == 25 + +def test_3d_points_on_a_plane(): + for default_filtration_value in [True, False]: + for precision in ['fast', 'safe', 'exact']: + _3d_points_on_a_plane(precision, default_filtration_value) + +def _3d_tetrahedrons(precision): + points = 10*np.random.rand(10, 3) + alpha = gd.AlphaComplex(points=points, precision = precision) + st_alpha = alpha.create_simplex_tree(default_filtration_value = False) + # New AlphaComplex for get_point to work + delaunay = gd.AlphaComplex(points=points, precision = precision) + st_delaunay = delaunay.create_simplex_tree(default_filtration_value = True) + + delaunay_tetra = [] + for sk in st_delaunay.get_skeleton(4): + if len(sk[0]) == 4: + tetra = [delaunay.get_point(sk[0][0]), + delaunay.get_point(sk[0][1]), + delaunay.get_point(sk[0][2]), + delaunay.get_point(sk[0][3]) ] + delaunay_tetra.append(sorted(tetra, key=lambda tup: tup[0])) + + alpha_tetra = [] + for sk in st_alpha.get_skeleton(4): + if len(sk[0]) == 4: + tetra = [alpha.get_point(sk[0][0]), + alpha.get_point(sk[0][1]), + alpha.get_point(sk[0][2]), + alpha.get_point(sk[0][3]) ] + alpha_tetra.append(sorted(tetra, key=lambda tup: tup[0])) + + # Check the tetrahedrons from one list are in the second one + assert len(alpha_tetra) == len(delaunay_tetra) + for tetra_from_del in delaunay_tetra: + assert tetra_from_del in alpha_tetra + +def test_3d_tetrahedrons(): + for precision in ['fast', 'safe', 'exact']: + _3d_tetrahedrons(precision) |