Point-based registration allows us to help the registration via pre-defined sets of corresponding points. The 'CorrespondingPointsEuclideanDistanceMetric' minimises the distance of between a points on the fixed image and corresponding points on the moving image. The metric can be used to help in a difficult registration task by taking into account positions are known to correspond. Think of it as a way of embedding expert knowledge in the registration procedure. We can manually or automatically select points via centroids of segmentations for example. Anything is possible. Point sets should always be given to elastix with their corresponding image.
import itk
import numpy as np
from itkwidgets import compare, checkerboard, view
In order for 3D registration to work with a point set, the 'CorrespondingPointsEuclideanDistanceMetric', should be set as metric. For the 3D case, this means that the metric should be a multimetric with the first metric of type AdvancedImageToImageMetric and the second the 'CorrespondingPointsEuclideanDistanceMetric'. The Registration parameter should therefore be set to 'MultiMetricMultiResolutionRegistration', to allow a multimetric parameter.
# Import Images
fixed_image = itk.imread('data/CT_3D_lung_fixed.mha', itk.F)
moving_image = itk.imread('data/CT_3D_lung_moving.mha', itk.F)
fixed_point_set = np.loadtxt('data/CT_3D_lung_fixed_point_set_corrected.txt', skiprows=2, delimiter=' ')
moving_point_set = np.loadtxt('data/CT_3D_lung_moving_point_set_corrected.txt', skiprows=2, delimiter=' ')
# Import and adjust Parameter Map
parameter_object = itk.ParameterObject.New()
parameter_map_rigid = parameter_object.GetDefaultParameterMap('rigid')
parameter_map_rigid['Registration'] = [
'MultiMetricMultiResolutionRegistration']
original_metric = parameter_map_rigid['Metric']
parameter_map_rigid['Metric'] = [original_metric[0],
'CorrespondingPointsEuclideanDistanceMetric']
parameter_object.AddParameterMap(parameter_map_rigid)
The point sets can be visualized and analysed with the itkwidgets package. Currently the point set visualization layer is put on top of the 3D image, which makes the slice view show all 'passed' points as well, this will change in future versions of the itkwidgets. (Changing the point set size and the gradient opacity of the image aids the visualization of the point sets in the 3D image)
view(fixed_image, point_sets=[fixed_point_set])
Viewer(geometries=[], gradient_opacity=0.22, point_set_colors=array([[0.8392157, 0. , 0. ]], dtype…
view(moving_image, point_sets=[moving_point_set])
Viewer(geometries=[], gradient_opacity=0.22, point_set_colors=array([[0.8392157, 0. , 0. ]], dtype…
With regards to the elastix function, the point sets do not need to be initialized first, so their file name + path can be given directly to elastix. Future version of ITKElastix will support initialization of point sets and passing these to elastix, just like the point set was initialized for the visualization with the itkwidgets.
Registration can either be done in one line with the registration function...
# Call registration function
result_image, result_transform_parameters = itk.elastix_registration_method(
fixed_image, moving_image,
fixed_point_set_file_name='data/CT_3D_lung_fixed_point_set_corrected.txt',
moving_point_set_file_name='data/CT_3D_lung_moving_point_set_corrected.txt',
log_to_console=False,
parameter_object=parameter_object)
.. 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.
elastix_object = itk.ElastixRegistrationMethod.New(fixed_image,moving_image)
elastix_object.SetFixedPointSetFileName('data/CT_3D_lung_fixed_point_set_corrected.txt')
elastix_object.SetMovingPointSetFileName(
'data/CT_3D_lung_moving_point_set_corrected.txt')
elastix_object.SetParameterObject(parameter_object)
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()