[Python] Image registration with SimpleITK: how to get transformed image

I recently learned image registration using SimpleITK (sitk).

My task is to subtract two CT images that were acquired at two different time points.

First of all, I need to register these images. I did that using the codes from this notebook. Below is what I implemented from there (I blindly used those codes. I don't quite get sitk documentations). What I understand from the code below is that the final_transform contains the transformation needed to register these images. Can someone explain:

1. what does this final_transform return?
2. how can I use final_transform so that I can do pixel-wise operations between two CT images, i.e. image subtraction?
3. how can I set one of the images as the reference image, and then transform the other so that it aligns with the reference?

reader = sitk.ImageSeriesReader()

dir_name_1 = "Directory of 1st CT image"
dir_name_2 = "Directory of 2nd CT image"

fixed_image = sitk.ReadImage(reader.GetGDCMSeriesFileNames(dir_name_1), sitk.sitkFloat32)
moving_image = sitk.ReadImage(reader.GetGDCMSeriesFileNames(dir_name_2), sitk.sitkFloat32)

# Metric evaluate method
# Dictionary with all the orientations we will try. We omit the identity (x=0, y=0, z=0) as we always use it. This
# set of rotations is arbitrary. For a complete grid coverage we would naively have 64 entries
# (0, pi/2, pi, 1.5pi for each angle), but we know better, there are only 24 unique rotation matrices defined by
# these parameter value combinations.
all_orientations = {
"x=0, y=0, z=180": (0.0, 0.0, np.pi),
"x=0, y=180, z=0": (0.0, np.pi, 0.0),
"x=0, y=180, z=180": (0.0, np.pi, np.pi),
}

# Registration framework setup.
registration_method = sitk.ImageRegistrationMethod()
registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
registration_method.SetMetricSamplingPercentage(0.01)
registration_method.SetInterpolator(sitk.sitkLinear)

# Evaluate the similarity metric using the rotation parameter space sampling, translation remains the same for all.
initial_transform = sitk.Euler3DTransform(
sitk.CenteredTransformInitializer(
fixed_image,
moving_image,
sitk.Euler3DTransform(),
sitk.CenteredTransformInitializerFilter.GEOMETRY,
)
)
registration_method.SetInitialTransform(initial_transform, inPlace=False)
best_orientation = (0.0, 0.0, 0.0)
best_similarity_value = registration_method.MetricEvaluate(fixed_image, moving_image)

# Iterate over all other rotation parameter settings.
for key, orientation in all_orientations.items():
initial_transform.SetRotation(*orientation)
registration_method.SetInitialTransform(initial_transform)
current_similarity_value = registration_method.MetricEvaluate(
fixed_image, moving_image
)
if current_similarity_value < best_similarity_value:
best_similarity_value = current_similarity_value
best_orientation = orientation
print("best orientation is: " + str(best_orientation))


from multiprocessing.pool import ThreadPool
from functools import partial


# This function evaluates the metric value in a thread safe manner
def evaluate_metric(current_rotation, tx, f_image, m_image):
registration_method = sitk.ImageRegistrationMethod()
registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
registration_method.SetMetricSamplingPercentage(0.01)
registration_method.SetInterpolator(sitk.sitkLinear)
current_transform = sitk.Euler3DTransform(tx)
current_transform.SetRotation(*current_rotation)
registration_method.SetInitialTransform(current_transform)
res = registration_method.MetricEvaluate(f_image, m_image)
return res


p = ThreadPool(len(all_orientations) + 1)
orientations_list = [(0, 0, 0)] + list(all_orientations.values())
all_metric_values = p.map(
partial(
evaluate_metric, tx=initial_transform, f_image=fixed_image, m_image=moving_image
),
orientations_list,
)
best_orientation = orientations_list[np.argmin(all_metric_values)]
print("best orientation is: " + str(best_orientation))


initial_transform.SetRotation(*best_orientation)
final_transform, _ = multires_registration(fixed_image, moving_image, initial_transform)


I am not sure how to proceed further with final_transform. How can I use it to subtract the two CT images?

Continue reading...