[MergeTree] saddle-saddle pairs processing

core/vtk/ttkMergeTree/ttkMergeTreeUtils.h

@@ -87,11 +87,29 @@ namespace ttk {
       return mergeTree;
     }
 
-    // Returns a branch decomposition tree
+    /**
+     * @brief Create a MergeTree (as a branch decomposition tree) object given a
+     * vtkUnstructuredGrid representing a persistence diagram (in the TTK
+     * format).
+     *
+     * @param[in] persistenceDiagram vtk object representing the persistence
+     * diagram.
+     * @param[in] useSecondPairsType boolean to choose which pair type to use
+     * depending on the value of diagramPairTypes, by default the choice is
+     * between min-saddle and saddle-max, true for the second type and false for
+     * the first type.
+     * @param[in] diagramPairTypes 0 to choose between min-saddle and
+     * saddle-max, 1 to choose between min-saddle and saddle-saddle, 2 to choose
+     * between saddle-saddle and saddle-max.
+     *
+     * @return a MergeTree object corresponding to the input persistence
+     * diagram.
+     */
     template <class dataType>
     MergeTree<dataType>
       makeBDTreeFromPDGrid(vtkUnstructuredGrid *persistenceDiagram,
-                           bool useSadMaxPairs = true) {
+                           bool useSecondPairsType = true,
+                           int diagramPairTypes = 0) {
       auto birthArray
         = persistenceDiagram->GetCellData()->GetArray(PersistenceBirthName);
       auto persArray
@@ -126,6 +144,18 @@ namespace ttk {
       // Init critical type enum values
       auto locMin = static_cast<int>(CriticalType::Local_minimum);
       auto locMax = static_cast<int>(CriticalType::Local_maximum);
+      auto locSad1 = static_cast<int>(CriticalType::Saddle1);
+      auto locSad2 = static_cast<int>(CriticalType::Saddle2);
+
+      // 0 : min-saddle ; 1 : saddle-saddle : 2 : saddle-max
+      int pairsType
+        = (useSecondPairsType
+               and (diagramPairTypes == 0 or diagramPairTypes == 2)
+             ? 2
+             : (not useSecondPairsType
+                    and (diagramPairTypes == 0 or diagramPairTypes == 1)
+                  ? 0
+                  : 1));
 
       // Get Min-Max pair index
       int minMaxPairIndex = -1;
@@ -151,8 +181,11 @@ namespace ttk {
         auto ct1 = criticalTypeArray->GetTuple1(pts[0]);
         auto ct2 = criticalTypeArray->GetTuple1(pts[1]);
         if((pairType == -1
-            or (useSadMaxPairs and ct1 != locMax and ct2 != locMax)
-            or (not useSadMaxPairs and ct1 != locMin and ct2 != locMin))
+            or (pairsType == 2 and ct1 != locMax and ct2 != locMax)
+            or (pairsType == 1
+                and not((ct1 == locSad1 and ct2 == locSad2)
+                        or (ct1 == locSad2 and ct2 == locSad1)))
+            or (pairsType == 0 and ct1 != locMin and ct2 != locMin))
            and i != minMaxPairIndex)
           continue;
         int const index1
@@ -218,7 +251,8 @@ namespace ttk {
       std::vector<vtkUnstructuredGrid *> &treesNodes,
       std::vector<vtkUnstructuredGrid *> &treesArcs,
       std::vector<vtkDataSet *> &treesSegmentation,
-      std::vector<bool> useSadMaxPairs) {
+      const std::vector<bool> &useSecondPairsTypeVec,
+      int diagramPairTypes = 0) {
       bool isPersistenceDiagram = false;
       const int numInputs = inputTrees.size();
       intermediateTrees.resize(numInputs);
@@ -242,7 +276,7 @@ namespace ttk {
           vtkUnstructuredGrid *persistenceDiagram
             = vtkUnstructuredGrid::SafeDownCast(inputTrees[i]->GetBlock(0));
           intermediateTrees[i] = makeBDTreeFromPDGrid<dataType>(
-            persistenceDiagram, useSadMaxPairs[i]);
+            persistenceDiagram, useSecondPairsTypeVec[i], diagramPairTypes);
           isPersistenceDiagram = true;
         }
       }
@@ -256,23 +290,27 @@ namespace ttk {
       std::vector<vtkUnstructuredGrid *> &treesNodes,
       std::vector<vtkUnstructuredGrid *> &treesArcs,
       std::vector<vtkDataSet *> &treesSegmentation,
-      bool useSadMaxPairs = true) {
-      std::vector<bool> const useSadMaxPairsVec(
-        inputTrees.size(), useSadMaxPairs);
+      bool useSecondPairsType = true,
+      int diagramPairTypes = 0) {
+      const std::vector<bool> useSecondPairsTypeVec(
+        inputTrees.size(), useSecondPairsType);
       return constructTrees(inputTrees, intermediateTrees, treesNodes,
-                            treesArcs, treesSegmentation, useSadMaxPairsVec);
+                            treesArcs, treesSegmentation, useSecondPairsTypeVec,
+                            diagramPairTypes);
     }
 
     template <class dataType>
     bool constructTrees(
       std::vector<vtkSmartPointer<vtkMultiBlockDataSet>> &inputTrees,
       std::vector<MergeTree<dataType>> &intermediateTrees,
-      bool useSadMaxPairs = true) {
+      bool useSecondPairsType = true,
+      int diagramPairTypes = 0) {
       std::vector<vtkUnstructuredGrid *> treesNodes;
       std::vector<vtkUnstructuredGrid *> treesArcs;
       std::vector<vtkDataSet *> treesSegmentation;
       return constructTrees(inputTrees, intermediateTrees, treesNodes,
-                            treesArcs, treesSegmentation, useSadMaxPairs);
+                            treesArcs, treesSegmentation, useSecondPairsType,
+                            diagramPairTypes);
     }
   } // namespace ftm
 } // namespace ttk

