infinite PH0 pair

Author: MClemot

core/base/delaunayRipsPersistenceDiagram/DelaunayRipsPersistenceDiagram.cpp

@@ -37,12 +37,6 @@ int ttk::DelaunayRipsPersistenceDiagram::execute(
     gph::tryDimensions(points, ph, getThreadNumber());
     ph[0].emplace_back(
       FiltratedSimplex{{-1}, 0.}, FiltratedSimplex{{-1}, inf}); // infinite pair
-    for(auto &diag : ph) {
-      for(auto &[b, d] : diag) {
-        b.first = {-1};
-        d.first = {-1};
-      }
-    }
   }
   return 0;
 #else

core/base/delaunayRipsPersistenceDiagram/geoPHd.h

@@ -249,7 +249,7 @@ namespace ttk::gph {
     }
     for(FiltratedDSimplex<DYN_DIM> const &e : MSA)
       ph[0].emplace_back(
-        FiltratedSimplex{{}, 0.}, FiltratedSimplex{e.s, sqrt(e.d)});
+        FiltratedSimplex{{-1}, 0.}, FiltratedSimplex{e.s, sqrt(e.d)});
   }
 
   template <unsigned DIM>
@@ -265,7 +265,7 @@ namespace ttk::gph {
       recurse(ph, nextMSA);
     } else if constexpr(D == 1) { // this is the 0-dimensional homology
       for(FiltratedDSimplex<1> const &e : MSA)
-        ph[0].emplace_back(FiltratedSimplex{{}, 0.},
+        ph[0].emplace_back(FiltratedSimplex{{-1}, 0.},
                            FiltratedSimplex{{e.s[0], e.s[1]}, sqrt(e.d)});
     }
   }

core/vtk/ttkRipsPersistenceDiagram/ttkRipsPersistenceDiagram.cpp

@@ -20,8 +20,17 @@ void DiagramToVTU(vtkUnstructuredGrid *vtu,
   int n_pairs = 0;
   for(auto const &diagram_d : diagram)
     n_pairs += diagram_d.size();
-  if(SimplexMaximumDiameter == ttk::rpd::inf)
-    n_pairs--;
+
+  if(SimplexMaximumDiameter == ttk::rpd::inf) {
+    double maxFiniteValue = 0.;
+    for(auto const &diag : diagram) {
+      for(auto const &[b, d] : diag) {
+        if(d.second < ttk::rpd::inf)
+          maxFiniteValue = std::max(maxFiniteValue, d.second);
+      }
+    }
+    SimplexMaximumDiameter = 1.5 * maxFiniteValue;
+  }
 
   // point data arrays
   vtkNew<ttkSimplexIdTypeArray> vertsId{};
@@ -74,10 +83,6 @@ void DiagramToVTU(vtkUnstructuredGrid *vtu,
   double birth_max = 0.;
   for(unsigned d = 0; d < diagram.size(); ++d) {
     for(auto const &pair : diagram[d]) {
-      if(d == 0 && pair.second.second == ttk::rpd::inf
-         && SimplexMaximumDiameter == ttk::rpd::inf)
-        continue;
-
       const unsigned i0 = 2 * i, i1 = 2 * i + 1;
       pairsId->SetTuple1(i, i);
       pairsDim->SetTuple1(i, d);