BUG: Jacobians not in physical space

The gradient is calculated in a neighborhood defined in index space.

DisplacementFieldTransform accounts for this by multiplying rows of the
gradient by the direction matrix. DisplacementFieldJacobianDeterminantFilter did not,
hence its determinants were different for the same physical transformation under different
voxel index ordering.

DisplacementFieldJacobianDeterminantFilter is now fixed to use the direction matrix.

The test is updated to check that the determinant for a given displacement field is consistent
at the same point in physical space, regardless of how the index space is oriented.

More details in issue #5358.

Author: cookpa

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

@@ -195,13 +195,13 @@ TestDisplacementJacobianDeterminantValue()
   auto dispacementfieldFlip = flip->GetOutput();
 
   // Adjust direction to compensate flip
-  itk::Matrix<double, 2, 2> newDirFlip = dispacementfield->GetDirection();
-  newDirFlip[0][0] *= -1.0;
-  newDirFlip[0][1] *= -1.0;
+  itk::Matrix<double, 2, 2> flipMat;
+  flipMat.SetIdentity();
+  flipMat[0][0] = -1.0;
+
+  itk::Matrix<double, 2, 2> newDirFlip = flipMat * dispacementfield->GetDirection();
 
   dispacementfieldFlip->SetDirection(newDirFlip);
-  dispacementfieldFlip->SetOrigin(dispacementfield->GetOrigin());
-  dispacementfieldFlip->SetSpacing(dispacementfield->GetSpacing());
 
   const float detFlipped = GetDeterminantAtPoint(dispacementfieldFlip, physPt);