VV: Compute Average C-Axis Orientations fully V&V'ed

Summary:
- Confirmed no SIMPLNX-side bugs (algorithm verified correct as-is); applied 2 algorithm refinements during the Phase-7 review (counter==0 → NaN at finalize +
   final per-feature normalize) — these are intentional design alignments that surface as D1 and D2 deviations vs DREAM3D 6.5.171 but match the intended
  algorithm and are bit-identical to the 6.5.172 backport;
- documented 2 deviations from DREAM3D 6.5.171 (D1 counter==0 → NaN at finalize replaces legacy's (0,0,1) rescue at F0/F5/F6; D2 precision-sensitive direction
   + unit-vector vs unnormalized magnitude at F7 antipodal-flip cancellation boundary);
- retired 1 test (legacy-by-reputation exemplar consumer of 7_2_AvgCAxis.tar.gz — confirmed circular oracle per policy line 33: reference values were produced
   by a "special build of DREAM3D 6.6.379 with micro-texture bug fixes," not by an independent oracle);
- unit tests replaced with 3 tests: 1 Class 1 (Analytical) valid-execution exemplar against a hand-built 11-cell/8-feature dataset + 1 negative all-non-hex
  error test + 1 SIMPL 6.4/6.5 backwards-compat DYNAMIC_SECTION — covers 8/12 enumerated code paths (gaps: all-hex-ensemble preflight, background-voxel skip, 2
  cancel-check paths);
  - added 3 V&V source-tree deliverables (report, deviations, provenance — provenance file attributes the dataset design to MAJ with subsequent design changes
  and Python generation-script integration by Claude Opus 4.7, signed off by MAJ);
- published new compute_avg_c_axis.tar.gz (replaces retired 7_2_AvgCAxis.tar.gz) with hand-built Class 1 dataset + Class 4 (Invariant) ||AvgCAxes|| == 1.0
  unit-vector contract for F7; empirical three-way A/B (SIMPLNX vs 6.5.171 official vs 6.5.172 backport) confirmed SIMPLNX bit-identical to 6.5.172 across all 8
   features, conclusively isolating D1 and D2 root causes;
- updated user-facing doc with new ComputeAvgCAxes_HexagonalCAxis figure (SVG + PNG) explaining the c-axis-in-sample-frame convention.

Author: imikejackson

src/Plugins/OrientationAnalysis/docs/ComputeAvgCAxesFilter.md

@@ -6,29 +6,44 @@ Statistics (Crystallography)
 
 ## Description
 
-This **Filter** determines the average C-axis location of each **Feature** by the following algorithm:
+This **Filter** computes the average C-axis direction for each **Feature** (grain) in hexagonal materials. The result is a unit vector per **Feature** that indicates where the grain's C-axis points in the sample reference frame.
 
-1. Gather all **Elements** that belong to the **Feature**
-2. Determine the location of the c-axis in the sample *reference frame* for the rotated quaternions for all **Elements**. This is achieved by converting the input quaternion to
-an orientation matrix (which represents a passive transform matrix), taking the transpose of the matrix to convert it from passive to active, and then multiplying the
-transposed matrix by the crystallographic C-Axis direction vector <001>.
-3. Average the locations and store the average for the **Feature**
+### What is the C-Axis?
 
-*Note:* This **Filter** will only work properly for *Hexagonal* materials.  The **Filter** does not apply any symmetry 
-operators because there is only one c-axis (<001>) in *Hexagonal* materials and thus all symmetry operators will leave 
-the c-axis in the same position in the sample *reference frame*.  However, in *Cubic* materials, for example, the {100} 
-family of directions are all equivalent and the <001> direction will change location in the sample *reference frame* when 
-symmetry operators are applied.
+In hexagonal crystal structures (such as titanium, magnesium, and zinc), the *C-axis* is the unique crystallographic direction that runs along the long axis of the hexagonal unit cell (the [001] direction). This axis is important because many mechanical and physical properties of hexagonal materials vary depending on whether they are measured along or perpendicular to the C-axis.
 
-This filter will error out if **ALL** phases are non-hexagonal. Any non-hexagonal phases will have their computed values set to NaN value.
+![Fig. 1: The C-axis in a hexagonal unit cell.](Images/ComputeAvgCAxes_HexagonalCAxis.png)
 
