BUG: Jacobian determinant filter results depended on index space orientation

Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.h

@@ -151,6 +151,10 @@ class ITK_TEMPLATE_EXPORT DisplacementFieldJacobianDeterminantFilter
   using RealVectorType = Vector<TRealType, InputPixelType::Dimension>;
   using RealVectorImageType = Image<RealVectorType, TInputImage::ImageDimension>;
 
+  /** Matrix type used to store the input direction matrix and the gradient matrix from
+      which the Jacobian determinant is computed. */
+  using InternalMatrixType = vnl_matrix_fixed<TRealType, ImageDimension, ImageDimension>;
+
   /** Type of the iterator that will be used to move through the image.  Also
       the type which will be passed to the evaluate function */
   using ConstNeighborhoodIteratorType = ConstNeighborhoodIterator<RealVectorImageType>;
@@ -271,6 +275,8 @@ class ITK_TEMPLATE_EXPORT DisplacementFieldJacobianDeterminantFilter
   typename ImageBaseType::ConstPointer m_RealValuedInputImage{};
 
   RadiusType m_NeighborhoodRadius{};
+
+  InternalMatrixType m_InputDirection{};
 };
 } // end namespace itk
 

Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.hxx

@@ -138,6 +138,16 @@ DisplacementFieldJacobianDeterminantFilter<TInputImage, TRealType, TOutputImage>
 {
   Superclass::BeforeThreadedGenerateData();
 
+  // Cache the direction matrix from the input
+  const auto & dir_d = this->GetInput()->GetDirection().GetVnlMatrix();
+  for (unsigned int i = 0; i < ImageDimension; ++i)
+  {
+    for (unsigned int j = 0; j < ImageDimension; ++j)
+    {
+      m_InputDirection(i, j) = static_cast<TRealType>(dir_d(i, j));
+    }
+  }
+
   // Set the weights on the derivatives.
   // Are we using image spacing in the calculations?  If so we must update now
   // in case our input image has changed.
@@ -208,18 +218,40 @@ TRealType
 DisplacementFieldJacobianDeterminantFilter<TInputImage, TRealType, TOutputImage>::EvaluateAtNeighborhood(
   const ConstNeighborhoodIteratorType & it) const
 {
-  vnl_matrix_fixed<TRealType, ImageDimension, VectorDimension> J;
-  for (unsigned int i = 0; i < ImageDimension; ++i)
+  InternalMatrixType localGrad;
+  InternalMatrixType physicalGrad;
+
+  for (unsigned int col = 0; col < ImageDimension; ++col)
+  {
+    // Take finite difference along index axis col
+    for (unsigned int row = 0; row < ImageDimension; ++row)
+    {
+      TRealType localComponentGrad = m_HalfDerivativeWeights[col] * (it.GetNext(col)[row] - it.GetPrevious(col)[row]);
+
+      localGrad[row][col] = localComponentGrad;
+    }
+  }
+
+  for (unsigned int row = 0; row < ImageDimension; ++row)
   {
-    for (unsigned int j = 0; j < VectorDimension; ++j)
+    vnl_vector_fixed<TRealType, ImageDimension> localComponentGrad_d;
+    for (unsigned int col = 0; col < ImageDimension; ++col)
+    {
+      localComponentGrad_d[col] = localGrad[row][col];
+    }
+
+    vnl_vector_fixed<TRealType, ImageDimension> physicalComponentGrad = m_InputDirection * localComponentGrad_d;
+
+    for (unsigned int col = 0; col < ImageDimension; ++col)
     {
-      J[i][j] = m_HalfDerivativeWeights[i] * (it.GetNext(i)[j] - it.GetPrevious(i)[j]);
+      physicalGrad[row][col] = physicalComponentGrad[col];
     }
-    // add one on the diagonal to consider the warping and not only the
-    // deformation field
-    J[i][i] += 1.0;
+
+    physicalGrad[row][row] += itk::NumericTraits<TRealType>::OneValue();
   }
-  return vnl_det(J);
+
+  // Return determinant of physical Jacobian:
+  return vnl_determinant(physicalGrad);
 }
 
 template <typename TInputImage, typename TRealType, typename TOutputImage>

Modules/Filtering/DisplacementField/test/itkDisplacementFieldJacobianDeterminantFilterTest.cxx

@@ -18,7 +18,9 @@
 
 #include <iostream>
 #include "itkDisplacementFieldJacobianDeterminantFilter.h"
+#include "itkFlipImageFilter.h"
 #include "itkNullImageToImageFilterDriver.hxx"
+#include "itkPermuteAxesImageFilter.h"
 #include "itkStdStreamStateSave.h"
 #include "itkTestingMacros.h"
 
@@ -36,23 +38,23 @@ TestDisplacementJacobianDeterminantValue()
   using VectorImageType = FieldType;
 
 
-  std::cout << "Create the dispacementfield image pattern." << std::endl;
+  std::cout << "Create the displacementField image pattern." << std::endl;
   VectorImageType::RegionType region;
   // NOTE:  Making the image size much larger than necessary in order to get
   //       some meaningful time measurements.  Simulate a 256x256x256 image.
   constexpr VectorImageType::SizeType size = { { 4096, 4096 } };
   region.SetSize(size);
 
-  auto dispacementfield = VectorImageType::New();
-  dispacementfield->SetLargestPossibleRegion(region);
-  dispacementfield->SetBufferedRegion(region);
-  dispacementfield->Allocate();
+  auto displacementField = VectorImageType::New();
+  displacementField->SetLargestPossibleRegion(region);
+  displacementField->SetBufferedRegion(region);
+  displacementField->Allocate();
 
   VectorType values;
   values[0] = 0;
   values[1] = 0;
   using Iterator = itk::ImageRegionIteratorWithIndex<VectorImageType>;
-  for (Iterator inIter(dispacementfield, region); !inIter.IsAtEnd(); ++inIter)
+  for (Iterator inIter(displacementField, region); !inIter.IsAtEnd(); ++inIter)
   {
     const unsigned int i = inIter.GetIndex()[0];
     const unsigned int j = inIter.GetIndex()[1];
@@ -94,7 +96,7 @@ TestDisplacementJacobianDeterminantValue()
 #endif
   ITK_TEST_SET_GET_BOOLEAN(filter, UseImageSpacing, useImageSpacing);
 
-  filter->SetInput(dispacementfield);
+  filter->SetInput(displacementField);
   filter->Update();
 
 
@@ -118,7 +120,113 @@ TestDisplacementJacobianDeterminantValue()
   }
   else
   {
-    std::cout << "Test passed." << std::endl;
+    std::cout << "Determinant value test passed" << std::endl;
+  }
+
+  // Now test that the determinant is consistent if the voxel space is different, but the physical space is the same.
+
+  // Define physical point to test
+  itk::Point<double, 2> physPt;
+  displacementField->TransformIndexToPhysicalPoint(index, physPt);
+
+  // Use this function to get the determinant at a specified physical point (it will be a different index between tests)
+  auto GetDeterminantAtPoint = [](FieldType::Pointer image, const itk::Point<double, 2> & pt) -> float {
+    auto filter = FilterType::New();
+    filter->SetInput(image);
+    filter->SetUseImageSpacing(true);
+    filter->Update();
+
+    itk::Index<2> mappedIdx = image->TransformPhysicalPointToIndex(pt);
+
+    return filter->GetOutput()->GetPixel(mappedIdx);
+  };
+
+  const float detOriginal = GetDeterminantAtPoint(displacementField, physPt);
+
+  // Check that this is the same as above, just to be sure the function above works
+  if (itk::Math::abs(detOriginal - expectedJacobianDeterminant) > epsilon)
+  {
+    std::cerr.precision(static_cast<int>(itk::Math::abs(std::log10(epsilon))));
+    std::cerr << "Test failed " << std::endl;
+    std::cerr << "Error in pixel value at physical point " << physPt << std::endl;
+    std::cerr << "Expected value " << detOriginal << std::endl;
+    std::cerr << " differs from " << expectedJacobianDeterminant;
+    std::cerr << " by more than " << epsilon << std::endl;
+    testPassed = false;
+  }
+
+  using PermuteFilterType = itk::PermuteAxesImageFilter<FieldType>;
+  auto permute = PermuteFilterType::New();
+
+  PermuteFilterType::PermuteOrderArrayType order;
+  order[0] = 1;
+  order[1] = 0; // swap X <-> Y
+  permute->SetOrder(order);
+  permute->SetInput(displacementField);
+  permute->Update();
+
+  auto displacementFieldPermuted = permute->GetOutput();
+
+  // Adjust direction to match permutation
+  itk::Matrix<double, 2, 2> origDir = displacementField->GetDirection();
+  itk::Matrix<double, 2, 2> newDirPermute;
+  newDirPermute[0][0] = origDir[1][0];
+  newDirPermute[0][1] = origDir[1][1];
+  newDirPermute[1][0] = origDir[0][0];
+  newDirPermute[1][1] = origDir[0][1];
+
+  displacementFieldPermuted->SetDirection(newDirPermute);
+  displacementFieldPermuted->SetOrigin(displacementField->GetOrigin());
+  displacementFieldPermuted->SetSpacing(displacementField->GetSpacing());
+
+  const float detPermuted = GetDeterminantAtPoint(displacementFieldPermuted, physPt);
+
+  using FlipFilterType = itk::FlipImageFilter<FieldType>;
+  auto flip = FlipFilterType::New();
+
+  FlipFilterType::FlipAxesArrayType flipAxes;
+  flipAxes[0] = true;  // Flip X
+  flipAxes[1] = false; // Keep Y
+  flip->SetFlipAxes(flipAxes);
+  flip->SetInput(displacementField);
+  flip->FlipAboutOriginOff();
+  flip->Update();
+
+  auto displacementFieldFlip = flip->GetOutput();
+
+  // Adjust direction to compensate flip
+  itk::Matrix<double, 2, 2> flipMat;
+  flipMat.SetIdentity();
+  flipMat[0][0] = -1.0;
+
+  itk::Matrix<double, 2, 2> newDirFlip = flipMat * displacementField->GetDirection();
+
+  displacementFieldFlip->SetDirection(newDirFlip);
+
+  const float detFlipped = GetDeterminantAtPoint(displacementFieldFlip, physPt);
+
+  std::cout << "Determinant at point " << physPt << ":" << std::endl;
+  std::cout << "  Original: " << detOriginal << std::endl;
+  std::cout << "  Permuted: " << detPermuted << std::endl;
+  std::cout << "  Flipped:  " << detFlipped << std::endl;
+
+  constexpr double delta = 1e-13;
+
+  if (itk::Math::abs(detPermuted - detOriginal) > delta)
+  {
+    std::cerr << "Test failed: determinant differs after Permute." << std::endl;
+    testPassed = false;
+  }
+
+  if (itk::Math::abs(detFlipped - detOriginal) > delta)
+  {
+    std::cerr << "Test failed: determinant differs after Flip." << std::endl;
+    testPassed = false;
+  }
+
+  if (testPassed)
+  {
+    std::cout << "Test passed: determinant consistent after Permute and Flip." << std::endl;
   }
 
   return testPassed;