VV: Compute Feature Face Misorientations has been V&V'ed

src/Plugins/OrientationAnalysis/docs/BadDataNeighborOrientationCheckFilter.md

@@ -35,12 +35,18 @@ changing voxels.
 ### 2D Versus 3D Note
 
 If the user is processing a 2D data set, **none** of the voxels can have 6 neighbors
-since there are no neighbors is the +-Z directions.
+since there are no neighbors in the +/-Z directions.
 
 ### Warning - Data Modification
 
 Only the *Mask* value defining the cell as *good* or *bad* is changed. No other cell level array is modified.
 
+### Memory Considerations
+
+The filter allocates a temporary `int32` neighbor-count array sized to the total voxel count
+(4 bytes per voxel). For a 1-billion-voxel dataset, that is approximately 4 GB of additional
+working memory during execution. This memory is released when the filter finishes.
+
 ## Example Data
 
 | Example Input Image                                                       | Example Output Image                                                                                                                          |
@@ -55,8 +61,12 @@ From the above before and after images you can see that this filter can help mod
 ## Example Pipelines
 
 + (02) Small IN100 Full Reconstruction
-+ INL Export
-+ 04_Steiner Compact
+
+## Related Filters
+
+- [Fill Bad Data](../SimplnxCore/FillBadDataFilter.md) — fills voxels still marked bad after this filter runs (or as a standalone alternative when no orientation data is available).
+- [Multi-Threshold Objects](../SimplnxCore/MultiThresholdObjectsFilter.md) — typical upstream filter that generates the initial *Mask* array (e.g., from `Confidence Index` and `Image Quality`).
+- [Replace Element Attributes with Neighbor Values](../SimplnxCore/ReplaceElementAttributesWithNeighborValuesFilter.md) — alternative cleanup approach that copies attribute values from neighboring cells rather than flipping a mask.
 
 ## License & Copyright
 

src/Plugins/OrientationAnalysis/docs/ComputeFeatureFaceMisorientationFilter.md

@@ -6,7 +6,13 @@ Processing (Crystallography)
 
 ## Description
 
-This **Filter** generates a 3 component vector for each **Triangle** in a **Triangle Geometry** that is the axis-angle of the misorientation associated with the **Features** that lie on either side of the **Triangle**.  The axis is normalized, so if the magnitude of the vector is viewed, it will be the *misorientation angle* in degrees.
+This **Filter** generates a 3 component vector for each **Triangle** in a **Triangle Geometry** that is the axis-angle of the misorientation associated with the **Features** that lie on either side of the **Triangle**. The axis is normalized, so if the magnitude of the vector is viewed, it will be the *misorientation angle* in degrees. This filter works on all valid crystal-structures/laue classes.
+
+If you want to get the "raw" (un-normalized) axis-angle of the misorientation, enable "Store Full Axis Angle" and then (optionally) modify the "Axis Angle Array Name" parameter. The actual axis-angle is stored in a 4 component DataArray with the format as follows:
+
+```text
+{{x_axis, y_axis, z_axis, angle}, ...}
+```
 
 % Auto generated parameter table will be inserted here
 

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

@@ -33,10 +33,15 @@ const std::atomic_bool& BadDataNeighborOrientationCheck::getCancel()
 // -----------------------------------------------------------------------------
 Result<> BadDataNeighborOrientationCheck::operator()()
 {
-  const float misorientationTolerance = m_InputValues->MisorientationTolerance * numbers::pi_v<float> / 180.0f;
+  // Compute the tolerance in double precision: numbers::pi_v<float> is the closest float to true pi, which is
+  // slightly *larger* than true pi; converting via float makes the radian tolerance ~5e-9 rad larger than the
+  // mathematically true k*pi/180. For boundary-exact misorientations (e.g., test fixtures landing on exactly the
+  // user-supplied tolerance), the float-converted tolerance can incorrectly include cases that should fail strict <.
+  // Using double-pi makes the conversion faithful and the strict < tolerance comparison match the analytical oracle.
+  const double misorientationTolerance = static_cast<double>(m_InputValues->MisorientationTolerance) * numbers::pi_v<double> / 180.0;
 
-  const auto* imageGeomPtr = m_DataStructure.getDataAs<ImageGeom>(m_InputValues->ImageGeomPath);
-  SizeVec3 udims = imageGeomPtr->getDimensions();
+  const auto& imageGeom = m_DataStructure.getDataRefAs<ImageGeom>(m_InputValues->ImageGeomPath);
+  SizeVec3 udims = imageGeom.getDimensions();
   const auto& cellPhases = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->CellPhasesArrayPath);
   const auto& quats = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->QuatsArrayPath);
   const auto& crystalStructures = m_DataStructure.getDataRefAs<UInt32Array>(m_InputValues->CrystalStructuresArrayPath);