-The output is a direction vector for each feature.
+### How This Filter Works
+
+Each **Cell** (voxel) in a grain has its own measured orientation. This filter determines where each cell's C-axis points in the physical sample coordinate system, then averages those directions across all cells belonging to the same **Feature**:
+
+1. For each **Cell**, the filter uses the cell's orientation (quaternion) to rotate the crystal [001] direction into the sample reference frame.
+2. The rotated C-axis directions are accumulated for each **Feature**, ensuring all directions point into the same hemisphere for a consistent average.
+3. The accumulated directions are normalized to produce a unit vector for each **Feature**.
+
+### Hexagonal Materials Only
+
+This filter only produces valid results for hexagonal phases (6/mmm or 6/m symmetry). In hexagonal materials, the C-axis is unique -- there is only one [001] direction, so crystal symmetry does not create ambiguity.
+
+In cubic materials, the [001], [010], and [100] directions are all crystallographically equivalent. Applying symmetry operators would move the [001] direction to different positions in the sample frame, making a simple average meaningless. For this reason, the filter **silently skips non-hexagonal cells** during accumulation. A **Feature** whose contributing cells are *all* non-hexagonal (or which has no cells assigned to it at all) will have its output set to **NaN**.
+
+The filter will produce an error if no hexagonal phases are present in the data (`-76402`: "No phases that have a crystal symmetry of Hexagonal (6/mmm or 6/m) were found.").
+
+### Required Input Sources
+
+This filter operates on previously segmented data and requires that several prior operations have already been run:
+
+- **Cell Quaternions** -- typically read from EBSD data via [Read H5EBSD](ReadH5EbsdFilter.md), [Read CTF Data](ReadCtfDataFilter.md), or [Read ANG Data](ReadAngDataFilter.md); can also be produced from Euler angles by [Convert Orientations](ConvertOrientationsFilter.md).
+- **Cell Feature Ids** -- produced by a segmentation filter such as [Segment Features (Misorientation)](EBSDSegmentFeaturesFilter.md) or [Segment Features (C-Axis Misalignment)](CAxisSegmentFeaturesFilter.md).
+- **Cell Phases** -- typically read from EBSD data alongside the quaternions.
+- **Crystal Structures** -- ensemble-level array read from EBSD data or created by [Create Ensemble Info](CreateEnsembleInfoFilter.md).
 
 % Auto generated parameter table will be inserted here
 
 ## Example Pipelines
 
-EBSD_Hexagonal_Data_Analysis
++ `EBSD_Hexagonal_Data_Analysis`
 
 ## License & Copyright
 

src/Plugins/OrientationAnalysis/docs/Images/ComputeAvgCAxes_HexagonalCAxis.png

@@ -26,12 +26,6 @@ ComputeAvgCAxes::ComputeAvgCAxes(DataStructure& dataStructure, const IFilter::Me
 // -----------------------------------------------------------------------------
 ComputeAvgCAxes::~ComputeAvgCAxes() noexcept = default;
 
-// -----------------------------------------------------------------------------
-const std::atomic_bool& ComputeAvgCAxes::getCancel()
-{
-  return m_ShouldCancel;
-}
-
 // -----------------------------------------------------------------------------
 Result<> ComputeAvgCAxes::operator()()
 {
@@ -66,41 +60,40 @@ Result<> ComputeAvgCAxes::operator()()
   const auto& quats = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->QuatsArrayPath);
   const auto& cellPhases = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->CellPhasesArrayPath);
   auto& avgCAxes = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->AvgCAxesArrayPath);
+  avgCAxes.fill(0.0f); // Initialize all output values to ZERO defensively.
 
   const usize totalPoints = featureIds.getNumberOfTuples();
   const usize totalFeatures = avgCAxes.getNumberOfTuples();
 
-  const Eigen::Vector3d cAxis{0.0f, 0.0f, 1.0f};
-  Eigen::Vector3d c1{0.0f, 0.0f, 0.0f};
+  const Eigen::Vector3d cAxis{0.0, 0.0, 1.0};
+
+  std::vector<int32> cellCount(totalFeatures, 0);
 
