ENH: Segmentation Filters OoC Optimization

CMakeLists.txt

@@ -537,6 +537,7 @@ set(SIMPLNX_HDRS
   ${SIMPLNX_SOURCE_DIR}/Utilities/DataGroupUtilities.hpp
   ${SIMPLNX_SOURCE_DIR}/Utilities/DataObjectUtilities.hpp
   ${SIMPLNX_SOURCE_DIR}/Utilities/DataStoreUtilities.hpp
+  ${SIMPLNX_SOURCE_DIR}/Utilities/AlgorithmDispatch.hpp
   ${SIMPLNX_SOURCE_DIR}/Utilities/FilePathGenerator.hpp
   ${SIMPLNX_SOURCE_DIR}/Utilities/ColorTableUtilities.hpp
   ${SIMPLNX_SOURCE_DIR}/Utilities/FileUtilities.hpp
@@ -558,6 +559,7 @@ set(SIMPLNX_HDRS
   ${SIMPLNX_SOURCE_DIR}/Utilities/ParallelTaskAlgorithm.hpp
   ${SIMPLNX_SOURCE_DIR}/Utilities/SamplingUtils.hpp
   ${SIMPLNX_SOURCE_DIR}/Utilities/SegmentFeatures.hpp
+  ${SIMPLNX_SOURCE_DIR}/Utilities/UnionFind.hpp
   ${SIMPLNX_SOURCE_DIR}/Utilities/TimeUtilities.hpp
   ${SIMPLNX_SOURCE_DIR}/Utilities/TooltipGenerator.hpp
   ${SIMPLNX_SOURCE_DIR}/Utilities/TooltipRowItem.hpp

src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/CAxisSegmentFeatures.cpp

@@ -5,6 +5,7 @@
 #include "simplnx/Common/Constants.hpp"
 #include "simplnx/DataStructure/DataArray.hpp"
 #include "simplnx/DataStructure/Geometry/ImageGeom.hpp"
+#include "simplnx/Utilities/AlgorithmDispatch.hpp"
 #include "simplnx/Utilities/ClusteringUtilities.hpp"
 
 #include <EbsdLib/Orientation/OrientationFwd.hpp>
@@ -68,8 +69,16 @@ Result<> CAxisSegmentFeatures::operator()()
   auto* active = m_DataStructure.getDataAs<UInt8Array>(m_InputValues->ActiveArrayPath);
   active->fill(1);
 
-  // Run the segmentation algorithm
-  execute(imageGeometry);
+  // Dispatch between DFS (in-core) and CCL (OOC) algorithms
+  if(IsOutOfCore(*m_FeatureIdsArray) || ForceOocAlgorithm())
+  {
+    auto& featureIdsStore = m_FeatureIdsArray->getDataStoreRef();
+    executeCCL(imageGeometry, featureIdsStore);
+  }
+  else
+  {
+    execute(imageGeometry);
+  }
   // Sanity check the result.
   if(this->m_FoundFeatures < 1)
   {
@@ -127,10 +136,7 @@ int64 CAxisSegmentFeatures::getSeed(int32 gnum, int64 nextSeed) const
   }
   if(seed >= 0)
   {
-    auto& cellFeatureAM = m_DataStructure.getDataRefAs<AttributeMatrix>(m_InputValues->CellFeatureAttributeMatrixPath);
     featureIds[static_cast<usize>(seed)] = gnum;
-    const ShapeType tDims = {static_cast<usize>(gnum) + 1};
-    cellFeatureAM.resizeTuples(tDims); // This will resize the active array
   }
   return seed;
 }
@@ -182,3 +188,59 @@ bool CAxisSegmentFeatures::determineGrouping(int64 referencepoint, int64 neighbo
   }
   return group;
 }
