ENH: Add option to keep input geometry origin for RotateSampleRefFrame Filter (#1355)

* Also fix unit test so the test is actually testing the code
* Add another unit test section to test the new option.
* Documentation updated

Signed-off-by: Michael Jackson <[email protected]>

Author: imikejackson

src/Plugins/SimplnxCore/docs/RotateSampleRefFrameFilter.md

@@ -6,7 +6,7 @@ Sampling (Rotating/Transforming)
 
 ## Description
 
-**NOTE: As of July 2023, this filter is only verified to work with a rotation angle of 180 degrees, a rotation axis of (010), and a (0, 0, 0) origin.**
+**NOTE: As of July 2023, this filter is only verified to work with a rotation angle of 90 or 180 degrees, a rotation axis of (010) || (100) || (001) The origin must also be (0, 0, 0).**
 
 This **Filter** will rotate the *spatial reference frame* around a user defined axis, by a user defined angle.  The **Filter** will modify the (X, Y, Z) positions of each **Cell** to correctly represent where the **Cell** sits in the newly defined reference frame. For example, if a user selected a *rotation angle* of 90<sup>o</sup> and a *rotation axis* of (001), then a **Cell** sitting at (10, 0, 0) would be transformed to (0, -10, 0), since the new *reference frame* would have x'=y and y'=-x.
 
@@ -22,6 +22,11 @@ The equivalent rotation matrix for the above rotation would be the following:
 
 When importing EBSD data from EDAX typically the user will need to rotate the sample reference frame about the <010> (Y) axis. This results in the image comparison below. Note that in the original image the origin of the data is at (0, 0) microns but after rotation the origin now becomes (-189, 0) microns. If you need to reset the origin back to (0,0) then the filter "Set Origin & Spacing" can be run.
 
+## Notes
+
+The transformation will most likely create an origin that is different from the input geometry's origin. If you wish to still keep the input geometry's origin
+then there is an option to allow you to do that. By default, the option is OFF so that the origin generated from the transformation is used.
+
 ![Imported EBSD Data Rotated about the <010> axis](Images/RotateSampleRefFrame_1.png)
 
 % Auto generated parameter table will be inserted here

src/Plugins/SimplnxCore/src/SimplnxCore/Filters/RotateSampleRefFrameFilter.cpp

@@ -10,6 +10,7 @@
 #include "simplnx/Filter/Actions/CreateImageGeometryAction.hpp"
 #include "simplnx/Filter/Actions/DeleteDataAction.hpp"
 #include "simplnx/Filter/Actions/RenameDataAction.hpp"
+#include "simplnx/Filter/Actions/UpdateImageGeomAction.hpp"
 #include "simplnx/Parameters/BoolParameter.hpp"
 #include "simplnx/Parameters/ChoicesParameter.hpp"
 #include "simplnx/Parameters/DataGroupCreationParameter.hpp"
@@ -23,18 +24,14 @@
 #include "simplnx/Utilities/ImageRotationUtilities.hpp"
 #include "simplnx/Utilities/ParallelAlgorithmUtilities.hpp"
 #include "simplnx/Utilities/ParallelTaskAlgorithm.hpp"
+#include "simplnx/Utilities/SIMPLConversion.hpp"
 #include "simplnx/Utilities/StringUtilities.hpp"
 
 #include <Eigen/Dense>
 
 #include <fmt/core.h>
 
 #include <algorithm>
-#include <array>
-
-#include "simplnx/Utilities/SIMPLConversion.hpp"
-
-#include <stdexcept>
 
 using namespace nx::core;
 using namespace nx::core::ImageRotationUtilities;
@@ -45,238 +42,6 @@ const std::string k_TempGeometryName = ".rotated_image_geometry";
 
 using RotationRepresentationType = RotateSampleRefFrameFilter::RotationRepresentation;
 
-constexpr float32 k_Threshold = 0.01f;
-
-const Eigen::Vector3f k_XAxis = Eigen::Vector3f::UnitX();
-const Eigen::Vector3f k_YAxis = Eigen::Vector3f::UnitY();
-const Eigen::Vector3f k_ZAxis = Eigen::Vector3f::UnitZ();
-#if 0
-void DetermineMinMax(const Matrix3fR& rotationMatrix, const FloatVec3& spacing, usize col, usize row, usize plane, float32& xMin, float32& xMax, float32& yMin, float32& yMax, float32& zMin,
-                     float32& zMax)
-{
-  Eigen::Vector3f coords(static_cast<float32>(col) * spacing[0], static_cast<float32>(row) * spacing[1], static_cast<float32>(plane) * spacing[2]);
-
-  Eigen::Vector3f newCoords = rotationMatrix * coords;
-
-  xMin = std::min(newCoords[0], xMin);
-  xMax = std::max(newCoords[0], xMax);
-
-  yMin = std::min(newCoords[1], yMin);
-  yMax = std::max(newCoords[1], yMax);
-
-  zMin = std::min(newCoords[2], zMin);
-  zMax = std::max(newCoords[2], zMax);
-}
-
-float32 CosBetweenVectors(const Eigen::Vector3f& vectorA, const Eigen::Vector3f& vectorB)
-{
-  float32 normA = vectorA.norm();
-  float32 normB = vectorB.norm();
-
-  if(normA == 0.0f || normB == 0.0f)
-  {
-    return 1.0f;
-  }
-
-  return vectorA.dot(vectorB) / (normA * normB);
-}
-
-float32 DetermineSpacing(const FloatVec3& spacing, const Eigen::Vector3f& axisNew)
-{
-  float32 xAngle = std::abs(CosBetweenVectors(k_XAxis, axisNew));
-  float32 yAngle = std::abs(CosBetweenVectors(k_YAxis, axisNew));
-  float32 zAngle = std::abs(CosBetweenVectors(k_ZAxis, axisNew));
-
-  std::array<float32, 3> axes = {xAngle, yAngle, zAngle};
-
-  auto result = std::max_element(axes.cbegin(), axes.cend());
-
-  usize index = std::distance(axes.cbegin(), result);
-
-  return spacing[index];
-}
-
-RotateArgs CreateRotateArgs(const ImageGeom& imageGeom, const Matrix3fR& rotationMatrix)
-{
-  const SizeVec3 origDims = imageGeom.getDimensions();
-  const FloatVec3 spacing = imageGeom.getSpacing();
-
-  float32 xMin = std::numeric_limits<float32>::max();
-  float32 xMax = std::numeric_limits<float32>::min();
-  float32 yMin = std::numeric_limits<float32>::max();
-  float32 yMax = std::numeric_limits<float32>::min();
-  float32 zMin = std::numeric_limits<float32>::max();
-  float32 zMax = std::numeric_limits<float32>::min();
-
-  const std::vector<std::array<usize, 3>> coords{{0, 0, 0},
-                                                 {origDims[0] - 1, 0, 0},
-                                                 {0, origDims[1] - 1, 0},
-                                                 {origDims[0] - 1, origDims[1] - 1, 0},
-                                                 {0, 0, origDims[2] - 1},
-                                                 {origDims[0] - 1, 0, origDims[2] - 1},
-                                                 {0, origDims[1] - 1, origDims[2] - 1},
-                                                 {origDims[0] - 1, origDims[1] - 1, origDims[2] - 1}};
-
-  for(const auto& item : coords)
-  {
-    DetermineMinMax(rotationMatrix, spacing, item[0], item[1], item[2], xMin, xMax, yMin, yMax, zMin, zMax);
-  }
-
-  Eigen::Vector3f xAxisNew = rotationMatrix * k_XAxis;
-  Eigen::Vector3f yAxisNew = rotationMatrix * k_YAxis;
-  Eigen::Vector3f zAxisNew = rotationMatrix * k_ZAxis;
-
-  float32 xResNew = DetermineSpacing(spacing, xAxisNew);
-  float32 yResNew = DetermineSpacing(spacing, yAxisNew);
-  float32 zResNew = DetermineSpacing(spacing, zAxisNew);
-
-  ImageGeom::MeshIndexType xpNew = static_cast<int64>(std::nearbyint((xMax - xMin) / xResNew) + 1);
-  ImageGeom::MeshIndexType ypNew = static_cast<int64>(std::nearbyint((yMax - yMin) / yResNew) + 1);
-  ImageGeom::MeshIndexType zpNew = static_cast<int64>(std::nearbyint((zMax - zMin) / zResNew) + 1);
-
-  RotateArgs params;
-
-  params.xp = static_cast<int64>(origDims[0]);
-  params.xRes = spacing[0];
-  params.yp = static_cast<int64>(origDims[1]);
-  params.yRes = spacing[1];
-  params.zp = static_cast<int64>(origDims[2]);
-  params.zRes = spacing[2];
-
-  params.xpNew = static_cast<int64>(xpNew);
-  params.xResNew = xResNew;
-  params.xMinNew = xMin;
-  params.ypNew = static_cast<int64>(ypNew);
-  params.yResNew = yResNew;
-  params.yMinNew = yMin;
-  params.zpNew = static_cast<int64>(zpNew);
-  params.zResNew = zResNew;
-  params.zMinNew = zMin;
-
-  return params;
-}
-
-template <typename K>
-bool closeEnough(const K& a, const K& b, const K& epsilon = std::numeric_limits<K>::epsilon())
-{
-  return (epsilon > fabs(a - b));
-}
-
-constexpr RotationRepresentationType CastIndexToRotationRepresentation(uint64 index)
-{
-  switch(index)
-  {
-  case to_underlying(RotationRepresentationType::AxisAngle): {
-    return RotationRepresentationType::AxisAngle;
-  }
-  case to_underlying(RotationRepresentationType::RotationMatrix): {
-    return RotationRepresentationType::RotationMatrix;
-  }
-  default: {
-    throw std::runtime_error(fmt::format("RotateSampleRefFrameFilter: Failed to cast index {} to RotationRepresentation", index));
-  }
-  }
-}
-
-Result<Matrix3fR> ConvertAxisAngleToRotationMatrix(const std::vector<float32>& pRotationValue)
-{
-  Matrix3fR transformationMatrix;
-
-  // Convert Degrees to Radians for the last element
-  const float rotAngle = pRotationValue[3] * Constants::k_PiOver180F;
-  // Ensure the axis part is normalized
-  FloatVec3 normalizedAxis(pRotationValue[0], pRotationValue[1], pRotationValue[2]);
-  MatrixMath::Normalize3x1<float32>(normalizedAxis.data());
-
-  const float cosTheta = cos(rotAngle);
-  const float oneMinusCosTheta = 1 - cosTheta;
-  const float sinTheta = sin(rotAngle);
-  const float l = normalizedAxis[0];
-  const float m = normalizedAxis[1];
-  const float n = normalizedAxis[2];
-
-  // First Row:
-  transformationMatrix(0, 0) = l * l * (oneMinusCosTheta) + cosTheta;
-  transformationMatrix(0, 1) = m * l * (oneMinusCosTheta) - (n * sinTheta);
-  transformationMatrix(0, 2) = n * l * (oneMinusCosTheta) + (m * sinTheta);
-
-  // Second Row:
-  transformationMatrix(1, 0) = l * m * (oneMinusCosTheta) + (n * sinTheta);
-  transformationMatrix(1, 1) = m * m * (oneMinusCosTheta) + cosTheta;
-  transformationMatrix(1, 2) = n * m * (oneMinusCosTheta) - (l * sinTheta);
-
-  // Third Row:
-  transformationMatrix(2, 0) = l * n * (oneMinusCosTheta) - (m * sinTheta);
-  transformationMatrix(2, 1) = m * n * (oneMinusCosTheta) + (l * sinTheta);
-  transformationMatrix(2, 2) = n * n * (oneMinusCosTheta) + cosTheta;
-
-  return {transformationMatrix};
-}
-
-Result<Matrix3fR> ConvertRotationTableToRotationMatrix(const std::vector<std::vector<float64>>& rotationMatrixTable)
-{
-  if(rotationMatrixTable.size() != 3)
-  {
-    return MakeErrorResult<Matrix3fR>(-45004, "Rotation Matrix must be 3 x 3");
-  }
-
-  for(const auto& row : rotationMatrixTable)
-  {
-    if(row.size() != 3)
-    {
-      return MakeErrorResult<Matrix3fR>(-45005, "Rotation Matrix must be 3 x 3");
-    }
-  }
-  Matrix3fR rotationMatrix;
-  const usize numTableRows = rotationMatrixTable.size();
-  const usize numTableCols = rotationMatrixTable[0].size();
-  for(size_t rowIndex = 0; rowIndex < numTableRows; rowIndex++)
-  {
-    std::vector<double> row = rotationMatrixTable[rowIndex];
-    for(size_t colIndex = 0; colIndex < numTableCols; colIndex++)
-    {
-      rotationMatrix(rowIndex, colIndex) = static_cast<float>(row[colIndex]);
-    }
-  }
-
-  float32 determinant = rotationMatrix.determinant();
-
-  if(!closeEnough(determinant, 1.0f, k_Threshold))
-  {
-    return MakeErrorResult<Matrix3fR>(-45006, fmt::format("Rotation Matrix must have a determinant of 1 (is {})", determinant));
-  }
-
-  Matrix3fR transpose = rotationMatrix.transpose();
-  Matrix3fR inverse = rotationMatrix.inverse();
-
-  if(!transpose.isApprox(inverse, k_Threshold))
-  {
-    return MakeErrorResult<Matrix3fR>(-45007, "Rotation Matrix's inverse and transpose must be equal");
-  }
-
-  return {rotationMatrix};
-}
-#endif
-
-// Result<Matrix3fR> ComputeRotationMatrix(const Arguments& args)
-//{
-//   auto rotationRepresentationIndex = args.value<uint64>(RotateSampleRefFrameFilter::k_RotationRepresentation_Key);
-//
-//   RotationRepresentationType rotationRepresentation = CastIndexToRotationRepresentation(rotationRepresentationIndex);
-//
-//   switch(rotationRepresentation)
-//   {
-//   case RotationRepresentationType::AxisAngle: {
-//     auto pRotationValue = args.value<std::vector<float32>>(RotateSampleRefFrameFilter::k_RotationAxisAngle_Key);
-//     return ConvertAxisAngleToRotationMatrix(pRotationValue);
-//   }
-//   case RotationRepresentationType::RotationMatrix: {
-//     auto rotationMatrixTable = args.value<DynamicTableParameter::ValueType>(RotateSampleRefFrameFilter::k_RotationMatrix_Key);
-//     return ConvertRotationTableToRotationMatrix(rotationMatrixTable);
-//   }
-//   }
-// }
-
 } // namespace
 
 namespace nx::core
