diff options
author | ROUVREAU Vincent <vincent.rouvreau@inria.fr> | 2021-04-03 10:27:08 +0200 |
---|---|---|
committer | ROUVREAU Vincent <vincent.rouvreau@inria.fr> | 2021-04-03 10:27:08 +0200 |
commit | e4381a3e2ad79d3150cd03704bef3fc006e7c54b (patch) | |
tree | 755aa40eb646798ed87152d4940b9f4b8379ea96 | |
parent | 98cc06acb5f4b7caf4c23645614a472f7f9b5f3a (diff) |
Python alpha complex specific 3d with weighted version and functor to get points
-rw-r--r-- | src/python/CMakeLists.txt | 2 | ||||
-rw-r--r-- | src/python/gudhi/alpha_complex_3d.pyx | 129 | ||||
-rw-r--r-- | src/python/include/Alpha_complex_factory.h | 48 | ||||
-rw-r--r-- | src/python/include/Alpha_complex_interface.h | 59 | ||||
-rw-r--r-- | src/python/include/Alpha_complex_interface_3d.h | 71 |
5 files changed, 261 insertions, 48 deletions
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 73303a24..307181b7 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -61,6 +61,7 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'subsampling', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'tangential_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'alpha_complex', ") + set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'alpha_complex_3d', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'euclidean_witness_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'euclidean_strong_witness_complex', ") # Modules that should not be auto-imported in __init__.py @@ -156,6 +157,7 @@ if(PYTHONINTERP_FOUND) endif () if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ") + set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex_3d', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'subsampling', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'tangential_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'euclidean_witness_complex', ") diff --git a/src/python/gudhi/alpha_complex_3d.pyx b/src/python/gudhi/alpha_complex_3d.pyx new file mode 100644 index 00000000..3959004a --- /dev/null +++ b/src/python/gudhi/alpha_complex_3d.pyx @@ -0,0 +1,129 @@ +# 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) 2021 Inria +# +# Modification(s): +# - YYYY/MM Author: Description of the modification + +from __future__ import print_function +from cython cimport numeric +from libcpp.vector cimport vector +from libcpp.utility cimport pair +from libcpp.string cimport string +from libcpp cimport bool +from libc.stdint cimport intptr_t +import errno +import os +import warnings + +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) 2021 Inria" +__license__ = "GPL v3" + +cdef extern from "Alpha_complex_interface_3d.h" namespace "Gudhi": + cdef cppclass Alpha_complex_interface_3d "Gudhi::alpha_complex::Alpha_complex_interface_3d": + Alpha_complex_interface_3d(vector[vector[double]] points, vector[double] weights, 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) nogil except + + +# AlphaComplex3D python interface +cdef class AlphaComplex3D: + """AlphaComplex3D is a simplicial complex constructed from the finite cells + of a Delaunay Triangulation. + + The filtration value of each simplex is computed as the square of the + circumradius of the simplex if the circumsphere is empty (the simplex is + then said to be Gabriel), and as the minimum of the filtration values of + the codimension 1 cofaces that make it not Gabriel otherwise. + + All simplices that have a filtration value strictly greater than a given + alpha squared value are not inserted into the complex. + + .. note:: + + When AlphaComplex3D is constructed with an infinite value of alpha, the + complex is a Delaunay complex. + + """ + + cdef Alpha_complex_interface_3d * this_ptr + + # Fake constructor that does nothing but documenting the constructor + def __init__(self, points=[], weights=[], precision='safe'): + """AlphaComplex3D constructor. + + :param points: A list of points in d-Dimension. + :type points: Iterable[Iterable[float]] + + :param weights: A list of weights. If set, the number of weights must correspond to the + number of points. + :type weights: Iterable[float] + + :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is + 'safe'. + :type precision: string + + :raises ValueError: In case of inconsistency between the number of points and weights. + """ + + # The real cython constructor + def __cinit__(self, points = [], weights=[], precision = 'safe'): + assert precision in ['fast', 'safe', 'exact'], \ + "Alpha complex precision can only be 'fast', 'safe' or 'exact'" + cdef bool fast = precision == 'fast' + cdef bool exact = precision == 'exact' + + # weights are set but is inconsistent with the number of points + if len(weights) != 0 and len(weights) != len(points): + raise ValueError("Inconsistency between the number of points and weights") + + # need to copy the points to use them without the gil + cdef vector[vector[double]] pts + cdef vector[double] wgts + pts = points + wgts = weights + with nogil: + self.this_ptr = new Alpha_complex_interface_3d(pts, wgts, fast, exact) + + def __dealloc__(self): + if self.this_ptr != NULL: + del self.this_ptr + + def __is_defined(self): + """Returns true if AlphaComplex3D pointer is not NULL. + """ + return self.this_ptr != NULL + + def get_point(self, vertex): + """This function returns the point corresponding to a given vertex from the :class:`~gudhi.SimplexTree`. + + :param vertex: The vertex. + :type vertex: int + :rtype: list of float + :returns: the point. + """ + return self.this_ptr.get_point(vertex) + + def create_simplex_tree(self, max_alpha_square = float('inf')): + """ + :param max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to + infinity, and there is very little point using anything else since it does not save time. + :type max_alpha_square: float + :returns: A simplex tree created from the Delaunay Triangulation. + :rtype: SimplexTree + """ + stree = SimplexTree() + cdef double mas = max_alpha_square + cdef intptr_t stree_int_ptr=stree.thisptr + with nogil: + self.this_ptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr, + mas) + return stree diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h index 36e98615..5d3bfb65 100644 --- a/src/python/include/Alpha_complex_factory.h +++ b/src/python/include/Alpha_complex_factory.h @@ -31,6 +31,34 @@ namespace Gudhi { namespace alpha_complex { +template<typename CgalPointType, bool Weighted> +struct Point_cgal_to_cython; + +template<typename CgalPointType> +struct Point_cgal_to_cython<CgalPointType, false> { + std::vector<double> operator()(CgalPointType const& point) const + { + 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> +struct Point_cgal_to_cython<CgalPointType, true> { + std::vector<double> operator()(CgalPointType const& weighted_point) const + { + auto point = weighted_point.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> std::vector<double> pt_cgal_to_cython(CgalPointType const& point) { std::vector<double> vd; @@ -159,13 +187,14 @@ class Inexact_weighted_alpha_complex_dD final : public Abstract_alpha_complex { Alpha_complex<Kernel, true> alpha_complex_; }; -template <complexity Complexity> +template <complexity Complexity, bool Weighted = false> class Alpha_complex_3D final : public Abstract_alpha_complex { private: - using Point = typename Alpha_complex_3d<Complexity, false, false>::Bare_point_3; + using Bare_point = typename Alpha_complex_3d<Complexity, Weighted, false>::Bare_point_3; + using Point = typename Alpha_complex_3d<Complexity, Weighted, false>::Point_3; - static Point pt_cython_to_cgal_3(std::vector<double> const& vec) { - return Point(vec[0], vec[1], vec[2]); + static Bare_point pt_cython_to_cgal_3(std::vector<double> const& vec) { + return Bare_point(vec[0], vec[1], vec[2]); } public: @@ -173,18 +202,23 @@ class Alpha_complex_3D final : public Abstract_alpha_complex { : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3)) { } + Alpha_complex_3D(const std::vector<std::vector<double>>& points, const std::vector<double>& weights) + : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3), weights) { + } + virtual std::vector<double> get_point(int vh) override { Point const& point = alpha_complex_.get_point(vh); - return pt_cgal_to_cython(point); + return Point_cgal_to_cython<Point, Weighted>()(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); + alpha_complex_.create_complex(*simplex_tree, max_alpha_square); + return true; } private: - Alpha_complex_3d<Complexity, false, false> alpha_complex_; + Alpha_complex_3d<Complexity, Weighted, false> alpha_complex_; }; diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index 43c96b2f..31a8147b 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -30,10 +30,20 @@ class Alpha_complex_interface { Alpha_complex_interface(const std::vector<std::vector<double>>& points, const std::vector<double>& weights, bool fast_version, bool exact_version) - : points_(points), - weights_(weights), - fast_version_(fast_version), - exact_version_(exact_version) { + : empty_point_set_(points.size() == 0) { + if (fast_version) { + if (weights.size() == 0) { + alpha_ptr_ = std::make_unique<Inexact_alpha_complex_dD>(points, exact_version); + } else { + alpha_ptr_ = std::make_unique<Inexact_weighted_alpha_complex_dD>(points, weights, exact_version); + } + } else { + if (weights.size() == 0) { + alpha_ptr_ = std::make_unique<Exact_alpha_complex_dD>(points, exact_version); + } else { + alpha_ptr_ = std::make_unique<Exact_weighted_alpha_complex_dD>(points, weights, exact_version); + } + } } std::vector<double> get_point(int vh) { @@ -42,47 +52,14 @@ class Alpha_complex_interface { void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) { - if (points_.size() > 0) { - std::size_t dimension = points_[0].size(); - if (dimension == 3 && weights_.size() == 0 && !default_filtration_value) { - if (fast_version_) - alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::FAST>>(points_); - else if (exact_version_) - alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::EXACT>>(points_); - else - alpha_ptr_ = std::make_unique<Alpha_complex_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 || weights_.size() != 0 || default_filtration_value) { - if (fast_version_) { - if (weights_.size() == 0) { - alpha_ptr_ = std::make_unique<Inexact_alpha_complex_dD>(points_, exact_version_); - } else { - alpha_ptr_ = std::make_unique<Inexact_weighted_alpha_complex_dD>(points_, weights_, exact_version_); - } - } else { - if (weights_.size() == 0) { - alpha_ptr_ = std::make_unique<Exact_alpha_complex_dD>(points_, exact_version_); - } else { - alpha_ptr_ = std::make_unique<Exact_weighted_alpha_complex_dD>(points_, weights_, exact_version_); - } - } - alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value); - } - } + // Nothing to be done in case of an empty point set + if (!empty_point_set_) + 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_; - std::vector<double> weights_; - bool fast_version_; - bool exact_version_; + bool empty_point_set_; }; } // namespace alpha_complex diff --git a/src/python/include/Alpha_complex_interface_3d.h b/src/python/include/Alpha_complex_interface_3d.h new file mode 100644 index 00000000..bb66b8e1 --- /dev/null +++ b/src/python/include/Alpha_complex_interface_3d.h @@ -0,0 +1,71 @@ +/* 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) 2021 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef INCLUDE_ALPHA_COMPLEX_INTERFACE_3D_H_ +#define INCLUDE_ALPHA_COMPLEX_INTERFACE_3D_H_ + +#include "Alpha_complex_factory.h" +#include <gudhi/Alpha_complex_options.h> + +#include "Simplex_tree_interface.h" + +#include <iostream> +#include <vector> +#include <string> +#include <memory> // for std::unique_ptr + +namespace Gudhi { + +namespace alpha_complex { + +class Alpha_complex_interface_3d { + public: + Alpha_complex_interface_3d(const std::vector<std::vector<double>>& points, + const std::vector<double>& weights, + bool fast_version, bool exact_version) + : empty_point_set_(points.size() == 0) { + const bool weighted = (weights.size() > 0); + if (fast_version) + if (weighted) + alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::FAST, true>>(points, weights); + else + alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::FAST>>(points); + else if (exact_version) + if (weighted) + alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::EXACT, true>>(points, weights); + else + alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::EXACT>>(points); + else + if (weighted) + alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::SAFE, true>>(points, weights); + else + alpha_ptr_ = std::make_unique<Alpha_complex_3D<Gudhi::alpha_complex::complexity::SAFE>>(points); + } + + std::vector<double> get_point(int vh) { + return alpha_ptr_->get_point(vh); + } + + void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square) { + // Nothing to be done in case of an empty point set + if (!empty_point_set_) + alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, false); + } + + private: + std::unique_ptr<Abstract_alpha_complex> alpha_ptr_; + bool empty_point_set_; +}; + +} // namespace alpha_complex + +} // namespace Gudhi + +#endif // INCLUDE_ALPHA_COMPLEX_INTERFACE_3D_H_ |