simpleITK - Registration - SimpleITKv4配准非刚性配准：自由变形

最新推荐文章于 2025-06-12 10:53:54 发布

原创最新推荐文章于 2025-06-12 10:53:54 发布 · 498 阅读 ·

大模型引用 1 次

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

收录于

图像配准

代码可运行

文章目录

非刚性配准：自由变形

非刚性配准：自由变形

本笔记本说明了 SimpleITK 中基于自由变形 (FFD) 的非刚性配准算法的使用。

我们使用的数据是 4D（3D+时间）胸腹部 CT，即基于点验证的像素的呼吸胸部模型 (POPI)。该数据由一组时间 CT 体积、一组将每个 CT 分割为空气/身体/肺部的掩模以及一组跨 CT 体积的对应点组成。

POPI 模型由法国里昂的 Léon Bérard 癌症中心和 CREATIS 实验室提供。相关出版物是：

J. Vandemeulebroucke、D. Sarrut、P. Clarysse，“POPI 模型，基于点验证的像素的呼吸胸部模型”，
Proc.第十五届计算机在放射治疗中的应用国际会议 (ICCR)，加拿大多伦多，2007 年。

可以从 CREATIS 实验室此处获取 POPI 数据以及带有参考点的其他 4D CT 数据集。

import SimpleITK as sitk
import registration_utilities as ru
import registration_callbacks as rc

%matplotlib inline
import matplotlib.pyplot as plt

from ipywidgets import interact, fixed

# utility method that either downloads data from the Girder repository or
# if already downloaded returns the file name for reading from disk (cached data)
%run update_path_to_download_script
from downloaddata import fetch_data as fdata

实用程序

加载特定于 POPI 数据的实用程序、用于加载地面实况数据、显示和掩码标签的函数。

%run popi_utilities_setup.py

加载数据

将所有图像、掩码和点数据加载到相应的列表中。如果本地没有数据，则会从原始远程存储库下载。

查看图像。根据 POPI 网站上的文档，体积一对应于最终吸气（最大空气量）。

images = []
masks = []
points = []
for i in range(0, 10):
    image_file_name = f"POPI/meta/{i}0-P.mhd"
    mask_file_name = f"POPI/masks/{i}0-air-body-lungs.mhd"
    points_file_name = f"POPI/landmarks/{i}0-Landmarks.pts"
    images.append(
        sitk.ReadImage(fdata(image_file_name), sitk.sitkFloat32)
    )  # read and cast to format required for registration
    masks.append(sitk.ReadImage(fdata(mask_file_name)))
    points.append(read_POPI_points(fdata(points_file_name)))

interact(
    display_coronal_with_overlay,
    temporal_slice=(0, len(images) - 1),
    coronal_slice=(0, images[0].GetSize()[1] - 1),
    images=fixed(images),
    masks=fixed(masks),
    label=fixed(lung_label),
    window_min=fixed(-1024),
    window_max=fixed(976),
);

Fetching POPI/meta/00-P.mhd
Fetching POPI/masks/00-air-body-lungs.mhd
Fetching POPI/landmarks/00-Landmarks.pts
Fetching POPI/meta/10-P.mhd
Fetching POPI/masks/10-air-body-lungs.mhd
Fetching POPI/landmarks/10-Landmarks.pts
Fetching POPI/meta/20-P.mhd
Fetching POPI/masks/20-air-body-lungs.mhd
Fetching POPI/landmarks/20-Landmarks.pts
Fetching POPI/meta/30-P.mhd
Fetching POPI/masks/30-air-body-lungs.mhd
Fetching POPI/landmarks/30-Landmarks.pts
Fetching POPI/meta/40-P.mhd
Fetching POPI/masks/40-air-body-lungs.mhd
Fetching POPI/landmarks/40-Landmarks.pts
Fetching POPI/meta/50-P.mhd
Fetching POPI/masks/50-air-body-lungs.mhd
Fetching POPI/landmarks/50-Landmarks.pts
Fetching POPI/meta/60-P.mhd
Fetching POPI/masks/60-air-body-lungs.mhd
Fetching POPI/landmarks/60-Landmarks.pts
Fetching POPI/meta/70-P.mhd
Fetching POPI/masks/70-air-body-lungs.mhd
Fetching POPI/landmarks/70-Landmarks.pts
Fetching POPI/meta/80-P.mhd
Fetching POPI/masks/80-air-body-lungs.mhd
Fetching POPI/landmarks/80-Landmarks.pts
Fetching POPI/meta/90-P.mhd
Fetching POPI/masks/90-air-body-lungs.mhd
Fetching POPI/landmarks/90-Landmarks.pts

