Unsafe Cech with Epick_d?

[mglisse]

I do not have an example at hand, but I think the computation done in Cech_complex and assign_MEB_filtration may be unsafe when the kernel is Epick_d. In the plane, with 4 points that are cocyclic (or almost), rounding errors may tell us that no point is in the MEB of the other 3. In this case, we will ask CGAL for the circumcenter of 4 points, which in the best case will give an assertion failure. If this comment is true, we should either

1. make this more robust with Epick_d. For the specific issue described above, we could just check if we have exceeded the ambient dimension, and in that case use the max of the faces no matter what the inclusion test might say (it could even count as an optimization, although I don't know how often this function will be used on simplices of dimension higher than ambient). But we need to think what other issues could happen.
2. document how risky Epick_d is for these constructions (unlike for Alpha_complex).

PoC (doesn't handle the empty case, probably not the best way to access the dimension)

@@ -45,6 +45,7 @@ void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRan
   std::vector<Sphere> cache_;
   std::vector<Point_d> pts;
   CGAL::NT_converter<FT, Filtration_value> cvt;
+  int ambient_dim = std::begin(points)->dimension();

   auto fun = [&](Simplex_handle sh, int dim){
     if (dim == 0) complex.assign_filtration(sh, 0);
@@ -62,6 +63,13 @@ void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRan
       complex.assign_key(sh, cache_.size());
       complex.assign_filtration(sh, std::max(cvt(r), Filtration_value(0)));
       cache_.emplace_back(std::move(m), std::move(r));
+    } else if (dim > ambient_dim) {
+      Filtration_value maxf = 0; // max filtration of the faces
+      for (auto face : complex.boundary_simplex_range(sh)) {
+        using std::max;
+        maxf = max(maxf, complex.filtration(face));
+      }
+      complex.assign_filtration(sh, maxf);
     } else {
       Filtration_value maxf = 0; // max filtration of the faces
       bool found = false;

While I am here, there is a property that the circumradius cannot exceed the max of the faces by a factor larger than $1/\sqrt{1-1/d^2}$, we could use it as an upper bound to ensure that the values we return never explode:

@@ -87,7 +95,12 @@ void assign_MEB_filtration(Kernel&&k, SimplicialComplexForMEB& complex, PointRan
         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);
+        int d2 = dim * dim;
+        Filtration_value max_sanity = maxf * d2 / (d2 - 1);
+        using std::min;
+        using std::max;
         maxf = max(maxf, cvt(r)); // maxf = cvt(r) except for rounding errors
+        maxf = min(maxf, max_sanity);
         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));

(of course it would only be done for Epick_d, not Epeck_d) I don't know if it is interesting though.