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?