Output of itk::DisplacementFieldJacobianDeterminantFilter depends on direction matrix

[cookpa]

Description

Jacobian determinants from itk::DisplacementFieldJacobianDeterminantFilter differ from those derived from itkDisplacementFieldTransform.hxx if the displacement field image has non-identity direction matrix.

Steps to Reproduce

Here's a link to ITK code that computes Jacobian two ways: by using the itk::DisplacementFieldJacobianDeterminantFilter, and by iterating over each voxel calling transform->ComputeJacobianWithRespectToPosition on a displacement field transform object.

https://github.com/cookpa/antsJacobianExample/tree/master/itk_jacobian

The warp fields I tested on were from ANTs, and are available here

https://upenn.box.com/s/voe2agpd30serkj4kn0fvt9bit8pg5nx

By integrating the Jacobian over each label in the fixed space, I computed the expected volume of the deformed moving image. For an identity direction matrix, the two methods are similar.

But for the case where the direction matrix is diag(-1,-1,1), the determinant filter is very different.

Expected behavior

Determinants are similar, if not exactly the same. They should both approximate the volume change within a region due to the warp.

Actual behavior

Differences of a few percent for oblique direction matrices, large differences for D = diag(-1,-1, 1).

Reproducibility

Consistent on my test data.

Versions

v5.4.3

Environment

Mac OS Intel
Python 3.11
Apple clang version 16.0.0 (clang-1600.0.26.6)

Additional Information

Code at https://github.com/cookpa/antsJacobianExample

The box link above contains data and code.

[github-actions]

Thank you for contributing an issue! 🙏

Welcome to the ITK community! 🤗👋☀️

We are glad you are here and appreciate your contribution. Please keep in mind our community participation guidelines. 📜
Also, please check existing open issues and consider discussion on the ITK Discourse. 📖

This is an automatic message. Allow for time for the ITK community to be able to read the issue and comment on it.

[cookpa]

I think I phrased this incorrectly - the dependency is only indirectly on the direction matrix, what matters is the layout of index space. If you change the voxel ordering, you get a different answer for the determinant at a particular point in physical space. This is because the filter evaluates the gradient of physical displacement vectors in an index-based neighborhood, and does not transform the result into physical space the way DisplacementFieldTransform does.

Working on a PR now.