Latest Page Update: 23-09-2026
Image RegistrationΒΆ
Start by importing some useful functions:
import numpy as np
import matplotlib.pyplot as plt
import SimpleITK as sitk
from IPython.display import clear_output
from skimage.util import img_as_ubyte
and defining some useful functions:
def imshow_orthogonal_view(sitkImage, origin = None, title=None):
"""
Display the orthogonal views of a 3D volume from the middle of the volume.
Parameters
----------
sitkImage : SimpleITK image
Image to display.
origin : array_like, optional
Origin of the orthogonal views, represented by a point [x,y,z].
If None, the middle of the volume is used.
title : str, optional
Super title of the figure.
Note:
On the axial and coronal views, patient's left is on the right
On the sagittal view, patient's anterior is on the left
"""
data = sitk.GetArrayFromImage(sitkImage)
if origin is None:
origin = np.array(data.shape) // 2
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
data = img_as_ubyte(data/np.max(data))
axes[0].imshow(data[origin[0], ::-1, ::-1], cmap='gray')
axes[0].set_title('Axial')
axes[1].imshow(data[::-1, origin[1], ::-1], cmap='gray')
axes[1].set_title('Coronal')
axes[2].imshow(data[::-1, ::-1, origin[2]], cmap='gray')
axes[2].set_title('Sagittal')
[ax.set_axis_off() for ax in axes]
if title is not None:
fig.suptitle(title, fontsize=16)
def overlay_slices(sitkImage0, sitkImage1, origin = None, title=None):
"""
Overlay the orthogonal views of a two 3D volume from the middle of the volume.
The two volumes must have the same shape. The first volume is displayed in red,
the second in green.
Parameters
----------
sitkImage0 : SimpleITK image
Image to display in red.
sitkImage1 : SimpleITK image
Image to display in green.
origin : array_like, optional
Origin of the orthogonal views, represented by a point [x,y,z].
If None, the middle of the volume is used.
title : str, optional
Super title of the figure.
Note:
On the axial and coronal views, patient's left is on the right
On the sagittal view, patient's anterior is on the left
"""
vol0 = sitk.GetArrayFromImage(sitkImage0)
vol1 = sitk.GetArrayFromImage(sitkImage1)
if vol0.shape != vol1.shape:
raise ValueError('The two volumes must have the same shape.')
if np.min(vol0) < 0 or np.min(vol1) < 0: # Remove negative values - Relevant for the noisy images
vol0[vol0 < 0] = 0
vol1[vol1 < 0] = 0
if origin is None:
origin = np.array(vol0.shape) // 2
sh = vol0.shape
R = img_as_ubyte(vol0/np.max(vol0))
G = img_as_ubyte(vol1/np.max(vol1))
vol_rgb = np.zeros(shape=(sh[0], sh[1], sh[2], 3), dtype=np.uint8)
vol_rgb[:, :, :, 0] = R
vol_rgb[:, :, :, 1] = G
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
axes[0].imshow(vol_rgb[origin[0], ::-1, ::-1, :])
axes[0].set_title('Axial')
axes[1].imshow(vol_rgb[::-1, origin[1], ::-1, :])
axes[1].set_title('Coronal')
axes[2].imshow(vol_rgb[::-1, ::-1, origin[2], :])
axes[2].set_title('Sagittal')
[ax.set_axis_off() for ax in axes]
if title is not None:
fig.suptitle(title, fontsize=16)
def composite2affine(composite_transform, result_center=None):
"""
Combine all of the composite transformation's contents to form an equivalent affine transformation.
Args:
composite_transform (SimpleITK.CompositeTransform): Input composite transform which contains only
global transformations, possibly nested.
result_center (tuple,list): The desired center parameter for the resulting affine transformation.
If None, then set to [0,...]. This can be any arbitrary value, as it is
possible to change the transform center without changing the transformation
effect.
Returns:
SimpleITK.AffineTransform: Affine transformation that has the same effect as the input composite_transform.
Source:
https://github.com/InsightSoftwareConsortium/SimpleITK-Notebooks/blob/master/Python/22_Transforms.ipynb
"""
# Flatten the copy of the composite transform, so no nested composites.
flattened_composite_transform = sitk.CompositeTransform(composite_transform)
flattened_composite_transform.FlattenTransform()
tx_dim = flattened_composite_transform.GetDimension()
A = np.eye(tx_dim)
c = np.zeros(tx_dim) if result_center is None else result_center
t = np.zeros(tx_dim)
for i in range(flattened_composite_transform.GetNumberOfTransforms() - 1, -1, -1):
curr_tx = flattened_composite_transform.GetNthTransform(i).Downcast()
# The TranslationTransform interface is different from other
# global transformations.
if curr_tx.GetTransformEnum() == sitk.sitkTranslation:
A_curr = np.eye(tx_dim)
t_curr = np.asarray(curr_tx.GetOffset())
c_curr = np.zeros(tx_dim)
else:
A_curr = np.asarray(curr_tx.GetMatrix()).reshape(tx_dim, tx_dim)
c_curr = np.asarray(curr_tx.GetCenter())
# Some global transformations do not have a translation
# (e.g. ScaleTransform, VersorTransform)
get_translation = getattr(curr_tx, "GetTranslation", None)
if get_translation is not None:
t_curr = np.asarray(get_translation())
else:
t_curr = np.zeros(tx_dim)
A = np.dot(A_curr, A)
t = np.dot(A_curr, t + c - c_curr) + t_curr + c_curr - c
return sitk.AffineTransform(A.flatten(), t, c)
# Callback invoked when the StartEvent happens, sets up our new data.
def start_plot():
global metric_values, multires_iterations
metric_values = []
multires_iterations = []
# Callback invoked when the EndEvent happens, do cleanup of data and figure.
def end_plot():
global metric_values, multires_iterations
del metric_values
del multires_iterations
# Close figure, we don't want to get a duplicate of the plot latter on.
plt.close()
# Callback invoked when the IterationEvent happens, update our data and display new figure.
def plot_values(registration_method):
global metric_values, multires_iterations
metric_values.append(registration_method.GetMetricValue())
# Clear the output area (wait=True, to reduce flickering), and plot current data
clear_output(wait=True)
# Plot the similarity metric values
plt.plot(metric_values, 'r')
plt.plot(multires_iterations, [metric_values[index] for index in multires_iterations], 'b*')
plt.xlabel('Iteration Number',fontsize=12)
plt.ylabel('Metric Value',fontsize=12)
plt.show()
# Callback invoked when the sitkMultiResolutionIterationEvent happens, update the index into the
# metric_values list.
def update_multires_iterations():
global metric_values, multires_iterations
multires_iterations.append(len(metric_values))
def command_iteration(method):
print(
f"{method.GetOptimizerIteration():3} "
+ f"= {method.GetMetricValue():10.5f} "
+ f": {method.GetOptimizerPosition()}"
)