ENH: Wrap StructureTensor and adds GetRotationMatrix

The type of the output is defaulted to VariableSizeMatrix<double> that is already wrapped in ITK. It wraps the needed ImageToImageFilter and ImageSource Also an utility function GetRotationMatrixFromOutputMatrix is added and tested.
InsightSoftwareConsortium · Jan 14, 2019 · 4b481f0 · 4b481f0
1 parent 18bc327
commit 4b481f0
Show file tree

Hide file tree

Showing 5 changed files with 116 additions and 4 deletions.
diff --git a/include/itkStructureTensor.h b/include/itkStructureTensor.h
@@ -77,7 +77,7 @@ and
  */
 template<typename TInputImage,
          typename TOutputImage = itk::Image<
-           itk::VariableSizeMatrix<typename TInputImage::PixelType>,
+           itk::VariableSizeMatrix<double>,
            TInputImage::ImageDimension> >
 class StructureTensor:
   public ImageToImageFilter< TInputImage, TOutputImage >
@@ -185,6 +185,23 @@ class StructureTensor:
    */
   InputImagePointer ComputeCoherencyImage() const;
 
+  /**
+   * From the output matrix at each index, computes the rotation matrix,
+   * which is the transpose of the eigenvector matrix computed in the output of
+   * this image filter.
+   * If reOrderLargestEigenvectorInFirstRow is true the first row of the rotation
+   * matrix will have the largest eigenvector. If false (the default), the largest
+   * eigenvector will be in the last row.
+   *
+   * @param outputMatrix matrix at any index stored in the output of this filter.
+   * @param reOrderLargestEigenvectorInFirstRow flag to reorder eigenvectors.
+   *
+   * @return rotationMatrix
+   */
+  EigenMatrixType GetRotationMatrixFromOutputMatrix(
+      const EigenMatrixType & outputMatrix,
+      bool reOrderLargestEigenvectorInFirstRow = false) const;
+
 protected:
   StructureTensor();
   ~StructureTensor() override {}

diff --git a/include/itkStructureTensor.hxx b/include/itkStructureTensor.hxx
@@ -230,6 +230,36 @@ StructureTensor< TInputImage, TOutputImage >
     } // end outIt
 }
 
+template< typename TInputImage, typename TOutputImage >
+typename StructureTensor< TInputImage, TOutputImage >::EigenMatrixType
+StructureTensor< TInputImage, TOutputImage >
+::GetRotationMatrixFromOutputMatrix(
+  const EigenMatrixType & outputMatrix,
+  bool reOrderLargestEigenvectorInFirstRow) const
+{
+  unsigned int nInputs = this->GetNumberOfInputs();
+  // Pre condition (output matrix is populated)
+  assert(outputMatrix.Rows() == nInputs);
+  EigenMatrixType rotationMatrix(nInputs, nInputs);
+  // transpose at copy
+  if(!reOrderLargestEigenvectorInFirstRow)
+    {
+    for ( unsigned int n = 0; n < nInputs; ++n )
+      {
+      rotationMatrix.GetVnlMatrix().set_row(n, outputMatrix.GetVnlMatrix().get_column(n));
+      }
+    }
+  // reOrder and transpose at copy
+  else
+    {
+    for ( unsigned int n = 0; n < nInputs; ++n )
+      {
+      rotationMatrix.GetVnlMatrix().set_row(nInputs - 1 - n, outputMatrix.GetVnlMatrix().get_column(n));
+      }
+    }
+  return rotationMatrix;
+}
+
 template< typename TInputImage, typename TOutputImage >
 typename StructureTensor< TInputImage, TOutputImage >::InputImagePointer
 StructureTensor< TInputImage, TOutputImage >