在这里插入图片描述

了解您的数据

虽然 POPI 网站指出图像 1 是吸气末期，并且目视检查似乎表明这是正确的，但我们可能应该查看一下肺容量，以确保我们预期的情况确实正在发生。

哪张图像是吸气末期，哪张是呼气末期？

label_shape_statistics_filter = sitk.LabelShapeStatisticsImageFilter()

for i, mask in enumerate(masks):
    label_shape_statistics_filter.Execute(mask)
    print(
        f"Lung volume in image {i} is {0.000001*label_shape_statistics_filter.GetPhysicalSize(lung_label)} liters."
    )

Lung volume in image 0 is 5.455734929689075 liters.
Lung volume in image 1 is 5.527017905172941 liters.
Lung volume in image 2 is 5.554014379499289 liters.
Lung volume in image 3 is 5.451784830895104 liters.
Lung volume in image 4 is 5.301415956621658 liters.
Lung volume in image 5 is 5.191841246039885 liters.
Lung volume in image 6 is 5.09130732168864 liters.
Lung volume in image 7 is 5.11128669632256 liters.
Lung volume in image 8 is 5.195802788867059 liters.
Lung volume in image 9 is 5.335378032491021 liters.

自由形式变形

此函数将使用 FFD 对齐固定和移动图像。如果给定掩码，将使用掩码内采样的点来评估相似度度量。如果给定固定和移动点，则在配准期间将显示相似度度量值和目标配准误差。

由于此笔记本执行模态内配准，我们使用 MeanSquares 相似度度量（计算简单且适合该任务）。

def bspline_intra_modal_registration(
    fixed_image,
    moving_image,
    fixed_image_mask=None,
    fixed_points=None,
    moving_points=None,
):
    registration_method = sitk.ImageRegistrationMethod()

    # Determine the number of BSpline control points using the physical spacing we want for the control grid.
    grid_physical_spacing = [50.0, 50.0, 50.0]  # A control point every 50mm
    image_physical_size = [
        size * spacing
        for size, spacing in zip(fixed_image.GetSize(), fixed_image.GetSpacing())
    ]
    mesh_size = [
        int(image_size / grid_spacing + 0.5)
        for image_size, grid_spacing in zip(image_physical_size, grid_physical_spacing)
    ]

    initial_transform = sitk.BSplineTransformInitializer(
        image1=fixed_image, transformDomainMeshSize=mesh_size, order=3
    )
    registration_method.SetInitialTransform(initial_transform)

    registration_method.SetMetricAsMeanSquares()
    # Settings for metric sampling, usage of a mask is optional. When given a mask the sample points will be
    # generated inside that region. Also, this implicitly speeds things up as the mask is smaller than the
    # whole image.
    registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
    registration_method.SetMetricSamplingPercentage(0.01)
    if fixed_image_mask:
        registration_method.SetMetricFixedMask(fixed_image_mask)

    # Multi-resolution framework.
    registration_method.SetShrinkFactorsPerLevel(shrinkFactors=[4, 2, 1])
    registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2, 1, 0])
    registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()

    registration_method.SetInterpolator(sitk.sitkLinear)
    registration_method.SetOptimizerAsLBFGSB(
        gradientConvergenceTolerance=1e-5, numberOfIterations=100
    )

    # If corresponding points in the fixed and moving image are given then we display the similarity metric
    # and the TRE during the registration.
    if fixed_points and moving_points:
        registration_method.AddCommand(
            sitk.sitkStartEvent, rc.metric_and_reference_start_plot
        )
        registration_method.AddCommand(
            sitk.sitkEndEvent, rc.metric_and_reference_end_plot
        )
        registration_method.AddCommand(
            sitk.sitkIterationEvent,
            lambda: rc.metric_and_reference_plot_values(
                registration_method, fixed_points, moving_points
            ),
        )

    return registration_method.Execute(fixed_image, moving_image)

执行配准

以下单元格允许您选择用于配准的图像，运行配准，然后计算统计数据，比较配准前后的目标配准误差并显示 TRE 的直方图。

要对配准进行计时，请取消注释 timeit magic。
注意：这会为单元格创建一个单独的范围。单元格内设置的变量（特别是 tx）将成为局部变量，因此它们的值在其他单元格中不可用。

# %%timeit -r1 -n1

# Select the fixed and moving images, valid entries are in [0,9].
fixed_image_index = 0
moving_image_index = 7