core/vtk/ttkMergeTreeAutoencoder/ttkMergeTreeAutoencoder.cpp

@@ -187,10 +187,10 @@ int ttkMergeTreeAutoencoder::runCompute(
   std::vector<ttk::ftm::MergeTree<float>> intermediateMTrees,
     intermediateMTrees2;
 
-  bool useSadMaxPairs = (mixtureCoefficient_ == 0);
+  bool useSecondPairsType = (mixtureCoefficient_ == 0);
   isPersistenceDiagram_ = ttk::ftm::constructTrees<float>(
     inputTrees, intermediateMTrees, treesNodes, treesArcs, treesSegmentation,
-    useSadMaxPairs);
+    useSecondPairsType, DiagramPairTypes);
   // If merge trees are provided in input and normalization is not asked
   convertToDiagram_
     = (not isPersistenceDiagram_ and not normalizedWasserstein_);
@@ -200,7 +200,7 @@ int ttkMergeTreeAutoencoder::runCompute(
       = (not isPersistenceDiagram_ ? inputTrees2 : inputTrees);
     ttk::ftm::constructTrees<float>(inputTrees2ToUse, intermediateMTrees2,
                                     treesNodes2, treesArcs2, treesSegmentation2,
-                                    !useSadMaxPairs);
+                                    !useSecondPairsType, DiagramPairTypes);
   }
   isPersistenceDiagram_ |= (not normalizedWasserstein_);
 
@@ -529,6 +529,11 @@ int ttkMergeTreeAutoencoder::runOutput(
   arrayActivateFunction->InsertNextTuple1(activationFunction_);
   output_coef->GetFieldData()->AddArray(arrayActivateFunction);
 
+  vtkNew<vtkIntArray> diagramPairTypesArray{};
+  diagramPairTypesArray->SetName("DiagramPairTypes");
+  diagramPairTypesArray->InsertNextTuple1(DiagramPairTypes);
+  output_coef->GetFieldData()->AddArray(diagramPairTypesArray);
+
   // ------------------------------------------
   // --- Axes Vectors
   // ------------------------------------------

core/vtk/ttkMergeTreeAutoencoder/ttkMergeTreeAutoencoder.h

@@ -68,6 +68,7 @@ class TTKMERGETREEAUTOENCODER_EXPORT ttkMergeTreeAutoencoder
    */
   // Input options
   double oldEpsilonTree1;
+  int DiagramPairTypes = 0;
 
   // ----------------------
   // Data for visualization
@@ -442,6 +443,13 @@ class TTKMERGETREEAUTOENCODER_EXPORT ttkMergeTreeAutoencoder
     return mixtureCoefficient_;
   }
 
