Groupwise registration methods try to mitigate uncertainties associated with any one image by simultaneously registering all images in a population. This incorporates all image information in registration process and eliminates bias towards a chosen reference frame. The method described here uses a 3D (2D+time) and 4D (3D+time) free-form B-spline deformation model and a similarity metric that minimizes variance of intensities under the constraint that the average deformation over images is zero. This constraint defines a true mean frame of reference that lie in the center of the population without having to calculate it explicitly.
The method can take into account temporal smoothness of the deformations and a cyclic transform in the time dimension. This may be appropriate if it is known a priori that the anatomical motion has a cyclic nature e.g. in cases of cardiac or respiratory motion.
# First import is currently necessary to run ITKElastix on MacOs
from itk import itkElastixRegistrationMethodPython
import itk
import numpy as np
import matplotlib.pyplot as plt
For this example image generators are used to generate 2D and 3D vector of images for registration.
# Image generator functions
def image_generator_2D(x1, x2, y1, y2):
image = np.zeros([100, 100], np.float32)
image[y1:y2, x1:x2] = 1
return image
def image_generator_3D(x1, x2, y1, y2, z1, z2):
image = np.zeros([100, 100, 100], np.float32)
image[z1:z2, y1:y2, x1:x2] = 1
return image
# Create a vector of images for a 2D+time example in numpy
vector_of_images = np.zeros([6, 100, 100], np.float32)
i = 0
for x in range(0, 30, 5):
image = image_generator_2D(x, x+50, x, x+50)
vector_of_images[i] = image
i += 1
vector_itk = itk.image_view_from_array(vector_of_images)
# Create a vector of images for a 3D+time example in numpy
vector_of_images = np.zeros([6, 100, 100, 100], np.float32)
i = 0
for x in range(0, 30, 5):
image = image_generator_3D(x, x+50, x, x+50, x, x+50)
vector_of_images[i] = image
i += 1
####### ----- Error for 3D+time case, itk package not build for 4D images,
####### ----- options to rebuild with different flag possible----- ######
# vector_itk = itk.image_view_from_array(vector_of_images)
# Or Load real imagesS
vector_of_images = []
for i in range(10):
image = itk.imread('data/00/00-slice00%s.dcm'%i, itk.F)
vector_of_images.append(image)
# Create Groupwise Parameter Object
parameter_object = itk.ParameterObject.New()
groupwise_parameter_map = parameter_object.GetDefaultParameterMap('groupwise')
parameter_object.AddParameterMap(groupwise_parameter_map)
Registration can either be done in one line with the registration function...
# Call registration function
# both fixed and moving image should be set with the vector_itk to prevent elastix from throwing errors
####### ----- Error for 3D+time case, itk package not build for 4D images,
####### ----- options to rebuild with different flag possible----- ######
# result_image, result_transform_parameters = itk.elastix_registration_method(
# vector_of_images, vector_of_images,
# parameter_object=parameter_object,
# log_to_console=True)
# Call registration function
# both fixed and moving image should be set with the vector_itk to prevent elastix from throwing errors
####### ----- Error for 3D+time case----- ######
result_image, result_transform_parameters = itk.elastix_registration_method(
vector_itk, vector_itk,
parameter_object=parameter_object,
log_to_console=True)
.. or by initiating an elastix image filter object.
# Load Elastix Image Filter Object
# Fixed and moving image should be given to the Elastix method to ensure that
# the correct 3D class is initialized.
# Both fixed and moving image should be set with the vector_itk to prevent elastix from throwing errors
####### ----- Error for 3D+time case, 2D+time case shown here----- ######
elastix_object = itk.ElastixRegistrationMethod.New(vector_itk, vector_itk)
elastix_object.SetParameterObject(parameter_object)
# Set additional options
elastix_object.SetLogToConsole(False)
# Update filter object (required)
elastix_object.UpdateLargestPossibleRegion()
# Results of Registration
result_image = elastix_object.GetOutput()
result_transform_parameters = elastix_object.GetTransformParameterObject()
The results of the groupwise image registration can be visualized in viewers such as the Napari viewer. For Napari the images have to be converted to n-dimensional numpy arrays.
# import napari
# Cast result images to numpy arrays for Napari Viewer
# result_image_np = np.asarray(result_image).astype(np.float32)
# View output with Napari Viewer
# with napari.gui_qt():
# viewer = napari.Viewer()
# viewer.add_image(result_image_np)
# viewer.add_image(vector_of_images_2D[0])
# viewer.add_image(vector_of_images_2D[1])
# viewer.add_image(vector_of_images_2D[2])
# viewer.add_image(vector_of_images_2D[3])
# viewer.add_image(vector_of_images_2D[4])
# viewer.add_image(vector_of_images_2D[5])