ENH: Microtexture related filter cleanup (#1438)

* API: Update to latest EbsdLib API changes

* BUG: CAxisSegmentFeatures - Ensure all voxels that are segmented have a hexagonal phase type.
* BUG: ComputeAvgOrientations fix incorrect computation.
* BUG: ConvertOrientations * make sure the quaternion is properly formed before conversion
* BUG: Edax,Bruker,Oxford H5 EBSD Readers should check for null PatternData pointer
* BUG: FindFeatureReferenceCAxisOrientation - Use doubles to accumulate the StdDev values. Output values now agree with DREAM.3D > 6.5.172
* BUG: Fixes ComputeFeatureNeighborCAxisAlignments crash if "Find Avg Misalignments" was not enabled.
* BUG: Fixes incorrect progress message in ComputeNeighborhoods filter
* BUG: ReadH5Ebsd: EulerAngles must be read to perform the Euler Ref frame transform
* BUG: ReadH5OimData-Fix index calculation when copying data into final arrays
* Compute Average C-Axis was updated and documented with comments

* DOC: Add Point Group and Rotation Point Group to the Laue class table
* DOC: Updates CreateEnsembleInfo documentation to add in the Rotation Point Groups

* ENH: Adds function to create a randomized "index" array that is used elsewhere
* ENH: Fixes SegmentFeatures algorithm to use a throttled messenger for progress updates
* ENH: Fixes small bugs related to running Microtexture pipelines
* ENH: Remove warning about hex phases for CAxisSegmentFeatures

* REV: Compute Feature Neighbor Code Review
* REV: Compute Feature Reference CAxis Misorientations has been reviewed.
* REV: ComputeFeatureReferenceCAxisMisorientations: Document and small code refactorings

* STY: Use proper include syntax for EbsdLib since it is an external library

* TEST: Update ComputeAvgOrientations Unit Test. This is really ONLY testing Cubic High Laue group.
* TEST: Update ComputeFZQuaternions Unit Test to allow for a tolerance for results

Author: imikejackson

src/Plugins/OrientationAnalysis/docs/ComputeAvgCAxesFilter.md

@@ -10,14 +10,20 @@ This **Filter** determines the average C-axis location of each **Feature** by th
 
 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
-and orientation matrix (which represents a passive transform matrix), taking the transpose of the matrix to convert it from passive to active, and thne multiplying the
+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 as the average for the **Feature**
+3. Average the locations and store the average for the **Feature**
 
-*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.
+*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.
 
 This filter will error out if **ALL** phases are non-hexagonal. Any non-hexagonal phases will have their computed values set to NaN value.
 
+The output is a direction vector for each feature.
+
 % Auto generated parameter table will be inserted here
 
 ## Example Pipelines

src/Plugins/OrientationAnalysis/docs/CreateEnsembleInfoFilter.md

@@ -8,22 +8,22 @@ Processing (Generation)
 
 This **Filter** allows the user to enter basic crystallographic information about each phase. The Laue class, Phase Type, and Phase Name can all be entered by the user. The information is stored in an EnsembleAttributeMatrix. These values are needed to allow the calculation of certain kinds of crystallographic statistics on the volume, if they have not already been provided by some other means. Each row in the table lists the **Crystal Structure**, **Phase Type**, and **Phase Name**. The proper values for the crystal structure and phase type come from internal constants within DREAM3D-NX and are listed here:
 
-### Crystal Structure
-
-| String Name | Internal Value | Laue Name |
-| ------------|----------------|----------|
-| Hexagonal_High | 0 |  Hexagonal-High 6/mmm |
-| Cubic_High | 1 |  Cubic Cubic-High m3m |
-| Hexagonal_Low | 2 |  Hexagonal-Low 6/m |
-| Cubic_Low | 3 |  Cubic Cubic-Low m3 (Tetrahedral) |
-| Triclinic | 4 |  Triclinic -1 |
-| Monoclinic | 5 |  Monoclinic 2/m |
-| OrthoRhombic | 6 |  Orthorhombic mmm |
-| Tetragonal_Low | 7 |  Tetragonal-Low 4/m |
-| Tetragonal_High | 8 |  Tetragonal-High 4/mmm |
-| Trigonal_Low | 9 |  Trigonal-Low -3 |
-| Trigonal_High | 10 |  Trigonal-High -3m |
-| UnknownCrystalStructure | 999 |  Undefined Crystal Structure |
+## DREAM3D-NX Laue Group to Point Group Table
+
+| Internal Value | String Name              | HM Sym | Point Group | Rotation Point Group |
+|----------------|--------------------------|--------|-------------|----------------------|
+| 0              | Hexagonal_High           | 6/mmm  | 27          | 622                  |
+| 1              | Cubic_High               | m-3m   | 32          | 432                  |
+| 2              | Hexagonal_Low            | 6/m    | 23          | 6                    |
+| 3              | Cubic_Low                | m-3    | 29          | 23                   |
+| 4              | Triclinic                | -1     | 2           | 1                    |
+| 5              | Monoclinic               | 2/m    | 5           | 2                    |
+| 6              | OrthoRhombic             | mmm    | 8           | 222                  |
+| 7              | Tetragonal_Low           | 4/m    | 11          | 4                    |
+| 8              | Tetragonal_High          | 4/mmm  | 15          | 422                  |
+| 9              | Trigonal_Low             | -3     | 17          | 3                    |
+| 10             | Trigonal_High            | -3m    | 20          | 32                   |
+| 999            | UnknownCrystalStructure  |        |             |                      |
 
 ### Phase Type
 

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

@@ -6,7 +6,7 @@
 #include "simplnx/Utilities/FilterUtilities.hpp"
 #include "simplnx/Utilities/MaskCompareUtilities.hpp"
 
-#include "EbsdLib/LaueOps/LaueOps.h"
+#include <EbsdLib/LaueOps/LaueOps.h>
 
 #include <iostream>
 

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

@@ -8,7 +8,7 @@
 #include "simplnx/Utilities/FilterUtilities.hpp"
 #include "simplnx/Utilities/StringUtilities.hpp"
 
-#include "EbsdLib/LaueOps/LaueOps.h"
+#include <EbsdLib/LaueOps/LaueOps.h>
 
 #include <vector>
 

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

@@ -8,7 +8,7 @@
 #include "simplnx/Utilities/MessageHelper.hpp"
 
 #include "EbsdLib/Core/Orientation.hpp"
-#include "EbsdLib/LaueOps/LaueOps.h"
+#include <EbsdLib/LaueOps/LaueOps.h>
 
 using namespace nx::core;
 

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

@@ -8,7 +8,7 @@
 #include "simplnx/Utilities/ClusteringUtilities.hpp"
 #include "simplnx/Utilities/Math/MatrixMath.hpp"
 
-#include "EbsdLib/LaueOps/LaueOps.h"
+#include <EbsdLib/Core/OrientationTransformation.hpp>
 
 using namespace nx::core;
 using namespace nx::core::OrientationUtilities;
@@ -42,10 +42,16 @@ Result<> CAxisSegmentFeatures::operator()()
       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
   const auto& crystalStructures = m_DataStructure.getDataRefAs<UInt32Array>(m_InputValues->CrystalStructuresArrayPath);
-  for(usize i = 1; i < crystalStructures.size(); ++i)
+  usize numCells = m_CellPhases->getNumberOfTuples();
+  for(usize cellIdx = 0; cellIdx < numCells; ++cellIdx)
   {
-    const auto crystalStructureType = crystalStructures[i];
+    int32 currentPhaseIdx = m_CellPhases->getValue(cellIdx);
+    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 "

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

@@ -4,6 +4,7 @@
 
 #include "simplnx/DataStructure/DataArray.hpp"
 #include "simplnx/Utilities/ImageRotationUtilities.hpp"
+#include "simplnx/Utilities/Math/GeometryMath.hpp"
 #include "simplnx/Utilities/Math/MatrixMath.hpp"
 
 #include "EbsdLib/Core/Orientation.hpp"
@@ -34,6 +35,8 @@ const std::atomic_bool& ComputeAvgCAxes::getCancel()
 // -----------------------------------------------------------------------------
 Result<> ComputeAvgCAxes::operator()()
 {
+
+  // Figure out if all phases are either Hexagonal-Low 6/m or Hexagonal-High 6/mmm Laue Phases
   const auto& crystalStructures = m_DataStructure.getDataRefAs<UInt32Array>(m_InputValues->CrystalStructuresArrayPath);
   bool allPhasesHexagonal = true;
   bool noPhasesHexagonal = true;
@@ -45,17 +48,18 @@ Result<> ComputeAvgCAxes::operator()()
     noPhasesHexagonal = noPhasesHexagonal && !isHex;
   }
 
+  // If NONE of the phases are hexagonal then bail out now with an error
   if(noPhasesHexagonal)
   {
-    return MakeErrorResult(-6402, "Finding the average c-axes requires at least one phase to be Hexagonal-Low 6/m or Hexagonal-High 6/mmm type crystal structures but none were found.");
+    return MakeErrorResult(-76402, "No phases that have a crystal symmetry of Hexagonal (6/mmm or 6/m) were found.");
   }
 
   Result<> result;
+
+  // Throw a warning for any NON-Hex Laue Phases
   if(!allPhasesHexagonal)
   {
-    result.warnings().push_back(
-        {-6403,
-         "Finding the average c-axes requires Hexagonal-Low 6/m or Hexagonal-High 6/mmm type crystal structures. All calculations for non Hexagonal phases will be skipped and a NaN value inserted."});
+    result.warnings().push_back({-76403, "Non Hexagonal phases were found. All calculations for non Hexagonal phases will be skipped and a NaN value inserted."});
   }
 
   const auto& featureIds = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->FeatureIdsArrayPath);
@@ -66,77 +70,87 @@ Result<> ComputeAvgCAxes::operator()()
   const usize totalPoints = featureIds.getNumberOfTuples();
   const usize totalFeatures = avgCAxes.getNumberOfTuples();
 
-  Matrix3fR g1T;
+  Matrix3dR g1T;
   g1T.fill(0.0f);
-  const Eigen::Vector3f cAxis{0.0f, 0.0f, 1.0f};
-  Eigen::Vector3f c1{0.0f, 0.0f, 0.0f};
+  const Eigen::Vector3d cAxis{0.0f, 0.0f, 1.0f};
+  Eigen::Vector3d c1{0.0f, 0.0f, 0.0f};
 
   std::vector<int32> counter(totalFeatures, 0);
 
-  Eigen::Vector3f curCAxis{0.0f, 0.0f, 0.0f};
-  usize cAxesIndex = 0;
-  float32 w = 0.0f;
-
+  // Loop over each cell
   for(usize i = 0; i < totalPoints; i++)
   {
-    if(featureIds[i] > 0)
+    int32 currentFeatureId = featureIds[i];
+    // If the featureId for a given cell is valid ( > 0) then analyze that value
+    if(currentFeatureId > 0)
     {
-      cAxesIndex = 3 * featureIds[i];
-      const auto crystalStructureType = crystalStructures[cellPhases[i]];
-      if(crystalStructureType == EbsdLib::CrystalStructure::Hexagonal_High || crystalStructureType == EbsdLib::CrystalStructure::Hexagonal_Low)
-      {
-        const usize quatIndex = i * 4;
-
-        // Create the 3x3 Orientation Matrix from the Quaternion. This represents a passive rotation matrix
-        OrientationF oMatrix = OrientationTransformation::qu2om<QuatF, OrientationF>({quats[quatIndex], quats[quatIndex + 1], quats[quatIndex + 2], quats[quatIndex + 3]});
-
-        // Convert the passive rotation matrix to an active rotation matrix
-        g1T = OrientationMatrixToGMatrixTranspose(oMatrix);
-
-        // 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 = g1T * cAxis;
-
-        // normalize so that the magnitude is 1
-        c1.normalize();
-
-        // Compute the running average c-axis and normalize the result
-        curCAxis[0] = avgCAxes[cAxesIndex] / counter[featureIds[i]];
-        curCAxis[1] = avgCAxes[cAxesIndex + 1] / counter[featureIds[i]];
-        curCAxis[2] = avgCAxes[cAxesIndex + 2] / counter[featureIds[i]];
-        curCAxis.normalize();
-
-        // Ensure that angle between the current point's sample reference frame C-Axis
-        // and the running average sample C-Axis is positive
-        w = ImageRotationUtilities::CosBetweenVectors(c1, curCAxis);
-        if(w < 0)
-        {
-          c1 *= -1.0f;
-        }
-        // Continue summing up the rotations
-        counter[featureIds[i]]++;
-        avgCAxes[cAxesIndex] += c1[0];
-        avgCAxes[cAxesIndex + 1] += c1[1];
-        avgCAxes[cAxesIndex + 2] += c1[2];
-      }
-      else
+      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 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)
       {
         avgCAxes[cAxesIndex] = NAN;
         avgCAxes[cAxesIndex + 1] = NAN;
         avgCAxes[cAxesIndex + 2] = NAN;
+        continue;
+      }
+
+      counter[currentFeatureId]++; // Increment the count
+      const usize quatIndex = i * 4;
+
+      // Create the 3x3 Orientation Matrix from the Quaternion. This represents a passive rotation matrix
+      OrientationD oMatrix = OrientationTransformation::qu2om<QuatD, OrientationD>({quats[quatIndex], quats[quatIndex + 1], quats[quatIndex + 2], quats[quatIndex + 3]});
+
+      // Convert the passive rotation matrix to an active rotation matrix by taking the transpose
+      g1T = OrientationMatrixToGMatrixTranspose(oMatrix);
+
+      // 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 = g1T * cAxis;
+
+      // normalize so that the magnitude is 1
+      c1.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();
+
+      // 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)
+      {
+        c1 *= -1.0f;
       }
+
+      // Continue summing up the rotations
+      float value = avgCAxes[cAxesIndex] + c1[0];
+      avgCAxes[cAxesIndex] = value;
+
+      value = avgCAxes[cAxesIndex + 1] + c1[1];
+      avgCAxes[cAxesIndex + 1] = value;
+
+      value = avgCAxes[cAxesIndex + 2] + c1[2];
+      avgCAxes[cAxesIndex + 2] = value;
     }
   }
 
   for(size_t i = 1; i < totalFeatures; i++)
   {
     const usize tupleIndex = i * 3;
-    if(std::isnan(avgCAxes[tupleIndex]))
+    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)
     {
       avgCAxes[tupleIndex] = 0;
@@ -145,19 +159,18 @@ Result<> ComputeAvgCAxes::operator()()
     }
     else
     {
-      float32 a = avgCAxes.getValue(tupleIndex);
-      float32 b = avgCAxes.getValue(tupleIndex + 1);
-      float32 c = avgCAxes.getValue(tupleIndex + 2);
-
-      a /= counter[i];
-      b /= counter[i];
-      c /= counter[i];
-
-      MatrixMath::Normalize3x1(a, b, c);
-
-      avgCAxes.setValue(tupleIndex, a);
-      avgCAxes.setValue(tupleIndex + 1, b);
-      avgCAxes.setValue(tupleIndex + 2, c);
+      // 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;
     }
   }
   return result;