+  void SetDiagramPairTypes(int diagramPairTypes) {
+    DiagramPairTypes = diagramPairTypes;
+    Modified();
+    resetDataVisualization();
+  }
+  vtkGetMacro(DiagramPairTypes, int);
+
   void SetEpsilon1UseFarthestSaddle(bool epsilon1UseFarthestSaddle) {
     epsilon1UseFarthestSaddle_ = epsilon1UseFarthestSaddle;
     Modified();

core/vtk/ttkMergeTreeAutoencoderDecoding/ttkMergeTreeAutoencoderDecoding.cpp

@@ -141,6 +141,10 @@ int ttkMergeTreeAutoencoderDecoding::RequestData(
   else
     printMsg("Computation without normalized Wasserstein.");
 
+  auto diagramPairTypesArray = fd->GetArray("DiagramPairTypes");
+  if(diagramPairTypesArray)
+    DiagramPairTypes = diagramPairTypesArray->GetTuple1(0);
+
   // -----------------
   // Origins
   // -----------------
@@ -156,13 +160,13 @@ int ttkMergeTreeAutoencoderDecoding::RequestData(
   std::vector<vtkDataSet *> originsTreeSegmentations,
     originsPrimeTreeSegmentations;
 
-  bool useSadMaxPairs = (mixtureCoefficient_ == 0);
+  bool useSecondPairsType = (mixtureCoefficient_ == 0);
   isPersistenceDiagram_ = ttk::ftm::constructTrees<float>(
     origins, originsTrees, originsTreeNodes, originsTreeArcs,
-    originsTreeSegmentations, useSadMaxPairs);
-  ttk::ftm::constructTrees<float>(originsPrime, originsPrimeTrees,
-                                  originsTreeNodes, originsTreeArcs,
-                                  originsTreeSegmentations, useSadMaxPairs);
+    originsTreeSegmentations, useSecondPairsType, DiagramPairTypes);
+  ttk::ftm::constructTrees<float>(
+    originsPrime, originsPrimeTrees, originsTreeNodes, originsTreeArcs,
+    originsTreeSegmentations, useSecondPairsType, DiagramPairTypes);
   // If merge trees are provided in input and normalization is not asked
   convertToDiagram_
     = (not isPersistenceDiagram_ and not normalizedWasserstein_);

core/vtk/ttkMergeTreeAutoencoderDecoding/ttkMergeTreeAutoencoderDecoding.h

@@ -59,6 +59,7 @@ class TTKMERGETREEAUTOENCODERDECODING_EXPORT ttkMergeTreeAutoencoderDecoding
    * Add all filter parameters only as private member variables and
    * initialize them here.
    */
+  int DiagramPairTypes = 0;
 
 public:
   /**

core/vtk/ttkMergeTreeClustering/ttkMergeTreeClustering.cpp

@@ -185,17 +185,18 @@ int ttkMergeTreeClustering::runCompute(
   std::vector<FTMTree_MT *> intermediateTrees(numInputs),
     intermediateTrees2(numInputs2);
 
-  bool const useSadMaxPairs = (JoinSplitMixtureCoefficient == 0);
-  IsPersistenceDiagram
-    = constructTrees<dataType>(inputTrees, intermediateMTrees, treesNodes,
-                               treesArcs, treesSegmentation, useSadMaxPairs);
+  bool const useSecondPairsType = (JoinSplitMixtureCoefficient == 0);
+  IsPersistenceDiagram = constructTrees<dataType>(
+    inputTrees, intermediateMTrees, treesNodes, treesArcs, treesSegmentation,
+    useSecondPairsType, DiagramPairTypes);
   if(not IsPersistenceDiagram
      or (JoinSplitMixtureCoefficient != 0
          and JoinSplitMixtureCoefficient != 1)) {
     auto &inputTrees2ToUse
       = (not IsPersistenceDiagram ? inputTrees2 : inputTrees);
     constructTrees<dataType>(inputTrees2ToUse, intermediateMTrees2, treesNodes2,
-                             treesArcs2, treesSegmentation2, !useSadMaxPairs);
+                             treesArcs2, treesSegmentation2,
+                             !useSecondPairsType, DiagramPairTypes);
   }
 
   mergeTreeToFTMTree<dataType>(intermediateMTrees, intermediateTrees);

core/vtk/ttkMergeTreeClustering/ttkMergeTreeClustering.h

@@ -92,6 +92,7 @@ class TTKMERGETREECLUSTERING_EXPORT ttkMergeTreeClustering
   bool oldNW = NormalizedWasserstein;
   bool oldKS = KeepSubtree;
   double JoinSplitMixtureCoefficient = 0.5;
+  int DiagramPairTypes = 0;
   bool ComputeBarycenter = false;
   unsigned int NumberOfBarycenters = 1;
   double BarycenterSizeLimitPercent = 0.0;
@@ -314,6 +315,13 @@ class TTKMERGETREECLUSTERING_EXPORT ttkMergeTreeClustering
   }
   vtkGetMacro(JoinSplitMixtureCoefficient, double);
 
+  void SetDiagramPairTypes(int diagramPairTypes) {
+    DiagramPairTypes = diagramPairTypes;
+    Modified();
+    resetDataVisualization();
+  }
+  vtkGetMacro(DiagramPairTypes, int);
+
   void SetComputeBarycenter(bool computeBarycenter) {
     ComputeBarycenter = computeBarycenter;
     Modified();

core/vtk/ttkMergeTreeDistanceMatrix/ttkMergeTreeDistanceMatrix.cpp

@@ -102,14 +102,16 @@ int ttkMergeTreeDistanceMatrix::run(
   // Construct trees
   const int numInputs = inputTrees.size();
   std::vector<MergeTree<dataType>> intermediateTrees, intermediateTrees2;
-  bool const useSadMaxPairs = (mixtureCoefficient_ == 0); // only for PD support
-  isPersistenceDiagram_
-    = constructTrees(inputTrees, intermediateTrees, useSadMaxPairs);
+  bool const useSecondPairsType
+    = (mixtureCoefficient_ == 0); // only for PD support
+  isPersistenceDiagram_ = constructTrees(
+    inputTrees, intermediateTrees, useSecondPairsType, DiagramPairTypes);
   if(not isPersistenceDiagram_
      or (mixtureCoefficient_ != 0 and mixtureCoefficient_ != 1)) {
     auto &inputTrees2ToUse
       = (not isPersistenceDiagram_ ? inputTrees2 : inputTrees);
-    constructTrees(inputTrees2ToUse, intermediateTrees2, !useSadMaxPairs);
+    constructTrees(inputTrees2ToUse, intermediateTrees2, !useSecondPairsType,
+                   DiagramPairTypes);
   }
 
   // Verify parameters

core/vtk/ttkMergeTreeDistanceMatrix/ttkMergeTreeDistanceMatrix.h

@@ -62,6 +62,7 @@ class TTKMERGETREEDISTANCEMATRIX_EXPORT ttkMergeTreeDistanceMatrix
   bool oldKS = keepSubtree_;
 
   bool UseFieldDataParameters = false;
+  int DiagramPairTypes = 0;
 
 public:
   /**
@@ -190,6 +191,9 @@ class TTKMERGETREEDISTANCEMATRIX_EXPORT ttkMergeTreeDistanceMatrix
   vtkSetMacro(mixtureCoefficient_, double);
   vtkGetMacro(mixtureCoefficient_, double);
 
+  vtkSetMacro(DiagramPairTypes, int);
+  vtkGetMacro(DiagramPairTypes, int);
+
   /**
    * This static method and the macro below are VTK conventions on how to
    * instantiate VTK objects. You don't have to modify this.

core/vtk/ttkMergeTreePrincipalGeodesics/ttkMergeTreePrincipalGeodesics.cpp

@@ -154,10 +154,10 @@ int ttkMergeTreePrincipalGeodesics::runCompute(
   std::vector<ttk::ftm::MergeTree<dataType>> intermediateMTrees,
     intermediateMTrees2;
 
-  bool const useSadMaxPairs = (mixtureCoefficient_ == 0);
+  bool const useSecondPairsType = (mixtureCoefficient_ == 0);
   isPersistenceDiagram_ = ttk::ftm::constructTrees<dataType>(
     inputTrees, intermediateMTrees, treesNodes, treesArcs, treesSegmentation,
-    useSadMaxPairs);
+    useSecondPairsType, DiagramPairTypes);
   // If merge trees are provided in input and normalization is not asked
   convertToDiagram_
     = (not isPersistenceDiagram_ and not normalizedWasserstein_);
@@ -169,9 +169,9 @@ int ttkMergeTreePrincipalGeodesics::runCompute(
      or (mixtureCoefficient_ != 0 and mixtureCoefficient_ != 1)) {
     auto &inputTrees2ToUse
       = (not isPersistenceDiagram_ ? inputTrees2 : inputTrees);
-    ttk::ftm::constructTrees<dataType>(inputTrees2ToUse, intermediateMTrees2,
-                                       treesNodes2, treesArcs2,
-                                       treesSegmentation2, !useSadMaxPairs);
+    ttk::ftm::constructTrees<dataType>(
+      inputTrees2ToUse, intermediateMTrees2, treesNodes2, treesArcs2,
+      treesSegmentation2, !useSecondPairsType, DiagramPairTypes);
   }
   isPersistenceDiagram_ |= (not normalizedWasserstein_);
 
@@ -329,6 +329,10 @@ int ttkMergeTreePrincipalGeodesics::runOutput(
     array->InsertNextTuple1(getParamValueFromName(paramName));
     output_coef->GetFieldData()->AddArray(array);
   }
+  vtkNew<vtkIntArray> diagramPairTypesArray{};
+  diagramPairTypesArray->SetName("DiagramPairTypes");
+  diagramPairTypesArray->InsertNextTuple1(DiagramPairTypes);
+  output_coef->GetFieldData()->AddArray(diagramPairTypesArray);
 
   // ------------------------------------------
   // --- Geodesics Vectors

core/vtk/ttkMergeTreePrincipalGeodesics/ttkMergeTreePrincipalGeodesics.h

@@ -71,6 +71,7 @@ class TTKMERGETREEPRINCIPALGEODESICS_EXPORT ttkMergeTreePrincipalGeodesics
    */
   // Input options
   double oldEpsilonTree1;
+  int DiagramPairTypes = 0;
   // Output options
 
   // ----------------------
@@ -169,6 +170,13 @@ class TTKMERGETREEPRINCIPALGEODESICS_EXPORT ttkMergeTreePrincipalGeodesics
     return mixtureCoefficient_;
   }
 
+  void SetDiagramPairTypes(int diagramPairTypes) {
+    DiagramPairTypes = diagramPairTypes;
+    Modified();
+    resetDataVisualization();
+  }
+  vtkGetMacro(DiagramPairTypes, int);
+
   void SetKeepState(bool keepState) {
     keepState_ = keepState;
     Modified();

core/vtk/ttkMergeTreePrincipalGeodesicsDecoding/ttkMergeTreePrincipalGeodesicsDecoding.cpp

@@ -154,6 +154,11 @@ int ttkMergeTreePrincipalGeodesicsDecoding::RequestData(
   else
     printMsg("Computation without normalized Wasserstein.");
 
+  auto diagramPairTypesArray
+    = tableCoefficients->GetFieldData()->GetArray("DiagramPairTypes");
+  if(diagramPairTypesArray)
+    DiagramPairTypes = diagramPairTypesArray->GetTuple1(0);
+
   // ------------------------------------------------------------------------------------
   // --- Load tables
   // ------------------------------------------------------------------------------------
@@ -316,21 +321,21 @@ int ttkMergeTreePrincipalGeodesicsDecoding::runCompute(
 
   std::vector<ttk::ftm::MergeTree<dataType>> baryDTree, inputDTrees;
 
-  std::vector<bool> useSadMaxPairsVec{false, true};
+  std::vector<bool> useSecondPairsTypeVec{false, true};
   if(not useDoubleInput_ and mixtureCoefficient_ == 0)
-    useSadMaxPairsVec.erase(useSadMaxPairsVec.begin()); // {true}
+    useSecondPairsTypeVec.erase(useSecondPairsTypeVec.begin()); // {true}
   ttk::ftm::constructTrees<dataType>(inputBary, baryDTree, baryTreeNodes,
                                      baryTreeArcs, baryTreeSegmentation,
-                                     useSadMaxPairsVec);
+                                     useSecondPairsTypeVec, DiagramPairTypes);
 
   if(OutputInputTrees
      or (ReconstructInputTrees
          and (computeReconstructionError_ or transferInputTreesInformation_))) {
-    bool const useSadMaxPairs
+    bool const useSecondPairsType
       = (useDoubleInput_ and not processFirstInput) or mixtureCoefficient_ == 0;
     bool const isInputPD = ttk::ftm::constructTrees<dataType>(
       inputTrees, inputDTrees, inputTreesNodes, inputTreesArcs,
-      inputTreesSegmentation, useSadMaxPairs);
+      inputTreesSegmentation, useSecondPairsType, DiagramPairTypes);
     if(not isInputPD and isPersistenceDiagram_)
       mtsFlattening(inputDTrees);
   }