Jacobian determinant

Beginner Questions
1

Dear all,
I’m trying to compute the Jacobian Determinant of a deformation field.
The transformation and its deformation field were not computed with simple ITK. I simply applied sitk.DisplacementFieldJacobianDeterminant to the Image corresponding to the deformation field. However, the output contains both negative and positive values, which I did not expect. Is it because the output corresponds to the log of the determinant? I found no mention of it in the documentation. Or am I misusing the function?

Code:

#deformed_points : 3D array describing the deformation field.
sitk_displacement_field = sitk.GetImageFromArray(deformed_points, isVector=True)
jacobian_det_volume = sitk.DisplacementFieldJacobianDeterminant(sitk_displacement_field)

Thanks in advance!

2

Hello @Fleur,

Your expectation is correct, output is the magnitude of the Jacobian’s determinant. Likely what you are seeing is your usage of unit sized voxels, as this is the default voxel size when creating an image from a numpy array. Before computing the Jacobian determinant, you should add a line sitk_displacement_field.SetSpacing(correct_voxel_spacing_here).

To illustrate this, run the code below as is and then set the use_image_spacing to False (equivalent to unit image spacing):

import SimpleITK as sitk
import numpy as np

image = sitk.Image((128,128), sitk.sitkUInt8)
image.SetSpacing((0.5, 3))

# tx corresponds to a volume contraction
scales = (0.4, 0.2)
tx = sitk.ScaleTransform(2, scales)
tx.SetCenter(image.TransformContinuousIndexToPhysicalPoint([(i-1)/2.0 for i in image.GetSize()]))

# Convert the transformation into deformation fields (tx results in volume contraction, inverse in expansion)
f = sitk.TransformToDisplacementFieldFilter()
f.SetReferenceImage(image)

use_image_spacing = True
jac_det_contracting_field = sitk.DisplacementFieldJacobianDeterminant(f.Execute(tx), use_image_spacing)
jac_det_expanding_field = sitk.DisplacementFieldJacobianDeterminant(f.Execute(tx.GetInverse()), use_image_spacing)
jac_det_identity_field = sitk.DisplacementFieldJacobianDeterminant(f.Execute(sitk.Euler2DTransform()), use_image_spacing)

contracting_arr = sitk.GetArrayViewFromImage(jac_det_contracting_field)
print(f'Volume contracting: min {np.min(contracting_arr)}, max {np.max(contracting_arr)}.') 
expanding_arr = sitk.GetArrayViewFromImage(jac_det_expanding_field)
print(f'Volume expanding: min {np.min(expanding_arr)}, max {np.max(expanding_arr)}.') 
identity_arr = sitk.GetArrayViewFromImage(jac_det_identity_field)
print(f'Volume no change: min {np.min(identity_arr)}, max {np.max(identity_arr)}.')
1 Like
3

Thank you for your help!
I understand what you are saying but I don’t think the issue comes from incorrect voxel spacing. My original images are made of unit sized voxels…

I actually think that my displacement field is wrong. As I did not compute the transformation using sitk, I could not compute the displacement field with
sitk.TransformToDisplacementFieldFilter().

In your example, when looking at the array corresponding to the f.Execute(tx) Image, values of the displacement field are comprise approximately between -30 and 30.
In my displacement field array, values are comprised between -0.5 and 120…
I’m actually wondering if there is a cleaner way to generate a displacement field from a transformation that did not come from sitk: for example, perhaps is it possible to convert this transformation (which is a diffeomorphism) to a simpleITK Transform? Then it would allow me to apply sitk.TransformToDisplacementFieldFilter() to the transformation.

4

Hello @Fleur,

The way to create the displacement field image is exactly what you did (the DisplacementFieldJacobianDeterminantFilter expects an image not a transform). What I was doing was a “trick”, creating a simple transformation with the desired contracting/expanding volume changes and then obtaining the displacement field from that.

Not sure what exactly is going on with your displacement field given that it is originally with unit spacing in all dimensions. The filter should give you positive determinant values.

1 Like
5

I think I know what I did wrong.
The 3D array I used to compute the Jacobian described the new position of the points after the transformation, not their displacement from their old position to the new one.
Thanks for your help zivy!

1 Like