VV: Rotate Euler Reference Frame fully V&V'ed

src/Plugins/OrientationAnalysis/docs/RotateEulerRefFrameFilter.md

@@ -6,16 +6,31 @@ Processing (Conversion)
 
 ## Description
 
-This **Filter** performs a passive rotation (Right hand rule) of the Euler Angles about a user defined axis. The *reference frame* is being rotated and thus the *Euler Angles* necessary to represent the same orientation must change to account for the new *reference frame*.  The user can set an *angle* and an *axis* to define the rotation of the *reference frame*.
+This **Filter** performs a passive rotation (Right hand rule) of the Euler Angles about a user defined axis. The *reference frame* is being rotated and thus the *Euler Angles* necessary to represent the same orientation must change to account for the new *reference frame*. The user can set an *angle* and an *axis* to define the rotation of the *reference frame*.
+
+For each Euler angle triplet the filter computes `g' = g · R(n, ω)`, where `g` is the orientation matrix of the input Euler angles (Bunge ZXZ convention), `R(n, ω)` is the rotation matrix for the user-supplied axis `n` and angle `ω`, and `g'` is converted back to Euler angles. For a rotation about the sample Z axis this reduces to `φ1' = φ1 − ω (mod 2π)` with `Φ` and `φ2` unchanged.
+
+### Units & Conventions
+
+- The rotation **angle parameter is in degrees**; the **Euler angle data is in radians** (Bunge ZXZ convention).
+- The rotation axis does not need to be normalized — the filter normalizes it internally. A zero-length axis (0,0,0) is rejected during preflight.
+- The selected Euler angles array is modified **in place** (no new array is created).
+- Output Euler angles are canonicalized to `φ1, φ2 ∈ [0, 2π)` and `Φ ∈ [0, π]` by the double-precision conversion, then stored as `float32`. Because a `double` just below `2π` rounds up to the `float32` value `6.2831855`, a stored `φ1`/`φ2` can be numerically equal to (or a hair above) `2π`; treat the upper bound as inclusive at float32 precision.
 
 % Auto generated parameter table will be inserted here
 
 ## Example Pipelines
 
