VV: Segment Features (C-Axis Misalignment) fully V&V'ed

docs/vv_templates/vv_policy.md

@@ -93,4 +93,4 @@ When SIMPLNX and 6.5.171 differ, each Deviation entry names one root cause (or a
 - **Precision** — different floating-point width or intermediate-math type
 - **Order of operations** — associativity differences in parallel reductions, different loop order, different accumulation
 - **Library** — Eigen vs. hand-rolled, different EbsdLib version, different HDF5 versions
-- **Algorithmic choice** — deliberate change in method (used only for Rewrite relationships)
+- **Algorithmic choice** — a deliberate difference in method or behavior (an added validation guard, a changed default, exposed randomization, a replacement algorithm). Any Algorithm Relationship may carry algorithmic-choice deviations — a **Port** or **Minor changes** filter that adds a guard legacy lacked is the common case. A **Rewrite** consists of algorithmic choices by definition, and a Rewrite whose outputs diverge must defend the shared UUID (see above)

src/Plugins/OrientationAnalysis/docs/CAxisSegmentFeaturesFilter.md

@@ -6,17 +6,27 @@ Reconstruction (Segmentation)
 
 ## Description
 
-This **Filter** segments the **Features** by grouping neighboring **Cells** that satisfy the *C-axis misalignment tolerance*, i.e., have misalignment angle less than the value set by the user. The *C-axis misalignment* refers to the angle between the <001> directions (C-axis in the hexagonal system) that is present between neighboring **Cells**.  The process by which the **Features** are identified is given below and is a standard *burn algorithm*.
+This **Filter** segments the **Features** by grouping neighboring **Cells** that satisfy the *c-axis misalignment tolerance*, i.e., have misalignment angle less than the value set by the user. The *c-axis misalignment* refers to the angle between the [0001] directions (the c-axis in the hexagonal system) that is present between neighboring **Cells**. Because the c-axis is a direction (not a vector), the misalignment is folded into [0°, 90°]: two nearly antiparallel c-axes are considered aligned. The process by which the **Features** are identified is given below and is a standard *burn algorithm*.
 
-1. Randomly select a **Cell**, add it to an empty list and set its *Feature Id* to the current **Feature**
-2. Compare the **Cell** to each of its six (6) face-sharing neighbors (i.e., calculate the c-axis misalignment with each neighbor)
+1. Select the next unassigned **Cell** in row-major order that is eligible to seed a **Feature** (not excluded by the mask and with a phase value greater than 0), add it to an empty list and set its *Feature Id* to the current **Feature**
+2. Compare the **Cell** to each of its neighbors as selected by the *Neighbor Scheme* parameter (i.e., calculate the c-axis misalignment with each neighbor)
 3. Add each neighbor **Cell** that has a C-axis misalignment below the user defined tolerance to the list created in 1. and set the *Feature Id* of the neighbor **Cell** to the current **Feature**
 4. Remove the current **Cell** from the list and move to the next **Cell** and repeat 2. and 3.; continue until no **Cells** are left in the list
-5. Increment the current **Feature** counter and repeat steps 1. through 4.; continue until no **Cells** remain unassigned in the dataset
+5. Increment the current **Feature** counter and repeat steps 1. through 4.; continue until no eligible **Cells** remain unassigned in the dataset
 
-The user has the option to *Use Mask Array*, which allows the user to set a boolean array for the **Cells** that remove **Cells** with a value of *false* from consideration in the above algorithm. This option is useful if the user has an array that either specifies the domain of the "sample" in the "image" or specifies if the orientation on the **Cell** is trusted/correct.
+The user has the option to *Use Mask Array*, which allows the user to set a boolean (or uint8) array for the **Cells** that removes **Cells** with a value of *false* from consideration in the above algorithm. This option is useful if the user has an array that either specifies the domain of the "sample" in the "image" or specifies if the orientation on the **Cell** is trusted/correct. Masked-out **Cells** and unindexed **Cells** (phase 0) never join a **Feature** and keep a *Feature Id* of 0.
 
-After all the **Features** have been identified, a **Feature Attribute Matrix** is created for the **Features** and each **Feature** is flagged as *Active* in a boolean array in the matrix.
+After all the **Features** have been identified, a **Feature Attribute Matrix** is created for the **Features** and each **Feature** is flagged as *Active* in the matrix (index 0 is reserved and never active).
+
+The input geometry may be either an **Image Geometry** or a **RectGrid Geometry**.
+
+### Hexagonal Crystal Structures Required
+
+The c-axis is only a physically meaningful unique axis for hexagonal Laue classes. Every **Cell** that can participate in the segmentation (phase > 0 and not excluded by the mask) must belong to an **Ensemble** whose crystal structure is *Hexagonal-Low (6/m)* or *Hexagonal-High (6/mmm)*; otherwise the filter fails with error `-8363`. A phase value with no corresponding entry in the *Crystal Structures* array produces error `-8364`. Unindexed **Cells** (phase 0) and masked-out **Cells** are exempt from this requirement, so datasets with unindexed points — or with a non-hexagonal phase that has been masked out — process normally.
+
+### Randomize Feature Ids
+
+When *Randomize Feature Ids* is enabled the final *Feature Ids* are relabeled with a deterministic (fixed-seed) random permutation, which improves visual contrast between neighboring **Features** when coloring by *Feature Id*. The segmentation itself is unchanged, and repeated runs produce identical output. (Legacy DREAM.3D 6.x always randomized with a clock-derived seed, so its labeling differed on every run; see the migration deviation notes in `vv/deviations/CAxisSegmentFeaturesFilter.md` of the simplnx source tree.)
 
 ### Neighbor Scheme
 

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

@@ -4,14 +4,16 @@
 
 #include "simplnx/Common/Constants.hpp"
 #include "simplnx/DataStructure/DataArray.hpp"
-#include "simplnx/DataStructure/Geometry/ImageGeom.hpp"
+#include "simplnx/DataStructure/Geometry/IGridGeometry.hpp"
 #include "simplnx/Utilities/ClusteringUtilities.hpp"
 
 #include <EbsdLib/Orientation/OrientationFwd.hpp>
 #include <EbsdLib/Orientation/OrientationMatrix.hpp>
 #include <EbsdLib/Orientation/Quaternion.hpp>
 
+#include <algorithm>
 #include <cmath>
+#include <vector>
 
 using namespace nx::core;
 using namespace nx::core::OrientationUtilities;
@@ -29,8 +31,9 @@ CAxisSegmentFeatures::~CAxisSegmentFeatures() noexcept = default;
 // -----------------------------------------------------------------------------
 Result<> CAxisSegmentFeatures::operator()()
 {
-  this->m_NeighborScheme = m_InputValues->NeighborScheme;
-  auto* imageGeometry = m_DataStructure.getDataAs<ImageGeom>(m_InputValues->ImageGeometryPath);
+  m_NeighborScheme = m_InputValues->NeighborScheme;
+  // The geometry parameter accepts Image AND RectGrid geometries; both derive from IGridGeometry.
+  auto* gridGeom = m_DataStructure.getDataAs<IGridGeometry>(m_InputValues->ImageGeometryPath);
   m_QuatsArray = m_DataStructure.getDataAs<Float32Array>(m_InputValues->QuatsArrayPath);
   m_CellPhases = m_DataStructure.getDataAs<Int32Array>(m_InputValues->CellPhasesArrayPath);
   if(m_InputValues->UseMask)
@@ -42,67 +45,92 @@ Result<> CAxisSegmentFeatures::operator()()
     {
       // This really should NOT be happening as the path was verified during preflight BUT we may be calling this from
       // somewhere else that is NOT going through the normal nx::core::IFilter API of Preflight and Execute
-      return MakeErrorResult(-8362, fmt::format("Mask Array DataPath does not exist or is not of the correct type (Bool | UInt8) {}", m_InputValues->MaskArrayPath.toString()));
+      return MakeErrorResult(-8362, fmt::format("Mask Array DataPath '{}' does not exist or is not of the correct type (Bool | UInt8).", m_InputValues->MaskArrayPath.toString()));
     }
   }
 
   // Loop through all the "Phase" cell values and validate that any phase found is
   // a hexagonal phase. This guards against there being multiple phases defined in
-  // and EBSD file but the non-hexagonal phases were actually never found
+  // an EBSD file where the non-hexagonal phases were never actually indexed to a cell.
   const auto& crystalStructures = m_DataStructure.getDataRefAs<UInt32Array>(m_InputValues->CrystalStructuresArrayPath);
-  usize numCells = m_CellPhases->getNumberOfTuples();
+  const usize numCells = m_CellPhases->getNumberOfTuples();
+  const usize numEnsembles = crystalStructures.getNumberOfTuples();
+  // Each phase only needs to be validated once. Caching validated phases keeps the loop from
+  // paying a per-cell mask lookup (expensive per-element access on out-of-core stores) for every
+  // cell of an already-validated phase — in the common all-hexagonal case the mask is consulted
+  // only for the first cell of each phase.
+  std::vector<uint8> phaseValidated(numEnsembles, 0);
   for(usize cellIdx = 0; cellIdx < numCells; ++cellIdx)
   {
-    int32 currentPhaseIdx = m_CellPhases->getValue(cellIdx);
-    // Only consider valid ebsd phase values
-    if(currentPhaseIdx < 1)
+    const int32 currentPhaseIdx = m_CellPhases->getValue(cellIdx);
+    // Cells that can never seed or join a feature are exempt from the crystal-structure
+    // requirement: phase 0 is the conventional "unindexed" phase (its CrystalStructures entry is
+    // the 999 "unknown" sentinel), and masked-out cells are excluded from segmentation entirely.
+    if(currentPhaseIdx <= 0)
     {
       continue;
     }
-    // Only consider valid mask values
-    if(m_InputValues->UseMask && !m_GoodVoxelsArray->isTrue(cellIdx))
+    if(static_cast<usize>(currentPhaseIdx) < numEnsembles && phaseValidated[currentPhaseIdx] != 0)
     {
       continue;
     }
+    if(m_GoodVoxelsArray != nullptr && !m_GoodVoxelsArray->isTrue(cellIdx))
+    {
+      continue;
+    }
+    if(static_cast<usize>(currentPhaseIdx) >= numEnsembles)
+    {
+      return MakeErrorResult(-8364, fmt::format("Cell {} has a phase value of {} but the Crystal Structures array '{}' only has {} entries.", cellIdx, currentPhaseIdx,
+                                                m_InputValues->CrystalStructuresArrayPath.toString(), numEnsembles));
+    }
     const auto crystalStructureType = crystalStructures[currentPhaseIdx];
 
     if(crystalStructureType != ebsdlib::CrystalStructure::Hexagonal_High && crystalStructureType != ebsdlib::CrystalStructure::Hexagonal_Low)
     {
-      return MakeErrorResult(-8363, fmt::format("Input data is using {} type crystal structures but segmenting features via c-axis mis orientation requires all phases to be either Hexagonal-Low 6/m "
-                                                "or Hexagonal-High 6/mmm type crystal structures.",
+      return MakeErrorResult(-8363, fmt::format("Input data is using {} type crystal structures but segmenting features via c-axis misorientation requires every phase that participates in the "
+                                                "segmentation to be either Hexagonal-Low 6/m or Hexagonal-High 6/mmm type crystal structures.",
                                                 CrystalStructureEnumToString(crystalStructureType)));
     }
+    phaseValidated[currentPhaseIdx] = 1;
   }
 
   m_FeatureIdsArray = m_DataStructure.getDataAs<Int32Array>(m_InputValues->FeatureIdsArrayPath);
   m_FeatureIdsArray->fill(0);
-  auto* active = m_DataStructure.getDataAs<UInt8Array>(m_InputValues->ActiveArrayPath);
-  active->fill(1);
 
   // Run the segmentation algorithm
-  execute(imageGeometry);
+  Result<> segmentResult = execute(gridGeom);
+  if(segmentResult.invalid())
+  {
+    return segmentResult;
+  }
+  // A canceled run returns successfully with no feature count; do not misreport it below as
+  // "no Features were detected".
+  if(m_ShouldCancel)
+  {
+    return {};
+  }
   // Sanity check the result.
-  if(this->m_FoundFeatures < 1)
+  if(m_FoundFeatures < 1)
   {
-    return MakeErrorResult(-87000, fmt::format("The number of Features is '{}' which means no Features were detected. A threshold value may be set incorrectly", this->m_FoundFeatures));
+    return MakeErrorResult(-87000, "No Features were detected: no Cell was eligible to seed a Feature. Every Cell is either excluded by the Mask or has a Phase value of 0 (unindexed).");
   }
 
   // Resize the Feature Attribute Matrix
-  ShapeType tDims = {static_cast<usize>(this->m_FoundFeatures + 1)};
-  auto& cellFeaturesAM = m_DataStructure.getDataRefAs<AttributeMatrix>(m_InputValues->CellFeatureAttributeMatrixPath);
-  cellFeaturesAM.resizeTuples(tDims); // This will resize the active array
+  ShapeType tDims = {static_cast<usize>(m_FoundFeatures + 1)};
+  auto& cellFeatureAM = m_DataStructure.getDataRefAs<AttributeMatrix>(m_InputValues->CellFeatureAttributeMatrixPath);
+  cellFeatureAM.resizeTuples(tDims); // This will resize the active array
 
   // make sure all values are initialized and "re-reserve" index 0
   auto* activeArray = m_DataStructure.getDataAs<UInt8Array>(m_InputValues->ActiveArrayPath);
   activeArray->getDataStore()->fill(1);
   (*activeArray)[0] = 0;
 
-  // Randomize the feature Ids for purely visual clarify. Having random Feature Ids
+  // Randomize the feature Ids purely for visual clarity. Having random Feature Ids
   // allows users visualizing the data to better discern each grain otherwise the coloring
   // would look like a smooth gradient. This is a user input parameter
   if(m_InputValues->RandomizeFeatureIds)
   {
-    ClusterUtilities::RandomizeFeatureIds(m_FeatureIdsArray->getDataStoreRef(), this->m_FoundFeatures + 1);
+    ClusterUtilities::RandomizeFeatureIds(m_FeatureIdsArray->getDataStoreRef(), m_FoundFeatures + 1);
   }
 
   return {};
@@ -111,63 +139,63 @@ Result<> CAxisSegmentFeatures::operator()()
 // -----------------------------------------------------------------------------
 int64 CAxisSegmentFeatures::getSeed(int32 gnum, int64 nextSeed) const
 {
-  DataArray<int32>::store_type& featureIds = m_FeatureIdsArray->getDataStoreRef();
-  const usize totalPoints = featureIds.getNumberOfTuples();
-  AbstractDataStore<int32>& cellPhases = m_CellPhases->getDataStoreRef();
+  AbstractDataStore<int32>& featureIdsRef = m_FeatureIdsArray->getDataStoreRef();
+  const usize totalPoints = featureIdsRef.getNumberOfTuples();
+  const AbstractDataStore<int32>& cellPhasesRef = m_CellPhases->getDataStoreRef();
 
-  // start with the next voxel after the last seed
-  auto randPoint = static_cast<usize>(nextSeed);
+  // Linearly scan for the next eligible voxel, starting just after the last seed
+  auto candidatePoint = static_cast<usize>(nextSeed);
   int64 seed = -1;
-  while(seed == -1 && randPoint < totalPoints)
+  while(seed == -1 && candidatePoint < totalPoints)
   {
-    if(featureIds[randPoint] == 0) // If the GrainId of the voxel is ZERO then we can use this as a seed point
+    if(featureIdsRef[candidatePoint] == 0) // If the GrainId of the voxel is ZERO then we can use this as a seed point
     {
-      if((!m_InputValues->UseMask || m_GoodVoxelsArray->isTrue(randPoint)) && cellPhases[randPoint] > 0)
+      if((!m_InputValues->UseMask || m_GoodVoxelsArray->isTrue(candidatePoint)) && cellPhasesRef[candidatePoint] > 0)
       {
-        seed = static_cast<int64>(randPoint);
+        seed = static_cast<int64>(candidatePoint);
       }
       else
       {
-        randPoint += 1;
+        candidatePoint += 1;
       }
     }
     else
     {
-      randPoint += 1;
+      candidatePoint += 1;
     }
   }
   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
+    // Stamp the seed only; the Feature AttributeMatrix is resized once in operator() after the
+    // segmentation completes (matching the sibling EBSD/Scalar algorithms). Nothing reads the
+    // feature-level arrays while the flood fill runs.
+    featureIdsRef[static_cast<usize>(seed)] = gnum;
   }
   return seed;
 }
 
 // -----------------------------------------------------------------------------
-bool CAxisSegmentFeatures::determineGrouping(int64 referencepoint, int64 neighborpoint, int32 gnum) const
+bool CAxisSegmentFeatures::determineGrouping(int64 referencePoint, int64 neighborPoint, int32 gnum) const
 {
   bool group = false;
 
   const Eigen::Vector3f cAxis{0.0f, 0.0f, 1.0f};
-  Float32Array& currentQuat = *m_QuatsArray;
+  const Float32Array& quats = *m_QuatsArray;
   Int32Array& featureIds = *m_FeatureIdsArray;
-  Int32Array& cellPhases = *m_CellPhases;
+  const Int32Array& cellPhases = *m_CellPhases;
 
   bool neighborPointIsGood = false;
   if(m_GoodVoxelsArray != nullptr)
   {
-    neighborPointIsGood = m_GoodVoxelsArray->isTrue(neighborpoint);
+    neighborPointIsGood = m_GoodVoxelsArray->isTrue(neighborPoint);
   }
 
-  if(featureIds[neighborpoint] == 0 && (!m_InputValues->UseMask || neighborPointIsGood))
+  if(featureIds[neighborPoint] == 0 && (!m_InputValues->UseMask || neighborPointIsGood))
   {
-    if(cellPhases[referencepoint] == cellPhases[neighborpoint])
+    if(cellPhases[referencePoint] == cellPhases[neighborPoint])
     {
-      const ebsdlib::QuatF q1(currentQuat[referencepoint * 4], currentQuat[referencepoint * 4 + 1], currentQuat[referencepoint * 4 + 2], currentQuat[referencepoint * 4 + 3]);
-      const ebsdlib::QuatF q2(currentQuat[neighborpoint * 4 + 0], currentQuat[neighborpoint * 4 + 1], currentQuat[neighborpoint * 4 + 2], currentQuat[neighborpoint * 4 + 3]);
+      const ebsdlib::QuatF q1(quats[referencePoint * 4], quats[referencePoint * 4 + 1], quats[referencePoint * 4 + 2], quats[referencePoint * 4 + 3]);
+      const ebsdlib::QuatF q2(quats[neighborPoint * 4 + 0], quats[neighborPoint * 4 + 1], quats[neighborPoint * 4 + 2], quats[neighborPoint * 4 + 3]);
 
       const ebsdlib::OrientationMatrixFType oMatrix1 = q1.toOrientationMatrix();
       const ebsdlib::OrientationMatrixFType oMatrix2 = q2.toOrientationMatrix();
@@ -187,7 +215,7 @@ bool CAxisSegmentFeatures::determineGrouping(int64 referencepoint, int64 neighbo
       if(w <= m_InputValues->MisorientationTolerance || (Constants::k_PiD - w) <= m_InputValues->MisorientationTolerance)
       {
         group = true;
-        featureIds[neighborpoint] = gnum;
+        featureIds[neighborPoint] = gnum;
       }
     }
   }

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

@@ -29,9 +29,9 @@ struct ORIENTATIONANALYSIS_EXPORT CAxisSegmentFeaturesInputValues
 
 /**
  * @class CAxisSegmentFeatures
- * @brief This filter segments the Features by grouping neighboring Cells that satisfy the C-axis misalignment tolerance, i.e., have misalignment angle less than the value set by the user.
+ * @brief Segmentation algorithm that groups neighboring Cells whose c-axis misalignment angle is within the user tolerance. Provides the seeding and grouping hooks for the shared
+ * SegmentFeatures flood-fill driver.
  */
-
 class ORIENTATIONANALYSIS_EXPORT CAxisSegmentFeatures : public SegmentFeatures
 {
 public:

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

@@ -44,11 +44,21 @@ Result<> EBSDSegmentFeatures::operator()()
   m_FeatureIdsArray->fill(0); // initialize the output array with zeros
 
   // Run the segmentation algorithm
-  execute(gridGeom);
+  Result<> segmentResult = execute(gridGeom);
+  if(segmentResult.invalid())
+  {
+    return segmentResult;
+  }
+  // A canceled run returns successfully with no feature count; do not misreport it below as
+  // "no Features were detected".
+  if(m_ShouldCancel)
+  {
+    return {};
+  }
   // Sanity check the result.
   if(this->m_FoundFeatures < 1)
   {
-    return MakeErrorResult(-87000, fmt::format("The number of Features is '{}' which means no Features were detected. A threshold value may be set incorrectly", this->m_FoundFeatures));
+    return MakeErrorResult(-87000, "No Features were detected: no Cell was eligible to seed a Feature. Every Cell is either excluded by the Mask or has a Phase value of 0 (unindexed).");
   }
 
   // Resize the Feature Attribute Matrix