+
+// -----------------------------------------------------------------------------
+bool CAxisSegmentFeatures::isValidVoxel(int64 point) const
+{
+  // Check mask
+  if(m_InputValues->UseMask && !m_GoodVoxelsArray->isTrue(point))
+  {
+    return false;
+  }
+  // Check that the voxel has a valid phase (> 0)
+  Int32Array& cellPhases = *m_CellPhases;
+  if(cellPhases[point] <= 0)
+  {
+    return false;
+  }
+  return true;
+}
+
+// -----------------------------------------------------------------------------
+bool CAxisSegmentFeatures::areNeighborsSimilar(int64 point1, int64 point2) const
+{
+  // The neighbor must also be valid
+  if(!isValidVoxel(point2))
+  {
+    return false;
+  }
+
+  Int32Array& cellPhases = *m_CellPhases;
+
+  // Must be same phase
+  if(cellPhases[point1] != cellPhases[point2])
+  {
+    return false;
+  }
+
+  // Calculate c-axis misalignment
+  const Eigen::Vector3f cAxis{0.0f, 0.0f, 1.0f};
+  Float32Array& quats = *m_QuatsArray;
+
+  const ebsdlib::QuatF q1(quats[point1 * 4], quats[point1 * 4 + 1], quats[point1 * 4 + 2], quats[point1 * 4 + 3]);
+  const ebsdlib::QuatF q2(quats[point2 * 4], quats[point2 * 4 + 1], quats[point2 * 4 + 2], quats[point2 * 4 + 3]);
+
+  const ebsdlib::OrientationMatrixFType oMatrix1 = q1.toOrientationMatrix();
+  const ebsdlib::OrientationMatrixFType oMatrix2 = q2.toOrientationMatrix();
+
+  Eigen::Vector3f c1 = oMatrix1.transpose() * cAxis;
+  Eigen::Vector3f c2 = oMatrix2.transpose() * cAxis;
+
+  c1.normalize();
+  c2.normalize();
+
+  float32 w = std::clamp(((c1[0] * c2[0]) + (c1[1] * c2[1]) + (c1[2] * c2[2])), -1.0F, 1.0F);
+  w = std::acos(w);
+
+  return w <= m_InputValues->MisorientationTolerance || (Constants::k_PiD - w) <= m_InputValues->MisorientationTolerance;
+}

src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/CAxisSegmentFeatures.hpp

@@ -49,6 +49,21 @@ class ORIENTATIONANALYSIS_EXPORT CAxisSegmentFeatures : public SegmentFeatures
   int64 getSeed(int32 gnum, int64 nextSeed) const override;
   bool determineGrouping(int64 referencePoint, int64 neighborPoint, int32 gnum) const override;
 
+  /**
+   * @brief Checks whether a voxel can participate in C-axis segmentation based on mask and phase.
+   * @param point Linear voxel index.
+   * @return true if the voxel passes mask and phase checks.
+   */
+  bool isValidVoxel(int64 point) const override;
+
+  /**
+   * @brief Determines whether two neighboring voxels belong to the same C-axis segment.
+   * @param point1 First voxel index.
+   * @param point2 Second (neighbor) voxel index.
+   * @return true if both voxels share the same phase and their C-axis misalignment is within tolerance.
+   */
+  bool areNeighborsSimilar(int64 point1, int64 point2) const override;
+
 private:
   const CAxisSegmentFeaturesInputValues* m_InputValues = nullptr;
 

src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/EBSDSegmentFeatures.cpp

@@ -2,6 +2,7 @@
 
 #include "simplnx/DataStructure/DataStore.hpp"
 #include "simplnx/DataStructure/Geometry/IGridGeometry.hpp"
+#include "simplnx/Utilities/AlgorithmDispatch.hpp"
 
 using namespace nx::core;
 
@@ -43,8 +44,16 @@ Result<> EBSDSegmentFeatures::operator()()
   m_FeatureIdsArray = m_DataStructure.getDataAs<Int32Array>(m_InputValues->FeatureIdsArrayPath);
   m_FeatureIdsArray->fill(0); // initialize the output array with zeros
 