-+ INL Export
-+ Export Small IN100 ODF Data (StatsGenerator)
-+ Edax IPF Colors
-+ Confidence Index Histogram
++ EBSD_File_Processing/Read_EDAX_Ang_File
++ EBSD_File_Processing/Edax_IPF_Colors
++ EBSD_File_Processing/CI_Histogram
++ EBSD_File_Processing/EBSD_Hexagonal_Data_Analysis
++ EBSD_File_Processing/ImportEdaxOIMData
+
+## References
+
+- D. Rowenhorst, A. D. Rollett, G. S. Rohrer, M. Groeber, M. Jackson, P. J. Konijnenberg, and M. De Graef, "Consistent representations of and conversions between 3D rotations," *Modelling and Simulation in Materials Science and Engineering*, vol. 23, no. 8, p. 083501, 2015. DOI: [10.1088/0965-0393/23/8/083501](https://doi.org/10.1088/0965-0393/23/8/083501)
+
 
 ## License & Copyright
 

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

@@ -12,6 +12,9 @@
 #include <EbsdLib/Orientation/OrientationFwd.hpp>
 #include <EbsdLib/Orientation/OrientationMatrix.hpp>
 
+#include <cmath>
+#include <limits>
+
 using namespace nx::core;
 
 namespace
@@ -27,25 +30,27 @@ class RotateEulerRefFrameImpl
 public:
   RotateEulerRefFrameImpl(Float32Array& data, const FloatVec3& rotAxis, float angle, const std::atomic_bool& shouldCancel, ProgressMessageHelper& progressMessageHelper)
   : m_CellEulerAngles(data)
-  , m_AxisAngle(rotAxis)
+  , m_RotationAxis(rotAxis)
   , m_Angle(angle)
   , m_ShouldCancel(shouldCancel)
   , m_ProgressMessageHelper(progressMessageHelper)
   {
   }
-  virtual ~RotateEulerRefFrameImpl() = default;
+  ~RotateEulerRefFrameImpl() = default;
 
   void convert(size_t start, size_t end) const
   {
-    ebsdlib::OrientationMatrixDType om = ebsdlib::AxisAngleDType(m_AxisAngle[0], m_AxisAngle[1], m_AxisAngle[2], m_Angle * nx::core::numbers::pi / 180.0).toOrientationMatrix();
+    // m_Angle arrives in degrees (user-facing parameter) while the Euler angle data is in
+    // radians. The axis-angle pair produces the active rotation matrix R, so gNew = g * R
+    // implements a passive rotation of the sample reference frame by +angle (right-hand rule).
+    ebsdlib::OrientationMatrixDType om = ebsdlib::AxisAngleDType(m_RotationAxis[0], m_RotationAxis[1], m_RotationAxis[2], m_Angle * nx::core::numbers::pi / 180.0).toOrientationMatrix();
 
     OrientationUtilities::Matrix3dR rotMat = om.toEigenGMatrix();
 
     ProgressMessenger progressMessenger = m_ProgressMessageHelper.createProgressMessenger();
 
     usize counter = 0;
     usize counterIncrement = (end - start) / 100;
-    // float ea1 = 0, ea2 = 0, ea3 = 0;
     for(size_t i = start; i < end; i++)
     {
       if(m_ShouldCancel)
@@ -77,7 +82,7 @@ class RotateEulerRefFrameImpl
 
 private:
   Float32Array& m_CellEulerAngles;
-  FloatVec3 m_AxisAngle;
+  FloatVec3 m_RotationAxis;
   float m_Angle = 0.0F;
   const std::atomic_bool& m_ShouldCancel;
   ProgressMessageHelper& m_ProgressMessageHelper;
@@ -109,16 +114,29 @@ Result<> RotateEulerRefFrame::operator()()
   size_t totalElements = eulerAngles.getNumberOfTuples();
 
   nx::core::FloatVec3 axis = {m_InputValues->rotationAxis[0], m_InputValues->rotationAxis[1], m_InputValues->rotationAxis[2]};
+  // The filter's preflight rejects a zero-length axis, but guard here as well so that any direct reuse of
+  // this Algorithm class cannot silently NaN-corrupt the data (normalize() of a zero vector is NaN).
+  const float32 axisMagnitude = std::sqrt(axis[0] * axis[0] + axis[1] * axis[1] + axis[2] * axis[2]);
+  if(axisMagnitude < std::numeric_limits<float32>::epsilon())
+  {
+    return MakeErrorResult(-67050, "The rotation axis has zero length; a rotation axis must be a non-zero vector.");
+  }
   axis = axis.normalize();
 
   MessageHelper messageHelper(m_MessageHandler);
   ProgressMessageHelper progressMessageHelper = messageHelper.createProgressMessageHelper();
   progressMessageHelper.setMaxProgresss(totalElements);
   progressMessageHelper.setProgressMessageTemplate("RotateEulerRefFrame: {:.2f}% complete");
 
-  // Allow data-based parallelization
+  // Data-based parallelization: each worker reads and writes only its own disjoint tuple range of the
+  // in-place Euler array. Per the project thread-safety policy, concurrent DataStore access is unsafe for
+  // out-of-core stores, so requireArraysInMemory disables parallelization unless the array is resident in
+  // memory (the codebase-sanctioned pattern; see ConvertOrientations / PartitionGeometry).
   ParallelDataAlgorithm dataAlg;
   dataAlg.setRange(0, totalElements);
+  IParallelAlgorithm::AlgorithmArrays algArrays;
+  algArrays.push_back(&eulerAngles);
+  dataAlg.requireArraysInMemory(algArrays);
   dataAlg.execute(RotateEulerRefFrameImpl(eulerAngles, axis, m_InputValues->rotationAxis[3], m_ShouldCancel, progressMessageHelper));
   return {};
 }

src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/RotateEulerRefFrameFilter.cpp

@@ -9,6 +9,9 @@
 
 #include "simplnx/Parameters/VectorParameter.hpp"
 
+#include <cmath>
+#include <limits>
+
 using namespace nx::core;
 
 namespace nx::core
@@ -83,6 +86,15 @@ IFilter::PreflightResult RotateEulerRefFrameFilter::preflightImpl(const DataStru
 
   std::vector<PreflightValue> preflightUpdatedValues;
 
+  // A zero-length axis cannot be normalized and would silently fill the Euler angles with NaN
+  const float32 axisMagnitude =
+      std::sqrt(pRotationAxisAngleValue[0] * pRotationAxisAngleValue[0] + pRotationAxisAngleValue[1] * pRotationAxisAngleValue[1] + pRotationAxisAngleValue[2] * pRotationAxisAngleValue[2]);
+  if(axisMagnitude < std::numeric_limits<float32>::epsilon())
+  {
+    return {MakeErrorResult<OutputActions>(-96200, fmt::format("Rotation axis <{}, {}, {}> has zero length. The ijk components must define a non-zero direction.", pRotationAxisAngleValue[0],
+                                                               pRotationAxisAngleValue[1], pRotationAxisAngleValue[2]))};
+  }
+
   resultOutputActions.value().modifiedActions.emplace_back(
       DataObjectModification{pCellEulerAnglesArrayPathValue, DataObjectModification::ModifiedType::Modified, dataStructure.getData(pCellEulerAnglesArrayPathValue)->getDataObjectType()});
 

src/Plugins/OrientationAnalysis/test/RotateEulerRefFrameTest.cpp

@@ -14,13 +14,87 @@
 #include "simplnx/Pipeline/PipelineFilter.hpp"
 #include "simplnx/UnitTest/UnitTestCommon.hpp"
 
+#include "simplnx/DataStructure/AttributeMatrix.hpp"
+#include "simplnx/DataStructure/DataArray.hpp"
+
+#include <array>
 #include <filesystem>
 #include <fstream>
 
 namespace fs = std::filesystem;
 
 using namespace nx::core;
 
+namespace AnalyticalFixtures
+{
+// Class 1 (Analytical) oracle fixtures for RotateEulerRefFrameFilter.
+//
+// The filter rotates the sample reference frame by angle w (degrees, right-hand rule)
+// about sample-frame axis n. Derivation (independent of EbsdLib/SIMPLNX/legacy code):
+//   - Bunge passive convention: v_crystal = g . v_sample
+//   - Rotating the reference frame by +w about unit axis n gives new-frame coordinates
+//     u_new = P(w) . u_old with P(w) the passive rotation matrix, so
+//     g' = g . P(w)^T = g . R_active(n, w)   (R_active = Rodrigues rotation matrix)
+//   - Euler <-> matrix conversions per Rowenhorst et al 2015, Modelling Simul. Mater.
+//     Sci. Eng. 23 083501, Eq. A.5 (eu2om) and Eq. A.9 (om2eu), with outputs
+//     canonicalized to phi1, phi2 in [0, 2pi), Phi in [0, pi].
+//
+// Closed forms used below (all verified by the independent numpy script archived in
+// vv/provenance — see the V&V report):
+//   - Z-axis rotation: phi1' = phi1 - w (mod 2pi), Phi and phi2 unchanged.
+//   - Identity orientation + axis-angle (n, w): g' = R_active(n, w), om2eu by hand.
+struct FixtureData
+{
+  std::string name;
+  std::array<float32, 3> eulerIn;   // radians
+  std::array<float32, 4> axisAngle; // i, j, k, w (degrees)
+  std::array<float64, 3> expected;  // radians, canonical ranges
+};
+
+const std::vector<FixtureData> k_Fixtures = {
+    // F1: g = I, rotate frame +90 about z. phi1' = 0 - pi/2 mod 2pi = 3pi/2.
+    {"F1 identity + z90", {0.0F, 0.0F, 0.0F}, {0.0F, 0.0F, 1.0F, 90.0F}, {4.712388980384709, 0.0, 0.0}},
+    // F2: general orientation, z-axis: phi1' = 1.0 - pi/6; Phi, phi2 unchanged.
+    {"F2 general + z30", {1.0F, 0.5F, 0.3F}, {0.0F, 0.0F, 1.0F, 30.0F}, {0.476401224401701, 0.5, 0.3}},
+    // F3: wrap-around: phi1' = 0.1 - pi/2 + 2pi.
+    {"F3 wraparound + z90", {0.1F, 0.5F, 0.3F}, {0.0F, 0.0F, 1.0F, 90.0F}, {4.812388980384732, 0.5, 0.3}},
+    // F4: g = I, +90 about x: g' = R_active(x, 90) => om2eu = (pi, pi/2, pi) by hand.
+    {"F4 identity + x90", {0.0F, 0.0F, 0.0F}, {1.0F, 0.0F, 0.0F, 90.0F}, {3.141592653589793, 1.570796326794897, 3.141592653589793}},
+    // F5: axis (2,0,0) must normalize to (1,0,0); expected identical to F4.
+    {"F5 unnormalized axis x90", {0.0F, 0.0F, 0.0F}, {2.0F, 0.0F, 0.0F, 90.0F}, {3.141592653589793, 1.570796326794897, 3.141592653589793}},
+    // F6: w = 0 is a no-op for already-canonical inputs.
+    {"F6 zero angle no-op", {5.0F, 1.0F, 2.0F}, {0.0F, 0.0F, 1.0F, 0.0F}, {5.0, 1.0, 2.0}},
+    // F7: g = I, +120 about (1,1,1): R_active maps x->y->z->x; om2eu = (pi, pi/2, pi/2) by hand.
+    {"F7 identity + 111 120", {0.0F, 0.0F, 0.0F}, {1.0F, 1.0F, 1.0F, 120.0F}, {3.141592653589793, 1.570796326794896, 1.570796326794897}},
+    // F8: general orientation, y-axis 45 deg. No closed form; value from the numpy oracle script.
+    {"F8 general + y45", {0.8F, 1.1F, 0.4F}, {0.0F, 1.0F, 0.0F, 45.0F}, {0.208427151998385, 0.687549596714496, 1.288701818256954}},
+};
+
+// Creates a DataStructure holding a cell AttributeMatrix with one float32 3-component
+// Euler angles array of numTuples tuples, returning the array pointer.
+inline Float32Array* CreateEulerArray(DataStructure& dataStructure, usize numTuples)
+{
+  AttributeMatrix* cellData = AttributeMatrix::Create(dataStructure, "CellData", ShapeType{numTuples});
+  return UnitTest::CreateTestDataArray<float32>(dataStructure, "EulerAngles", {numTuples}, {3}, cellData->getId());
+}
+
+// Runs RotateEulerRefFrameFilter in place on the given euler array.
+inline void RunRotateFilter(DataStructure& dataStructure, const DataPath& eulerPath, const std::array<float32, 4>& axisAngle)
+{
+  RotateEulerRefFrameFilter filter;
+  Arguments args;
+  args.insertOrAssign(RotateEulerRefFrameFilter::k_RotationAxisAngle_Key,
+                      std::make_any<VectorFloat32Parameter::ValueType>(std::vector<float32>{axisAngle[0], axisAngle[1], axisAngle[2], axisAngle[3]}));
+  args.insertOrAssign(RotateEulerRefFrameFilter::k_EulerAnglesArrayPath_Key, std::make_any<DataPath>(eulerPath));
+
+  auto preflightResult = filter.preflight(dataStructure, args);
+  SIMPLNX_RESULT_REQUIRE_VALID(preflightResult.outputActions);
+
+  auto executeResult = filter.execute(dataStructure, args);
+  SIMPLNX_RESULT_REQUIRE_VALID(executeResult.result);
+}
+} // namespace AnalyticalFixtures
+
 TEST_CASE("OrientationAnalysis::RotateEulerRefFrame", "[OrientationAnalysis][RotateEulerRefFrameFilter]")
 {
   const nx::core::UnitTest::TestFileSentinel testDataSentinel(nx::core::unit_test::k_TestFilesDir, "ASCIIData.tar.gz", "ASCIIData");
@@ -168,3 +242,162 @@ TEST_CASE("OrientationAnalysis::RotateEulerRefFrameFilter: SIMPL Backwards Compa
     }
   }
 }
+
+TEST_CASE("OrientationAnalysis::RotateEulerRefFrameFilter: Class 1 Analytical Fixtures", "[OrientationAnalysis][RotateEulerRefFrameFilter]")
+{
+  UnitTest::LoadPlugins();
+
+  // Single-pass double-precision math on float32 storage: float32 rounding of the decimal
+  // inputs perturbs the output by O(1e-7); 1e-5 rad is a comfortable-but-tight bound.
+  constexpr float64 k_Tolerance = 1.0e-5;
+
+  const DataPath k_EulerPath({"CellData", "EulerAngles"});
+
+  for(const auto& fixture : AnalyticalFixtures::k_Fixtures)
+  {
+    DYNAMIC_SECTION(fixture.name)
+    {
+      DataStructure dataStructure;
+      Float32Array* eulerArray = AnalyticalFixtures::CreateEulerArray(dataStructure, 1);
+      for(usize c = 0; c < 3; c++)
+      {
+        (*eulerArray)[c] = fixture.eulerIn[c];
+      }
+
+      AnalyticalFixtures::RunRotateFilter(dataStructure, k_EulerPath, fixture.axisAngle);
+
+      REQUIRE_NOTHROW(dataStructure.getDataRefAs<Float32Array>(k_EulerPath));
+      const auto& resultRef = dataStructure.getDataRefAs<Float32Array>(k_EulerPath);
+      for(usize c = 0; c < 3; c++)
+      {
+        const float64 absDif = std::fabs(static_cast<float64>(resultRef[c]) - fixture.expected[c]);
+        INFO(fmt::format("{}: component {}: actual {} expected {}", fixture.name, c, resultRef[c], fixture.expected[c]));
+        REQUIRE(absDif < k_Tolerance);
+      }
+
+      UnitTest::CheckArraysInheritTupleDims(dataStructure);
+    }
+  }
+}
+
+TEST_CASE("OrientationAnalysis::RotateEulerRefFrameFilter: Zero-Length Axis Fails Preflight", "[OrientationAnalysis][RotateEulerRefFrameFilter]")
+{
+  UnitTest::LoadPlugins();
+
+  const DataPath k_EulerPath({"CellData", "EulerAngles"});
+
+  DataStructure dataStructure;
+  AnalyticalFixtures::CreateEulerArray(dataStructure, 1);
+
+  RotateEulerRefFrameFilter filter;
+  Arguments args;
+  args.insertOrAssign(RotateEulerRefFrameFilter::k_RotationAxisAngle_Key, std::make_any<VectorFloat32Parameter::ValueType>(std::vector<float32>{0.0F, 0.0F, 0.0F, 90.0F}));
+  args.insertOrAssign(RotateEulerRefFrameFilter::k_EulerAnglesArrayPath_Key, std::make_any<DataPath>(k_EulerPath));
+
+  auto preflightResult = filter.preflight(dataStructure, args);
+  SIMPLNX_RESULT_REQUIRE_INVALID(preflightResult.outputActions);
+}
+
+TEST_CASE("OrientationAnalysis::RotateEulerRefFrameFilter: Class 4 Invariants", "[OrientationAnalysis][RotateEulerRefFrameFilter]")
+{
+  UnitTest::LoadPlugins();
+
+  const DataPath k_EulerPath({"CellData", "EulerAngles"});
+
+  // A batch of canonical-range orientations (phi1, phi2 in [0, 2pi), Phi in [0, pi]).
+  const std::vector<std::array<float32, 3>> k_BatchEulers = {
+      {0.0F, 0.0F, 0.0F}, {1.0F, 0.5F, 0.3F}, {0.7F, 0.9F, 4.2F}, {2.0F, 2.5F, 1.0F}, {5.9F, 0.1F, 3.3F}, {3.1F, 1.6F, 6.0F},
+  };
+
+  constexpr float32 k_TwoPi = 6.2831855F;
+  constexpr float32 k_Pi = 3.1415927F;
+
+  SECTION("Output ranges: phi1, phi2 in [0, 2pi], Phi in [0, pi]")
+  {
+    DataStructure dataStructure;
+    Float32Array* eulerArray = AnalyticalFixtures::CreateEulerArray(dataStructure, k_BatchEulers.size());
+    for(usize t = 0; t < k_BatchEulers.size(); t++)
+    {
+      for(usize c = 0; c < 3; c++)
+      {
+        (*eulerArray)[t * 3 + c] = k_BatchEulers[t][c];
+      }
+    }
+
+    AnalyticalFixtures::RunRotateFilter(dataStructure, k_EulerPath, {0.3F, -0.5F, 0.81F, 37.0F});
+
+    const auto& resultRef = dataStructure.getDataRefAs<Float32Array>(k_EulerPath);
+    for(usize t = 0; t < k_BatchEulers.size(); t++)
+    {
+      INFO(fmt::format("tuple {}", t));
+      REQUIRE(resultRef[t * 3 + 0] >= 0.0F);
+      REQUIRE(resultRef[t * 3 + 0] <= k_TwoPi);
+      REQUIRE(resultRef[t * 3 + 1] >= 0.0F);
+      REQUIRE(resultRef[t * 3 + 1] <= k_Pi);
+      REQUIRE(resultRef[t * 3 + 2] >= 0.0F);
+      REQUIRE(resultRef[t * 3 + 2] <= k_TwoPi);
+    }
+
+    UnitTest::CheckArraysInheritTupleDims(dataStructure);
+  }
+
+  SECTION("Round-trip: rotate by (n, w) then (-n, w) recovers the input")
+  {
+    DataStructure dataStructure;
+    Float32Array* eulerArray = AnalyticalFixtures::CreateEulerArray(dataStructure, k_BatchEulers.size());
+    for(usize t = 0; t < k_BatchEulers.size(); t++)
+    {
+      for(usize c = 0; c < 3; c++)
+      {
+        (*eulerArray)[t * 3 + c] = k_BatchEulers[t][c];
+      }
+    }
+
+    AnalyticalFixtures::RunRotateFilter(dataStructure, k_EulerPath, {0.3F, -0.5F, 0.81F, 37.0F});
+    AnalyticalFixtures::RunRotateFilter(dataStructure, k_EulerPath, {-0.3F, 0.5F, -0.81F, 37.0F});
+
+    const auto& resultRef = dataStructure.getDataRefAs<Float32Array>(k_EulerPath);
+    constexpr float32 k_Tolerance = 1.0e-4F; // two float32 store/load passes
+    for(usize i = 0; i < k_BatchEulers.size() * 3; i++)
+    {
+      const float32 absDif = std::fabs(resultRef[i] - k_BatchEulers[i / 3][i % 3]);
+      INFO(fmt::format("flat index {}", i));
+      REQUIRE(absDif < k_Tolerance);
+    }
+
+    UnitTest::CheckArraysInheritTupleDims(dataStructure);
+  }
+
+  SECTION("Composability: 45 deg twice equals 90 deg once")
+  {
+    DataStructure twiceDataStructure;
+    Float32Array* twiceArray = AnalyticalFixtures::CreateEulerArray(twiceDataStructure, k_BatchEulers.size());
+    DataStructure onceDataStructure;
+    Float32Array* onceArray = AnalyticalFixtures::CreateEulerArray(onceDataStructure, k_BatchEulers.size());
+    for(usize t = 0; t < k_BatchEulers.size(); t++)
+    {
+      for(usize c = 0; c < 3; c++)
+      {
+        (*twiceArray)[t * 3 + c] = k_BatchEulers[t][c];
+        (*onceArray)[t * 3 + c] = k_BatchEulers[t][c];
+      }
+    }
+
+    AnalyticalFixtures::RunRotateFilter(twiceDataStructure, k_EulerPath, {0.0F, 0.0F, 1.0F, 45.0F});
+    AnalyticalFixtures::RunRotateFilter(twiceDataStructure, k_EulerPath, {0.0F, 0.0F, 1.0F, 45.0F});
+    AnalyticalFixtures::RunRotateFilter(onceDataStructure, k_EulerPath, {0.0F, 0.0F, 1.0F, 90.0F});
+
+    const auto& twiceRef = twiceDataStructure.getDataRefAs<Float32Array>(k_EulerPath);
+    const auto& onceRef = onceDataStructure.getDataRefAs<Float32Array>(k_EulerPath);
+    constexpr float32 k_Tolerance = 1.0e-4F;
+    for(usize i = 0; i < k_BatchEulers.size() * 3; i++)
+    {
+      const float32 absDif = std::fabs(twiceRef[i] - onceRef[i]);
+      INFO(fmt::format("flat index {}", i));
+      REQUIRE(absDif < k_Tolerance);
+    }
+
+    UnitTest::CheckArraysInheritTupleDims(twiceDataStructure);
+    UnitTest::CheckArraysInheritTupleDims(onceDataStructure);
+  }
+}