参考：https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch3.html#x26-1740003.18
3.8 Multi-Stage Registration
多阶段，多分辨率的方法来运行一个多模态配准过程使用两个线性itk::TranslationTransform和itk::AffineTransform。此外，它还展示了在使用仿射变换时使用Scale Estimators来微调优化器的尺度参数。该itk:::RegistrationParameterScalesFromPhysicalShift过滤器用于自动估计的参数规模。
在多阶段场景中，每个阶段都需要itk::ImageRegistrationMethodv4的单独实例化，因此每个阶段可能有不同的transform、不同的optimizer和不同的metric度量，可以在多个级别执行。配准方法在每个阶段的配置都严格遵循前一节中的程序。
在早期阶段，我们可以使用更简单的transform和更积极的optimization参数，向最优值迈进一大步。然后，在最后阶段，我们可以有一个更复杂的转换来做最终参数的微调。
一种可能的方案是对初始粗配准级别使用简单的平移变换，然后在较细的配准级别升级为仿射变换。由于我们有两种不同类型的转换，我们可以使用当前示例中所示的多级配准方法。
仿射变换中可优化的参数集具有不同的动态范围。通常，与matrix相关的参数的值在[-1:1]附近，但它们不受此区间的限制。另一方面，与平移相关的参数往往具有更高的值，通常在10.0到100.0之间。这种动态范围的差异对梯度下降优化器的性能有负面影响。ITK提供了一些机制来补偿传递给优化器的参数之间值的这种差异。
第一种机制包括向优化器提供一组scale factors规模因子。这些因素在使用梯度组件在当前迭代中计算优化器的步骤之前对梯度组件进行重新规格化re-normalize。如本章前面的例子所示，这些scales是由用户直观地估计出来的。在我们的特殊情况下，尺度参数的常见选择是将所有与matrix系数相关的参数设置为1.0，即前N乘以N个因子。然后，我们将其余的比例因子设置为一个小的值。
这里，仿射变换由矩阵M和向量T表示。点P到P'的转换表示为
矩阵M负责缩放，旋转和剪切，而矩阵T负责平移。
向量T (T_x, T_y)的平移参数比矩阵M (M_11, M_12, M_21, M_22)的平移参数需要更小的尺度。然而，在处理较大的参数空间时，要直观地估计所有参数尺度是不容易的。
幸运的是，ITKv4为自动参数scaling提供了一个框架。RegistrationParameterScalesEstimator极大地降低了针对不同transform/metric组合调整参数的难度。通过分析小参数更新对形变引起的物理空间变形幅度变化的结果，估计了参数尺度。
参数单位变化的影响可以通过多种方式定义，如索引空间或物理空间中体素的最大位移，或变换雅可比矩阵的平均范数。过滤器itk::RegistrationParameterScalesFromPhysicalShift和itk::RegistrationParameterScalesFromIndexShift使用第一个定义估计规模，而itk:: registrationparameterscalesfrom雅可比滤波器估计规模基于后一个定义。在所有的方法中，目标是缩放变换参数，使每个缩放参数的单位变化对变形产生相同的影响。
步长必须与搜索空间中参数的期望值成比例。由于matrix系数的期望值在1.0左右，因此优化的初始步骤应该小于1.0。作为指导原则，将矩阵系数考虑为cos(θ)和sin(θ)的组合是有用的。这就导致使用接近于以弧度测量的预期旋转的值。例如，一个1.0度的旋转大约是0.017弧度。
重要的是要注意，一旦图像在亚像素级配准，任何进一步的改进配准严重依赖于插值器的质量。它可能然后是合理的使用一个粗糙和快速的插值器在较低的分辨率水平和切换到一个高质量但缓慢的插值器在最终的分辨率水平。然而，在这个例子中，我们使用了一个线性插值为所有阶段和不同的配准级别，因为它是如此之快。
#include "itkImageRegistrationMethodv4.h"//头文件
#include "itkMattesMutualInformationImageToImageMetricv4.h"
#include "itkRegularStepGradientDescentOptimizerv4.h"
#include "itkConjugateGradientLineSearchOptimizerv4.h"
#include "itkTranslationTransform.h"
#include "itkAffineTransform.h"
#include "itkCompositeTransform.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkCheckerBoardImageFilter.h"
#include "itkCommand.h"
template <typename TRegistration>
class RegistrationInterfaceCommand : public itk::Command
{
public:
using Self = RegistrationInterfaceCommand;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
protected:
RegistrationInterfaceCommand() = default;
public:
using RegistrationType = TRegistration;
void
Execute(itk::Object * object, const itk::EventObject & event) override
{
Execute((const itk::Object *)object, event);
}
void
Execute(const itk::Object * object, const itk::EventObject & event) override
{
if (!(itk::MultiResolutionIterationEvent().CheckEvent(&event)))
{
return;
}
std::cout << "\nObserving from class " << object->GetNameOfClass();
if (!object->GetObjectName().empty())
{
std::cout << " \"" << object->GetObjectName() << "\"" << std::endl;
}
const auto * registration = static_cast<const RegistrationType *>(object);
unsigned int currentLevel = registration->GetCurrentLevel();
typename RegistrationType::ShrinkFactorsPerDimensionContainerType
shrinkFactors =
registration->GetShrinkFactorsPerDimension(currentLevel);
typename RegistrationType::SmoothingSigmasArrayType smoothingSigmas =
registration->GetSmoothingSigmasPerLevel();
std::cout << "-------------------------------------" << std::endl;
std::cout << " Current multi-resolution level = " << currentLevel
<< std::endl;
std::cout << " shrink factor = " << shrinkFactors << std::endl;
std::cout << " smoothing sigma = " << smoothingSigmas[currentLevel]
<< std::endl;
std::cout << std::endl;
}
};
class CommandIterationUpdate : public itk::Command
{
public:
using Self = CommandIterationUpdate;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
protected:
CommandIterationUpdate() = default;
public:
using OptimizerType = itk::GradientDescentOptimizerv4Template<double>;
using OptimizerPointer = const OptimizerType *;
void
Execute(itk::Object * caller, const itk::EventObject & event) override
{
Execute((const itk::Object *)caller, event);
}
void
Execute(const itk::Object * object, const itk::EventObject & event) override
{
auto optimizer = static_cast<OptimizerPointer>(object);
if (!(itk::IterationEvent().CheckEvent(&event)))
{
return;
}
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue() << " ";
std::cout << optimizer->GetCurrentPosition() << " "
<< m_CumulativeIterationIndex++ << std::endl;
}
private:
unsigned int m_CumulativeIterationIndex{ 0 };
};
int
main(int argc, char * argv[])
{
if (argc < 4)
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << " outputImagefile [backgroundGrayLevel]";
std::cerr << " [checkerboardbefore] [CheckerBoardAfter]";
std::cerr << " [numberOfBins] " << std::endl;
return EXIT_FAILURE;
}
constexpr unsigned int Dimension = 2;
using PixelType = float;
using FixedImageType = itk::Image<PixelType, Dimension>;
using MovingImageType = itk::Image<PixelType, Dimension>;
using TTransformType = itk::TranslationTransform<double, Dimension>;//首先，我们需要配置初始阶段的配准组件。transform类型的实例化仅需要空间的维数和用于表示空间坐标的类型。
using TOptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;//这里定义了其他配准组件的类型。使用itk::RegularStepGradientDescentOptimizerv4作为第一阶段的优化器。
using MetricType =
itk::MattesMutualInformationImageToImageMetricv4<FixedImageType,
MovingImageType>;//此外，我们使用itk::MattesMutualInformationImageToImageMetricv4作为度量，因为它适合多模式配准。
using TRegistrationType = itk::ImageRegistrationMethodv4<FixedImageType,
MovingImageType,
TTransformType>;
//所有组件都使用它们的New（）方法实例化并连接到配准对象，
TOptimizerType::Pointer transOptimizer = TOptimizerType::New();
MetricType::Pointer transMetric = MetricType::New();
TRegistrationType::Pointer transRegistration = TRegistrationType::New();
transRegistration->SetOptimizer(transOptimizer);
transRegistration->SetMetric(transMetric);
TTransformType::Pointer movingInitTx = TTransformType::New();//配准过程的输出transform将在配准过滤器内部构造，因为相关的TransformType已经作为模板参数传递给注册方法。但是，如果需要，我们应该为配准方法提供一个初始移动变换。
using ParametersType = TOptimizerType::ParametersType;
ParametersType initialParameters(movingInitTx->GetNumberOfParameters());
initialParameters[0] = 3.0; // Initial offset in mm along X
initialParameters[1] = 5.0; // Initial offset in mm along Y
movingInitTx->SetParameters(initialParameters);
transRegistration->SetMovingInitialTransform(movingInitTx);//设置初始参数后，可以通过SetMovingInitialTransform（）方法将初始变换传递到配准过滤器。
using CompositeTransformType = itk::CompositeTransform<double, Dimension>;//我们可以使用itk::CompositeTransform来堆叠来自多个阶段的所有输出transforms。这个复合转换还应该包含移动的初始转换(如果存在的话)，因为正如3.6.1节中解释的那样，每个配准阶段的输出不包括到该阶段的输入初始转换。
CompositeTransformType::Pointer compositeTransform =
CompositeTransformType::New();
compositeTransform->AddTransform(movingInitTx);
using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
fixedImageReader->SetFileName(argv[1]);
movingImageReader->SetFileName(argv[2]);
transRegistration->SetFixedImage(fixedImageReader->GetOutput());
transRegistration->SetMovingImage(movingImageReader->GetOutput());
transRegistration->SetObjectName("TranslationRegistration");
constexpr unsigned int numberOfLevels1 = 1;//在此简单示例的情况下，第一阶段仅以粗略分辨率在一个配准级别中运行。
TRegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel1;
shrinkFactorsPerLevel1.SetSize(numberOfLevels1);
shrinkFactorsPerLevel1[0] = 3;
TRegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel1;
smoothingSigmasPerLevel1.SetSize(numberOfLevels1);
smoothingSigmasPerLevel1[0] = 2;
transRegistration->SetNumberOfLevels(numberOfLevels1);
transRegistration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel1);
transRegistration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel1);
transMetric->SetNumberOfHistogramBins(24);
if (argc > 7)
{
transMetric->SetNumberOfHistogramBins(std::stoi(argv[7]));
}
transOptimizer->SetNumberOfIterations(200);//同样，对于这个初始阶段，我们可以通过采用较大的步长和放宽停止条件来为优化器使用更具积极的参数集。
transOptimizer->SetRelaxationFactor(0.5);
transOptimizer->SetLearningRate(16);
transOptimizer->SetMinimumStepLength(1.5);
CommandIterationUpdate::Pointer observer1 = CommandIterationUpdate::New();
transOptimizer->AddObserver(itk::IterationEvent(), observer1);
using TranslationCommandType =
RegistrationInterfaceCommand<TRegistrationType>;
TranslationCommandType::Pointer command1 = TranslationCommandType::New();
transRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
command1);
try
{
transRegistration->Update();//一旦所有配准组件就位，我们通过调用Update()来触发配准过程，并将结果输出transform添加到最终的复合变换，因此该复合变换可用于初始化下一个配准阶段。
std::cout
<< "Optimizer stop condition: "
<< transRegistration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch (const itk::ExceptionObject & err)
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return EXIT_FAILURE;
}
compositeTransform->AddTransform(
transRegistration->GetModifiableTransform());
using ATransformType = itk::AffineTransform<double, Dimension>;//现在我们可以升级到一个仿射变换作为配准过程的第二阶段。仿射变换是将直线映射为直线的线性变换。它可以用来表示平移、旋转、各向异性缩放、剪切或它们的任何组合。关于仿射变换的细节可以在3.9.16节中看到。转换类型的实例化只需要空间的维度和用于表示空间坐标的类型。
using AOptimizerType =
itk::ConjugateGradientLineSearchOptimizerv4Template<double>;//第二阶段的配置中使用了不同的优化程序，同时该度量标准保持与以前相同。
using ARegistrationType = itk::ImageRegistrationMethodv4<FixedImageType,
MovingImageType,
ATransformType>;
//同样，所有组件都使用它们的New（）方法实例化并连接到配准对象，就像之前的阶段一样。
AOptimizerType::Pointer affineOptimizer = AOptimizerType::New();
MetricType::Pointer affineMetric = MetricType::New();
ARegistrationType::Pointer affineRegistration = ARegistrationType::New();
affineRegistration->SetOptimizer(affineOptimizer);
affineRegistration->SetMetric(affineMetric);
//可以使用配准的初始变换和前一阶段的结果变换来初始化当前阶段，以便将两者都合并为复合变换。
affineRegistration->SetMovingInitialTransform(compositeTransform);
affineRegistration->SetFixedImage(fixedImageReader->GetOutput());
affineRegistration->SetMovingImage(movingImageReader->GetOutput());
affineRegistration->SetObjectName("AffineRegistration");
affineMetric->SetNumberOfHistogramBins(24);
if (argc > 7)
{
affineMetric->SetNumberOfHistogramBins(std::stoi(argv[7]));
}
ATransformType::Pointer affineTx = ATransformType::New();//用户必须在配准方法之外设置配准变换的fixed参数，因此首先我们实例化输出transform类型的对象。
using SpacingType = FixedImageType::SpacingType;
using OriginType = FixedImageType::PointType;
using RegionType = FixedImageType::RegionType;
using SizeType = FixedImageType::SizeType;
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
const SpacingType fixedSpacing = fixedImage->GetSpacing();
const OriginType fixedOrigin = fixedImage->GetOrigin();
const RegionType fixedRegion = fixedImage->GetLargestPossibleRegion();
const SizeType fixedSize = fixedRegion.GetSize();
ATransformType::InputPointType centerFixed;//在仿射变换中，旋转中心由fixed参数集定义，默认设置为[0,0]。但是，考虑这样一种情况，即运行配准的virtual空间的原点远离零原点。在这种情况下，将旋转中心作为默认值会使优化过程不稳定。因此，我们一直热衷于将旋转中心设置在virtual空间的中心，virtual空间通常是fixed的图像空间。
//注意，任何一个重心或几何中心都可以用作旋转中心。在这个例子中，旋转中心被设置为fixed图像的几何中心。我们也可以使用itk::ImageMomentsCalculator 过滤器来计算质心。
//然后，我们计算fixed图像的物理中心，并将其设置为输出仿射变换的中心。
centerFixed[0] = fixedOrigin[0] + fixedSpacing[0] * fixedSize[0] / 2.0;
centerFixed[1] = fixedOrigin[1] + fixedSpacing[1] * fixedSize[1] / 2.0;
const unsigned int numberOfFixedParameters =
affineTx->GetFixedParameters().Size();
ATransformType::ParametersType fixedParameters(numberOfFixedParameters);
for (unsigned int i = 0; i < numberOfFixedParameters; ++i)
{
fixedParameters[i] = centerFixed[i];
}
affineTx->SetFixedParameters(fixedParameters);
affineRegistration->SetInitialTransform(affineTx);//然后，应使用SetInitialTransform（）方法将初始化后的输出transform连接到配准对象。重要的是要区分SetInitialTransform()和SetMovingInitialTransform()，后者用于根据前几个阶段的结果初始化配准阶段。您可以假设第一个用于直接操作当前配准过程中的可优化变换。
using ScalesEstimatorType =
itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;
ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();//选择physical滤波器以估计参数scales。scales estimator然后将传递给optimizer。
scalesEstimator->SetMetric(affineMetric);
scalesEstimator->SetTransformForward(true);
//然而，我们不需要担心上述考虑。由于ScalesEstimator，初始步长也可以自动估计，可以在每次迭代或仅在第一次迭代。在这个例子中，我们选择在注册过程开始时估计一次学习率。
affineOptimizer->SetScalesEstimator(scalesEstimator);
affineOptimizer->SetDoEstimateLearningRateOnce(true);
affineOptimizer->SetDoEstimateLearningRateAtEachIteration(false);
affineOptimizer->SetLowerLimit(0);
affineOptimizer->SetUpperLimit(2);
affineOptimizer->SetEpsilon(0.2);
affineOptimizer->SetNumberOfIterations(200);
affineOptimizer->SetMinimumConvergenceValue(1e-6);
affineOptimizer->SetConvergenceWindowSize(5);
CommandIterationUpdate::Pointer observer2 = CommandIterationUpdate::New();
affineOptimizer->AddObserver(itk::IterationEvent(), observer2);
constexpr unsigned int numberOfLevels2 = 2;//我们运行两个配准级别，其中第二个级别以全分辨率运行，在此我们对输出参数进行最终调整。
ARegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel2;
shrinkFactorsPerLevel2.SetSize(numberOfLevels2);
shrinkFactorsPerLevel2[0] = 2;
shrinkFactorsPerLevel2[1] = 1;
ARegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel2;
smoothingSigmasPerLevel2.SetSize(numberOfLevels2);
smoothingSigmasPerLevel2[0] = 1;
smoothingSigmasPerLevel2[1] = 0;
affineRegistration->SetNumberOfLevels(numberOfLevels2);
affineRegistration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel2);
affineRegistration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel2);
using AffineCommandType = RegistrationInterfaceCommand<ARegistrationType>;
AffineCommandType::Pointer command2 = AffineCommandType::New();
affineRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
command2);
try
{
//最后，我们通过调用Update()触发配准过程，并将最后阶段的输出转换添加到复合转换中。此复合变换将被视为此多级配准过程的最终变换，并将被重采样器用于将moving图像重采样到virtual域空间(如果没有fixed的初始变换，则为fixed的图像空间)。
affineRegistration->Update();
std::cout
<< "Optimizer stop condition: "
<< affineRegistration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch (const itk::ExceptionObject & err)
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return EXIT_FAILURE;
}
compositeTransform->AddTransform(
affineRegistration->GetModifiableTransform());
std::cout << "\nInitial parameters of the registration process: "
<< std::endl
<< movingInitTx->GetParameters() << std::endl;
std::cout << "\nTranslation parameters after registration: " << std::endl
<< transOptimizer->GetCurrentPosition() << std::endl
<< " Last LearningRate: "
<< transOptimizer->GetCurrentStepLength() << std::endl;
std::cout << "\nAffine parameters after registration: " << std::endl
<< affineOptimizer->GetCurrentPosition() << std::endl
<< " Last LearningRate: " << affineOptimizer->GetLearningRate()
<< std::endl;
using ResampleFilterType =
itk::ResampleImageFilter<MovingImageType, FixedImageType>;
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(compositeTransform);
resample->SetInput(movingImageReader->GetOutput());
PixelType backgroundGrayLevel = 100;
if (argc > 4)
{
backgroundGrayLevel = std::stoi(argv[4]);
}
resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resample->SetOutputOrigin(fixedImage->GetOrigin());
resample->SetOutputSpacing(fixedImage->GetSpacing());
resample->SetOutputDirection(fixedImage->GetDirection());
resample->SetDefaultPixelValue(backgroundGrayLevel);
using OutputPixelType = unsigned char;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
using CastFilterType =
itk::CastImageFilter<FixedImageType, OutputImageType>;
using WriterType = itk::ImageFileWriter<OutputImageType>;
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName(argv[3]);
caster->SetInput(resample->GetOutput());
writer->SetInput(caster->GetOutput());
writer->Update();
using CheckerBoardFilterType = itk::CheckerBoardImageFilter<FixedImageType>;
CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();
checker->SetInput1(fixedImage);
checker->SetInput2(resample->GetOutput());
caster->SetInput(checker->GetOutput());
writer->SetInput(caster->GetOutput());
resample->SetDefaultPixelValue(0);
using TransformType = itk::IdentityTransform<double, Dimension>;
TransformType::Pointer identityTransform;
try
{
identityTransform = TransformType::New();
}
catch (const itk::ExceptionObject & err)
{
err.Print(std::cerr);
return EXIT_FAILURE;
}
identityTransform->SetIdentity();
resample->SetTransform(identityTransform);
if (argc > 5)
{
writer->SetFileName(argv[5]);
writer->Update();
}
resample->SetTransform(compositeTransform);
if (argc > 6)
{
writer->SetFileName(argv[6]);
writer->Update();
}
return EXIT_SUCCESS;
}3.8.2 Cascaded Multistage Registration
多阶段配准，第一阶段的输出直接链接到第二阶段，整个配准过程只在最后阶段之后通过调用Update()触发一次。
#include "itkImageRegistrationMethodv4.h"
#include "itkMattesMutualInformationImageToImageMetricv4.h"
#include "itkRegularStepGradientDescentOptimizerv4.h"
#include "itkConjugateGradientLineSearchOptimizerv4.h"
#include "itkTranslationTransform.h"
#include "itkAffineTransform.h"
#include "itkCompositeTransform.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImageMomentsCalculator.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkCheckerBoardImageFilter.h"
#include "itkCommand.h"
template <typename TRegistration>
class RegistrationInterfaceCommand : public itk::Command
{
public:
using Self = RegistrationInterfaceCommand;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
protected:
RegistrationInterfaceCommand() = default;
public:
using RegistrationType = TRegistration;
void
Execute(itk::Object * object, const itk::EventObject & event) override
{
Execute((const itk::Object *)object, event);
}
void
Execute(const itk::Object * object, const itk::EventObject & event) override
{
if (!(itk::MultiResolutionIterationEvent().CheckEvent(&event)))
{
return;
}
std::cout << "\nObserving from class " << object->GetNameOfClass();
if (!object->GetObjectName().empty())
{
std::cout << " \"" << object->GetObjectName() << "\"" << std::endl;
}
const auto * registration = static_cast<const RegistrationType *>(object);
if (registration == nullptr)
{
itkExceptionMacro(<< "Dynamic cast failed, object of type "
<< object->GetNameOfClass());
}
unsigned int currentLevel = registration->GetCurrentLevel();
typename RegistrationType::ShrinkFactorsPerDimensionContainerType
shrinkFactors =
registration->GetShrinkFactorsPerDimension(currentLevel);
typename RegistrationType::SmoothingSigmasArrayType smoothingSigmas =
registration->GetSmoothingSigmasPerLevel();
std::cout << "-------------------------------------" << std::endl;
std::cout << " Current multi-resolution level = " << currentLevel
<< std::endl;
std::cout << " shrink factor = " << shrinkFactors << std::endl;
std::cout << " smoothing sigma = " << smoothingSigmas[currentLevel]
<< std::endl;
std::cout << std::endl;
}
};
class CommandIterationUpdate : public itk::Command
{
public:
using Self = CommandIterationUpdate;
using Superclass = itk::Command;
using Pointer = itk::SmartPointer<Self>;
itkNewMacro(Self);
protected:
CommandIterationUpdate() = default;
public:
using OptimizerType = itk::GradientDescentOptimizerv4Template<double>;
using OptimizerPointer = const OptimizerType *;
void
Execute(itk::Object * caller, const itk::EventObject & event) override
{
Execute((const itk::Object *)caller, event);
}
void
Execute(const itk::Object * object, const itk::EventObject & event) override
{
auto optimizer = static_cast<OptimizerPointer>(object);
if (optimizer == nullptr)
{
return;
}
if (!(itk::IterationEvent().CheckEvent(&event)))
{
return;
}
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue() << " ";
std::cout << optimizer->GetCurrentPosition() << " "
<< m_CumulativeIterationIndex++ << std::endl;
}
private:
unsigned int m_CumulativeIterationIndex{ 0 };
};
int
main(int argc, char * argv[])
{
if (argc < 4)
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile ";
std::cerr << " outputImagefile [backgroundGrayLevel]";
std::cerr << " [checkerboardbefore] [CheckerBoardAfter]";
std::cerr << " [numberOfBins] " << std::endl;
return EXIT_FAILURE;
}
constexpr unsigned int Dimension = 2;
using PixelType = float;
using FixedImageType = itk::Image<PixelType, Dimension>;
using MovingImageType = itk::Image<PixelType, Dimension>;
using TTransformType = itk::TranslationTransform<double, Dimension>;//定义第一阶段的不同类型。
using TOptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
using MetricType =
itk::MattesMutualInformationImageToImageMetricv4<FixedImageType,
MovingImageType>;
using TRegistrationType =
itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;//多阶配准中，类型定义与前面的示例相同，但有一个重要的细微变化:转换类型不再作为模板参数传递给配准方法。在这种情况下，配准过滤器将把转换基类itk:: transform作为输出转换的类型。
TOptimizerType::Pointer transOptimizer = TOptimizerType::New();
MetricType::Pointer transMetric = MetricType::New();
TRegistrationType::Pointer transRegistration = TRegistrationType::New();
transRegistration->SetOptimizer(transOptimizer);
transRegistration->SetMetric(transMetric);
TTransformType::Pointer translationTx = TTransformType::New();//我们没有传递transform类型，而是在配准过滤器外部创建transform对象的显式实例化，并使用SetInitialTransform()方法将其连接到配准对象。另外，通过调用InPlaceOn()方法，这个转换对象将成为配准过滤器的输出转换，或者将被嫁接到输出。
transRegistration->SetInitialTransform(translationTx);
transRegistration->InPlaceOn();
using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
fixedImageReader->SetFileName(argv[1]);
movingImageReader->SetFileName(argv[2]);
transRegistration->SetFixedImage(fixedImageReader->GetOutput());
transRegistration->SetMovingImage(movingImageReader->GetOutput());
transRegistration->SetObjectName("TranslationRegistration");
constexpr unsigned int numberOfLevels1 = 1;
TRegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel1;
shrinkFactorsPerLevel1.SetSize(numberOfLevels1);
shrinkFactorsPerLevel1[0] = 3;
TRegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel1;
smoothingSigmasPerLevel1.SetSize(numberOfLevels1);
smoothingSigmasPerLevel1[0] = 2;
transRegistration->SetNumberOfLevels(numberOfLevels1);
transRegistration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel1);
transRegistration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel1);
transMetric->SetNumberOfHistogramBins(24);
if (argc > 7)
{
transMetric->SetNumberOfHistogramBins(std::stoi(argv[7]));
}
transOptimizer->SetNumberOfIterations(200);
transOptimizer->SetRelaxationFactor(0.5);
transOptimizer->SetLearningRate(16);
transOptimizer->SetMinimumStepLength(1.5);
CommandIterationUpdate::Pointer observer1 = CommandIterationUpdate::New();
transOptimizer->AddObserver(itk::IterationEvent(), observer1);
using TranslationCommandType =
RegistrationInterfaceCommand<TRegistrationType>;
TranslationCommandType::Pointer command1 = TranslationCommandType::New();
transRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
command1);
//在粗分辨率级别上，第一个阶段仅使用一个配准级别运行。但是，请注意，在这一步我们不需要更新平移配准过滤器，因为这一阶段的输出将直接连接到下一阶段的初始输入。由于ITK的管道结构，当我们在最后阶段调用Update()时，第一个阶段也将被更新。现在我们升级到一个仿射变换作为配准过程的第二阶段，和之前一样，我们最初定义并实例化当前配准阶段的不同组件。我们使用了新的优化器，但在新的配置中使用了相同的度量。
//再次注意，TransformType没有传递给配准过滤器的类型定义。这一点很重要，因为当配准过滤器将转换基类itk::transform作为其输出变换的类型时，它可以防止在两个阶段级联时出现类型不匹配。
using ATransformType = itk::AffineTransform<double, Dimension>;
using AOptimizerType =
itk::ConjugateGradientLineSearchOptimizerv4Template<double>;
using ARegistrationType =
itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;
AOptimizerType::Pointer affineOptimizer = AOptimizerType::New();
MetricType::Pointer affineMetric = MetricType::New();
ARegistrationType::Pointer affineRegistration = ARegistrationType::New();
affineRegistration->SetOptimizer(affineOptimizer);
affineRegistration->SetMetric(affineMetric);
affineMetric->SetNumberOfHistogramBins(24);
if (argc > 7)
{
affineMetric->SetNumberOfHistogramBins(std::stoi(argv[7]));
}
fixedImageReader->Update();
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
using FixedImageCalculatorType =
itk::ImageMomentsCalculator<FixedImageType>;
FixedImageCalculatorType::Pointer fixedCalculator =
FixedImageCalculatorType::New();//然后，使用它们的New()方法实例化所有组件，并在转换类型中连接到配准对象。除了前面的例子，这里我们使用fixed图像的质心来初始化仿射变换的fixed参数。ImageMomentsCalculator过滤器用于此目的。
fixedCalculator->SetImage(fixedImage);
fixedCalculator->Compute();
FixedImageCalculatorType::VectorType fixedCenter =
fixedCalculator->GetCenterOfGravity();
ATransformType::Pointer affineTx = ATransformType::New();//然后，我们在Affine变换中初始化fixed参数（旋转中心），并将其连接到配准对象。
const unsigned int numberOfFixedParameters =
affineTx->GetFixedParameters().Size();
ATransformType::ParametersType fixedParameters(numberOfFixedParameters);
for (unsigned int i = 0; i < numberOfFixedParameters; ++i)
{
fixedParameters[i] = fixedCenter[i];
}
affineTx->SetFixedParameters(fixedParameters);
affineRegistration->SetInitialTransform(affineTx);
affineRegistration->InPlaceOn();
affineRegistration->SetFixedImage(fixedImageReader->GetOutput());
affineRegistration->SetMovingImage(movingImageReader->GetOutput());
affineRegistration->SetObjectName("AffineRegistration");
affineRegistration->SetMovingInitialTransformInput(
transRegistration->GetTransformOutput());//现在，第一阶段的输出通过itk::DataObjectDecorator进行包装，并通过SetMovingInitialTransformInput()方法作为moving初始变换传递给第二阶段的输入。注意，这个API在另一个初始化方法SetMovingInitialTransform()的名称后面附加了一个“Input”字，该方法在前面的示例中已经使用过了。此扩展意味着以下API需要数据对象装饰器类型。
using ScalesEstimatorType =
itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;
ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();
scalesEstimator->SetMetric(affineMetric);
scalesEstimator->SetTransformForward(true);
affineOptimizer->SetScalesEstimator(scalesEstimator);
affineOptimizer->SetDoEstimateLearningRateOnce(true);
affineOptimizer->SetDoEstimateLearningRateAtEachIteration(false);
affineOptimizer->SetLowerLimit(0);
affineOptimizer->SetUpperLimit(2);
affineOptimizer->SetEpsilon(0.2);
affineOptimizer->SetNumberOfIterations(200);
affineOptimizer->SetMinimumConvergenceValue(1e-6);
affineOptimizer->SetConvergenceWindowSize(10);
CommandIterationUpdate::Pointer observer2 = CommandIterationUpdate::New();
affineOptimizer->AddObserver(itk::IterationEvent(), observer2);
constexpr unsigned int numberOfLevels2 = 2;
ARegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel2;
shrinkFactorsPerLevel2.SetSize(numberOfLevels2);
shrinkFactorsPerLevel2[0] = 2;
shrinkFactorsPerLevel2[1] = 1;
ARegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel2;
smoothingSigmasPerLevel2.SetSize(numberOfLevels2);
smoothingSigmasPerLevel2[0] = 1;
smoothingSigmasPerLevel2[1] = 0;
affineRegistration->SetNumberOfLevels(numberOfLevels2);
affineRegistration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel2);
affineRegistration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel2);
using AffineCommandType = RegistrationInterfaceCommand<ARegistrationType>;
AffineCommandType::Pointer command2 = AffineCommandType::New();
affineRegistration->AddObserver(itk::MultiResolutionIterationEvent(),
command2);
try
{
//一旦所有配准组件都准备就绪，最后，我们通过在最后一个阶段的配准过滤器上调用Update（）来触发整个配准过程，包括两个级联的配准阶段。
affineRegistration->Update();
std::cout
<< "Optimizer stop condition: "
<< affineRegistration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch (const itk::ExceptionObject & err)
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return EXIT_FAILURE;
}
using CompositeTransformType = itk::CompositeTransform<double, Dimension>;//最后,复合变换用于所有阶段的结果连接在一起,这将被视为最终的输出的多级过程,将传递给重新取样重新取样图像moving到virtual域空间(如果没有fxied的初始fixed图像空间变换)。
CompositeTransformType::Pointer compositeTransform =
CompositeTransformType::New();
compositeTransform->AddTransform(translationTx);
compositeTransform->AddTransform(affineTx);
std::cout << " Translation transform parameters after registration: "
<< std::endl
<< transOptimizer->GetCurrentPosition() << std::endl
<< " Last LearningRate: "
<< transOptimizer->GetCurrentStepLength() << std::endl;
std::cout << " Affine transform parameters after registration: "
<< std::endl
<< affineOptimizer->GetCurrentPosition() << std::endl
<< " Last LearningRate: " << affineOptimizer->GetLearningRate()
<< std::endl;
using ResampleFilterType =
itk::ResampleImageFilter<MovingImageType, FixedImageType>;
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(compositeTransform);
resample->SetInput(movingImageReader->GetOutput());
PixelType backgroundGrayLevel = 100;
if (argc > 4)
{
backgroundGrayLevel = std::stoi(argv[4]);
}
resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resample->SetOutputOrigin(fixedImage->GetOrigin());
resample->SetOutputSpacing(fixedImage->GetSpacing());
resample->SetOutputDirection(fixedImage->GetDirection());
resample->SetDefaultPixelValue(backgroundGrayLevel);
using OutputPixelType = unsigned char;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
using CastFilterType =
itk::CastImageFilter<FixedImageType, OutputImageType>;
using WriterType = itk::ImageFileWriter<OutputImageType>;
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName(argv[3]);
caster->SetInput(resample->GetOutput());
writer->SetInput(caster->GetOutput());
writer->Update();
using CheckerBoardFilterType = itk::CheckerBoardImageFilter<FixedImageType>;
CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();
checker->SetInput1(fixedImage);
checker->SetInput2(resample->GetOutput());
caster->SetInput(checker->GetOutput());
writer->SetInput(caster->GetOutput());
resample->SetDefaultPixelValue(0);
using TransformType = itk::IdentityTransform<double, Dimension>;
TransformType::Pointer identityTransform;
try
{
identityTransform = TransformType::New();
}
catch (const itk::ExceptionObject & err)
{
err.Print(std::cerr);
return EXIT_FAILURE;
}
identityTransform->SetIdentity();
resample->SetTransform(identityTransform);
if (argc > 5)
{
writer->SetFileName(argv[5]);
writer->Update();
}
resample->SetTransform(compositeTransform);
if (argc > 6)
{
writer->SetFileName(argv[6]);
writer->Update();
}
return EXIT_SUCCESS;
}
463
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