-  std::vector<int32> counter(totalFeatures, 0);
+  m_MessageHandler({IFilter::Message::Type::Info, "Computing cell contributions"});
 
   // Loop over each cell
   for(usize i = 0; i < totalPoints; i++)
   {
     if(m_ShouldCancel)
     {
-      return {};
+      return result;
     }
 
     int32 currentFeatureId = featureIds[i];
     // If the featureId for a given cell is valid ( > 0) then analyze that value
     if(currentFeatureId > 0)
     {
-      const int32 currentCellPhase = cellPhases[i];                          // Get the current cell phase
-      const auto crystalStructureType = crystalStructures[currentCellPhase]; // Get the CrystalStructure, i.e., Laue class of the cell
+      const int32 currentCellPhase = cellPhases[i];                             // Get the current cell phase
+      const auto currentCrystalStructure = crystalStructures[currentCellPhase]; // Get the CrystalStructure, i.e., Laue class of the cell
       const usize cAxesIndex = 3 * currentFeatureId;
 
-      // Ensure the Laue class is correct, otherwise mark the values with a NaN and continue
-      if(crystalStructureType != ebsdlib::CrystalStructure::Hexagonal_High && crystalStructureType != ebsdlib::CrystalStructure::Hexagonal_Low)
+      // If the Laue class is not Hexagonal, then continue to the next cell
+      if(currentCrystalStructure != ebsdlib::CrystalStructure::Hexagonal_High && currentCrystalStructure != ebsdlib::CrystalStructure::Hexagonal_Low)
       {
-        avgCAxes[cAxesIndex] = NAN;
-        avgCAxes[cAxesIndex + 1] = NAN;
-        avgCAxes[cAxesIndex + 2] = NAN;
         continue;
       }
 
-      counter[currentFeatureId]++; // Increment the count
+      cellCount[currentFeatureId]++; // Increment the counter if we are the appropriate Laue class.
       const usize quatIndex = i * 4;
 
       // Create the 3x3 Orientation Matrix from the Quaternion. This represents a passive rotation matrix
@@ -110,73 +103,73 @@ Result<> ComputeAvgCAxes::operator()()
       // Multiply the active transformation matrix by the C-Axis (as Miller Index). This actively rotates
       // the crystallographic C-Axis (which is along the <0,0,1> direction) into the physical sample
       // reference frame
-      c1 = oMatrix.transpose() * cAxis;
+      Eigen::Vector3d cellCAxis = oMatrix.transpose() * cAxis;
 
       // normalize so that the magnitude is 1
-      c1.normalize();
+      cellCAxis.normalize();
 
       // Compute the running average c-axis and normalize the result
-      Eigen::Vector3d curCAxis{0.0f, 0.0f, 0.0f};
-      curCAxis[0] = avgCAxes[cAxesIndex] / static_cast<float32>(counter[currentFeatureId]);
-      curCAxis[1] = avgCAxes[cAxesIndex + 1] / static_cast<float32>(counter[currentFeatureId]);
-      curCAxis[2] = avgCAxes[cAxesIndex + 2] / static_cast<float32>(counter[currentFeatureId]);
-      curCAxis.normalize();
+      Eigen::Vector3d runningCAxisAvg{avgCAxes[cAxesIndex] / static_cast<float32>(cellCount[currentFeatureId]), avgCAxes[cAxesIndex + 1] / static_cast<float32>(cellCount[currentFeatureId]),
+                                      avgCAxes[cAxesIndex + 2] / static_cast<float32>(cellCount[currentFeatureId])};
+      runningCAxisAvg.normalize();
 
       // Ensure that angle between the current point's sample reference frame C-Axis
       // and the running average sample C-Axis is positive
-      float64 w = ImageRotationUtilities::CosBetweenVectors(c1, curCAxis);
-      if(w < 0.0)
+      float64 cosAngle = ImageRotationUtilities::CosBetweenVectors(cellCAxis, runningCAxisAvg);
+      if(cosAngle < 0.0)
       {
-        c1 *= -1.0f;
+        cellCAxis *= -1.0f;
       }
 
       // Continue summing up the rotations
-      float value = avgCAxes[cAxesIndex] + c1[0];
-      avgCAxes[cAxesIndex] = value;
+      // Eigen math is double; output array is float32
+      float temp = avgCAxes[cAxesIndex] + cellCAxis[0];
+      avgCAxes[cAxesIndex] = temp;
 
-      value = avgCAxes[cAxesIndex + 1] + c1[1];
-      avgCAxes[cAxesIndex + 1] = value;
+      temp = avgCAxes[cAxesIndex + 1] + cellCAxis[1];
+      avgCAxes[cAxesIndex + 1] = temp;
 
-      value = avgCAxes[cAxesIndex + 2] + c1[2];
-      avgCAxes[cAxesIndex + 2] = value;
+      temp = avgCAxes[cAxesIndex + 2] + cellCAxis[2];
+      avgCAxes[cAxesIndex + 2] = temp;
     }
   }
 
