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 const& 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 // This block is only needed to get ambient_dim
50 if(std::begin(points) == std::end(points)) {
51 // assert(complex.is_empty());
52 return;
53 }
54 int ambient_dim = k.point_dimension_d_object()(*std::begin(points));
55
56 auto fun = [&](Simplex_handle sh, int dim){
57 using std::max;
58 if (dim == 0) complex.assign_filtration(sh, 0);
59 else if (dim == 1) {
60 // For a Simplex_tree, this would be a bit faster, but that's probably negligible
61 // Vertex_handle u = sh->first; Vertex_handle v = self_siblings(sh)->parent();
62 auto verts = complex.simplex_vertex_range(sh);
63 auto vert_it = verts.begin();
64 Vertex_handle u = *vert_it;
65 Vertex_handle v = *++vert_it;
66 auto&& pu = points[u];
67 Point_d m = k.midpoint_d_object()(pu, points[v]);
68 FT r = k.squared_distance_d_object()(m, pu);
69 if (exact) CGAL::exact(r);
70 complex.assign_key(sh, cache_.size());
71 complex.assign_filtration(sh, max(cvt(r), Filtration_value(0)));
72 cache_.emplace_back(std::move(m), std::move(r));
73 } else if (dim > ambient_dim) {
74 // The sphere is always defined by at most d+1 points
75 Filtration_value maxf = 0; // max filtration of the faces
76 for (auto face : complex.boundary_simplex_range(sh)) {
77 maxf = max(maxf, complex.filtration(face));
78 }
79 complex.assign_filtration(sh, maxf);
80 } else {
81 Filtration_value maxf = 0; // max filtration of the faces
82 bool found = false;
83 for (auto face_opposite_vertex : complex.boundary_opposite_vertex_simplex_range(sh)) {
84 maxf = max(maxf, complex.filtration(face_opposite_vertex.first));
85 if (!found) {
86 auto key = complex.key(face_opposite_vertex.first);
87 Sphere const& sph = cache_[key];
88 if (k.squared_distance_d_object()(sph.first, points[face_opposite_vertex.second]) > sph.second) continue;
89 found = true;
90 complex.assign_key(sh, key);
91 // With exact computations, we could stop here
92 // complex.assign_filtration(sh, complex.filtration(face_opposite_vertex.first)); return;
93 // but because of possible rounding errors, we continue with the equivalent of make_filtration_non_decreasing
94 }
95 }
96 if (!found) {
97 // None of the faces are good enough, MEB must be the circumsphere.
98 pts.clear();
99 for (auto vertex : complex.simplex_vertex_range(sh))
100 pts.push_back(points[vertex]);
101 Point_d c = k.construct_circumcenter_d_object()(pts.begin(), pts.end());
102 FT r = k.squared_distance_d_object()(c, pts.front());
103 if (exact) CGAL::exact(r);
104 // For Epick_d, if the circumcenter computation is too unstable, we could compute
105 // int d2 = dim * dim;
106 // Filtration_value max_sanity = maxf * d2 / (d2 - 1);
107 // and use min(max_sanity, ...), which would limit how bad numerical errors can be.
108 maxf = max(maxf, cvt(r)); // maxf = cvt(r) except for rounding errors
109 complex.assign_key(sh, cache_.size());
110 // We could check if the simplex is maximal and avoiding adding it to the cache in that case.
111 cache_.emplace_back(std::move(c), std::move(r));
112 }
113 complex.assign_filtration(sh, maxf);
114 }
115 };
116 complex.for_each_simplex(fun);
117
118 // We could avoid computing maxf, but when !exact rounding errors may cause
119 // the filtration values to be non-monotonous, so we would need to call
120 // if (!exact) complex.make_filtration_non_decreasing();
121 // which is way more costly than computing maxf. The exact case is already so
122 // costly that it isn't worth maintaining code without maxf just for it.
123 // Cech_complex has "free" access to the max of the faces, because
124 // expansion_with_blockers computes it before the callback.
125
126 // TODO: use a map if complex does not provide key?
127}
128} // namespace Gudhi::cech_complex
129
130#endif // MEB_FILTRATION_H_
void assign_MEB_filtration(Kernel &&k, SimplicialComplexForMEB &complex, PointRange const &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