BUG: Non-isotropic spacing on Jacobian computation.

Spacing was not correctly taking into account in ComputeJacobianWithRespectToPositionInternal. The contribution in each direction was wrongly assigned to jacobian column instead of row. Modified test to take into account the anisotropic spacing. The test does not pass with the provious version.
InsightSoftwareConsortium · Feb 8, 2022 · 0bb068d · 0bb068d
1 parent 46e2ab9
commit 0bb068d
Show file tree

Hide file tree

Showing 2 changed files with 7 additions and 1 deletion.
diff --git a/Modules/Filtering/DisplacementField/include/itkDisplacementFieldTransform.hxx b/Modules/Filtering/DisplacementField/include/itkDisplacementFieldTransform.hxx
@@ -303,7 +303,7 @@ DisplacementFieldTransform<TParametersValueType, NDimensions>::ComputeJacobianWi
       //      ( lpix - rpix )*space/(2.0*h); //2nd order centered difference
       for (unsigned int col = 0; col < NDimensions; ++col)
       {
-        TParametersValueType val = dPix[col] / spacing[col];
+        TParametersValueType val = dPix[col] / spacing[row];
         if (row == col)
         {
           val += 1.0;

diff --git a/Modules/Filtering/DisplacementField/test/itkDisplacementFieldTransformTest.cxx b/Modules/Filtering/DisplacementField/test/itkDisplacementFieldTransformTest.cxx
@@ -293,6 +293,12 @@ itkDisplacementFieldTransformTest(int argc, char * argv[])
   fieldJTruth(0, 1) = 0.02;
   fieldJTruth(1, 1) = 1.1;
 
+  /* Test the correctness of the Jacobian computation with anisotropic spacing */
+  FieldType::SpacingType anisotropicSpacing;
+  anisotropicSpacing[0] = 1.0;
+  anisotropicSpacing[1] = 0.8;
+  field->SetSpacing(anisotropicSpacing);
+
   itk::ImageRegionIteratorWithIndex<FieldType> it(field, field->GetLargestPossibleRegion());
   it.GoToBegin();