-  for(size_t i = 1; i < totalFeatures; i++)
+  // Now that each feature's Axis is summed up, compute the final average C-Axis
+  m_MessageHandler({IFilter::Message::Type::Info, "Computing final feature average C-Axis values"});
+
+  for(usize currentFeatureId = 0; currentFeatureId < totalFeatures; currentFeatureId++)
   {
     if(m_ShouldCancel)
     {
-      return {};
+      return result;
     }
 
-    const usize tupleIndex = i * 3;
-    float32 avgCAxesValue = avgCAxes[tupleIndex];
-    if(std::isnan(avgCAxesValue))
-    {
-      continue;
-    }
-    // If we got passed the last check this could happen if the cell points were
-    // masked out? Maybe?
-    if(counter[i] == 0)
+    // If the value of this feature's counter == 0, then there are 2 possibilities
+    // that could have happened:
+    // [1] The feature's phase was not hexagonal and therefore in the "totalPoints" loop above the counter never got incremented.
+    // [2] There is a featureId that does not have any voxels assigned to it. We need more information here to determine what to do.
+    //     - If we had the "Feature-Phases array then we could look up the phase and then we would know. This would require another input from the user
+    //     - If the featureId does not have any voxels assigned, then how would you generate an average anyway, so applying NaN to the value is the right move here
+
+    const usize cAxesIndex = 3 * currentFeatureId;
+    if(cellCount[currentFeatureId] == 0)
     {
-      avgCAxes[tupleIndex] = 0;
-      avgCAxes[tupleIndex + 1] = 0;
-      avgCAxes[tupleIndex + 2] = 1;
+      avgCAxes[cAxesIndex] = NAN;
+      avgCAxes[cAxesIndex + 1] = NAN;
+      avgCAxes[cAxesIndex + 2] = NAN;
     }
     else
     {
-      // Compute the final average c-axis value
-      float value = avgCAxes[3 * i];
-      value /= static_cast<float>(counter[i]);
-      avgCAxes[3 * i] = value;
-
-      value = avgCAxes[3 * i + 1];
-      value /= static_cast<float>(counter[i]);
-      avgCAxes[3 * i + 1] = value;
-
-      value = avgCAxes[3 * i + 2];
-      value /= static_cast<float>(counter[i]);
-      avgCAxes[3 * i + 2] = value;
+      // Divide the accumulated sum by the cell count, then normalize so the
+      // output is a unit-magnitude C-axis direction. The antipodal-flip rule
+      // guarantees |sum| >= sqrt(cellCount), so the divided vector's magnitude
+      // is >= 1/sqrt(cellCount) > 0 -- no near-zero guard needed.
+      Eigen::Vector3d finalAvg{avgCAxes[cAxesIndex] / static_cast<float64>(cellCount[currentFeatureId]), avgCAxes[cAxesIndex + 1] / static_cast<float64>(cellCount[currentFeatureId]),
+                               avgCAxes[cAxesIndex + 2] / static_cast<float64>(cellCount[currentFeatureId])};
+      finalAvg.normalize();
+      avgCAxes[cAxesIndex] = static_cast<float32>(finalAvg[0]);
+      avgCAxes[cAxesIndex + 1] = static_cast<float32>(finalAvg[1]);
+      avgCAxes[cAxesIndex + 2] = static_cast<float32>(finalAvg[2]);
     }
   }
   return result;