@@ -274,9 +304,10 @@ StructureTensor< TInputImage, TOutputImage >
     while ( !outIt.IsAtEndOfLine() )
       {
       InputImagePixelType value = 0;
+      auto outputMatrix = outIt.Get();
       for ( unsigned int r = 0; r < nInputs; r++ )
         {
-        value += outIt.Get()[r][eigen_number] * inputIts[r].Get();
+        value += outputMatrix[r][eigen_number] * inputIts[r].Get();
         ++inputIts[r];
         }
       projectIt.Set(value);

diff --git a/test/itkStructureTensorTest.cxx b/test/itkStructureTensorTest.cxx
@@ -132,14 +132,39 @@ runStructureTensorTest()
   tensor->Update();
 
   auto eigenImage = tensor->GetOutput();
-  unsigned eigenMatrixRows = eigenImage->GetPixel(start).Rows();
-  unsigned eigenMatrixCols = eigenImage->GetPixel(start).Cols();
+
+  auto outMatrixAtStart = eigenImage->GetPixel(start);
+  unsigned eigenMatrixRows = outMatrixAtStart.Rows();
+  unsigned eigenMatrixCols = outMatrixAtStart.Cols();
   if ( eigenMatrixRows != nInputs || eigenMatrixCols != nInputs + 1 )
     {
     testFailed = true;
     std::cout << "The resulting eigenMatrix size is wrong. Rows: " << eigenMatrixRows << " . Columns: "
               << eigenMatrixCols << std::endl;
     }
+  auto rotationMatrix = tensor->GetRotationMatrixFromOutputMatrix(outMatrixAtStart);
+  auto rotationMatrixReOrdered = tensor->GetRotationMatrixFromOutputMatrix(outMatrixAtStart, true);
+  for (unsigned int r = 0; r < nInputs; ++r)
+    {
+    for (unsigned int c = 0; c < nInputs; ++c)
+      {
+      if(rotationMatrix[r][c] != outMatrixAtStart[c][r])
+        {
+        testFailed = true;
+        std::cout << "GetRotationMatrixFromOutputMatrix fails. rotationMatrix[r][c] "
+          << rotationMatrix[r][c] << " !=  outMatrixAtStart[c][r] "
+          << outMatrixAtStart[c][r] << std::endl;
+        }
+      if(rotationMatrixReOrdered[r][c] != outMatrixAtStart[nInputs - 1 - c][r])
+        {
+        testFailed = true;
+        std::cout << "GetRotationMatrixFromOutputMatrix fails. rotationMatrixReOrdered[r][c] "
+          << rotationMatrixReOrdered[r][c] << " !=  outMatrixAtStart[nInputs - 1 - c][r] "
+          << outMatrixAtStart[nInputs - 1 - c][r] << std::endl;
+        }
+      }
+    }
+
 #ifdef ITK_VISUALIZE_TESTS
   itk::ViewImage<ImageType>::View( tensor->GetGaussianSource()->GetOutput(), "Gaussian" );
 #endif

diff --git a/wrapping/CMakeLists.txt b/wrapping/CMakeLists.txt
@@ -1,5 +1,6 @@
 itk_wrap_module(IsotropicWavelets)
 set(WRAPPER_SUBMODULE_ORDER
+  itkStructureTensor
   itkMonogenicSignalFrequencyImageFilter
   itkVectorInverseFFTImageFilter
   )

diff --git a/wrapping/itkStructureTensor.wrap b/wrapping/itkStructureTensor.wrap
@@ -0,0 +1,38 @@
+itk_wrap_include("itkVariableSizeMatrix.h")
+
+itk_wrap_class("itk::Image" POINTER)
+  foreach(d ${ITK_WRAP_IMAGE_DIMS})
+    foreach(t ${WRAP_ITK_REAL})
+      itk_wrap_template("VSM${ITKM_I${t}}${d}"
+        "itk::VariableSizeMatrix< ${ITKT_D} >, ${d}")
+    endforeach()
+  endforeach()
+itk_end_wrap_class()
+
+itk_wrap_class("itk::ImageSource" POINTER)
+  foreach(d ${ITK_WRAP_IMAGE_DIMS})
+    foreach(t ${WRAP_ITK_REAL})
+      itk_wrap_template("VSM${ITKM_I${ITKM_${t}}}${d}"
+        "itk::Image<itk::VariableSizeMatrix< ${ITKT_D} >, ${d} >")
+    endforeach()
+  endforeach()
+itk_end_wrap_class()
+
+itk_wrap_class("itk::ImageToImageFilter" POINTER)
+  foreach(d ${ITK_WRAP_IMAGE_DIMS})
+    foreach(t ${WRAP_ITK_REAL})
+      itk_wrap_template("${ITKM_I${t}${d}}VSM${ITKM_I${ITKM_${t}}}${d}"
+        "${ITKT_I${t}${d}}, itk::Image<itk::VariableSizeMatrix< ${ITKT_D} >, ${d} >")
+    endforeach()
+  endforeach()
+itk_end_wrap_class()
+
+itk_wrap_class("itk::StructureTensor" POINTER)
+  foreach(d ${ITK_WRAP_IMAGE_DIMS})
+    foreach(t ${WRAP_ITK_REAL})
+      itk_wrap_template("${ITKM_I${t}${d}}VSM${ITKM_I${ITKM_${t}}}${d}"
+        "${ITKT_I${t}${d}}, itk::Image<itk::VariableSizeMatrix< ${ITKT_D} >, ${d} >")
+    endforeach()
+  endforeach()
+itk_end_wrap_class()
+