@@ -59,32 +64,66 @@ Result<> BadDataNeighborOrientationCheck::operator()()
       static_cast<int64>(udims[2]),
   };
 
+  // VoxelNeighbors<Image3D>::k_FaceNeighborCount = 6 is the maximum possible face-neighbor count.
+  // computeValidFaceNeighbors() runtime-skips +/-Z neighbors when dims[2] == 1 (2D images), so this
+  // 3D-typed array correctly handles 2D images without any change here.
   constexpr FaceNeighborType k_NumFaceNeighbors = VoxelNeighbors<Image3D>::k_FaceNeighborCount;
   const std::array<int64, k_NumFaceNeighbors> neighborVoxelIndexOffsets = initializeFaceNeighborOffsets(dims);
   constexpr std::array<FaceNeighborType, k_NumFaceNeighbors> faceNeighborInternalIdx = initializeFaceNeighborInternalIdx();
 
   const std::vector<ebsdlib::LaueOps::Pointer> orientationOps = ebsdlib::LaueOps::GetAllOrientationOps();
 
+  // Validate that every entry in the CrystalStructures ensemble array is a valid Laue-group index
+  // (< orientationOps.size()). Catches malformed inputs such as a legacy CreateEnsembleInfo sentinel
+  // value (999) at ensemble index 0 before they cause an out-of-bounds dereference in the per-voxel
+  // loop below. The UnknownCrystalStructure value is explicitly allowed as a sentinel; voxels whose
+  // phase resolves to it will be skipped by the cellPhases > 0 guard. CrystalStructures is typically
+  // tiny (2-4 entries), so the cost is negligible.
+  const usize numOrientationOps = orientationOps.size();
+  for(usize i = 0; i < crystalStructures.getSize(); ++i)
+  {
+    if(crystalStructures[i] >= numOrientationOps && crystalStructures[i] != static_cast<uint32>(ebsdlib::CrystalStructure::UnknownCrystalStructure))
+    {
+      return MakeErrorResult(
+          -54901, fmt::format("Crystal structure at ensemble index {} has value {}, which is not a valid Laue-group index. Valid range is [0, {}).", i, crystalStructures[i], numOrientationOps));
+    }
+  }
+
+  // Per-voxel running count of within-tolerance face-neighbors. Allocated proportional to the
+  // input geometry size: 4 bytes per voxel (~4 GB for a 1B-voxel dataset). Cannot be in-place
+  // on the mask array because the algorithm needs to distinguish "newly flipped" from "still bad".
   std::vector<int32> neighborCount(totalPoints, 0);
 
   MessageHelper messageHelper(m_MessageHandler);
   ThrottledMessenger throttledMessenger = messageHelper.createThrottledMessenger();
   // Loop over every point finding the number of neighbors that fall within the
   // user defined angle tolerance.
