From 0f64fd5d394519fdfdd42ff7fbc6261eaba24bca Mon Sep 17 00:00:00 2001 From: Philip Cook Date: Sun, 8 Jun 2025 22:23:54 -0400 Subject: [PATCH 1/3] ENH: Test showing Jacobian differs under different voxel ordering --- ...mentFieldJacobianDeterminantFilterTest.cxx | 110 +++++++++++++++++- 1 file changed, 109 insertions(+), 1 deletion(-) diff --git a/Modules/Filtering/DisplacementField/test/itkDisplacementFieldJacobianDeterminantFilterTest.cxx b/Modules/Filtering/DisplacementField/test/itkDisplacementFieldJacobianDeterminantFilterTest.cxx index 7c86fd254fe..914d55a622d 100644 --- a/Modules/Filtering/DisplacementField/test/itkDisplacementFieldJacobianDeterminantFilterTest.cxx +++ b/Modules/Filtering/DisplacementField/test/itkDisplacementFieldJacobianDeterminantFilterTest.cxx @@ -18,7 +18,9 @@ #include #include "itkDisplacementFieldJacobianDeterminantFilter.h" +#include "itkFlipImageFilter.h" #include "itkNullImageToImageFilterDriver.hxx" +#include "itkPermuteAxesImageFilter.h" #include "itkStdStreamStateSave.h" #include "itkTestingMacros.h" @@ -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 physPt; + dispacementfield->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 & 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(dispacementfield, 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(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; + auto permute = PermuteFilterType::New(); + + PermuteFilterType::PermuteOrderArrayType order; + order[0] = 1; + order[1] = 0; // swap X <-> Y + permute->SetOrder(order); + permute->SetInput(dispacementfield); + permute->Update(); + + auto dispacementfieldPermuted = permute->GetOutput(); + + // Adjust direction to match permutation + itk::Matrix origDir = dispacementfield->GetDirection(); + itk::Matrix 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]; + + dispacementfieldPermuted->SetDirection(newDirPermute); + dispacementfieldPermuted->SetOrigin(dispacementfield->GetOrigin()); + dispacementfieldPermuted->SetSpacing(dispacementfield->GetSpacing()); + + const float detPermuted = GetDeterminantAtPoint(dispacementfieldPermuted, physPt); + + using FlipFilterType = itk::FlipImageFilter; + auto flip = FlipFilterType::New(); + + FlipFilterType::FlipAxesArrayType flipAxes; + flipAxes[0] = true; // Flip X + flipAxes[1] = false; // Keep Y + flip->SetFlipAxes(flipAxes); + flip->SetInput(dispacementfield); + flip->FlipAboutOriginOff(); + flip->Update(); + + auto dispacementfieldFlip = flip->GetOutput(); + + // Adjust direction to compensate flip + itk::Matrix newDirFlip = dispacementfield->GetDirection(); + newDirFlip[0][0] *= -1.0; + newDirFlip[0][1] *= -1.0; + + dispacementfieldFlip->SetDirection(newDirFlip); + dispacementfieldFlip->SetOrigin(dispacementfield->GetOrigin()); + dispacementfieldFlip->SetSpacing(dispacementfield->GetSpacing()); + + const float detFlipped = GetDeterminantAtPoint(dispacementfieldFlip, 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; From 64fb5730f593e5b8894c40df18651a87a5ced0c0 Mon Sep 17 00:00:00 2001 From: Philip Cook Date: Mon, 9 Jun 2025 16:45:28 -0400 Subject: [PATCH 2/3] 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. --- ...splacementFieldJacobianDeterminantFilter.h | 6 +++ ...lacementFieldJacobianDeterminantFilter.hxx | 48 +++++++++++++++---- ...mentFieldJacobianDeterminantFilterTest.cxx | 10 ++-- 3 files changed, 51 insertions(+), 13 deletions(-) diff --git a/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.h b/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.h index a02e6363680..e1d0cabb009 100644 --- a/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.h +++ b/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.h @@ -151,6 +151,10 @@ class ITK_TEMPLATE_EXPORT DisplacementFieldJacobianDeterminantFilter using RealVectorType = Vector; using RealVectorImageType = Image; + /** 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; + /** 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; @@ -271,6 +275,8 @@ class ITK_TEMPLATE_EXPORT DisplacementFieldJacobianDeterminantFilter typename ImageBaseType::ConstPointer m_RealValuedInputImage{}; RadiusType m_NeighborhoodRadius{}; + + InternalMatrixType m_InputDirection{}; }; } // end namespace itk diff --git a/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.hxx b/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.hxx index fce8877f2f7..e9aed8d88a3 100644 --- a/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.hxx +++ b/Modules/Filtering/DisplacementField/include/itkDisplacementFieldJacobianDeterminantFilter.hxx @@ -138,6 +138,16 @@ DisplacementFieldJacobianDeterminantFilter { 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(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::EvaluateAtNeighborhood( const ConstNeighborhoodIteratorType & it) const { - vnl_matrix_fixed 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 localComponentGrad_d; + for (unsigned int col = 0; col < ImageDimension; ++col) + { + localComponentGrad_d[col] = localGrad[row][col]; + } + + vnl_vector_fixed 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::OneValue(); } - return vnl_det(J); + + // Return determinant of physical Jacobian: + return vnl_determinant(physicalGrad); } template diff --git a/Modules/Filtering/DisplacementField/test/itkDisplacementFieldJacobianDeterminantFilterTest.cxx b/Modules/Filtering/DisplacementField/test/itkDisplacementFieldJacobianDeterminantFilterTest.cxx index 914d55a622d..7146cfb2e78 100644 --- a/Modules/Filtering/DisplacementField/test/itkDisplacementFieldJacobianDeterminantFilterTest.cxx +++ b/Modules/Filtering/DisplacementField/test/itkDisplacementFieldJacobianDeterminantFilterTest.cxx @@ -195,13 +195,13 @@ TestDisplacementJacobianDeterminantValue() auto dispacementfieldFlip = flip->GetOutput(); // Adjust direction to compensate flip - itk::Matrix newDirFlip = dispacementfield->GetDirection(); - newDirFlip[0][0] *= -1.0; - newDirFlip[0][1] *= -1.0; + itk::Matrix flipMat; + flipMat.SetIdentity(); + flipMat[0][0] = -1.0; + + itk::Matrix newDirFlip = flipMat * dispacementfield->GetDirection(); dispacementfieldFlip->SetDirection(newDirFlip); - dispacementfieldFlip->SetOrigin(dispacementfield->GetOrigin()); - dispacementfieldFlip->SetSpacing(dispacementfield->GetSpacing()); const float detFlipped = GetDeterminantAtPoint(dispacementfieldFlip, physPt); From b084fc00e3db58bc8b845a0ecae359fc4b642134 Mon Sep 17 00:00:00 2001 From: Philip Cook Date: Mon, 16 Jun 2025 11:40:06 -0400 Subject: [PATCH 3/3] STYLE: Fix variable name typo --- ...mentFieldJacobianDeterminantFilterTest.cxx | 42 +++++++++---------- 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/Modules/Filtering/DisplacementField/test/itkDisplacementFieldJacobianDeterminantFilterTest.cxx b/Modules/Filtering/DisplacementField/test/itkDisplacementFieldJacobianDeterminantFilterTest.cxx index 7146cfb2e78..5ddc6719f38 100644 --- a/Modules/Filtering/DisplacementField/test/itkDisplacementFieldJacobianDeterminantFilterTest.cxx +++ b/Modules/Filtering/DisplacementField/test/itkDisplacementFieldJacobianDeterminantFilterTest.cxx @@ -38,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; - 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]; @@ -96,7 +96,7 @@ TestDisplacementJacobianDeterminantValue() #endif ITK_TEST_SET_GET_BOOLEAN(filter, UseImageSpacing, useImageSpacing); - filter->SetInput(dispacementfield); + filter->SetInput(displacementField); filter->Update(); @@ -127,7 +127,7 @@ TestDisplacementJacobianDeterminantValue() // Define physical point to test itk::Point physPt; - dispacementfield->TransformIndexToPhysicalPoint(index, 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 & pt) -> float { @@ -141,7 +141,7 @@ TestDisplacementJacobianDeterminantValue() return filter->GetOutput()->GetPixel(mappedIdx); }; - const float detOriginal = GetDeterminantAtPoint(dispacementfield, physPt); + 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) @@ -162,24 +162,24 @@ TestDisplacementJacobianDeterminantValue() order[0] = 1; order[1] = 0; // swap X <-> Y permute->SetOrder(order); - permute->SetInput(dispacementfield); + permute->SetInput(displacementField); permute->Update(); - auto dispacementfieldPermuted = permute->GetOutput(); + auto displacementFieldPermuted = permute->GetOutput(); // Adjust direction to match permutation - itk::Matrix origDir = dispacementfield->GetDirection(); + itk::Matrix origDir = displacementField->GetDirection(); itk::Matrix 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]; - dispacementfieldPermuted->SetDirection(newDirPermute); - dispacementfieldPermuted->SetOrigin(dispacementfield->GetOrigin()); - dispacementfieldPermuted->SetSpacing(dispacementfield->GetSpacing()); + displacementFieldPermuted->SetDirection(newDirPermute); + displacementFieldPermuted->SetOrigin(displacementField->GetOrigin()); + displacementFieldPermuted->SetSpacing(displacementField->GetSpacing()); - const float detPermuted = GetDeterminantAtPoint(dispacementfieldPermuted, physPt); + const float detPermuted = GetDeterminantAtPoint(displacementFieldPermuted, physPt); using FlipFilterType = itk::FlipImageFilter; auto flip = FlipFilterType::New(); @@ -188,22 +188,22 @@ TestDisplacementJacobianDeterminantValue() flipAxes[0] = true; // Flip X flipAxes[1] = false; // Keep Y flip->SetFlipAxes(flipAxes); - flip->SetInput(dispacementfield); + flip->SetInput(displacementField); flip->FlipAboutOriginOff(); flip->Update(); - auto dispacementfieldFlip = flip->GetOutput(); + auto displacementFieldFlip = flip->GetOutput(); // Adjust direction to compensate flip itk::Matrix flipMat; flipMat.SetIdentity(); flipMat[0][0] = -1.0; - itk::Matrix newDirFlip = flipMat * dispacementfield->GetDirection(); + itk::Matrix newDirFlip = flipMat * displacementField->GetDirection(); - dispacementfieldFlip->SetDirection(newDirFlip); + displacementFieldFlip->SetDirection(newDirFlip); - const float detFlipped = GetDeterminantAtPoint(dispacementfieldFlip, physPt); + const float detFlipped = GetDeterminantAtPoint(displacementFieldFlip, physPt); std::cout << "Determinant at point " << physPt << ":" << std::endl; std::cout << " Original: " << detOriginal << std::endl;