Merge pull request #115 from phcerdan/set_riesz_rotation_matrix_to_complex

ENH: Set RieszRotationMatrix ValueType to std::complex

Author: phcerdan

include/itkRieszRotationMatrix.h

@@ -17,10 +17,12 @@
  *=========================================================================*/
 #ifndef itkRieszRotationMatrix_h
 #define itkRieszRotationMatrix_h
-#include "itkVariableSizeMatrix.h"
 #include <vector>
+#include "itkImage.h"
 #include "itkMatrix.h"
 #include "itkRieszUtilities.h"
+#include "itkVariableSizeMatrix.h"
+#include "itkVectorContainer.h"
 
 namespace itk
 {
@@ -35,26 +37,28 @@ namespace itk
  *
  * \f[ M := p(N,d) = \frac{(N+d-1)!}{(d-1)! N!} \f]
  *
- * The rotation matrix is a dxd matrix.
+ * The rotation matrix is a dxd matrix of real type.
  *
  * \sa RieszFrequencyFunction
  * \sa RieszFrequencyFilterBankGenerator
  *
  * \ingroup IsotropicWavelets
  */
 
-template <typename T = double, unsigned int VImageDimension = 3>
-class RieszRotationMatrix : public itk::VariableSizeMatrix<T>
+template <unsigned int VImageDimension>
+class RieszRotationMatrix : public itk::VariableSizeMatrix<std::complex<double>>
 {
 public:
   /** Standard type alias */
   using Self = RieszRotationMatrix;
-  using Superclass = itk::VariableSizeMatrix<T>;
+  using ValueType = std::complex<double>;
+  using Superclass = itk::VariableSizeMatrix<ValueType>;
 
   /** Component value type */
-  using ValueType = typename Superclass::ValueType;
+  using RealType = typename ValueType::value_type;
   using InternalMatrixType = typename Superclass::InternalMatrixType;
-  using SpatialRotationMatrixType = itk::Matrix<T, VImageDimension, VImageDimension>;
+  using SpatialRotationMatrixType = itk::Matrix<RealType, VImageDimension, VImageDimension>;
+  using ComplexImageType = itk::Image<ValueType, VImageDimension>;
 
   /** Matrix by std::vector<TImage> multiplication.
    * To perform the rotation with the output of
@@ -194,7 +198,7 @@ class RieszRotationMatrix : public itk::VariableSizeMatrix<T>
     return this->m_MaxAbsoluteDifferenceCloseToZero;
   }
   inline void
-  SetMaxAbsoluteDifferenceCloseToZero(const ValueType & maxAbsoluteDifference)
+  SetMaxAbsoluteDifferenceCloseToZero(const RealType & maxAbsoluteDifference)
   {
     this->m_MaxAbsoluteDifferenceCloseToZero = maxAbsoluteDifference;
   }
@@ -228,10 +232,11 @@ class RieszRotationMatrix : public itk::VariableSizeMatrix<T>
 #endif
 
 private:
+  using ResultValueType = std::complex<long double>;
   SpatialRotationMatrixType m_SpatialRotationMatrix;
   unsigned int              m_Order{ 0 };
   unsigned int              m_Components{ 0 };
-  ValueType                 m_MaxAbsoluteDifferenceCloseToZero;
+  RealType                  m_MaxAbsoluteDifferenceCloseToZero;
   bool                      m_Debug{ false };
 
 }; // end of class

include/itkRieszRotationMatrix.hxx

@@ -25,42 +25,42 @@
 
 namespace itk
 {
-template <typename T, unsigned int VImageDimension>
-RieszRotationMatrix<T, VImageDimension>::RieszRotationMatrix()
+template <unsigned int VImageDimension>
+RieszRotationMatrix<VImageDimension>::RieszRotationMatrix()
   : Superclass()
   , m_SpatialRotationMatrix()
   ,
 
-  m_MaxAbsoluteDifferenceCloseToZero(1 * itk::NumericTraits<ValueType>::epsilon())
+  m_MaxAbsoluteDifferenceCloseToZero(1 * itk::NumericTraits<RealType>::epsilon())
 
 {}
 
-template <typename T, unsigned int VImageDimension>
-RieszRotationMatrix<T, VImageDimension>::RieszRotationMatrix(const Self & rieszMatrix)
+template <unsigned int VImageDimension>
+RieszRotationMatrix<VImageDimension>::RieszRotationMatrix(const Self & rieszMatrix)
   : Superclass(rieszMatrix)
   , m_SpatialRotationMatrix(rieszMatrix.GetSpatialRotationMatrix())
   , m_Order(rieszMatrix.GetOrder())
   , m_Components(rieszMatrix.GetComponents())
-  , m_MaxAbsoluteDifferenceCloseToZero(1 * itk::NumericTraits<ValueType>::epsilon())
+  , m_MaxAbsoluteDifferenceCloseToZero(1 * itk::NumericTraits<RealType>::epsilon())
 
 {}
 
-template <typename T, unsigned int VImageDimension>
-RieszRotationMatrix<T, VImageDimension>::RieszRotationMatrix(const SpatialRotationMatrixType & spatialRotationMatrix,
-                                                             const unsigned int &              order)
+template <unsigned int VImageDimension>
+RieszRotationMatrix<VImageDimension>::RieszRotationMatrix(const SpatialRotationMatrixType & spatialRotationMatrix,
+                                                          const unsigned int &              order)
   : Superclass()
   , m_SpatialRotationMatrix(spatialRotationMatrix)
-  , m_MaxAbsoluteDifferenceCloseToZero(1 * itk::NumericTraits<ValueType>::epsilon())
+  , m_MaxAbsoluteDifferenceCloseToZero(1 * itk::NumericTraits<RealType>::epsilon())
 
 {
   this->SetOrder(order);
   this->ComputeSteerableMatrix();
 }
 
-template <typename T, unsigned int VImageDimension>
+template <unsigned int VImageDimension>
 template <typename TInputValue>
 std::vector<TInputValue>
-RieszRotationMatrix<T, VImageDimension>::MultiplyWithVector(const std::vector<TInputValue> & vect) const
+RieszRotationMatrix<VImageDimension>::MultiplyWithVector(const std::vector<TInputValue> & vect) const
 {
   unsigned int rows = this->Rows();
   unsigned int cols = this->Cols();
@@ -76,10 +76,10 @@ RieszRotationMatrix<T, VImageDimension>::MultiplyWithVector(const std::vector<TI
   return resultVector;
 }
 
-template <typename T, unsigned int VImageDimension>
+template <unsigned int VImageDimension>
 template <typename TInputValue>
 itk::VariableSizeMatrix<TInputValue>
-RieszRotationMatrix<T, VImageDimension>::MultiplyWithColumnMatrix(
+RieszRotationMatrix<VImageDimension>::MultiplyWithColumnMatrix(
   const itk::VariableSizeMatrix<TInputValue> & inputColumn) const
 {
   unsigned int rows = this->Rows();
@@ -100,10 +100,10 @@ RieszRotationMatrix<T, VImageDimension>::MultiplyWithColumnMatrix(
   return columnMatrix;
 }
 
-template <typename T, unsigned int VImageDimension>
+template <unsigned int VImageDimension>
 template <typename TImage>
 std::vector<typename TImage::Pointer>
-RieszRotationMatrix<T, VImageDimension>::MultiplyWithVectorOfImages(
+RieszRotationMatrix<VImageDimension>::MultiplyWithVectorOfImages(
   const std::vector<typename TImage::Pointer> & vect) const
 {
   unsigned int rows = this->Rows();
@@ -148,9 +148,9 @@ RieszRotationMatrix<T, VImageDimension>::MultiplyWithVectorOfImages(
   return result;
 }
 
-template <typename T, unsigned int VImageDimension>
-typename RieszRotationMatrix<T, VImageDimension>::IndicesMatrix
-RieszRotationMatrix<T, VImageDimension>::GenerateIndicesMatrix()
+template <unsigned int VImageDimension>
+typename RieszRotationMatrix<VImageDimension>::IndicesMatrix
+RieszRotationMatrix<VImageDimension>::GenerateIndicesMatrix()
 {
   using LocalIndicesArrayType = std::vector<unsigned int>;
   using LocalIndicesVector = std::vector<LocalIndicesArrayType>;
@@ -185,9 +185,9 @@ RieszRotationMatrix<T, VImageDimension>::GenerateIndicesMatrix()
   return allIndicesPairs;
 }
 
-template <typename T, unsigned int VImageDimension>
-const typename RieszRotationMatrix<T, VImageDimension>::InternalMatrixType &
-RieszRotationMatrix<T, VImageDimension>::ComputeSteerableMatrix()
+template <unsigned int VImageDimension>
+const typename RieszRotationMatrix<VImageDimension>::InternalMatrixType &
+RieszRotationMatrix<VImageDimension>::ComputeSteerableMatrix()
 {
   // precondition
   if (this->m_Order == 0)
@@ -198,7 +198,13 @@ RieszRotationMatrix<T, VImageDimension>::ComputeSteerableMatrix()
   InternalMatrixType & S = this->GetVnlMatrix();
   if (this->m_Order == 1)
   {
-    S = this->GetSpatialRotationMatrix().GetVnlMatrix().as_matrix();
+    for (unsigned int i = 0; i < this->m_Components; ++i)
+    {
+      for (unsigned int j = 0; j < this->m_Components; ++j)
+      {
+        S.put(i, j, static_cast<ValueType>(m_SpatialRotationMatrix.GetVnlMatrix().get(i, j)));
+      }
+    }
     return this->GetVnlMatrix();
   }
 
@@ -291,19 +297,19 @@ RieszRotationMatrix<T, VImageDimension>::ComputeSteerableMatrix()
       // matrix R = [r1,...,rd], where r_i are columns of the rotation matrix.
       // we sum and normalize them.
       S[i][j] = 0;
-      long double result = 0;
-      long        nFactorial = 1;
-      long        mFactorial = 1;
+      ResultValueType result = 0;
+      long            nFactorial = 1;
+      long            mFactorial = 1;
       for (unsigned int dim = 0; dim < VImageDimension; ++dim)
       {
         nFactorial *= itk::utils::Factorial(n[dim]);
         mFactorial *= itk::utils::Factorial(m[dim]);
       }
-      auto nFactorialReal = static_cast<double>(nFactorial);
+      auto nFactorialReal = static_cast<RealType>(nFactorial);
       for (auto & kValidIndex : kValidIndices)
       {
-        long double rotationFactor = 1;
-        long        kFactorialMultiplication = 1;
+        ValueType rotationFactor = 1;
+        long      kFactorialMultiplication = 1;
         // There are always VImageDimension indices. (k1,k2,...,kd)
         for (unsigned int kIndex = 0; kIndex < VImageDimension; ++kIndex)
         {
@@ -335,12 +341,19 @@ RieszRotationMatrix<T, VImageDimension>::ComputeSteerableMatrix()
       result *= sqrt(mFactorial / nFactorialReal);
       S[i][j] = static_cast<ValueType>(result);
       // Try to fix close to zero float errors
-      if (itk::Math::FloatAlmostEqual(S[i][j],
-                                      static_cast<ValueType>(0),
+      if (itk::Math::FloatAlmostEqual(S[i][j].real(),
+                                      static_cast<typename ValueType::value_type>(0),
+                                      4, // default maxULPS from Math::AlmostFloatEqual
+                                      this->m_MaxAbsoluteDifferenceCloseToZero))
+      {
+        S[i][j].real(0);
+      }
+      if (itk::Math::FloatAlmostEqual(S[i][j].imag(),
+                                      static_cast<typename ValueType::value_type>(0),
                                       4, // default maxULPS from Math::AlmostFloatEqual
                                       this->m_MaxAbsoluteDifferenceCloseToZero))
       {
-        S[i][j] = 0;
+        S[i][j].imag(0);
       }
 
       if (this->GetDebug())

test/itkRieszRotationMatrixTest.cxx

@@ -16,26 +16,27 @@
  *
  *=========================================================================*/
 
-#include "itkRieszRotationMatrix.h"
 #include <complex>
+#include "itkComplexToRealImageFilter.h"
 #include "itkImage.h"
 #include "itkImageDuplicator.h"
 #include "itkMath.h"
+#include "itkRieszRotationMatrix.h"
 #include "itkTestingComparisonImageFilter.h"
 
 template <typename TImage>
 int
-compareImagesAndReport(typename TImage::Pointer validImage, typename TImage::Pointer testImage)
+compareRealImagesAndReport(typename TImage::Pointer validImage, typename TImage::Pointer testImage)
 {
   using ImageType = TImage;
   using ComparisonType = itk::Testing::ComparisonImageFilter<ImageType, ImageType>;
   auto diff = ComparisonType::New();
   diff->SetValidInput(validImage);
   diff->SetTestInput(testImage);
   diff->UpdateLargestPossibleRegion();
-  bool                differenceFailed = false;
-  const double        averageIntensityDifference = diff->GetTotalDifference();
-  const unsigned long numberOfPixelsWithDifferences = diff->GetNumberOfPixelsWithDifferences();
+  bool                             differenceFailed = false;
+  const typename TImage::PixelType averageIntensityDifference = diff->GetTotalDifference();
+  const unsigned long              numberOfPixelsWithDifferences = diff->GetNumberOfPixelsWithDifferences();
   if (averageIntensityDifference > 0.0)
   {
     if (static_cast<int>(numberOfPixelsWithDifferences) > 0)
@@ -67,8 +68,8 @@ runRieszRotationMatrixInterfaceWithRieszFrequencyFilterBankGeneratorTest()
 {
   constexpr unsigned int Dimension = 2;
   // Create a rotation matrix
-  using ValueType = double;
-  using SteerableMatrix = itk::RieszRotationMatrix<ValueType, Dimension>;
+  using ValueType = std::complex<double>;
+  using SteerableMatrix = itk::RieszRotationMatrix<Dimension>;
   using SpatialRotationMatrix = SteerableMatrix::SpatialRotationMatrixType;
   SpatialRotationMatrix R;
   double                angle = itk::Math::pi_over_2;
@@ -94,9 +95,12 @@ runRieszRotationMatrixInterfaceWithRieszFrequencyFilterBankGeneratorTest()
   // emulating the output of RieszFrequencyFilterBankGenerator
   // of Order = 2, Dimension = 2
 
-  using PixelType = double;
+  using PixelType = std::complex<double>;
+  using RealPixelType = double;
   using ImageType = itk::Image<PixelType, Dimension>;
   using ImagePointer = ImageType::Pointer;
+  using RealImageType = itk::Image<RealPixelType, Dimension>;
+  using RealImagePointer = RealImageType::Pointer;
   ImagePointer         image = ImageType::New();
   ImageType::IndexType start;
   start.Fill(0);
@@ -137,7 +141,7 @@ runRieszRotationMatrixInterfaceWithRieszFrequencyFilterBankGeneratorTest()
   expectedMultiplyResult[2] = 1;
   for (unsigned int i = 0; i < M; ++i)
   {
-    if (!itk::Math::FloatAlmostEqual(vectorRotated[i], expectedMultiplyResult[i]))
+    if (!itk::Math::FloatAlmostEqual(vectorRotated[i].real(), expectedMultiplyResult[i].real()))
     {
       std::cout << "vectorRotated not equal!: ";
       std::cout << vectorRotated[i] << " != " << expectedMultiplyResult[i] << std::endl;
@@ -152,7 +156,7 @@ runRieszRotationMatrixInterfaceWithRieszFrequencyFilterBankGeneratorTest()
   auto columnMatrixRotated = S.MultiplyWithColumnMatrix(inputColumnMatrix);
   for (unsigned int i = 0; i < M; ++i)
   {
-    if (!itk::Math::FloatAlmostEqual(columnMatrixRotated.GetVnlMatrix()(i, 0), expectedMultiplyResult[i]))
+    if (!itk::Math::FloatAlmostEqual(columnMatrixRotated.GetVnlMatrix()(i, 0).real(), expectedMultiplyResult[i].real()))
     {
       std::cout << "columnMatrixRotated not Equal!: ";
       std::cout << columnMatrixRotated.GetVnlMatrix()(i, 0) << " != " << expectedMultiplyResult[i] << std::endl;
@@ -164,22 +168,50 @@ runRieszRotationMatrixInterfaceWithRieszFrequencyFilterBankGeneratorTest()
   std::cout << "Size: ";
   std::cout << imagesMultipliedByRieszRotationMatrix.size() << std::endl;
 
+  // Convert input images to real to perform comparison.
+  std::vector<RealImagePointer> realImages(M);
+  {
+    for (unsigned int i = 0; i < M; ++i)
+    {
+      using ComplexToRealFilter = itk::ComplexToRealImageFilter<ImageType, RealImageType>;
+      auto complexToRealFilter = ComplexToRealFilter::New();
+      complexToRealFilter->SetInput(images[i]);
+      complexToRealFilter->Update();
+      realImages[i] = complexToRealFilter->GetOutput();
+    }
+  }
+  std::vector<RealImagePointer> realImagesMultipliedByRieszRotationMatrix(M);
+  {
+    for (unsigned int i = 0; i < M; ++i)
+    {
+      using ComplexToRealFilter = itk::ComplexToRealImageFilter<ImageType, RealImageType>;
+      auto complexToRealFilter = ComplexToRealFilter::New();
+      complexToRealFilter->SetInput(imagesMultipliedByRieszRotationMatrix[i]);
+      complexToRealFilter->Update();
+      realImagesMultipliedByRieszRotationMatrix[i] = complexToRealFilter->GetOutput();
+    }
+  }
+
   // First
   // g(0, 0) = 0; g(0, 1) = 0;  g(0, 2) = 1;
   // result = 1.0 * images[2]
-  int firstComponentStatus = compareImagesAndReport<ImageType>(images[2], imagesMultipliedByRieszRotationMatrix[0]);
+  int firstComponentStatus =
+    compareRealImagesAndReport<RealImageType>(realImages[2], realImagesMultipliedByRieszRotationMatrix[0]);
   // Second
   // g(1, 0) = 0; g(1, 1) = -1; g(1, 2) = 0;
   // expectedresult = -1.0 * images[1]
-  typename MultiplyImageFilterType::Pointer multiplyImageFilterInvert = MultiplyImageFilterType::New();
-  multiplyImageFilterInvert->SetInput(imagesMultipliedByRieszRotationMatrix[1]);
+  using MultiplyRealImageFilterType = itk::MultiplyImageFilter<RealImageType>;
+  typename MultiplyRealImageFilterType::Pointer multiplyImageFilterInvert = MultiplyRealImageFilterType::New();
+  multiplyImageFilterInvert->SetInput(realImagesMultipliedByRieszRotationMatrix[1]);
   multiplyImageFilterInvert->SetConstant(-1.0);
   multiplyImageFilterInvert->Update();
-  int secondComponentStatus = compareImagesAndReport<ImageType>(images[1], multiplyImageFilterInvert->GetOutput());
+  int secondComponentStatus =
+    compareRealImagesAndReport<RealImageType>(realImages[1], multiplyImageFilterInvert->GetOutput());
   // Third
   // g(2, 0) = 1; g(2, 1) = 0;  g(2, 2) = 0;
   // result = 1.0 * images[0]
-  int thirdComponentStatus = compareImagesAndReport<ImageType>(images[0], imagesMultipliedByRieszRotationMatrix[2]);
+  int thirdComponentStatus =
+    compareRealImagesAndReport<RealImageType>(realImages[0], realImagesMultipliedByRieszRotationMatrix[2]);
 
   if (!multiplyWithSomethingPassed)
   {
@@ -195,8 +227,8 @@ runRieszRotationMatrixTest()
   bool               testPassed = true;
   const unsigned int Dimension = VDimension;
 
-  using ValueType = double;
-  using SteerableMatrix = itk::RieszRotationMatrix<ValueType, Dimension>;
+  using ValueType = std::complex<double>;
+  using SteerableMatrix = itk::RieszRotationMatrix<Dimension>;
   using SpatialRotationMatrix = typename SteerableMatrix::SpatialRotationMatrixType;
 
   // Define a spatial rotation matrix.

wrapping/CMakeLists.txt

@@ -1,5 +1,7 @@
 itk_wrap_module(IsotropicWavelets)
 set(WRAPPER_SUBMODULE_ORDER
+  itkMatrixComplex
+  itkVariableSizeMatrixComplex
   itkStructureTensor
   itkMonogenicSignalFrequencyImageFilter
   itkVectorInverseFFTImageFilter