With the transformix algorithm the spatial jacobian and the determinant of the spatial jacobian of the transformation can be calculated. Especially the determinant of the spatial Jacobian, which identifies the amount of local compression or expansion and can be useful, for example in lung ventilation studies.
import itk
import numpy as np
# Import Images
fixed_image = itk.imread('data/CT_2D_head_fixed.mha', itk.F)
moving_image = itk.imread('data/CT_2D_head_moving.mha', itk.F)
# Import Default Parameter Map
parameter_object = itk.ParameterObject.New()
parameter_map_rigid = parameter_object.GetDefaultParameterMap('rigid')
parameter_object.AddParameterMap(parameter_map_rigid)
Registration with the registration function. The output directory has to be specified, otherwise elastix will not save the transformparameter file as .txt file.
# Call registration function and specify output directory
result_image, result_transform_parameters = itk.elastix_registration_method(
fixed_image, moving_image,
parameter_object=parameter_object,
output_directory='exampleoutput/')
# Import Image to transform, transformix is transforming from moving -> fixed;
# for this example the exact same moving image is used, this however is normally not
# very usefull since the elastix algorithm already transformed this image.
moving_image_transformix = itk.imread('data/CT_2D_head_moving.mha', itk.F)
The calculation of the Jacobian matrix and it's determinant can be done in one line...
# Calculate Jacobian matrix and it's determinant in a tuple
jacobians = itk.transformix_jacobian(moving_image_transformix, result_transform_parameters)
# Casting tuple to two numpy matrices for further calculations.
spatial_jacobian = np.asarray(jacobians[0]).astype(np.float32)
det_spatial_jacobian = np.asarray(jacobians[1]).astype(np.float32)
.. or by initiating an transformix image filter object, calculating the jacobian and reading the jacobian from Disk IO.
# Load Transformix Object
transformix_object = itk.TransformixFilter.New(moving_image_transformix)
transformix_object.SetTransformParameterObject(result_transform_parameters)
# Set advanced options
transformix_object.SetComputeSpatialJacobian(True)
transformix_object.SetComputeDeterminantOfSpatialJacobian(True)
# Set output directory for spatial jacobian and its determinant,
# default directory is current directory.
transformix_object.SetOutputDirectory('exampleoutput/')
# Update object (required)
transformix_object.UpdateLargestPossibleRegion()
# Results of Transformation
result_image_transformix = transformix_object.GetOutput()
# Load Jacobian from Disk IO and cast to numpy matrices
spatial_jacobian = itk.imread('exampleoutput/fullSpatialJacobian.nii', itk.F)
spatial_jacobian = np.asarray(spatial_jacobian).astype(np.float32)
det_spatial_jacobian = itk.imread('exampleoutput/spatialJacobian.nii', itk.F)
det_spatial_jacobian = np.asarray(det_spatial_jacobian).astype(np.float32)
Inspect the deformation field by looking at the determinant of the Jacobian of Tµ(x). Values smaller than 1 indicate local compression, values larger than 1 indicate local expansion, and 1 means volume preservation. The measure is quantitative: a value of 1.1 means a 10% increase in volume. If this value deviates substantially from 1, you may be worried (but maybe not if this is what you expect for your application). In case it is negative you have “foldings” in your transformation, and you definitely should be worried. For more information see elastix manual.
print("Number of foldings in transformation:",np.sum(det_spatial_jacobian < 0))
Number of foldings in transformation: 0