使用 SimpleITK 进行图像配准：如何获取变换后的图像

Question

我最近学习了使用 SimpleITK (sitk) 进行图像配准。

我的任务是减去在两个不同时间点获取的两张 CT 图像。

首先，我需要注册这些图像。我使用此笔记本中的代码做到了这一点。下面是我从那里实现的（我盲目地使用了这些代码。我不太了解 satk 文档）。 我从下面的代码中了解到的是

final_transform

包含注册这些图像所需的转换。谁能解释一下：

1. 这个

   final_transform

   返回什么？
2. 如何使用

   final_transform

   以便我可以在两个CT图像之间进行逐像素操作，即图像相减？
3. 如何将其中一张图像设置为参考图像，然后变换另一张图像以使其与参考图像对齐？

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

。如何用它来减去两张 CT 图像？

Answer 1

配准方法是返回在固定图像和运动图像的两个空间之间映射的变换。 假设您想要对移动图像重新采样，使其与固定图像匹配。 为此，您可以使用 SimpleITK 的 Resample 函数。

要对移动图像重新采样，您可以执行以下操作：

resampled_moving = sitk.Resample(moving_image, fixed_image, final_transform)

这将为您提供与固定图像的方向、分辨率和间距相匹配的移动图像版本。 因此，您可以将这两个图像相减。