diff options
Diffstat (limited to 'src/Cech_complex/include/gudhi')
| -rw-r--r-- | src/Cech_complex/include/gudhi/Cech_complex.h | 92 | ||||
| -rw-r--r-- | src/Cech_complex/include/gudhi/Cech_complex_blocker.h | 90 | ||||
| -rw-r--r-- | src/Cech_complex/include/gudhi/Miniball.COPYRIGHT | 4 | ||||
| -rw-r--r-- | src/Cech_complex/include/gudhi/Miniball.README | 26 | ||||
| -rw-r--r-- | src/Cech_complex/include/gudhi/Miniball.hpp | 523 | ||||
| -rw-r--r-- | src/Cech_complex/include/gudhi/Sphere_circumradius.h | 78 |
6 files changed, 209 insertions, 604 deletions
diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index b0871e10..dbdf5e93 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -1,24 +1,24 @@ /* 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 + * Author(s): Vincent Rouvreau, Hind Montassif * * Copyright (C) 2018 Inria * * Modification(s): * - YYYY/MM Author: Description of the modification + * - 2022/02 Hind Montassif : Replace MiniBall with Sphere_circumradius */ #ifndef CECH_COMPLEX_H_ #define CECH_COMPLEX_H_ -#include <gudhi/distance_functions.h> // for Gudhi::Minimal_enclosing_ball_radius +#include <gudhi/Sphere_circumradius.h> // for Gudhi::cech_complex::Sphere_circumradius #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 { @@ -26,55 +26,55 @@ namespace cech_complex { /** * \class Cech_complex - * \brief Cech complex data structure. + * \brief Cech complex class. * * \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. + * Cech complex is a simplicial complex where the set of all simplices is filtered + * by the radius of their minimal enclosing ball and bounded by the given max_radius. * - * \tparam SimplicialComplexForProximityGraph furnishes `Vertex_handle` and `Filtration_value` type definition required - * by `Gudhi::Proximity_graph`. + * \tparam Kernel CGAL kernel: either Epick_d or Epeck_d. + * + * \tparam SimplicialComplexForCechComplex furnishes `Vertex_handle` and `Filtration_value` type definition required + * by `Gudhi::Proximity_graph` and Cech blocker. * - * \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> +template <typename Kernel, typename SimplicialComplexForCechComplex> 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. + using Vertex_handle = typename SimplicialComplexForCechComplex::Vertex_handle; + using Filtration_value = typename SimplicialComplexForCechComplex::Filtration_value; + using Proximity_graph = Gudhi::Proximity_graph<SimplicialComplexForCechComplex>; + + using cech_blocker = Cech_blocker<SimplicialComplexForCechComplex, Cech_complex, Kernel>; + + using Point_d = typename cech_blocker::Point_d; + using Point_cloud = std::vector<Point_d>; + + // Numeric type of coordinates in the kernel + using FT = typename cech_blocker::FT; + // Sphere is a pair of point and squared radius. + using Sphere = typename cech_blocker::Sphere; + + public: + /** \brief Cech_complex constructor from a range of points. * - * @param[in] points Range of points. + * @param[in] points Range of points where each point is defined as `kernel::Point_d`. * @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. + * @param[in] exact Exact filtration values computation. Not exact if `Kernel` is not <a target="_blank" + * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>. + * Default is false. * */ - 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)); + template<typename InputPointRange > + Cech_complex(const InputPointRange & points, Filtration_value max_radius, const bool exact = false) : max_radius_(max_radius), exact_(exact) { - cech_skeleton_graph_ = Gudhi::compute_proximity_graph<SimplicialComplexForProximityGraph>( - point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); + point_cloud_.assign(std::begin(points), std::end(points)); + + cech_skeleton_graph_ = Gudhi::compute_proximity_graph<SimplicialComplexForCechComplex>( + point_cloud_, max_radius_, Sphere_circumradius<Kernel, Filtration_value>(exact)); } /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal @@ -85,7 +85,6 @@ class Cech_complex { * @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")); @@ -93,8 +92,7 @@ class Cech_complex { // 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)); + complex.expansion_with_blockers(dim_max, cech_blocker(&complex, this)); } /** @return max_radius value given at construction. */ @@ -103,12 +101,24 @@ class Cech_complex { /** @param[in] vertex Point position in the range. * @return The point. */ - const Point& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; } + const Point_d& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; } + + /** + * @return Vector of cached spheres. + */ + std::vector<Sphere> & get_cache() { return cache_; } + + /** \brief Check exact option + * @return Exact option. + */ + const bool is_exact() { return exact_; } private: Proximity_graph cech_skeleton_graph_; Filtration_value max_radius_; Point_cloud point_cloud_; + std::vector<Sphere> cache_; + const bool exact_; }; } // namespace cech_complex diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 31b9aab5..e78e37b7 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -11,10 +11,12 @@ #ifndef CECH_COMPLEX_BLOCKER_H_ #define CECH_COMPLEX_BLOCKER_H_ -#include <gudhi/distance_functions.h> // for Gudhi::Minimal_enclosing_ball_radius +#include <CGAL/NT_converter.h> // for casting from FT to Filtration_value +#include <CGAL/Lazy_exact_nt.h> // for CGAL::exact #include <iostream> #include <vector> +#include <set> #include <cmath> // for std::sqrt namespace Gudhi { @@ -30,37 +32,104 @@ namespace 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, + * \tparam SimplicialComplexForCech 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. + * \tparam Cech_complex is required by the blocker. + * + * \tparam Kernel CGAL kernel: either Epick_d or Epeck_d. */ -template <typename SimplicialComplexForCech, typename Cech_complex> +template <typename SimplicialComplexForCech, typename Cech_complex, typename Kernel> class Cech_blocker { + + public: + + using Point_d = typename Kernel::Point_d; + // Numeric type of coordinates in the kernel + using FT = typename Kernel::FT; + // Sphere is a pair of point and squared radius. + using Sphere = typename std::pair<Point_d, FT>; + private: - using Point_cloud = typename Cech_complex::Point_cloud; using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; using Filtration_value = typename SimplicialComplexForCech::Filtration_value; + using Simplex_key = typename SimplicialComplexForCech::Simplex_key; + + template<class PointIterator> + Sphere get_sphere(PointIterator begin, PointIterator end) const { + Point_d c = kernel_.construct_circumcenter_d_object()(begin, end); + FT r = kernel_.squared_distance_d_object()(c, *begin); + return std::make_pair(std::move(c), std::move(r)); + } 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) { + using Point_cloud = std::vector<Point_d>; + Filtration_value radius = 0; + bool is_min_enclos_ball = false; Point_cloud points; - for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { - points.push_back(cc_ptr_->get_point(vertex)); + points.reserve(sc_ptr_->dimension(sh)+1); + + // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices + for (auto face_opposite_vertex : sc_ptr_->boundary_opposite_vertex_simplex_range(sh)) { + auto k = sc_ptr_->key(face_opposite_vertex.first); + Simplex_key sph_key; + if(k != sc_ptr_->null_key()) { + sph_key = k; + } + else { + for (auto vertex : sc_ptr_->simplex_vertex_range(face_opposite_vertex.first)) { + points.push_back(cc_ptr_->get_point(vertex)); +#ifdef DEBUG_TRACES + std::clog << "#(" << vertex << ")#"; +#endif // DEBUG_TRACES + } + // Put edge sphere in cache + sph_key = cc_ptr_->get_cache().size(); + sc_ptr_->assign_key(face_opposite_vertex.first, sph_key); + cc_ptr_->get_cache().push_back(get_sphere(points.cbegin(), points.cend())); + // Clear face points + points.clear(); + } + // Check if the minimal enclosing ball of current face contains the extra point/opposite vertex + Sphere const& sph = cc_ptr_->get_cache()[sph_key]; + if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(face_opposite_vertex.second)) <= sph.second) { + is_min_enclos_ball = true; + sc_ptr_->assign_key(sh, sph_key); + radius = sc_ptr_->filtration(face_opposite_vertex.first); #ifdef DEBUG_TRACES - std::clog << "#(" << vertex << ")#"; + std::clog << "center: " << sph.first << ", radius: " << radius << std::endl; #endif // DEBUG_TRACES + break; + } } - Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); + // Spheres of each face don't contain the whole simplex + if(!is_min_enclos_ball) { + for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { + points.push_back(cc_ptr_->get_point(vertex)); + } + Sphere sph = get_sphere(points.cbegin(), points.cend()); +#if CGAL_VERSION_NR >= 1050000000 + if(cc_ptr_->is_exact()) CGAL::exact(sph.second); +#endif + CGAL::NT_converter<FT, Filtration_value> cast_to_fv; + radius = std::sqrt(cast_to_fv(sph.second)); + + sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); + cc_ptr_->get_cache().push_back(std::move(sph)); + } + #ifdef DEBUG_TRACES if (radius > cc_ptr_->max_radius()) std::clog << "radius > max_radius => expansion is blocked\n"; #endif // DEBUG_TRACES - sc_ptr_->assign_filtration(sh, radius); + // Check that the filtration to be assigned (radius) would be valid + if (radius > sc_ptr_->filtration(sh)) sc_ptr_->assign_filtration(sh, radius); return (radius > cc_ptr_->max_radius()); } @@ -70,6 +139,7 @@ class Cech_blocker { private: SimplicialComplexForCech* sc_ptr_; Cech_complex* cc_ptr_; + Kernel kernel_; }; } // namespace cech_complex diff --git a/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT b/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT deleted file mode 100644 index dbe4c553..00000000 --- a/src/Cech_complex/include/gudhi/Miniball.COPYRIGHT +++ /dev/null @@ -1,4 +0,0 @@ -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 deleted file mode 100644 index 033d8953..00000000 --- a/src/Cech_complex/include/gudhi/Miniball.README +++ /dev/null @@ -1,26 +0,0 @@ -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 deleted file mode 100644 index ce6cbb5b..00000000 --- a/src/Cech_complex/include/gudhi/Miniball.hpp +++ /dev/null @@ -1,523 +0,0 @@ -// 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/include/gudhi/Sphere_circumradius.h b/src/Cech_complex/include/gudhi/Sphere_circumradius.h new file mode 100644 index 00000000..2f916c0a --- /dev/null +++ b/src/Cech_complex/include/gudhi/Sphere_circumradius.h @@ -0,0 +1,78 @@ +/* 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): Hind Montassif + * + * Copyright (C) 2021 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SPHERE_CIRCUMRADIUS_H_ +#define SPHERE_CIRCUMRADIUS_H_ + +#include <CGAL/Epick_d.h> // for #include <CGAL/NT_converter.h> which is not working/compiling alone +#include <CGAL/Lazy_exact_nt.h> // for CGAL::exact + +#include <cmath> // for std::sqrt +#include <vector> + +namespace Gudhi { + +namespace cech_complex { + +/** \private @brief Compute the circumradius of the sphere passing through points given by a range of coordinates. + * The points are assumed to have the same dimension. */ +template<typename Kernel, typename Filtration_value> +class Sphere_circumradius { + private: + Kernel kernel_; + const bool exact_; + public: + using FT = typename Kernel::FT; + using Point = typename Kernel::Point_d; + using Point_cloud = typename std::vector<Point>; + + CGAL::NT_converter<FT, Filtration_value> cast_to_fv; + + /** \brief Circumradius of sphere passing through two points using CGAL. + * + * @param[in] point_1 + * @param[in] point_2 + * @return Sphere circumradius passing through two points. + * \tparam Point must be a Kernel::Point_d from CGAL. + * + */ + Filtration_value operator()(const Point& point_1, const Point& point_2) const { + auto squared_dist_obj = kernel_.squared_distance_d_object()(point_1, point_2); + if(exact_) CGAL::exact(squared_dist_obj); + return std::sqrt(cast_to_fv(squared_dist_obj)) / 2.; + } + + /** \brief Circumradius of sphere passing through point cloud using CGAL. + * + * @param[in] point_cloud The points. + * @return Sphere circumradius passing through the points. + * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. + * + */ + Filtration_value operator()(const Point_cloud& point_cloud) const { + auto squared_radius_obj = kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()); + if(exact_) CGAL::exact(squared_radius_obj); + return std::sqrt(cast_to_fv(squared_radius_obj)); + } + + /** \brief Constructor + * @param[in] exact Option for exact filtration values computation. Not exact if `Kernel` is not <a target="_blank" + * href="https://doc.cgal.org/latest/Kernel_d/structCGAL_1_1Epeck__d.html">CGAL::Epeck_d</a>. + * Default is false. + */ + Sphere_circumradius(const bool exact = false) : exact_(exact) {} + +}; + +} // namespace cech_complex + +} // namespace Gudhi + +#endif // SPHERE_CIRCUMRADIUS_H_ |