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 7c86fd254fe..5ddc6719f38 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" @@ -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; - 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 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 { + 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(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(displacementField); + permute->Update(); + + auto displacementFieldPermuted = permute->GetOutput(); + + // Adjust direction to match permutation + 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]; + + displacementFieldPermuted->SetDirection(newDirPermute); + displacementFieldPermuted->SetOrigin(displacementField->GetOrigin()); + displacementFieldPermuted->SetSpacing(displacementField->GetSpacing()); + + const float detPermuted = GetDeterminantAtPoint(displacementFieldPermuted, 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(displacementField); + flip->FlipAboutOriginOff(); + flip->Update(); + + auto displacementFieldFlip = flip->GetOutput(); + + // Adjust direction to compensate flip + itk::Matrix flipMat; + flipMat.SetIdentity(); + flipMat[0][0] = -1.0; + + itk::Matrix 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;