-  for(int64 voxelIndex = 0; voxelIndex < totalPoints; voxelIndex++)
+  for(usize voxelIndex = 0; voxelIndex < totalPoints; voxelIndex++)
   {
+    if(m_ShouldCancel)
+    {
+      return {};
+    }
     throttledMessenger.sendThrottledMessage([&] { return fmt::format("Processing Data {:.2f}% completed", CalculatePercentComplete(voxelIndex, totalPoints)); });
     // If the mask was set to false, then we check this voxel
     if(!maskCompare->isTrue(voxelIndex))
     {
       // We precalculate the positive voxel quaternion and laue class here to prevent reading and recalculating it for each face below
       ebsdlib::QuatD quat1(quats[voxelIndex * 4], quats[voxelIndex * 4 + 1], quats[voxelIndex * 4 + 2], quats[voxelIndex * 4 + 3]);
       quat1.positiveOrientation();
-      const uint32 laueClass1 = crystalStructures[cellPhases[voxelIndex]];
+      const uint32 laueClassIndex = crystalStructures[cellPhases[voxelIndex]];
+      // Defensive: skip voxels whose phase resolves to an out-of-range Laue index (e.g., the
+      // UnknownCrystalStructure sentinel allowed by the validation above). Without this, the
+      // orientationOps[laueClassIndex] dereference below would be out-of-bounds.
+      if(laueClassIndex >= numOrientationOps)
+      {
+        continue;
+      }
 
-      int64 xIdx = voxelIndex % dims[0];
-      int64 yIdx = (voxelIndex / dims[0]) % dims[1];
-      int64 zIdx = voxelIndex / (dims[0] * dims[1]);
+      const int64 voxelIndexI64 = static_cast<int64>(voxelIndex);
+      int64 xIdx = voxelIndexI64 % dims[0];
+      int64 yIdx = (voxelIndexI64 / dims[0]) % dims[1];
+      int64 zIdx = voxelIndexI64 / (dims[0] * dims[1]);
 
       // Loop over the 6 face neighbors of the voxel
       const std::array<bool, k_NumFaceNeighbors> isValidFaceNeighbor = computeValidFaceNeighbors(xIdx, yIdx, zIdx, dims);
@@ -94,7 +133,7 @@ Result<> BadDataNeighborOrientationCheck::operator()()
         {
           continue;
         }
-        const int64 neighborPoint = voxelIndex + neighborVoxelIndexOffsets[faceIndex];
+        const int64 neighborPoint = voxelIndexI64 + neighborVoxelIndexOffsets[faceIndex];
 
         // Now compare the mask of the neighbor. If the mask is TRUE, i.e., that voxel
         // did not fail the threshold filter that most likely produced the mask array,
@@ -107,7 +146,7 @@ Result<> BadDataNeighborOrientationCheck::operator()()
             ebsdlib::QuatD quat2(quats[neighborPoint * 4], quats[neighborPoint * 4 + 1], quats[neighborPoint * 4 + 2], quats[neighborPoint * 4 + 3]);
             quat2.positiveOrientation();
             // Compute the Axis_Angle misorientation between those 2 quaternions
-            ebsdlib::AxisAngleDType axisAngle = orientationOps[laueClass1]->calculateMisorientation(quat1, quat2);
+            ebsdlib::AxisAngleDType axisAngle = orientationOps[laueClassIndex]->calculateMisorientation(quat1, quat2);
             // if the angle is less than our tolerance, then we increment the neighbor count
             // for this voxel
             if(axisAngle[3] < misorientationTolerance)
@@ -120,49 +159,66 @@ Result<> BadDataNeighborOrientationCheck::operator()()
     }
   }
 