-  // Run the segmentation algorithm
-  execute(gridGeom);
+  // Dispatch between DFS (in-core) and CCL (OOC) algorithms
+  if(IsOutOfCore(*m_FeatureIdsArray) || ForceOocAlgorithm())
+  {
+    auto& featureIdsStore = m_FeatureIdsArray->getDataStoreRef();
+    executeCCL(gridGeom, featureIdsStore);
+  }
+  else
+  {
+    execute(gridGeom);
+  }
   // Sanity check the result.
   if(this->m_FoundFeatures < 1)
   {
@@ -152,3 +161,55 @@ bool EBSDSegmentFeatures::determineGrouping(int64 referencePoint, int64 neighbor
 
   return group;
 }
+
+// -----------------------------------------------------------------------------
+bool EBSDSegmentFeatures::isValidVoxel(int64 point) const
+{
+  // Check mask
+  if(m_InputValues->UseMask && !m_GoodVoxelsArray->isTrue(point))
+  {
+    return false;
+  }
+  // Check that the voxel has a valid phase (> 0)
+  AbstractDataStore<int32>& cellPhases = m_CellPhases->getDataStoreRef();
+  if(cellPhases[point] <= 0)
+  {
+    return false;
+  }
+  return true;
+}
+
+// -----------------------------------------------------------------------------
+bool EBSDSegmentFeatures::areNeighborsSimilar(int64 point1, int64 point2) const
+{
+  // The neighbor must also be valid
+  if(!isValidVoxel(point2))
+  {
+    return false;
+  }
+
+  AbstractDataStore<int32>& cellPhases = m_CellPhases->getDataStoreRef();
+
+  // Must be same phase
+  if(cellPhases[point1] != cellPhases[point2])
+  {
+    return false;
+  }
+
+  // Check crystal structure validity
+  int32 laueClass = (*m_CrystalStructures)[cellPhases[point1]];
+  if(static_cast<usize>(laueClass) >= m_OrientationOps.size())
+  {
+    return false;
+  }
+
+  // Calculate misorientation
+  Float32Array& quats = *m_QuatsArray;
+  const ebsdlib::QuatD q1(quats[point1 * 4], quats[point1 * 4 + 1], quats[point1 * 4 + 2], quats[point1 * 4 + 3]);
+  const ebsdlib::QuatD q2(quats[point2 * 4], quats[point2 * 4 + 1], quats[point2 * 4 + 2], quats[point2 * 4 + 3]);
+
+  ebsdlib::AxisAngleDType axisAngle = m_OrientationOps[laueClass]->calculateMisorientation(q1, q2);
+  float w = static_cast<float>(axisAngle[3]);
+
+  return w < m_InputValues->MisorientationTolerance;
+}

src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/EBSDSegmentFeatures.hpp

@@ -56,26 +56,23 @@ class ORIENTATIONANALYSIS_EXPORT EBSDSegmentFeatures : public SegmentFeatures
   Result<> operator()();
 
 protected:
+  int64_t getSeed(int32 gnum, int64 nextSeed) const override;
+  bool determineGrouping(int64 referencePoint, int64 neighborPoint, int32 gnum) const override;
+
   /**
-   * @brief
-   * @param data
-   * @param args
-   * @param gnum
-   * @param nextSeed
-   * @return int64
+   * @brief Checks whether a voxel can participate in EBSD segmentation based on mask and phase.
+   * @param point Linear voxel index.
+   * @return true if the voxel passes mask and phase checks.
    */
-  int64_t getSeed(int32 gnum, int64 nextSeed) const override;
+  bool isValidVoxel(int64 point) const override;
 
   /**
-   * @brief
-   * @param data
-   * @param args
-   * @param referencepoint
-   * @param neighborpoint
-   * @param gnum
-   * @return bool
+   * @brief Determines whether two neighboring voxels belong to the same EBSD segment.
+   * @param point1 First voxel index.
+   * @param point2 Second (neighbor) voxel index.
+   * @return true if both voxels share the same phase and their misorientation is within tolerance.
    */
-  bool determineGrouping(int64 referencePoint, int64 neighborPoint, int32 gnum) const override;
+  bool areNeighborsSimilar(int64 point1, int64 point2) const override;
 
 private:
   const EBSDSegmentFeaturesInputValues* m_InputValues = nullptr;