文章目录
- Mutual Information Metric 互信息度量
- Retrieve fixed and moving images for registration 检索固定和移动图像进行注册
- Prepare images for registration 准备注册图像
- Plot the MutualInformationImageToImageMetric surface 绘制 MutualInformationImageToImageMetric 表面
- Follow gradient ascent 按照梯度上升
- Resample the moving image 重新采样运动图像
- Clean up 清理图像
- 完整代码
- itkExhaustiveOptimizer说明
原为url: https://examples.itk.org/src/core/transform/mutualinformationaffine/mutualinformationaffine?highlight=mattesmutualinformationimagetoimagemetricv4
Mutual Information Metric 互信息度量
MutualInformationImageToImageMetric 类计算两幅图像之间的互信息，即一幅图像中的信息内容对另一幅图像的依赖程度。此示例展示了如何使用 MutualInformationImageToImageMetric 映射仿射变换参数并使用梯度上升算法配准两幅图像。
import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from urllib.request import urlretrieve
import itk
from itkwidgets import view
Retrieve fixed and moving images for registration 检索固定和移动图像进行注册
我们的目标是配准两幅切片图像，其中一幅具有任意偏移和旋转。我们试图使用仿射变换适当旋转和平移移动图像，使其与固定图像配准。
fixed_image_path = "fixed.png"
moving_image_path = "moving.png"
if not os.path.exists(fixed_image_path):
url = "https://data.kitware.com/api/v1/file/602c10a22fa25629b97d2896/download"
urlretrieve(url, fixed_image_path)
if not os.path.exists(moving_image_path):
url = "https://data.kitware.com/api/v1/file/602c10a32fa25629b97d28a0/download"
urlretrieve(url, moving_image_path)
fixed_image = itk.imread(fixed_image_path, itk.F)
moving_image = itk.imread(moving_image_path, itk.F)
view(fixed_image)
view(moving_image)
Prepare images for registration 准备注册图像
ImageType = type(fixed_image)
fixed_normalized_image = itk.normalize_image_filter(fixed_image)
fixed_smoothed_image = itk.discrete_gaussian_image_filter(fixed_normalized_image, variance=2.0)
moving_normalized_image = itk.normalize_image_filter(moving_image)
moving_smoothed_image = itk.discrete_gaussian_image_filter(moving_normalized_image, variance=2.0)
view(fixed_smoothed_image)
Plot the MutualInformationImageToImageMetric surface 绘制 MutualInformationImageToImageMetric 表面
对于这个相对简单的例子，我们试图使用 TranslationTransform 仅调整移动图像的 x 和 y 偏移。我们可以使用 ExhaustiveOptimizer 比较两个图像在许多不同的可能偏移对中获取 MutualInformationImageToImageMetric 值，并使用 matplotlib 将此数据集可视化为表面。
仿射变换包含六个参数，表示仿射矩阵 A 中的每个元素，这将决定如何对移动图像进行采样。我们知道移动图像已经平移，因此我们将可视化两个平移参数，但我们可以设置 X_INDEX 和 Y_INDEX 来可视化任何一对参数。有关仿射变换的更多信息，请参阅 https://en.wikipedia.org/wiki/Affine_transformation#Image_transformation。
X_INDEX = 4 # Translation in the X direction
Y_INDEX = 5 # Translation in the Y direction
# Move at most 20 pixels away from the initial position
window_size = [0] * 6
window_size[X_INDEX] = 20 # Set lower if visualizing elements 0-3
window_size[Y_INDEX] = 20 # Set lower if visualizing elements 0-3
# Collect 50 steps of data along each axis
n_steps = [0] * 6
n_steps[X_INDEX] = 50
n_steps[Y_INDEX] = 50
dim = fixed_image.GetImageDimension()
TransformType = itk.AffineTransform[itk.D, dim]
transform = TransformType.New()
InterpolatorType = itk.LinearInterpolateImageFunction[ImageType, itk.D]
interpolator = InterpolatorType.New()
MetricType = itk.MutualInformationImageToImageMetric[ImageType, ImageType]
metric = MetricType.New()
metric.SetNumberOfSpatialSamples(100)
metric.SetFixedImageStandardDeviation(5.0)
metric.SetMovingImageStandardDeviation(5.0)
metric.ReinitializeSeed(121212)
ExhaustiveOptimizerType = itk.ExhaustiveOptimizer
optimizer = ExhaustiveOptimizerType.New()
# Map out [n_steps] in each direction
optimizer.SetNumberOfSteps(n_steps)
# Move [window_size / n_steps] units with every step
scales = optimizer.GetScales()
scales.SetSize(6)
for i in range(0, 6):
scales.SetElement(i, (window_size[i] / n_steps[i]) if n_steps[i] != 0 else 1)
optimizer.SetScales(scales)
# Collect data describing the parametric surface with an observer
surface = dict()
def print_iteration():
surface[tuple(optimizer.GetCurrentPosition())] = optimizer.GetCurrentValue()
optimizer.AddObserver(itk.IterationEvent(), print_iteration)
RegistrationType = itk.ImageRegistrationMethod[ImageType, ImageType]
registrar = RegistrationType.New()
registrar.SetFixedImage(fixed_smoothed_image)
registrar.SetMovingImage(moving_smoothed_image)
registrar.SetOptimizer(optimizer)
registrar.SetTransform(transform)
registrar.SetInterpolator(interpolator)
registrar.SetMetric(metric)
registrar.SetFixedImageRegion(fixed_image.GetBufferedRegion())
registrar.SetInitialTransformParameters(transform.GetParameters())
registrar.Update()
# Check the extreme positions within the observed window
max_position = list(optimizer.GetMaximumMetricValuePosition())
min_position = list(optimizer.GetMinimumMetricValuePosition())
max_val = optimizer.GetMaximumMetricValue()
min_val = optimizer.GetMinimumMetricValue()
print(max_position)
print(min_position)
# Set up values for the plot
x_vals = [
list(set([x[i] for x in surface.keys()])) for i in range(0, transform.GetNumberOfParameters())
]
for i in range(0, transform.GetNumberOfParameters()):
x_vals[i].sort()
X, Y = np.meshgrid(x_vals[X_INDEX], x_vals[Y_INDEX])
Z = np.array([[surface[(1, 0, 0, 1, x0, x1)] for x1 in x_vals[X_INDEX]] for x0 in x_vals[Y_INDEX]])
# Plot the surface as a 2D heat map
fig = plt.figure()
# Invert the y-axis to represent the image coordinate system
plt.gca().invert_yaxis()
ax = plt.gca()
surf = ax.scatter(X, Y, c=Z, cmap=cm.coolwarm)
# Mark extremes on the plot
ax.plot(max_position[X_INDEX], max_position[Y_INDEX], "k^")
ax.plot(min_position[X_INDEX], min_position[Y_INDEX], "kv")
# Plot the surface as a 3D scatter plot
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
Follow gradient ascent 按照梯度上升
一旦我们理解了参数曲面的形状，就更容易将梯度上升算法可视化。我们看到曲面有一些粗糙度，但它有一个明显的向上坡度。我们希望最大化两幅图像之间的互信息，以优化配准。梯度上升优化的结果可以叠加到 matplotlib 图上。
transform = TransformType.New()
interpolator = InterpolatorType.New()
metric = MetricType.New()
metric.SetNumberOfSpatialSamples(100)
metric.SetFixedImageStandardDeviation(5.0)
metric.SetMovingImageStandardDeviation(5.0)
metric.ReinitializeSeed(121212)
n_iterations = 200
optimizer = itk.GradientDescentOptimizer.New()
optimizer.SetLearningRate(1.0)
optimizer.SetNumberOfIterations(n_iterations)
optimizer.MaximizeOn()
# Set scales so that the optimizer can take
# large steps along translation parameters,
# moderate steps along rotational parameters, and
# small steps along scale parameters
optimizer.SetScales([100, 0.5, 0.5, 100, 0.0001, 0.0001])
descent_data = dict()
descent_data[0] = (1, 0, 0, 1, 0, 0)
def log_iteration():
descent_data[optimizer.GetCurrentIteration() + 1] = tuple(optimizer.GetCurrentPosition())
optimizer.AddObserver(itk.IterationEvent(), log_iteration)
registrar = RegistrationType.New()
registrar.SetFixedImage(fixed_smoothed_image)
registrar.SetMovingImage(moving_smoothed_image)
registrar.SetTransform(transform)
registrar.SetInterpolator(interpolator)
registrar.SetMetric(metric)
registrar.SetOptimizer(optimizer)
registrar.SetFixedImageRegion(fixed_image.GetBufferedRegion())
registrar.SetInitialTransformParameters(transform.GetParameters())
registrar.Update()
print(f"Its: {optimizer.GetCurrentIteration()}")
print(f"Final Value: {optimizer.GetValue()}")
print(f"Final Position: {list(registrar.GetLastTransformParameters())}")
x_vals = [descent_data[i][X_INDEX] for i in range(0, n_iterations)]
y_vals = [descent_data[i][Y_INDEX] for i in range(0, n_iterations)]
我们在图中看到，随着每次迭代更新变换参数，度量通常会得到改善，但最终位置可能与图上的最大位置不一致。这是一种难以在超维空间中可视化梯度上升的情况，其中优化器正在逐步遍历六个参数维度，但我们使用 ExhaustiveOptimizer 收集的 2D 图表示空间中的“切片”，其中 x[0:4] 固定在 (1,0,0,1)。在这里，直接比较配准后的两幅图像以评估适应度可能更有用。
fig = plt.figure()
# Note: We invert the y-axis to represent the image coordinate system
plt.gca().invert_yaxis()
ax = plt.gca()
surf = ax.scatter(X, Y, c=Z, cmap=cm.coolwarm)
for i in range(0, n_iterations - 1):
plt.plot(x_vals[i : i + 2], y_vals[i : i + 2], "wx-")
plt.plot(descent_data[0][X_INDEX], descent_data[0][Y_INDEX], "bo")
plt.plot(descent_data[n_iterations - 1][X_INDEX], descent_data[n_iterations - 1][Y_INDEX], "ro")
plt.plot(max_position[X_INDEX], max_position[Y_INDEX], "k^")
plt.plot(min_position[X_INDEX], min_position[Y_INDEX], "kv")
Resample the moving image 重新采样运动图像
为了应用梯度上升的结果，我们必须将移动图像重新采样到固定图像的域中。TranslationTransform 的参数已通过梯度上升选择，用于指示如何从固定图像域采样移动图像。我们可以使用 itkwidgets 比较这两幅图像以验证配准是否成功。
ResampleFilterType = itk.ResampleImageFilter[ImageType, ImageType]
resample = ResampleFilterType.New(
Transform=transform,
Input=moving_image,
Size=fixed_image.GetLargestPossibleRegion().GetSize(),
OutputOrigin=fixed_image.GetOrigin(),
OutputSpacing=fixed_image.GetSpacing(),
OutputDirection=fixed_image.GetDirection(),
DefaultPixelValue=100,
)
resample.Update()
view(resample.GetOutput())
图像比较显示，图像已成功平移重叠，但未完全旋转以完全对齐。如果我们要进一步探索，我们可以使用具有度量的其他优化器，例如 LBFGSBOptimizer 类，它可能更适合在粗糙的参数表面上进行优化。我们还可以探索不同的度量，例如 MattesMutualInformationImageToImageMetricv4 类，以利用 ITK v4+ 注册框架，与本例中用作 v3 框架一部分的 MutualInformationImageToImageMetric 形成对比。
Clean up 清理图像
os.remove(fixed_image_path)
os.remove(moving_image_path)os.remove(fixed_image_path)
os.remove(moving_image_path)
完整代码
import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import itk
fixed_image_path = r"D:\learn\itk\ITK\Examples\Data\BrainProtonDensitySliceBorder20.png"
moving_image_path = r"D:\learn\itk\ITK\Examples\Data\BrainProtonDensitySliceR10X13Y17.png"
fixed_image = itk.imread(fixed_image_path, itk.F)
moving_image = itk.imread(moving_image_path, itk.F)
fixed_image_array = itk.array_view_from_image(fixed_image)
moving_image_array = itk.array_view_from_image(moving_image)
print('fixed image shape:', fixed_image_array.shape)
print('moving image shape:', moving_image_array.shape)
ImageType = type(fixed_image)
fixed_normalized_image = itk.normalize_image_filter(fixed_image)
fixed_smoothed_image = itk.discrete_gaussian_image_filter(fixed_normalized_image, variance=2.0)
moving_normalized_image = itk.normalize_image_filter(moving_image)
moving_smoothed_image = itk.discrete_gaussian_image_filter(moving_normalized_image, variance=2.0)
fixed_image_smoothed_array = itk.array_view_from_image(fixed_smoothed_image)
moving_image_smoothed_array = itk.array_view_from_image(moving_smoothed_image)
plt.figure()
plt.subplot(221), plt.imshow(fixed_image_array), plt.colorbar(), plt.title('fixed_image_array'), plt.axis('off')
plt.subplot(222), plt.imshow(moving_image_array, cmap='gray'), plt.colorbar(), plt.title('moving_image_array'), plt.axis('off')
plt.subplot(223), plt.imshow(fixed_image_smoothed_array), plt.colorbar(), plt.title('fixed_image_smoothed_array'), plt.axis('off')
plt.subplot(224), plt.imshow(moving_image_smoothed_array, cmap='gray'), plt.colorbar(), plt.title('moving_image_smoothed_array'), plt.axis('off')
plt.show()
X_INDEX = 4 # Translation in the X direction
Y_INDEX = 5 # Translation in the Y direction
# Move at most 20 pixels away from the initial position
window_size = [0] * 6
window_size[X_INDEX] = 20 # Set lower if visualizing elements 0-3
window_size[Y_INDEX] = 20 # Set lower if visualizing elements 0-3
# Collect 50 steps of data along each axis
n_steps = [0] * 6
n_steps[X_INDEX] = 50
n_steps[Y_INDEX] = 50
dim = fixed_image.GetImageDimension()
TransformType = itk.AffineTransform[itk.D, dim]
transform = TransformType.New()
InterpolatorType = itk.LinearInterpolateImageFunction[ImageType, itk.D]
interpolator = InterpolatorType.New()
MetricType = itk.MutualInformationImageToImageMetric[ImageType, ImageType]
metric = MetricType.New()
metric.SetNumberOfSpatialSamples(100)
metric.SetFixedImageStandardDeviation(5.0)
metric.SetMovingImageStandardDeviation(5.0)
metric.ReinitializeSeed(121212)
ExhaustiveOptimizerType = itk.ExhaustiveOptimizer
optimizer = ExhaustiveOptimizerType.New()
# Map out [n_steps] in each direction
optimizer.SetNumberOfSteps(n_steps)
# Move [window_size / n_steps] units with every step
scales = optimizer.GetScales()
scales.SetSize(6)
for i in range(0, 6):
scales.SetElement(i, (window_size[i] / n_steps[i]) if n_steps[i] != 0 else 1)
optimizer.SetScales(scales)
# Collect data describing the parametric surface with an observer
surface = dict()
def print_iteration():
surface[tuple(optimizer.GetCurrentPosition())] = optimizer.GetCurrentValue()
optimizer.AddObserver(itk.IterationEvent(), print_iteration)
RegistrationType = itk.ImageRegistrationMethod[ImageType, ImageType]
registrar = RegistrationType.New()
registrar.SetFixedImage(fixed_smoothed_image)
registrar.SetMovingImage(moving_smoothed_image)
registrar.SetOptimizer(optimizer)
registrar.SetTransform(transform)
registrar.SetInterpolator(interpolator)
registrar.SetMetric(metric)
registrar.SetFixedImageRegion(fixed_image.GetBufferedRegion())
registrar.SetInitialTransformParameters(transform.GetParameters())
registrar.Update()
# Check the extreme positions within the observed window
max_position = list(optimizer.GetMaximumMetricValuePosition())
min_position = list(optimizer.GetMinimumMetricValuePosition())
max_val = optimizer.GetMaximumMetricValue()
min_val = optimizer.GetMinimumMetricValue()
print('max_position:', max_position, max_val)
print('min_position:', min_position, min_val)
# Set up values for the plot
x_vals = [
list(set([x[i] for x in surface.keys()])) for i in range(0, transform.GetNumberOfParameters())
]
for i in range(0, transform.GetNumberOfParameters()):
x_vals[i].sort()
X, Y = np.meshgrid(x_vals[X_INDEX], x_vals[Y_INDEX])
Z = np.array([[surface[(1, 0, 0, 1, x0, x1)] for x1 in x_vals[X_INDEX]] for x0 in x_vals[Y_INDEX]])
# Plot the surface as a 2D heat map
fig = plt.figure()
# Invert the y-axis to represent the image coordinate system
plt.gca().invert_yaxis()
ax = plt.gca()
surf = ax.scatter(X, Y, c=Z, cmap=cm.coolwarm)
# Mark extremes on the plot
ax.plot(max_position[X_INDEX], max_position[Y_INDEX], "k^")
ax.plot(min_position[X_INDEX], min_position[Y_INDEX], "kv")
plt.show()
# Plot the surface as a 3D scatter plot
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()
transform = TransformType.New()
interpolator = InterpolatorType.New()
metric = MetricType.New()
metric.SetNumberOfSpatialSamples(100)
metric.SetFixedImageStandardDeviation(5.0)
metric.SetMovingImageStandardDeviation(5.0)
metric.ReinitializeSeed(121212)
n_iterations = 200
optimizer = itk.GradientDescentOptimizer.New()
optimizer.SetLearningRate(1.0)
optimizer.SetNumberOfIterations(n_iterations)
optimizer.MaximizeOn()
# Set scales so that the optimizer can take
# large steps along translation parameters,
# moderate steps along rotational parameters, and
# small steps along scale parameters
optimizer.SetScales([100, 0.5, 0.5, 100, 0.0001, 0.0001])
descent_data = dict()
descent_data[0] = (1, 0, 0, 1, 0, 0)
def log_iteration():
descent_data[optimizer.GetCurrentIteration() + 1] = tuple(optimizer.GetCurrentPosition())
optimizer.AddObserver(itk.IterationEvent(), log_iteration)
registrar = RegistrationType.New()
registrar.SetFixedImage(fixed_smoothed_image)
registrar.SetMovingImage(moving_smoothed_image)
registrar.SetTransform(transform)
registrar.SetInterpolator(interpolator)
registrar.SetMetric(metric)
registrar.SetOptimizer(optimizer)
registrar.SetFixedImageRegion(fixed_image.GetBufferedRegion())
registrar.SetInitialTransformParameters(transform.GetParameters())
registrar.Update()
print(f"Its: {optimizer.GetCurrentIteration()}")
print(f"Final Value: {optimizer.GetValue()}")
print(f"Final Position: {list(registrar.GetLastTransformParameters())}")
x_vals = [descent_data[i][X_INDEX] for i in range(0, n_iterations)]
y_vals = [descent_data[i][Y_INDEX] for i in range(0, n_iterations)]
fig = plt.figure()
# Note: We invert the y-axis to represent the image coordinate system
plt.gca().invert_yaxis()
ax = plt.gca()
surf = ax.scatter(X, Y, c=Z, cmap=cm.coolwarm)
for i in range(0, n_iterations - 1):
plt.plot(x_vals[i : i + 2], y_vals[i : i + 2], "wx-")
plt.plot(descent_data[0][X_INDEX], descent_data[0][Y_INDEX], "bo")
plt.plot(descent_data[n_iterations - 1][X_INDEX], descent_data[n_iterations - 1][Y_INDEX], "ro")
plt.plot(max_position[X_INDEX], max_position[Y_INDEX], "k^")
plt.plot(min_position[X_INDEX], min_position[Y_INDEX], "kv")
plt.show()
ResampleFilterType = itk.ResampleImageFilter[ImageType, ImageType]
resample = ResampleFilterType.New(
Transform=transform,
Input=moving_image,
Size=fixed_image.GetLargestPossibleRegion().GetSize(),
OutputOrigin=fixed_image.GetOrigin(),
OutputSpacing=fixed_image.GetSpacing(),
OutputDirection=fixed_image.GetDirection(),
DefaultPixelValue=100,
)
resample.Update()
result_image_array = itk.array_view_from_image(resample.GetOutput())
plt.figure()
plt.subplot(121), plt.imshow(fixed_image_array), plt.colorbar(), plt.title('fixed_image_array'), plt.axis('off')
plt.subplot(122), plt.imshow(result_image_array), plt.colorbar(), plt.title('result_image_array'), plt.axis('off')
plt.show()
itkExhaustiveOptimizer说明
itkExhaustiveOptimizer 是 ITK（Insight Segmentation and Registration Toolkit）中的一个优化器，用于执行网格搜索（grid search），在给定的参数空间内寻找最佳参数组合。这种优化方法对于解决一些小规模的问题非常有用，尤其是在需要全局搜索而参数空间又不是非常大的情况下。
基本概念
- 参数空间：指的是所有可能的参数值的集合。对于优化问题，我们需要在给定的参数空间内找到最优解。
- 网格搜索：是指在指定的参数范围内，按照一定的步长对每个参数进行遍历，从而找到使得目标函数（如相似性度量）达到最优的参数组合。
使用场景
itkExhaustiveOptimizer 通常用于图像配准（registration）中的参数优化。在配准过程中，需要找到一组变换参数，使得一个图像（浮动图像）与另一个参考图像之间的相似性最大化。
参数设置
当使用 itkExhaustiveOptimizer 进行优化时，需要设置以下几个关键参数：
- StepsPerAxis：每个参数轴上的步数。这定义了网格的分辨率。
- MaximumStepLength：最大步长，即在参数空间中的最大移动距离。
- MinimumStepLength：最小步长，即在参数空间中的最小移动距离。
- NumberOfSamples：样本数量，在网格搜索中，这通常是指在参数空间中考虑的点的数量。
示例代码
下面是一个简单的示例代码，展示了如何在 ITK 中使用 itkExhaustiveOptimizer 进行图像配准：
import itk
# 创建优化器实例
optimizer = itk.ExhaustiveOptimizer.New()
# 设置参数
steps_per_axis = [10, 10] # 指定每个参数轴上的步数
optimizer.SetMaximumStepLength(5.0) # 最大步长
optimizer.SetMinimumStepLength(0.01) # 最小步长
optimizer.SetNumberOfSamples(steps_per_axis) # 设置样本数量
# 创建配准器
registration_method = itk.ImageRegistrationMethodv2.New()
# 设置优化器
registration_method.SetOptimizer(optimizer)
# 设置变换类型
translation_transform = itk.TranslationTransform[itk.D, 2].New()
registration_method.SetTransform(translation_transform)
# 设置固定图像和移动图像
fixed_image = ... # 加载固定图像
moving_image = ... # 加载移动图像
registration_method.SetFixedImage(fixed_image)
registration_method.SetMovingImage(moving_image)
# 设置相似性度量
similarity_metric = itk.MeanSquaresImageToImageMetricv2.New()
registration_method.SetMetric(similarity_metric)
# 执行配准
registration_method.Update()
# 获取结果变换参数
final_parameters = registration_method.GetOutput().Get().GetParameters()
print("Final Parameters:", final_parameters)
在上面的示例中，我们创建了一个 ExhaustiveOptimizer 实例，并设置了其参数。然后定义了一个配准方法，并将优化器、变换类型、相似性度量等组件添加到配准方法中。最后，执行配准过程并输出最终的变换参数。
请注意，此示例假设你已经有了两个图像（固定图像和移动图像），并且你已经定义了合适的变换类型和相似性度量。在实际应用中，这些组件的选择和配置可能会有所不同，取决于具体的配准任务。
扫一扫
76
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
