Регистрация изображения с помощью SimpleITK: как получить преобразованное изображение ⇐ Python

[Anonymous]

Недавно я научился регистрировать изображения с помощью SimpleITK (sitk).
Моя задача — вычесть два КТ-изображения, полученные в два разных момента времени.
Прежде всего мне нужно зарегистрировать эти изображения. Я сделал это, используя коды из этого блокнота. Ниже показано, что я реализовал оттуда (я слепо использовал эти коды. Я не совсем разбираюсь в документации по ситку).
Из приведенного ниже кода я понимаю, что Final_transform содержит преобразование, необходимое для регистрации этих изображения. Может ли кто-нибудь объяснить:

• что возвращает этот Final_transform?
• как я могу использовать Final_transform чтобы я мог выполнять попиксельные операции между двумя изображениями КТ, т. е. вычитание изображений?

Код: Выделить всё

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)

Я не знаю, как действовать дальше с Final_transform. Как я могу использовать его для вычитания двух изображений КТ?

Подробнее здесь: https://stackoverflow.com/questions/789 ... rmed-image