tx = bspline_intra_modal_registration(
    fixed_image=images[fixed_image_index],
    moving_image=images[moving_image_index],
    fixed_image_mask=(masks[fixed_image_index] == lung_label),
    fixed_points=points[fixed_image_index],
    moving_points=points[moving_image_index],
)
(
    initial_errors_mean,
    initial_errors_std,
    _,
    initial_errors_max,
    initial_errors,
) = ru.registration_errors(
    sitk.Euler3DTransform(), points[fixed_image_index], points[moving_image_index]
)
(
    final_errors_mean,
    final_errors_std,
    _,
    final_errors_max,
    final_errors,
) = ru.registration_errors(tx, points[fixed_image_index], points[moving_image_index])

plt.hist(initial_errors, bins=20, alpha=0.5, label="before registration", color="blue")
plt.hist(final_errors, bins=20, alpha=0.5, label="after registration", color="green")
plt.legend()
plt.title("TRE histogram")
print(
    f"Initial alignment errors in millimeters, mean(std): {initial_errors_mean:.2f}({initial_errors_std:.2f}), max: {initial_errors_max:.2f}"
)
print(
    f"Final alignment errors in millimeters, mean(std): {final_errors_mean:.2f}({final_errors_std:.2f}), max: {final_errors_max:.2f}"
)

在这里插入图片描述

Initial alignment errors in millimeters, mean(std): 5.07(2.67), max: 14.02
Final alignment errors in millimeters, mean(std): 1.50(0.71), max: 3.69

在这里插入图片描述

评估配准的另一种选择是使用分割。在这种情况下，我们将分割从一张图像转移到另一张图像，并从视觉和定量两个方面比较重叠。

注意：可以在分割评估笔记本中找到此处描述的方法的更详细版本。

# Transfer the segmentation via the estimated transformation. Use Nearest Neighbor interpolation to retain the labels.
transformed_labels = sitk.Resample(
    masks[moving_image_index],
    images[fixed_image_index],
    tx,
    sitk.sitkNearestNeighbor,
    0.0,
    masks[moving_image_index].GetPixelID(),
)

segmentations_before_and_after = [masks[moving_image_index], transformed_labels]
interact(
    display_coronal_with_label_maps_overlay,
    coronal_slice=(0, images[0].GetSize()[1] - 1),
    mask_index=(0, len(segmentations_before_and_after) - 1),
    image=fixed(images[fixed_image_index]),
    masks=fixed(segmentations_before_and_after),
    label=fixed(lung_label),
    window_min=fixed(-1024),
    window_max=fixed(976),
)

# Compute the Dice coefficient and Hausdorff distance between the segmentations before, and after registration.
ground_truth = masks[fixed_image_index] == lung_label
before_registration = masks[moving_image_index] == lung_label
after_registration = transformed_labels == lung_label

label_overlap_measures_filter = sitk.LabelOverlapMeasuresImageFilter()
label_overlap_measures_filter.Execute(ground_truth, before_registration)
print(
    f"Dice coefficient before registration: {label_overlap_measures_filter.GetDiceCoefficient():.2f}"
)
label_overlap_measures_filter.Execute(ground_truth, after_registration)
print(
    f"Dice coefficient after registration: {label_overlap_measures_filter.GetDiceCoefficient():.2f}"
)

hausdorff_distance_image_filter = sitk.HausdorffDistanceImageFilter()
hausdorff_distance_image_filter.Execute(ground_truth, before_registration)
print(
    f"Hausdorff distance before registration: {hausdorff_distance_image_filter.GetHausdorffDistance():.2f}"
)
hausdorff_distance_image_filter.Execute(ground_truth, after_registration)
print(
    f"Hausdorff distance after registration: {hausdorff_distance_image_filter.GetHausdorffDistance():.2f}"
)

在这里插入图片描述

Dice coefficient before registration: 0.94
Dice coefficient after registration: 0.97
Hausdorff distance before registration: 18.04
Hausdorff distance after registration: 12.51

多分辨率控制点网格

在上面的例子中，我们使用了标准图像配准框架。这意味着在所有图像分辨率下都使用相同的变换模型。对于全局变换（例如刚性、仿射……），变换参数的数量与变化的分辨率无关。对于 BSpline 变换，我们可以对频率较低的图像、图像金字塔的较高级别使用较少的控制点，随着金字塔的下降，控制点的数量会增加。使用标准框架，我们对所有金字塔级别使用相同数量的控制点。

要使用多分辨率控制点网格，我们有一个用于 BSpline 变换的特定初始化程序，SetInitialTransformAsBSpline。

