ENH: Python, wrap structure tensor

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 {}

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);

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

wrapping/CMakeLists.txt

@@ -1,5 +1,6 @@
 itk_wrap_module(IsotropicWavelets)
 set(WRAPPER_SUBMODULE_ORDER
+  itkStructureTensor
   itkMonogenicSignalFrequencyImageFilter
   itkVectorInverseFFTImageFilter
   )

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()
+