Merge pull request #116 from phcerdan/add_multiply_methods_to_riesz_rotation_matrix

ENH: Add multiplyWith methods to RieszRotationMatrix

Author: phcerdan

include/itkRieszRotationMatrix.h

@@ -64,6 +64,12 @@ class RieszRotationMatrix:
   template <typename TImage>
   std::vector< typename TImage::Pointer > MultiplyWithVectorOfImages(const std::vector< typename TImage::Pointer > & vect) const;
 
+  template <typename TInputValue>
+  std::vector< TInputValue > MultiplyWithVector(const std::vector< TInputValue > & vect) const;
+
+  template <typename TInputValue>
+  VariableSizeMatrix< TInputValue > MultiplyWithColumnMatrix(const VariableSizeMatrix< TInputValue > & vect) const;
+
   /**
    * Multi-index notation
    * S[n = (n1,...,nd)][m = (m1,...,md)]

include/itkRieszRotationMatrix.hxx

@@ -64,6 +64,50 @@ RieszRotationMatrix< T, VImageDimension >
   this->ComputeSteerableMatrix();
 }
 
+template< typename T, unsigned int VImageDimension >
+template <typename TInputValue>
+std::vector<TInputValue>
+RieszRotationMatrix< T, VImageDimension >
+::MultiplyWithVector(const std::vector<TInputValue> & vect) const
+{
+  unsigned int rows = this->Rows();
+  unsigned int cols = this->Cols();
+  auto resultVector = std::vector<TInputValue>(rows, NumericTraits<TInputValue>::ZeroValue());
+
+  for ( unsigned int r = 0; r < rows; r++ )
+    {
+      for ( unsigned int c = 0; c < cols; c++ )
+      {
+        resultVector[r] += this->GetVnlMatrix()(r, c) * vect[c];
+      }
+    }
+  return resultVector;
+}
+
+template< typename T, unsigned int VImageDimension >
+template <typename TInputValue>
+itk::VariableSizeMatrix<TInputValue>
+RieszRotationMatrix< T, VImageDimension >
+::MultiplyWithColumnMatrix(const itk::VariableSizeMatrix<TInputValue> & inputColumn) const
+{
+  unsigned int rows = this->Rows();
+  unsigned int cols = this->Cols();
+  using ColumnMatrix = VariableSizeMatrix<TInputValue>;
+  ColumnMatrix columnMatrix(rows, 1);
+  columnMatrix.Fill(NumericTraits<TInputValue>::ZeroValue());
+
+  for ( unsigned int r = 0; r < rows; r++ )
+    {
+      TInputValue sum = NumericTraits<TInputValue>::ZeroValue();
+      for ( unsigned int c = 0; c < cols; c++ )
+      {
+        sum += this->GetVnlMatrix()(r, c) * inputColumn.GetVnlMatrix()(c, 0);
+      }
+      columnMatrix.GetVnlMatrix()(r, 0) = sum;
+    }
+  return columnMatrix;
+}
+
 template< typename T, unsigned int VImageDimension >
 template <typename TImage>
 std::vector< typename TImage::Pointer >

test/itkRieszRotationMatrixTest.cxx

@@ -129,6 +129,28 @@ runRieszRotationMatrixInterfaceWithRieszFrequencyFilterBankGeneratorTest()
     std::cout << images[i]->GetLargestPossibleRegion() << std::endl;
   }
 
+  std::vector<PixelType> inputVector(M, 0);
+  inputVector[0] = 1;
+  inputVector[1] = 2;
+  inputVector[2] = 3;
+  auto vectorRotated = S.MultiplyWithVector(inputVector);
+
+  itk::VariableSizeMatrix<PixelType> inputColumnMatrix(M, 1);
+  inputColumnMatrix.GetVnlMatrix()(0,0) = 1;
+  inputColumnMatrix.GetVnlMatrix()(1,0) = 2;
+  inputColumnMatrix.GetVnlMatrix()(2,0) = 3;
+  auto columnMatrixRotated = S.MultiplyWithColumnMatrix(inputColumnMatrix);
+  bool multiplyWithSomethingPassed = true;
+  for (unsigned int i = 0; i < M; ++i)
+    {
+    if(!itk::Math::FloatAlmostEqual(inputVector[i], columnMatrixRotated.GetVnlMatrix()(i, 0)))
+      {
+      std::cout << "Not Equal!: ";
+      std::cout << inputVector[i] << " != " << columnMatrixRotated.GetVnlMatrix()(i,0) << std::endl;
+      multiplyWithSomethingPassed = false;
+      }
+    }
+
   auto imagesMultipliedByRieszRotationMatrix =
     S.MultiplyWithVectorOfImages<ImageType>(images);
   std::cout << "Size: ";
@@ -155,6 +177,10 @@ runRieszRotationMatrixInterfaceWithRieszFrequencyFilterBankGeneratorTest()
   int thirdComponentStatus = compareImagesAndReport<ImageType>(
       images[0], imagesMultipliedByRieszRotationMatrix[2]);
 
+  if(!multiplyWithSomethingPassed)
+    {
+    return EXIT_FAILURE;
+    }
   return firstComponentStatus && secondComponentStatus && thirdComponentStatus;
 
 }