@@ -324,6 +89,8 @@ Parameters RotateSampleRefFrameFilter::parameters() const
   params.insert(std::make_unique<VectorFloat32Parameter>(k_RotationAxisAngle_Key, "Rotation Axis-Angle [<ijk>w]", "Axis-Angle in sample reference frame to rotate about.",
                                                          VectorFloat32Parameter::ValueType{0.0f, 0.0f, 1.0f, 90.0F}, std::vector<std::string>{"i", "j", "k", "w (Deg)"}));
   params.insertLinkableParameter(std::make_unique<BoolParameter>(k_RemoveOriginalGeometry_Key, "Perform In-Place Rotation", "Performs the rotation in-place for the given Image Geometry", true));
+  params.insertLinkableParameter(
+      std::make_unique<BoolParameter>(k_KeepInputGeometryOrigin_Key, "Keep Input Geometry's Origin", "The input geometry's origin is kept instead of the origin resulting from the transform", false));
 
   DynamicTableInfo tableInfo;
   tableInfo.setColsInfo(DynamicTableInfo::StaticVectorInfo({"1", "2", "3", "4"}));
@@ -362,10 +129,10 @@ IFilter::UniquePointer RotateSampleRefFrameFilter::clone() const
 IFilter::PreflightResult RotateSampleRefFrameFilter::preflightImpl(const DataStructure& dataStructure, const Arguments& filterArgs, const MessageHandler& messageHandler,
                                                                    const std::atomic_bool& shouldCancel, const ExecutionContext& executionContext) const
 {
-
   auto srcImagePath = filterArgs.value<DataPath>(k_SelectedImageGeometryPath_Key);
   auto destImagePath = filterArgs.value<DataPath>(k_CreatedImageGeometryPath_Key);
   auto pRemoveOriginalGeometry = filterArgs.value<bool>(k_RemoveOriginalGeometry_Key);
+  auto keepInputGeometryOrigin = filterArgs.value<bool>(k_KeepInputGeometryOrigin_Key);
 
   nx::core::Result<OutputActions> resultOutputActions;
 
@@ -394,9 +161,12 @@ IFilter::PreflightResult RotateSampleRefFrameFilter::preflightImpl(const DataStr
   const std::vector<usize> dims = {static_cast<usize>(rotateArgs.xpNew), static_cast<usize>(rotateArgs.ypNew), static_cast<usize>(rotateArgs.zpNew)};
   const std::vector<float32> spacing = {rotateArgs.xResNew, rotateArgs.yResNew, rotateArgs.zResNew};
   auto origin = selectedImageGeom.getOrigin().toContainer<std::vector<float32>>();
-  origin[0] += rotateArgs.xMinNew;
-  origin[1] += rotateArgs.yMinNew;
-  origin[2] += rotateArgs.zMinNew;
+  if(!keepInputGeometryOrigin)
+  {
+    origin[0] += rotateArgs.xMinNew;
+    origin[1] += rotateArgs.yMinNew;
+    origin[2] += rotateArgs.zMinNew;
+  }
 
   std::vector<usize> dataArrayShape = {dims[2], dims[1], dims[0]}; // The DataArray shape goes slowest to fastest (ZYX)
 
@@ -487,6 +257,7 @@ IFilter::PreflightResult RotateSampleRefFrameFilter::preflightImpl(const DataStr
   {
     resultOutputActions.value().appendDeferredAction(std::make_unique<RenameDataAction>(destImagePath, srcImagePath.getTargetName()));
   }
+
   // Return both the resultOutputActions and the preflightUpdatedValues via std::move()
   return {std::move(resultOutputActions), std::move(preflightUpdatedValues)};
 }
@@ -498,9 +269,10 @@ Result<> RotateSampleRefFrameFilter::executeImpl(DataStructure& dataStructure, c
   auto destImagePath = filterArgs.value<DataPath>(k_CreatedImageGeometryPath_Key);
   auto sliceBySlice = filterArgs.value<bool>(k_RotateSliceBySlice_Key);
   auto removeOriginalGeometry = filterArgs.value<bool>(k_RemoveOriginalGeometry_Key);
+  auto keepInputGeometryOrigin = filterArgs.value<bool>(k_KeepInputGeometryOrigin_Key);
 
   auto& srcImageGeom = dataStructure.getDataRefAs<ImageGeom>(srcImagePath);
-
+  auto sourceImageGeomorigin = srcImageGeom.getOrigin();
   if(removeOriginalGeometry)
   {
     auto tempPathVector = srcImagePath.getPathVector();
@@ -561,6 +333,10 @@ Result<> RotateSampleRefFrameFilter::executeImpl(DataStructure& dataStructure, c
 
   taskRunner.wait(); // This will spill over if the number of DataArrays to process does not divide evenly by the number of threads.
 
+  if(keepInputGeometryOrigin)
+  {
+    destImageGeom.setOrigin(srcImageGeom.getOrigin());
+  }
   return {};
 }
 

src/Plugins/SimplnxCore/src/SimplnxCore/Filters/RotateSampleRefFrameFilter.hpp

@@ -38,6 +38,7 @@ class SIMPLNXCORE_EXPORT RotateSampleRefFrameFilter : public IFilter
   static inline constexpr StringLiteral k_CreatedImageGeometryPath_Key = "output_image_geometry_path";
   static inline constexpr StringLiteral k_RotateSliceBySlice_Key = "rotate_slice_by_slice";
   static inline constexpr StringLiteral k_RemoveOriginalGeometry_Key = "remove_original_geometry";
+  static inline constexpr StringLiteral k_KeepInputGeometryOrigin_Key = "keep_input_geometry_origin";
 
   /**
    * @brief Reads SIMPL json and converts it simplnx Arguments.

src/Plugins/SimplnxCore/test/CMakeLists.txt

@@ -240,7 +240,7 @@ if(EXISTS "${DREAM3D_DATA_DIR}" AND SIMPLNX_DOWNLOAD_TEST_FILES)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME ResampleImageGeom_Exemplar_2.tar.gz SHA512 0a5dc046410a9d93682a40581b169cf3ca5b516623d90eb081de8a8623c55de45a2b00427cd298d6f17fb84bb93dee061a96a4b4386392f925ecdc3fbfe6e5a6)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME ResampleImageGeom_Exemplar.tar.gz SHA512 464029b7354b96a943d75c495ef02bac0f834032e5a86576dde9afee51febff3fd6ffd7d4f8f1e9f8315d8cda36971df26601c7212e1876151109ca5428b8659)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME reverse_triangle_winding.tar.gz SHA512 63247d7b8a0deee2fdb737527506312827331861758cd6106974d3aa1b9cb9c1d3b85d4b135e3eda27ac98e891198b4ee0498077ab231127dbb3dd8e83a5ea14)
-  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME Rotate_Sample_Ref_Frame_Test.tar.gz SHA512 75c02a11fbb06e9df464df542be8cecc82ec73a7324abf7f12a055570ecc32ec3add6c662ab071868f6708005a6027b89065dbd70605b576520d31e16ef3c372)
+  download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME Rotate_Sample_Ref_Frame_Test_v2.tar.gz SHA512 b7c1a0dca46e133233ef931c2b99879da87ef0facdff77458fdeeb5db2dabd9a477e6699fd006c3676789ca052b0c74e52f7d1f36648f934841672620ec4ce4e)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME point_sample_edge_geometry.tar.gz SHA512 86ee348aac83da0b29c820851bf3cfb03d94d17c8249c8754b3471db3566467c814ce5868526fa35389d32832fbf2991d62c8a5a8de815594bac34176e17d3d6)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME Small_IN100_dream3d_v3.tar.gz SHA512 bfa9547e787b0f8e8122702da0eb4e519f48a48bf4f94aad020f72479d071d32dfc96a1425705874c68507a61ed391d28606d9c4f4acd559043ef0ace64fd33f)
   download_test_data(DREAM3D_DATA_DIR ${DREAM3D_DATA_DIR} ARCHIVE_NAME Small_IN100_h5ebsd.tar.gz SHA512 31e606285ea9e8235dcb5f608fd2b252a5ab1492abd975e5ec33a21d083aa9720fe16fb8f752742c140f40e963d692f1a46256b9d36e96b1b09796c1e4ea3db9)