Loading...
Searching...
No Matches
MEB_filtration.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): Marc Glisse
4 *
5 * Copyright (C) 2023 Inria
6 *
7 * Modification(s):
8 * - YYYY/MM Author: Description of the modification
9 */
10
11#ifndef MEB_FILTRATION_H_
12#define MEB_FILTRATION_H_
13
14namespace Gudhi::cech_complex {
15
35template<typename Kernel, typename SimplicialComplexForMEB, typename PointRange>
36void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRange points, bool exact = false) {
37 using Point_d = typename Kernel::Point_d;
38 using FT = typename Kernel::FT;
39 using Sphere = std::pair<Point_d, FT>;
40
42 using Simplex_handle = typename SimplicialComplexForMEB::Simplex_handle;
44
45 std::vector<Sphere> cache_;
46 std::vector<Point_d> pts;
47 CGAL::NT_converter<FT, Filtration_value> cvt;
48
49 auto fun = [&](Simplex_handle sh, int dim){
50 if (dim == 0) complex.assign_filtration(sh, 0);
51 else if (dim == 1) {
52 // For a Simplex_tree, this would be a bit faster, but that's probably negligible
53 // Vertex_handle u = sh->first; Vertex_handle v = self_siblings(sh)->parent();
54 auto verts = complex.simplex_vertex_range(sh);
55 auto vert_it = verts.begin();
56 Vertex_handle u = *vert_it;
57 Vertex_handle v = *++vert_it;
58 auto&& pu = points[u];
59 Point_d m = k.midpoint_d_object()(pu, points[v]);
60 FT r = k.squared_distance_d_object()(m, pu);
61 if (exact) CGAL::exact(r);
62 complex.assign_key(sh, cache_.size());
63 complex.assign_filtration(sh, std::max(cvt(r), Filtration_value(0)));
64 cache_.emplace_back(std::move(m), std::move(r));
65 } else {
66 Filtration_value maxf = 0; // max filtration of the faces
67 bool found = false;
68 using std::max;
69 for (auto face_opposite_vertex : complex.boundary_opposite_vertex_simplex_range(sh)) {
70 maxf = max(maxf, complex.filtration(face_opposite_vertex.first));
71 if (!found) {
72 auto key = complex.key(face_opposite_vertex.first);
73 Sphere const& sph = cache_[key];
74 if (k.squared_distance_d_object()(sph.first, points[face_opposite_vertex.second]) > sph.second) continue;
75 found = true;
76 complex.assign_key(sh, key);
77 // With exact computations, we could stop here
78 // complex.assign_filtration(sh, complex.filtration(face_opposite_vertex.first)); return;
79 // but because of possible rounding errors, we continue with the equivalent of make_filtration_non_decreasing
80 }
81 }
82 if (!found) {
83 // None of the faces are good enough, MEB must be the circumsphere.
84 pts.clear();
85 for (auto vertex : complex.simplex_vertex_range(sh))
86 pts.push_back(points[vertex]);
87 Point_d c = k.construct_circumcenter_d_object()(pts.begin(), pts.end());
88 FT r = k.squared_distance_d_object()(c, pts.front());
89 if (exact) CGAL::exact(r);
90 maxf = max(maxf, cvt(r)); // maxf = cvt(r) except for rounding errors
91 complex.assign_key(sh, cache_.size());
92 // We could check if the simplex is maximal and avoiding adding it to the cache in that case.
93 cache_.emplace_back(std::move(c), std::move(r));
94 }
95 complex.assign_filtration(sh, maxf);
96 }
97 };
98 complex.for_each_simplex(fun);
99
100 // We could avoid computing maxf, but when !exact rounding errors may cause
101 // the filtration values to be non-monotonous, so we would need to call
102 // if (!exact) complex.make_filtration_non_decreasing();
103 // which is way more costly than computing maxf. The exact case is already so
104 // costly that it isn't worth maintaining code without maxf just for it.
105 // Cech_complex has "free" access to the max of the faces, because
106 // expansion_with_blockers computes it before the callback.
107
108 // TODO: use a map if complex does not provide key?
109}
110} // namespace Gudhi::cech_complex
111
112#endif // MEB_FILTRATION_H_
void assign_MEB_filtration(Kernel &&k, SimplicialComplexForMEB &complex, PointRange points, bool exact=false)
Given a simplicial complex and an embedding of its vertices, this assigns to each simplex a filtratio...
Definition MEB_filtration.h:36
Value type for a filtration function on a cell complex.
Definition FiltrationValue.h:20
Definition SimplicialComplexForMEB.h:22
unspecified Simplex_handle
Handle for a simplex.
Definition SimplicialComplexForMEB.h:24
Simplex_key key(Simplex_handle simplex)
Returns the key assigned to the 'simplex' with assign_key().
Filtration_value filtration(Simplex_handle simplex)
Returns the filtration value to the 'simplex'.
Boundary_opposite_vertex_simplex_range boundary_opposite_vertex_simplex_range(Simplex_handle simplex)
Returns a range of the pairs (simplex, opposite vertex) of the boundary of the 'simplex'.
unspecified Filtration_value
Type of filtration values.
Definition SimplicialComplexForMEB.h:29
void assign_key(Simplex_handle simplex, Simplex_key key)
Assigns this 'key' to the 'simplex'.
Simplex_vertex_range simplex_vertex_range(Simplex_handle simplex)
Returns a range over vertices (as Vertex_handle) of a given simplex.
unspecified Vertex_handle
Handle for a vertex. Must be a non-negative integer, it is also used as an index into the input list ...
Definition SimplicialComplexForMEB.h:27
int assign_filtration(Simplex_handle simplex, Filtration_value filtration)
Assigns this 'filtration' value to the 'simplex'.
void for_each_simplex(auto callback)
Calls callback(simplex, dim) for every simplex of the complex, with the guarantee that faces are visi...
Handle type for the vertices of a cell complex.
Definition VertexHandle.h:15