-  constexpr int32 startLevel = 6;
+  // Convergence loop starts at the maximum possible face-neighbor count (6 in 3D; 2D images
+  // simply never reach the top levels because no voxel can have count > 4). Tying this to
+  // k_NumFaceNeighbors keeps the upper bound consistent if VoxelNeighbors ever changes.
+  constexpr int32 startLevel = static_cast<int32>(k_NumFaceNeighbors);
   int32 currentLevel = startLevel;
   int32 counter = 0;
 
   // Now we loop over all the points again, but this time we do it as many times
   // as the user has requested to iteratively flip voxels
   while(currentLevel >= m_InputValues->NumberOfNeighbors)
   {
+    if(m_ShouldCancel)
+    {
+      return {};
+    }
     counter = 1;
     int32 loopNumber = 0;
     while(counter > 0)
     {
+      if(m_ShouldCancel)
+      {
+        return {};
+      }
       counter = 0; // Set this while control variable to zero
       for(usize voxelIndex = 0; voxelIndex < totalPoints; voxelIndex++)
       {
+        if(m_ShouldCancel)
+        {
+          return {};
+        }
         throttledMessenger.sendThrottledMessage([&] {
           return fmt::format("Level '{}' of '{}' || Processing Data ('{}') {:.2f}% completed", (startLevel - currentLevel) + 1, startLevel - m_InputValues->NumberOfNeighbors, loopNumber,
                              CalculatePercentComplete(voxelIndex, totalPoints));
         });
 
-        // We are comparing the number-of-neighbors of the current voxel, and if it
-        // is > the current level and the mask is FALSE, then we drop into this
-        // conditional. The first thing that happens in the conditional is that
-        // the current voxel's mask value is set to TRUE.
+        // If the current voxel's neighbor count is >= the current level and the mask is FALSE,
+        // we flip the voxel to TRUE and recompute its (still-bad) neighbors' counts below.
         if(neighborCount[voxelIndex] >= currentLevel && !maskCompare->isTrue(voxelIndex))
         {
-          maskCompare->setValue(voxelIndex, true); // the current voxel's mask value is set to TRUE.
-          counter++;                               // Increment the `counter` to force the loop to iterate again
+          maskCompare->setValue(voxelIndex, true);
+          counter++; // Increment the `counter` to force the loop to iterate again
 
           // We precalculate the positive voxel quaternion and laue class here to prevent reading and recalculating it for each face below
           ebsdlib::QuatD quat1(quats[voxelIndex * 4], quats[voxelIndex * 4 + 1], quats[voxelIndex * 4 + 2], quats[voxelIndex * 4 + 3]);
           quat1.positiveOrientation();
-          const uint32 laueClass1 = crystalStructures[cellPhases[voxelIndex]];
-
-          // This whole section below is to now look at the neighbor voxels of the
-          // current voxel that just got flipped to true. This is needed because
-          // if any of those neighbor's mask was `false`, then its neighbor count
-          // is now not correct and will be off-by-one. So we run _almost_ the same
-          // loop code as above but checking the specific neighbors of the current
-          // voxel. This part should be termed the "Update Neighbor's Neighbor Count"
-          int64 xIdx = voxelIndex % dims[0];
-          int64 yIdx = (voxelIndex / dims[0]) % dims[1];
-          int64 zIdx = voxelIndex / (dims[0] * dims[1]);
+          const uint32 laueClassIndex = crystalStructures[cellPhases[voxelIndex]];
+          // Defensive: skip voxels with out-of-range Laue index. See matching guard in pass 1.
+          if(laueClassIndex >= numOrientationOps)
+          {
+            continue;
+          }
+
+          // "Update Neighbor's Neighbor Count" pass: now that the current voxel just flipped to
+          // true, every still-bad face neighbor must have its neighborCount incremented by 1 if
+          // its misorientation to the freshly-flipped voxel is within tolerance. Skipping this
+          // update would leave the neighbor counts stale and prevent valid cascade flips later.
+          const int64 voxelIndexI64 = static_cast<int64>(voxelIndex);
+          int64 xIdx = voxelIndexI64 % dims[0];
+          int64 yIdx = (voxelIndexI64 / dims[0]) % dims[1];
+          int64 zIdx = voxelIndexI64 / (dims[0] * dims[1]);
 
           // Loop over the 6 face neighbors of the voxel
           const std::array<bool, k_NumFaceNeighbors> isValidFaceNeighbor = computeValidFaceNeighbors(xIdx, yIdx, zIdx, dims);
@@ -173,9 +229,9 @@ Result<> BadDataNeighborOrientationCheck::operator()()
               continue;
             }
 
-            int64 neighborPoint = voxelIndex + neighborVoxelIndexOffsets[faceIndex];
+            const int64 neighborPoint = voxelIndexI64 + neighborVoxelIndexOffsets[faceIndex];
 
-            // If the neighbor voxel's mask is false, then ....
+            // If the neighbor voxel's mask is false, then compute misorientation angle
             if(!maskCompare->isTrue(neighborPoint))
             {
               // Make sure both cells phase values are identical and valid
@@ -184,7 +240,7 @@ Result<> BadDataNeighborOrientationCheck::operator()()
                 ebsdlib::QuatD quat2(quats[neighborPoint * 4], quats[neighborPoint * 4 + 1], quats[neighborPoint * 4 + 2], quats[neighborPoint * 4 + 3]);
                 quat2.positiveOrientation();
                 // Quaternion Math is not commutative so do not reorder
-                ebsdlib::AxisAngleDType axisAngle = orientationOps[laueClass1]->calculateMisorientation(quat1, quat2);
+                ebsdlib::AxisAngleDType axisAngle = orientationOps[laueClassIndex]->calculateMisorientation(quat1, quat2);
                 if(axisAngle[3] < misorientationTolerance)
                 {
                   neighborCount[neighborPoint]++;