Register Two Point Sets 注册两个点集

文章目录

- Register Two Point Sets 注册两个点集
- Visualize Gradient Descent 可视化梯度下降
- Hyperparameter Search 超参数搜索
- JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4类说明

原文url: https://examples.itk.org/src/registration/metricsv4/registertwopointsets/registertwopointsets

Register Two Point Sets 注册两个点集

与图像配准类似，可以对 n 维“移动”点集进行重新采样，以与“固定”点集对齐。可以使用 ITK 点集度量和 ITK 优化器来注册这两个集合。

在此示例中，我们创建两个具有任意偏移量的 itk.PointSet 表示，并选择参数以将它们与 TranslationTransform 对齐。我们使用 JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4 类来量化点集之间的差异，并使用 GradientDescentOptimizerv4 类通过修改变换参数来迭代减少这种差异。我们的示例还包括使用 matplotlib 和 itkwidgets 可视化参数表面的示例代码以及用于优化梯度下降性能的示例超参数搜索。

import os
import sys
import itertools
from math import pi, sin, cos, sqrtimport matplotlib.pyplot as plt
import numpy as npimport itk
from itkwidgets import view# 构造二个点集
# Generate two circles with a small offset
def make_circles(dimension: int = 2, offset: list = None):PointSetType = itk.PointSet[itk.F, dimension]RADIUS = 100if not offset or len(offset) != dimension:offset = [2.0] * dimensionfixed_points = PointSetType.New()moving_points = PointSetType.New()fixed_points.Initialize()moving_points.Initialize()count = 0step = 0.1for count in range(0, int(2 * pi / step) + 1):theta = count * stepfixed_point = list()fixed_point.append(RADIUS * cos(theta))for dim in range(1, dimension):fixed_point.append(RADIUS * sin(theta))fixed_points.SetPoint(count, fixed_point)moving_point = [fixed_point[dim] + offset[dim] for dim in range(0, dimension)]moving_points.SetPoint(count, moving_point)return fixed_points, moving_pointsPOINT_SET_OFFSET = [15.0, 15.0]
fixed_set, moving_set = make_circles(offset=POINT_SET_OFFSET)# Visualize point sets with matplotlibfig = plt.figure()
ax = plt.axes()n_points = fixed_set.GetNumberOfPoints()
ax.scatter([fixed_set.GetPoint(i)[0] for i in range(0, n_points)],[fixed_set.GetPoint(i)[1] for i in range(0, n_points)],
)
ax.scatter([moving_set.GetPoint(i)[0] for i in range(0, n_points)],[moving_set.GetPoint(i)[1] for i in range(0, n_points)],
)# 运行梯度下降优化
# 我们将使用 JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4 量化点集偏移，并在 10 次梯度下降迭代中最小化度量值。ExhaustiveOptimizerType = itk.ExhaustiveOptimizerv4[itk.D]
dim = 2
# Define translation parameters to update iteratively
TransformType = itk.TranslationTransform[itk.D, dim]
transform = TransformType.New()
transform.SetIdentity()PointSetType = type(fixed_set)
# Define a metric to reflect the difference between point sets
PointSetMetricType = itk.JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4[PointSetType]
metric = PointSetMetricType.New(FixedPointSet=fixed_set,MovingPointSet=moving_set,MovingTransform=transform,PointSetSigma=5.0,KernelSigma=10.0,UseAnisotropicCovariances=False,CovarianceKNeighborhood=5,EvaluationKNeighborhood=10,Alpha=1.1,
)
metric.Initialize()
# Define an estimator to help determine step sizes along each transform parameter2
ShiftScalesType = itk.RegistrationParameterScalesFromPhysicalShift[PointSetMetricType]
shift_scale_estimator = ShiftScalesType.New(Metric=metric, VirtualDomainPointSet=metric.GetVirtualTransformedPointSet(), TransformForward=True
)max_iterations = 10
# Define the gradient descent optimzer
OptimizerType = itk.GradientDescentOptimizerv4Template[itk.D]
optimizer = OptimizerType.New(Metric=metric,NumberOfIterations=max_iterations,ScalesEstimator=shift_scale_estimator,MaximumStepSizeInPhysicalUnits=8.0,MinimumConvergenceValue=-1,DoEstimateLearningRateAtEachIteration=False,DoEstimateLearningRateOnce=True,ReturnBestParametersAndValue=True,
)
iteration_data = dict()# Track gradient descent iterations with observers
def print_iteration():print(f"It: {optimizer.GetCurrentIteration()}"f" metric value: {optimizer.GetCurrentMetricValue():.6f} "# f' transform position: {list(optimizer.GetCurrentPosition())}'f" learning rate: {optimizer.GetLearningRate()}")def log_iteration():iteration_data[optimizer.GetCurrentIteration() + 1] = list(optimizer.GetCurrentPosition())optimizer.AddObserver(itk.AnyEvent(), print_iteration)
optimizer.AddObserver(itk.IterationEvent(), log_iteration)# Set first value to default transform position
iteration_data[0] = list(optimizer.GetCurrentPosition())
# Run optimization and print out results
optimizer.StartOptimization()print(f"Number of iterations: {optimizer.GetCurrentIteration() - 1}")
print(f"Moving-source final value: {optimizer.GetCurrentMetricValue()}")
print(f"Moving-source final position: {list(optimizer.GetCurrentPosition())}")
print(f"Optimizer scales: {list(optimizer.GetScales())}")
print(f"Optimizer learning rate: {optimizer.GetLearningRate()}")
print(f"Stop reason: {optimizer.GetStopConditionDescription()}")# Resample Moving Point Set 重采样移动点集
moving_inverse = metric.GetMovingTransform().GetInverseTransform()
fixed_inverse = metric.GetFixedTransform().GetInverseTransform()
transformed_fixed_set = PointSetType.New()
transformed_moving_set = PointSetType.New()for n in range(0, metric.GetNumberOfComponents()):transformed_moving_point = moving_inverse.TransformPoint(moving_set.GetPoint(n))transformed_moving_set.SetPoint(n, transformed_moving_point)transformed_fixed_point = fixed_inverse.TransformPoint(fixed_set.GetPoint(n))transformed_fixed_set.SetPoint(n, transformed_fixed_point)# Compare fixed point set with resampled moving point set to see alignment
fig = plt.figure()
ax = plt.axes()n_points = fixed_set.GetNumberOfPoints()
ax.scatter([fixed_set.GetPoint(i)[0] for i in range(0, n_points)],[fixed_set.GetPoint(i)[1] for i in range(0, n_points)],
)
ax.scatter([transformed_moving_set.GetPoint(i)[0] for i in range(0, n_points)],[transformed_moving_set.GetPoint(i)[1] for i in range(0, n_points)],
)

使用 matplotlib 可视化点集:
在这里插入图片描述

将固定点集与重采样移动点集进行比较以查看对齐情况:
在这里插入图片描述

Visualize Gradient Descent 可视化梯度下降

我们可以使用 ITK ExhaustiveOptimizerv4 类来查看优化器如何沿着变换参数和度量定义的表面移动。

# Set up the new optimizer# Create a new transform and metric for analysis
transform = TransformType.New()
transform.SetIdentity()metric = PointSetMetricType.New(FixedPointSet=fixed_set,MovingPointSet=moving_set,MovingTransform=transform,PointSetSigma=5,KernelSigma=10.0,UseAnisotropicCovariances=False,CovarianceKNeighborhood=5,EvaluationKNeighborhood=10,Alpha=1.1,
)
metric.Initialize()# Create a new observer to map out the parameter surface
optimizer.RemoveAllObservers()
optimizer = ExhaustiveOptimizerType.New(Metric=metric)# Use observers to collect points on the surface
param_space = dict()def log_exhaustive_iteration():param_space[tuple(optimizer.GetCurrentPosition())] = optimizer.GetCurrentValue()optimizer.AddObserver(itk.IterationEvent(), log_exhaustive_iteration)# Collect a moderate number of steps along each transform parameter
step_count = 25
optimizer.SetNumberOfSteps([step_count, step_count])# Step a reasonable distance along each transform parameter
scales = optimizer.GetScales()
scales.SetSize(2)scale_size = 1.0
scales.SetElement(0, scale_size)
scales.SetElement(1, scale_size)optimizer.SetScales(scales)optimizer.StartOptimization()
print(f"MinimumMetricValue: {optimizer.GetMinimumMetricValue():.4f}\t"f"MaximumMetricValue: {optimizer.GetMaximumMetricValue():.4f}\n"f"MinimumMetricValuePosition: {list(optimizer.GetMinimumMetricValuePosition())}\t"f"MaximumMetricValuePosition: {list(optimizer.GetMaximumMetricValuePosition())}\n"f"StopConditionDescription: {optimizer.GetStopConditionDescription()}\t"
)# Reformat gradient descent data to overlay on the plot
descent_x_vals = [iteration_data[i][0] for i in range(0, len(iteration_data))]
descent_y_vals = [iteration_data[i][1] for i in range(0, len(iteration_data))]# Plot the surface, extrema, and gradient descent data in a matplotlib scatter plot
fig = plt.figure()
ax = plt.axes()ax.scatter([x for (x, y) in param_space.keys()],[y for (x, y) in param_space.keys()],c=list(param_space.values()),cmap="Greens",zorder=1,
)
ax.plot(optimizer.GetMinimumMetricValuePosition()[0], optimizer.GetMinimumMetricValuePosition()[1], "kv"
)
ax.plot(optimizer.GetMaximumMetricValuePosition()[0], optimizer.GetMaximumMetricValuePosition()[1], "w^"
)for i in range(0, len(iteration_data)):ax.plot(descent_x_vals[i : i + 2], descent_y_vals[i : i + 2], "rx-")
ax.plot(descent_x_vals[0], descent_y_vals[0], "ro")
ax.plot(descent_x_vals[len(iteration_data) - 1], descent_y_vals[len(iteration_data) - 1], "bo")

在 matplotlib 散点图中绘制表面、极值和梯度下降数据:
在这里插入图片描述

我们还可以使用 itkwidgets 查看表面并将其导出为图像.

x_vals = list(set(x for (x, y) in param_space.keys()))
y_vals = list(set(y for (x, y) in param_space.keys()))x_vals.sort()
y_vals.sort(reverse=True)
array = np.array([[param_space[(x, y)] for x in x_vals] for y in y_vals])image_view = itk.GetImageViewFromArray(array)view(image_view)

Hyperparameter Search 超参数搜索

为了找到具有不同变换、度量和优化器的适当结果，通常需要比较超参数变化的结果。在此示例中，需要评估 JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4.PointSetSigma 参数和 GradientDescentOptimizerv4.DoEstimateLearningRate 参数的不同值的性能。梯度下降迭代数据沿二维散点图绘制，以比较并选择产生最佳性能的超参数组合。

# Index values for gradient descent logging
FINAL_OPT_INDEX = 0
DESCENT_DATA_INDEX = 1hyper_data = dict()
final_optimizers = dict()# sigma must be sufficiently large to avoid negative entropy results
point_set_sigmas = (1.0, 2.5, 5.0, 10.0, 20.0, 50.0)# Compare performance with repeated or one-time learning rate estimation
estimate_rates = [(False, False), (False, True), (True, False)]# Run gradient descent optimization for each combination of hyperparameters
for trial_values in itertools.product(point_set_sigmas, estimate_rates):hyper_data[trial_values] = dict()(point_set_sigma, est_rate) = trial_valuesfixed_set, moving_set = make_circles(offset=POINT_SET_OFFSET)transform = TransformType.New()transform.SetIdentity()metric = PointSetMetricType.New(FixedPointSet=fixed_set,MovingPointSet=moving_set,PointSetSigma=point_set_sigma,KernelSigma=10.0,UseAnisotropicCovariances=False,CovarianceKNeighborhood=5,EvaluationKNeighborhood=10,MovingTransform=transform,Alpha=1.1,)shift_scale_estimator = ShiftScalesType.New(Metric=metric, VirtualDomainPointSet=metric.GetVirtualTransformedPointSet())metric.Initialize()optimizer = OptimizerType.New(Metric=metric,NumberOfIterations=100,MaximumStepSizeInPhysicalUnits=3.0,MinimumConvergenceValue=-1,DoEstimateLearningRateOnce=est_rate[0],DoEstimateLearningRateAtEachIteration=est_rate[1],LearningRate=1e6,  # Ignored if either est_rate argument is TrueReturnBestParametersAndValue=False,)optimizer.SetScalesEstimator(shift_scale_estimator)def log_hyper_iteration_data():hyper_data[trial_values][optimizer.GetCurrentIteration()] = round(optimizer.GetCurrentMetricValue(), 8)optimizer.AddObserver(itk.IterationEvent(), log_hyper_iteration_data)optimizer.StartOptimization()final_optimizers[trial_values] = optimizer# Print results for each set of hyperparameters
print("PS_sigma\test once/each:\tfinal index\tfinal metric")
for trial_values in itertools.product(point_set_sigmas, estimate_rates):print(f"{trial_values[0]}\t\t{trial_values[1]}:\t\t"f"{final_optimizers[trial_values].GetCurrentIteration()}\t"f"{final_optimizers[trial_values].GetCurrentMetricValue():10.8f}")

我们可以使用 matplotlib 子图和条形图来比较每组超参数的梯度下降性能和最终度量值。在这个例子中，我们看到，一次估计学习率通常会随着时间的推移带来最佳性能，而每次迭代都估计学习率可能会阻止梯度下降有效收敛。提供最佳和最一致度量值的超参数集是 PointSetSigma=5.0 和 DoEstimateLearningRateOnce=True，这是我们在本笔记本中使用的值。

# Visualize metric over gradient descent iterations as matplotlib subplots.f, axn = plt.subplots(len(point_set_sigmas), len(estimate_rates), sharex=True)for (i, j) in [(i, j) for i in range(0, len(point_set_sigmas)) for j in range(0, len(estimate_rates))]:axn[i, j].scatter(x=list(hyper_data[point_set_sigmas[i], estimate_rates[j]].keys())[1:],y=list(hyper_data[point_set_sigmas[i], estimate_rates[j]].values())[1:],)axn[i, j].set_title(f"sigma={point_set_sigmas[i]},est={estimate_rates[j]}")axn[i, j].set_ylim(-0.08, 0)plt.subplots_adjust(top=5, bottom=1, right=5)
plt.show()# Compare final metric magnitudes
fig = plt.figure()
ax = fig.gca()labels = [f'{round(sig,0)}{"T" if est0 else "F"}{"T" if est1 else "F"}'for (sig, (est0, est1)) in itertools.product(point_set_sigmas, estimate_rates)
]vals = [final_optimizers[trial_values].GetCurrentMetricValue()for trial_values in itertools.product(point_set_sigmas, estimate_rates)
]ax.bar(labels, vals)
plt.show()