Move RotateSampleRefFrame to algorithm class

Author: nyoungbq

src/Plugins/SimplnxCore/CMakeLists.txt

@@ -231,6 +231,7 @@ set(AlgorithmList
   ResampleImageGeom
   ResampleRectGridToImageGeom
   ReshapeDataArray
+  RotateSampleRefFrame
   ScalarSegmentFeatures
   SharedFeatureFace
   Silhouette

src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/RotateSampleRefFrame.cpp

@@ -0,0 +1,329 @@
+#include "RotateSampleRefFrame.hpp"
+
+#include "simplnx/Utilities/DataGroupUtilities.hpp"
+#include "simplnx/Utilities/ImageRotationUtilities.hpp"
+#include "simplnx/Utilities/ParallelAlgorithmUtilities.hpp"
+#include "simplnx/Utilities/ParallelTaskAlgorithm.hpp"
+
+#include <Eigen/Dense>
+
+using namespace nx::core;
+
+namespace
+{
+using RotationRepresentationType = RotateSampleRefFrame::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
+
+// -----------------------------------------------------------------------------
+RotateSampleRefFrame::RotateSampleRefFrame(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, RotateSampleRefFrameInputValues* inputValues)
+: m_DataStructure(dataStructure)
+, m_InputValues(inputValues)
+, m_ShouldCancel(shouldCancel)
+, m_MessageHandler(mesgHandler)
+{
+}
+
+// -----------------------------------------------------------------------------
+RotateSampleRefFrame::~RotateSampleRefFrame() noexcept = default;
+
+// -----------------------------------------------------------------------------
+void RotateSampleRefFrame::updateProgress(const std::string& message)
+{
+  m_MessageHandler(IFilter::Message::Type::Info, message);
+}
+
+// -----------------------------------------------------------------------------
+const std::atomic_bool& RotateSampleRefFrame::getCancel()
+{
+  return m_ShouldCancel;
+}
+
+// -----------------------------------------------------------------------------
+Result<> RotateSampleRefFrame::operator()()
+{
+  auto& srcImageGeom = m_DataStructure.getDataRefAs<ImageGeom>(m_InputValues->SourceGeometryPath);
+  auto& destImageGeom = m_DataStructure.getDataRefAs<ImageGeom>(m_InputValues->DestGeometryPath);
+
+  ImageRotationUtilities::Matrix4fR rotationMatrix;
+
+  switch(static_cast<RotationRepresentationType>(m_InputValues->RotationRepresentationIndex))
+  {
+  case RotationRepresentationType::AxisAngle: {
+    rotationMatrix = ImageRotationUtilities::GenerateRotationTransformationMatrix(m_InputValues->RotationAxisAngle);
+    break;
+  }
+  case RotationRepresentationType::RotationMatrix: {
+    rotationMatrix = ImageRotationUtilities::GenerateManualTransformationMatrix(m_InputValues->RotationMatrixTable);
+    break;
+  }
+  }
+
+  ImageRotationUtilities::RotateArgs rotateArgs = ImageRotationUtilities::CreateRotationArgs(srcImageGeom, rotationMatrix);
+
+  auto selectedCellDataChildren = GetAllChildArrayDataPaths(m_DataStructure, srcImageGeom.getCellDataPath());
+  auto selectedCellArrays = selectedCellDataChildren.has_value() ? selectedCellDataChildren.value() : std::vector<DataPath>{};
+
+  ImageRotationUtilities::FilterProgressCallback filterProgressCallback(m_MessageHandler, m_ShouldCancel);
+
+  // The actual rotating of the dataStructure arrays is done in parallel where parallel here
+  // refers to the cropping of each DataArray being done on a separate thread.
+  ParallelTaskAlgorithm taskRunner;
+  taskRunner.setParallelizationEnabled(true);
+  const DataPath srcCelLDataAMPath = srcImageGeom.getCellDataPath();
+  const auto& srcCellDataAM = srcImageGeom.getCellDataRef();
+
+  const DataPath destCellDataAMPath = destImageGeom.getCellDataPath();
+
+  for(const auto& [dataId, srcDataObject] : srcCellDataAM)
+  {
+    if(m_ShouldCancel)
+    {
+      return {};
+    }
+
+    const auto* srcDataArray = m_DataStructure.getDataAs<IDataArray>(srcCelLDataAMPath.createChildPath(srcDataObject->getName()));
+    auto* destDataArray = m_DataStructure.getDataAs<IDataArray>(destCellDataAMPath.createChildPath(srcDataObject->getName()));
+    m_MessageHandler(fmt::format("Rotating Volume || Copying Data Array {}", srcDataObject->getName()));
+
+    ExecuteParallelFunction<ImageRotationUtilities::RotateImageGeometryWithNearestNeighbor>(srcDataArray->getDataType(), taskRunner, srcDataArray, destDataArray, rotateArgs, rotationMatrix,
+                                                                                            m_InputValues->SliceBySlice, &filterProgressCallback);
+  }
+
+  taskRunner.wait(); // This will spill over if the number of DataArrays to process does not divide evenly by the number of threads.
+
+  return {};
+}

src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/RotateSampleRefFrame.hpp

@@ -0,0 +1,56 @@
+#pragma once
+
+#include "SimplnxCore/SimplnxCore_export.hpp"
+
+#include "simplnx/DataStructure/DataPath.hpp"
+#include "simplnx/DataStructure/DataStructure.hpp"
+#include "simplnx/Filter/IFilter.hpp"
+#include "simplnx/Parameters/ChoicesParameter.hpp"
+#include "simplnx/Parameters/DynamicTableParameter.hpp"
+#include "simplnx/Parameters/VectorParameter.hpp"
+
+namespace nx::core
+{
+
+struct SIMPLNXCORE_EXPORT RotateSampleRefFrameInputValues
+{
+  bool SliceBySlice;
+  ChoicesParameter::ValueType RotationRepresentationIndex;
+  DataPath SourceGeometryPath;
+  DataPath DestGeometryPath;
+  VectorFloat32Parameter::ValueType RotationAxisAngle;
+  DynamicTableParameter::ValueType RotationMatrixTable;
+};
+
+/**
+ * @class
+ */
+class SIMPLNXCORE_EXPORT RotateSampleRefFrame
+{
+public:
+  RotateSampleRefFrame(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, RotateSampleRefFrameInputValues* inputValues);
+  ~RotateSampleRefFrame() noexcept;
+
+  RotateSampleRefFrame(const RotateSampleRefFrame&) = delete;
+  RotateSampleRefFrame(RotateSampleRefFrame&&) noexcept = delete;
+  RotateSampleRefFrame& operator=(const RotateSampleRefFrame&) = delete;
+  RotateSampleRefFrame& operator=(RotateSampleRefFrame&&) noexcept = delete;
+
+  enum class RotationRepresentation : uint64
+  {
+    AxisAngle = 0,
+    RotationMatrix = 1
+  };
+
+  Result<> operator()();
+  void updateProgress(const std::string& message);
+  const std::atomic_bool& getCancel();
+
+private:
+  DataStructure& m_DataStructure;
+  const RotateSampleRefFrameInputValues* m_InputValues = nullptr;
+  const std::atomic_bool& m_ShouldCancel;
+  const IFilter::MessageHandler& m_MessageHandler;
+};
+
+} // namespace nx::core