Loading...
Searching...
No Matches
MEB_filtration.h
/* 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): Marc Glisse
*
* Copyright (C) 2023 Inria
*
* Modification(s):
* - 2024/10 Vincent Rouvreau: Add Output_squared_values argument to enable/disable squared radii computation
* - YYYY/MM Author: Description of the modification
*/
#ifndef MEB_FILTRATION_H_
#define MEB_FILTRATION_H_
#include <CGAL/NT_converter.h>
#include <vector>
#include <utility> // for std::pair
#include <cmath> // for std::sqrt
namespace Gudhi::cech_complex {
template<bool Output_squared_values = true, typename Kernel, typename SimplicialComplexForMEB, typename PointRange>
void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRange const& points, bool exact = false) {
using Point_d = typename Kernel::Point_d;
using FT = typename Kernel::FT;
using Sphere = std::pair<Point_d, FT>;
using Vertex_handle = typename SimplicialComplexForMEB::Vertex_handle;
using Simplex_handle = typename SimplicialComplexForMEB::Simplex_handle;
using Filtration_value = typename SimplicialComplexForMEB::Filtration_value;
// For users to be able to define their own sqrt function on their desired Filtration_value type
using std::sqrt;
std::vector<Sphere> cache_;
std::vector<Point_d> pts;
CGAL::NT_converter<FT, Filtration_value> cvt;
// This block is only needed to get ambient_dim
if(std::begin(points) == std::end(points)) {
// assert(complex.is_empty());
return;
}
int ambient_dim = k.point_dimension_d_object()(*std::begin(points));
auto fun = [&](Simplex_handle sh, int dim){
using std::max;
if (dim == 0) complex.assign_filtration(sh, 0);
else if (dim == 1) {
// For a Simplex_tree, this would be a bit faster, but that's probably negligible
// Vertex_handle u = sh->first; Vertex_handle v = self_siblings(sh)->parent();
auto verts = complex.simplex_vertex_range(sh);
auto vert_it = verts.begin();
Vertex_handle u = *vert_it;
Vertex_handle v = *++vert_it;
auto&& pu = points[u];
Point_d m = k.midpoint_d_object()(pu, points[v]);
FT r = k.squared_distance_d_object()(m, pu);
if (exact) CGAL::exact(r);
complex.assign_key(sh, cache_.size());
Filtration_value filt{max(cvt(r), Filtration_value(0))};
if constexpr (!Output_squared_values)
filt = sqrt(filt);
complex.assign_filtration(sh, filt);
cache_.emplace_back(std::move(m), std::move(r));
} else if (dim > ambient_dim) {
// The sphere is always defined by at most d+1 points
Filtration_value maxf = 0; // max filtration of the faces
for (auto face : complex.boundary_simplex_range(sh)) {
maxf = max(maxf, complex.filtration(face));
}
complex.assign_filtration(sh, maxf);
} else {
Filtration_value maxf = 0; // max filtration of the faces
bool found = false;
for (auto face_opposite_vertex : complex.boundary_opposite_vertex_simplex_range(sh)) {
maxf = max(maxf, complex.filtration(face_opposite_vertex.first));
if (!found) {
auto key = complex.key(face_opposite_vertex.first);
Sphere const& sph = cache_[key];
if (k.squared_distance_d_object()(sph.first, points[face_opposite_vertex.second]) > sph.second) continue;
found = true;
complex.assign_key(sh, key);
// With exact computations, we could stop here
// complex.assign_filtration(sh, complex.filtration(face_opposite_vertex.first)); return;
// but because of possible rounding errors, we continue with the equivalent of make_filtration_non_decreasing
}
}
if (!found) {
// None of the faces are good enough, MEB must be the circumsphere.
pts.clear();
for (auto vertex : complex.simplex_vertex_range(sh))
pts.push_back(points[vertex]);
Point_d c = k.construct_circumcenter_d_object()(pts.begin(), pts.end());
FT r = k.squared_distance_d_object()(c, pts.front());
if (exact) CGAL::exact(r);
// For Epick_d, if the circumcenter computation is too unstable, we could compute
// int d2 = dim * dim;
// Filtration_value max_sanity = maxf * d2 / (d2 - 1);
// and use min(max_sanity, ...), which would limit how bad numerical errors can be.
Filtration_value filt{cvt(r)};
if constexpr (!Output_squared_values)
filt = sqrt(max(filt, Filtration_value(0)));
// maxf = filt except for rounding errors
maxf = max(maxf, filt);
complex.assign_key(sh, cache_.size());
// We could check if the simplex is maximal and avoiding adding it to the cache in that case.
cache_.emplace_back(std::move(c), std::move(r));
}
complex.assign_filtration(sh, maxf);
}
};
complex.for_each_simplex(fun);
// We could avoid computing maxf, but when !exact rounding errors may cause
// the filtration values to be non-monotonous, so we would need to call
// if (!exact) complex.make_filtration_non_decreasing();
// which is way more costly than computing maxf. The exact case is already so
// costly that it isn't worth maintaining code without maxf just for it.
// Cech_complex has "free" access to the max of the faces, because
// expansion_with_blockers computes it before the callback.
// TODO: use a map if complex does not provide key?
}
} // namespace Gudhi::cech_complex
#endif // MEB_FILTRATION_H_
Gudhi::cech_complex::assign_MEB_filtration
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:46
Gudhi::cech_complex::SimplicialComplexForMEB
Definition SimplicialComplexForMEB.h:22
Gudhi::cech_complex::SimplicialComplexForMEB::Simplex_handle
unspecified Simplex_handle
Handle for a simplex.
Definition SimplicialComplexForMEB.h:24
Gudhi::cech_complex::SimplicialComplexForMEB::key
Simplex_key key(Simplex_handle simplex)
Returns the key assigned to the 'simplex' with assign_key().
Gudhi::cech_complex::SimplicialComplexForMEB::filtration
Filtration_value filtration(Simplex_handle simplex)
Returns the filtration value to the 'simplex'.
Gudhi::cech_complex::SimplicialComplexForMEB::boundary_opposite_vertex_simplex_range
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'.
Gudhi::cech_complex::SimplicialComplexForMEB::Filtration_value
unspecified Filtration_value
Type of filtration values.
Definition SimplicialComplexForMEB.h:29
Gudhi::cech_complex::SimplicialComplexForMEB::assign_key
void assign_key(Simplex_handle simplex, Simplex_key key)
Assigns this 'key' to the 'simplex'.
Gudhi::cech_complex::SimplicialComplexForMEB::simplex_vertex_range
Simplex_vertex_range simplex_vertex_range(Simplex_handle simplex)
Returns a range over vertices (as Vertex_handle) of a given simplex.
Gudhi::cech_complex::SimplicialComplexForMEB::Vertex_handle
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
Gudhi::cech_complex::SimplicialComplexForMEB::assign_filtration
int assign_filtration(Simplex_handle simplex, Filtration_value filtration)
Assigns this 'filtration' value to the 'simplex'.
Gudhi::cech_complex::SimplicialComplexForMEB::for_each_simplex
void for_each_simplex(auto callback)
Calls callback(simplex, dim) for every simplex of the complex, with the guarantee that faces are visi...