Cech_complex_blocker.h
1 /* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2  * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3  * Author(s): Vincent Rouvreau
4  *
5  * Copyright (C) 2018 Inria
6  *
7  * Modification(s):
8  * - YYYY/MM Author: Description of the modification
9  */
10 
11 #ifndef CECH_COMPLEX_BLOCKER_H_
12 #define CECH_COMPLEX_BLOCKER_H_
13 
14 #include <CGAL/NT_converter.h> // for casting from FT to Filtration_value
15 #include <CGAL/Lazy_exact_nt.h> // for CGAL::exact
16 
17 #include <iostream>
18 #include <vector>
19 #include <set>
20 #include <cmath> // for std::sqrt
21 
22 namespace Gudhi {
23 
24 namespace cech_complex {
25 
42 template <typename SimplicialComplexForCech, typename Cech_complex, typename Kernel>
43 class Cech_blocker {
44 
45  public:
46 
47  using Point_d = typename Kernel::Point_d;
48  // Numeric type of coordinates in the kernel
49  using FT = typename Kernel::FT;
50  // Sphere is a pair of point and squared radius.
51  using Sphere = typename std::pair<Point_d, FT>;
52 
53  private:
54 
55  using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle;
57  using Simplex_key = typename SimplicialComplexForCech::Simplex_key;
58 
59  template<class PointIterator>
60  Sphere get_sphere(PointIterator begin, PointIterator end) const {
61  Point_d c = kernel_.construct_circumcenter_d_object()(begin, end);
62  FT r = kernel_.squared_distance_d_object()(c, *begin);
63  return std::make_pair(std::move(c), std::move(r));
64  }
65 
66  public:
67 
72  bool operator()(Simplex_handle sh) {
73  using Point_cloud = std::vector<Point_d>;
74  Filtration_value radius = 0;
75  bool is_min_enclos_ball = false;
76  Point_cloud points;
77  points.reserve(sc_ptr_->dimension(sh)+1);
78 
79  // 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
80  for (auto face_opposite_vertex : sc_ptr_->boundary_opposite_vertex_simplex_range(sh)) {
81  auto k = sc_ptr_->key(face_opposite_vertex.first);
82  Simplex_key sph_key;
83  if(k != sc_ptr_->null_key()) {
84  sph_key = k;
85  }
86  else {
87  for (auto vertex : sc_ptr_->simplex_vertex_range(face_opposite_vertex.first)) {
88  points.push_back(cc_ptr_->get_point(vertex));
89 #ifdef DEBUG_TRACES
90  std::clog << "#(" << vertex << ")#";
91 #endif // DEBUG_TRACES
92  }
93  // Put edge sphere in cache
94  sph_key = cc_ptr_->get_cache().size();
95  sc_ptr_->assign_key(face_opposite_vertex.first, sph_key);
96  cc_ptr_->get_cache().push_back(get_sphere(points.cbegin(), points.cend()));
97  // Clear face points
98  points.clear();
99  }
100  // Check if the minimal enclosing ball of current face contains the extra point/opposite vertex
101  Sphere const& sph = cc_ptr_->get_cache()[sph_key];
102  if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(face_opposite_vertex.second)) <= sph.second) {
103  is_min_enclos_ball = true;
104  sc_ptr_->assign_key(sh, sph_key);
105  radius = sc_ptr_->filtration(face_opposite_vertex.first);
106 #ifdef DEBUG_TRACES
107  std::clog << "center: " << sph.first << ", radius: " << radius << std::endl;
108 #endif // DEBUG_TRACES
109  break;
110  }
111  }
112  // Spheres of each face don't contain the whole simplex
113  if(!is_min_enclos_ball) {
114  for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) {
115  points.push_back(cc_ptr_->get_point(vertex));
116  }
117  Sphere sph = get_sphere(points.cbegin(), points.cend());
118 #if CGAL_VERSION_NR >= 1050000000
119  if(cc_ptr_->is_exact()) CGAL::exact(sph.second);
120 #endif
121  CGAL::NT_converter<FT, Filtration_value> cast_to_fv;
122  radius = std::sqrt(cast_to_fv(sph.second));
123 
124  sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size());
125  cc_ptr_->get_cache().push_back(std::move(sph));
126  }
127 
128 #ifdef DEBUG_TRACES
129  if (radius > cc_ptr_->max_radius()) std::clog << "radius > max_radius => expansion is blocked\n";
130 #endif // DEBUG_TRACES
131  // Check that the filtration to be assigned (radius) would be valid
132  if (radius > sc_ptr_->filtration(sh)) sc_ptr_->assign_filtration(sh, radius);
133  return (radius > cc_ptr_->max_radius());
134  }
135 
137  Cech_blocker(SimplicialComplexForCech* sc_ptr, Cech_complex* cc_ptr) : sc_ptr_(sc_ptr), cc_ptr_(cc_ptr) {}
138 
139  private:
140  SimplicialComplexForCech* sc_ptr_;
141  Cech_complex* cc_ptr_;
142  Kernel kernel_;
143 };
144 
145 } // namespace cech_complex
146 
147 } // namespace Gudhi
148 
149 #endif // CECH_COMPLEX_BLOCKER_H_
Cech complex class.
Definition: Cech_complex.h:43
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:20
unspecified Filtration_value
Definition: SimplicialComplexForCech.h:27
unspecified Simplex_handle
Definition: SimplicialComplexForCech.h:23