以下代码解决了与上述相同的注册任务，只是使用多分辨率控制点网格。

def bspline_intra_modal_registration2(
    fixed_image,
    moving_image,
    fixed_image_mask=None,
    fixed_points=None,
    moving_points=None,
):
    registration_method = sitk.ImageRegistrationMethod()

    # Determine the number of BSpline control points using the physical spacing we
    # want for the finest resolution control grid.
    grid_physical_spacing = [50.0, 50.0, 50.0]  # A control point every 50mm
    image_physical_size = [
        size * spacing
        for size, spacing in zip(fixed_image.GetSize(), fixed_image.GetSpacing())
    ]
    mesh_size = [
        int(image_size / grid_spacing + 0.5)
        for image_size, grid_spacing in zip(image_physical_size, grid_physical_spacing)
    ]

    # The starting mesh size will be 1/4 of the original, it will be refined by
    # the multi-resolution framework.
    mesh_size = [int(sz / 4 + 0.5) for sz in mesh_size]

    initial_transform = sitk.BSplineTransformInitializer(
        image1=fixed_image, transformDomainMeshSize=mesh_size, order=3
    )
    # Instead of the standard SetInitialTransform we use the BSpline specific method which also
    # accepts the scaleFactors parameter to refine the BSpline mesh. In this case we start with
    # the given mesh_size at the highest pyramid level then we double it in the next lower level and
    # in the full resolution image we use a mesh that is four times the original size.
    registration_method.SetInitialTransformAsBSpline(
        initial_transform, inPlace=True, scaleFactors=[1, 2, 4]
    )
    registration_method.SetMetricAsMeanSquares()
    # Settings for metric sampling, usage of a mask is optional. When given a mask the sample points will be
    # generated inside that region. Also, this implicitly speeds things up as the mask is smaller than the
    # whole image.
    registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
    registration_method.SetMetricSamplingPercentage(0.01)
    if fixed_image_mask:
        registration_method.SetMetricFixedMask(fixed_image_mask)

    # Multi-resolution framework.
    registration_method.SetShrinkFactorsPerLevel(shrinkFactors=[4, 2, 1])
    registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2, 1, 0])
    registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()

    registration_method.SetInterpolator(sitk.sitkLinear)
    # Use the LBFGS2 instead of LBFGS. The latter cannot adapt to the changing control grid resolution.
    registration_method.SetOptimizerAsLBFGS2(
        solutionAccuracy=1e-2, numberOfIterations=100, deltaConvergenceTolerance=0.01
    )

    # If corresponding points in the fixed and moving image are given then we display the similarity metric
    # and the TRE during the registration.
    if fixed_points and moving_points:
        registration_method.AddCommand(
            sitk.sitkStartEvent, rc.metric_and_reference_start_plot
        )
        registration_method.AddCommand(
            sitk.sitkEndEvent, rc.metric_and_reference_end_plot
        )
        registration_method.AddCommand(
            sitk.sitkIterationEvent,
            lambda: rc.metric_and_reference_plot_values(
                registration_method, fixed_points, moving_points
            ),
        )

    return registration_method.Execute(fixed_image, moving_image)

# %%timeit -r1 -n1

# Select the fixed and moving images, valid entries are in [0,9].
fixed_image_index = 0
moving_image_index = 7


tx = bspline_intra_modal_registration2(
    fixed_image=images[fixed_image_index],
    moving_image=images[moving_image_index],
    fixed_image_mask=(masks[fixed_image_index] == lung_label),
    fixed_points=points[fixed_image_index],
    moving_points=points[moving_image_index],
)
(
    initial_errors_mean,
    initial_errors_std,
    _,
    initial_errors_max,
    initial_errors,
) = ru.registration_errors(
    sitk.Euler3DTransform(), points[fixed_image_index], points[moving_image_index]
)
(
    final_errors_mean,
    final_errors_std,
    _,
    final_errors_max,
    final_errors,
) = ru.registration_errors(tx, points[fixed_image_index], points[moving_image_index])

plt.hist(initial_errors, bins=20, alpha=0.5, label="before registration", color="blue")
plt.hist(final_errors, bins=20, alpha=0.5, label="after registration", color="green")
plt.legend()
plt.title("TRE histogram")
print(
    f"Initial alignment errors in millimeters, mean(std): {initial_errors_mean:.2f}({initial_errors_std:.2f}), max: {initial_errors_max:.2f}"
)
print(
    f"Final alignment errors in millimeters, mean(std): {final_errors_mean:.2f}({final_errors_std:.2f}), max: {final_errors_max:.2f}"
)