参考：https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch3.html#x26-1740003.18
3.6中心初始化
ITK图像坐标原点通常位于图像的一个角上。当涉及到旋转和缩放时，这会导致反直觉的转换行为。用户倾向于假设旋转和缩放是围绕图像中心的一个固定点进行的。为了弥补预期解释中的这种差异，在工具包中引入了转换中心的概念。该参数通常是一个固定参数，在配准过程中没有进行优化，因此初始化对于获得高效准确的结果至关重要。下面的部分描述初始化转换中心的主要特征和效果。
3.6.1二维刚性配准
为了更好地underestanding配准过程,考虑一个例子，3阶段配准过程使用一个初始启动moving变换(Γ_mi)。通过itk::ImageRegistrationMethodv4类的多个实例，可以处理多个阶段。配准过滤器内的第一阶段,最初的moving变换被添加到一个内部复合变换以及一个可更新的单位变换(Γu)。尽管整个复合变换用于度量评估,只有Γu将被优化器在每次迭代更新。Γu将被视为当前阶段的输出变换融合优化过程。这意味着此阶段的输出不包含初始化参数，因此我们需要将输出和初始化转换连接到一个复合转换，作为第一个注册阶段的最终转换。
T 1(x) = Γmi(Γstage1(x))
上述变换是从视域(即不固定初始变换时的固定图像空间)到移动图像空间的映射。
然后，第一阶段的结果变换将用作配准过程第二阶段的初始moving变换，并且这种方法一直持续到配准过程的最后阶段。
在注册过程结束时，可以将Γmi和每个阶段的输出连接起来，形成最终的复合变换，该变换被认为是整个注册过程的最终输出。
I′m(x) = Im(Γmi(Γstage1(Γstage2(Γstage3(x)))))
如果每个阶段的特征是不同类型的变换，上述方法特别有用，例如，当我们运行一个严格的配准过程，该过程由仿射配准进行，而最后由b样条配准完成。
除了上述方法外，还有一种直接初始化方法，直接对初始变换进行优化。这样，初始变换在配准过程中会被修改，以便在配准过程完成时作为最终变换使用。这种直接的方法在概念上与在ITKv3注册中所发生的事情很接近。
当我们只有一个配准级别时，使用这种方法是非常简单和有效的，这就是本例中的情况。另外,一个很好的应用程序初始化方法在多级场景，是当两个顺向阶段具有相同的变换类型,或至少在初始参数可以很容易地推断来自前一个阶段的结果,例如当平移变换是紧随其后的是一个刚性变换。
当前示例显示了直接初始化方法，其中我们尝试直接初始化可优化变换（Γu）的参数。
#include "itkImageRegistrationMethodv4.h"
#include "itkMeanSquaresImageToImageMetricv4.h"
#include "itkRegularStepGradientDescentOptimizerv4.h"
#include "itkEuler2DTransform.h"//追加头文件
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkCommand.h"
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::RegularStepGradientDescentOptimizerv4<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() << std::endl;
}
};
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 [differenceAfterRegistration] ";
std::cerr << " [differenceBeforeRegistration] ";
std::cerr << " [initialStepLength] " << std::endl;
return EXIT_FAILURE;
}
//我们使用Fixed/Moving初始转换来初始化配准配置。这种方法很好地直观地了解了配准方法，特别是当我们打算运行一个多阶段配准过程时，每个阶段的输出可以用来初始化下一个配准阶段。
constexpr unsigned int Dimension = 2;
using PixelType = float;
using FixedImageType = itk::Image<PixelType, Dimension>;
using MovingImageType = itk::Image<PixelType, Dimension>;
using TransformType = itk::Euler2DTransform<double>;//这个类的唯一模板参数是空间坐标的表示类型。
using OptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
using MetricType =
itk::MeanSquaresImageToImageMetricv4<FixedImageType, MovingImageType>;
using RegistrationType = itk::
ImageRegistrationMethodv4<FixedImageType, MovingImageType, TransformType>;
MetricType::Pointer metric = MetricType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetMetric(metric);
registration->SetOptimizer(optimizer);
TransformType::Pointer initialTransform = TransformType::New();//为此，首先，在下面构造初始转换对象。此变换将被初始化，并在配准过程开始时使用其初始参数。
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]);
registration->SetFixedImage(fixedImageReader->GetOutput());
registration->SetMovingImage(movingImageReader->GetOutput());
fixedImageReader->Update();//在本例中，输入图像取自阅读器。下面的代码更新了读取器，以确保在初始化转换时图像参数(size、origin和spacing)是有效的。我们打算使用fixed图像的中心作为旋转中心，然后使用固定图像中心和移动图像中心之间的向量作为旋转后应用的初始平移。
movingImageReader->Update();
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();
TransformType::InputPointType centerFixed;
//使用fixed图像的origin，size和spacing计算旋转中心。
centerFixed[0] = fixedOrigin[0] + fixedSpacing[0] * fixedSize[0] / 2.0;
centerFixed[1] = fixedOrigin[1] + fixedSpacing[1] * fixedSize[1] / 2.0;
MovingImageType::Pointer movingImage = movingImageReader->GetOutput();
const SpacingType movingSpacing = movingImage->GetSpacing();
const OriginType movingOrigin = movingImage->GetOrigin();
const RegionType movingRegion = movingImage->GetLargestPossibleRegion();
const SizeType movingSize = movingRegion.GetSize();
TransformType::InputPointType centerMoving;
//计算moving图像的中心
centerMoving[0] = movingOrigin[0] + movingSpacing[0] * movingSize[0] / 2.0;
centerMoving[1] = movingOrigin[1] + movingSpacing[1] * movingSize[1] / 2.0;
//然后，我们使用SetCenter()方法将fixed图像的中心作为旋转中心传递给它，从而初始化转换。并将平移设置为移动图像中心与固定图像中心之间的向量。最后这个向量是通过方法SetTranslation()传递的。
initialTransform->SetCenter(centerFixed);
initialTransform->SetTranslation(centerMoving - centerFixed);
//零角度初始化旋转。
initialTransform->SetAngle(0.0);
//现在的当前参数初始变换将被设置为配准方法,所以他们可以直接分配给Γu。注意，您不应该将下面的函数与Hello World中使用的SetMoving(Fixed)InitialTransform()方法混淆。
registration->SetInitialTransform(initialTransform);
//请记住，旋转和平移的单位比例是完全不同的。例如,在这里我们知道的参数数组的第一个元素对应的角度以弧度,而其他参数对应的翻译可以用毫米,所以天真的应用梯度下降优化器不会产生平滑变化的参数,因为类似的改变δ的每个参数变换将产生不同程度的影响。因此，我们需要参数尺度来定制每个参数的学习率。我们可以利用优化器提供的缩放功能。
//平移相关的尺度中使用小因子。但是，对于具有较大参数集的变换，用户设置缩放不是很直观。幸运的是，ITKv4提供了一个自动估计参数尺度的框架。
using OptimizerScalesType = OptimizerType::ScalesType;
OptimizerScalesType optimizerScales(
initialTransform->GetNumberOfParameters());
const double translationScale = 1.0 / 1000.0;
optimizerScales[0] = 1.0;
optimizerScales[1] = translationScale;
optimizerScales[2] = translationScale;
optimizer->SetScales(optimizerScales);
double initialStepLength = 0.1;
if (argc > 6)
{
initialStepLength = std::stod(argv[6]);
}
//然后对优化方法的常规参数进行设置。在本例中，我们使用的是itk::RegularStepGradientDescentOptimizerv4。下面，我们定义优化参数如松弛因子、学习率(初始步长)、最小步长和迭代次数。最后两个作为优化的停止标准。
optimizer->SetRelaxationFactor(0.6);
optimizer->SetLearningRate(initialStepLength);
optimizer->SetMinimumStepLength(0.001);
optimizer->SetNumberOfIterations(200);
//optimizer->SetMaximumStepLength( 1.3 );//为了加速收敛，可以方便地使用较大的步长。
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver(itk::IterationEvent(), observer);
constexpr unsigned int numberOfLevels = 1;
RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
shrinkFactorsPerLevel.SetSize(1);
shrinkFactorsPerLevel[0] = 1;
RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
smoothingSigmasPerLevel.SetSize(1);
smoothingSigmasPerLevel[0] = 0;
registration->SetNumberOfLevels(numberOfLevels);
registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);
try
{
registration->Update();
std::cout << "Optimizer stop condition: "
<< registration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch (const itk::ExceptionObject & err)
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
const TransformType::ParametersType finalParameters =
registration->GetOutput()->Get()->GetParameters();
const double finalAngle = finalParameters[0];
const double finalTranslationX = finalParameters[1];
const double finalTranslationY = finalParameters[2];
const double rotationCenterX =
registration->GetOutput()->Get()->GetCenter()[0];
const double rotationCenterY =
registration->GetOutput()->Get()->GetCenter()[1];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
const double finalAngleInDegrees = finalAngle * 180.0 / itk::Math::pi;
std::cout << "Result = " << std::endl;
std::cout << " Angle (radians) = " << finalAngle << std::endl;
std::cout << " Angle (degrees) = " << finalAngleInDegrees << std::endl;
std::cout << " Translation X = " << finalTranslationX << std::endl;
std::cout << " Translation Y = " << finalTranslationY << std::endl;
std::cout << " Fixed Center X = " << rotationCenterX << std::endl;
std::cout << " Fixed Center Y = " << rotationCenterY << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
using ResampleFilterType =
itk::ResampleImageFilter<MovingImageType, FixedImageType>;
// TransformType::ConstPointer finalTransform = TransformType::New();
// TransformType::ConstPointer finalTransform =
// registration->GetTransform();
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(registration->GetTransform());
resample->SetInput(movingImageReader->GetOutput());
resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resample->SetOutputOrigin(fixedImage->GetOrigin());
resample->SetOutputSpacing(fixedImage->GetSpacing());
resample->SetOutputDirection(fixedImage->GetDirection());
resample->SetDefaultPixelValue(100);
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());
try
{
writer->Update();
}
catch (const itk::ExceptionObject & excp)
{
std::cerr << "ExceptionObject while writing the resampled image !"
<< std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
using DifferenceImageType = itk::Image<float, Dimension>;
using DifferenceFilterType = itk::
SubtractImageFilter<FixedImageType, FixedImageType, DifferenceImageType>;
DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
using RescalerType =
itk::RescaleIntensityImageFilter<DifferenceImageType, OutputImageType>;
RescalerType::Pointer intensityRescaler = RescalerType::New();
intensityRescaler->SetOutputMinimum(0);
intensityRescaler->SetOutputMaximum(255);
difference->SetInput1(fixedImageReader->GetOutput());
difference->SetInput2(resample->GetOutput());
resample->SetDefaultPixelValue(1);
intensityRescaler->SetInput(difference->GetOutput());
WriterType::Pointer writer2 = WriterType::New();
writer2->SetInput(intensityRescaler->GetOutput());
try
{
if (argc > 4)
{
writer2->SetFileName(argv[4]);
writer2->Update();
}
TransformType::Pointer identityTransform = TransformType::New();
identityTransform->SetIdentity();
resample->SetTransform(identityTransform);
if (argc > 5)
{
writer2->SetFileName(argv[5]);
writer2->Update();
}
}
catch (const itk::ExceptionObject & excp)
{
std::cerr << "Error while writing difference images" << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}3.6.2 图像矩初始化
CenteredTransformInitializer支持两种操作模式。在第一种模式下，利用图像的origin、size和spacing计算图像的中心作为空间坐标。将fixed图像的中心指定为变换的旋转中心，通过从fixed图像中心到moving图像中心的向量作为变换的初始平移。在第二种模式下，图像中心不是几何计算，而是利用灰度强度的矩。每个图像的质心是使用辅助类itk::ImageMomentsCalculator计算的。通过fixed图像的质心作为变换的旋转中心，通过从fixed图像质心到moving图像质心的矢量作为变换的初始平移。当感兴趣的解剖结构不在图像中心时，第二种手术方式非常方便。在这种情况下，质心对齐提供了比简单使用几何中心更好的粗略初始配准。当两幅图像以不同的成像方式获得时，初始配准的有效性应受到质疑。在这些情况下，一种成像模态的强度质心不一定与另一种成像模态的强度质心相匹配。
#include "itkImageRegistrationMethodv4.h"
#include "itkMeanSquaresImageToImageMetricv4.h"
#include "itkRegularStepGradientDescentOptimizerv4.h"
#include "itkEuler2DTransform.h"//头文件包含
#include "itkCenteredTransformInitializer.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkCommand.h"
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::RegularStepGradientDescentOptimizerv4<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() << std::endl;
}
};
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 [differenceBeforeRegistration] ";
std::cerr << " [differenceAfterRegistration] " << 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 TransformType = itk::Euler2DTransform<double>;//此类的唯一模板参数是空间坐标的表示类型。
using OptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
using MetricType =
itk::MeanSquaresImageToImageMetricv4<FixedImageType, MovingImageType>;
using RegistrationType =
itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType>;
MetricType::Pointer metric = MetricType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetMetric(metric);
registration->SetOptimizer(optimizer);
TransformType::Pointer transform = TransformType::New();
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]);///输入图像从reader。没有必要显式地在读取器上调用Update()，因为CenteredTransformInitializer类将在其初始化过程中调用它。
registration->SetFixedImage(fixedImageReader->GetOutput());
registration->SetMovingImage(movingImageReader->GetOutput());
using TransformInitializerType =
itk::CenteredTransformInitializer<TransformType,
FixedImageType,
MovingImageType>;//初始化器。该类在fixed和moving图像类型以及变换类型上模板化。然后通过调用New()方法并将结果赋值给itk::SmartPointer来构造初始化器。
TransformInitializerType::Pointer initializer =
TransformInitializerType::New();
initializer->SetTransform(transform);//initializer连接transform、fixed、moving图像
initializer->SetFixedImage(fixedImageReader->GetOutput());
initializer->SetMovingImage(movingImageReader->GetOutput());
initializer->MomentsOn();//使用质心配准，如果是几何中心的话就是GeometryOn()
initializer->InitializeTransform();//使用InitializeTransform计算center和平移，结果会直接传递给transform
transform->SetAngle(0.0);//初始化transform的其余参数。
registration->SetInitialTransform(transform);//将已初始化的transform对象将设置到registration方法，并且配准的起点由其初始参数定义。
registration->InPlaceOn();//如果调用InPlaceOn（）方法，则此initialized transform将是输出transform对象或“grafted”到输出。否则，此“ InitialTransform”将被深复制或“clone”到输出。因此将transform对象移植到输出并通过registration方法进行更新。
using OptimizerScalesType = OptimizerType::ScalesType;
OptimizerScalesType optimizerScales(transform->GetNumberOfParameters());
const double translationScale = 1.0 / 1000.0;
optimizerScales[0] = 1.0;
optimizerScales[1] = translationScale;
optimizerScales[2] = translationScale;
optimizer->SetScales(optimizerScales);
optimizer->SetLearningRate(0.1);
optimizer->SetMinimumStepLength(0.001);
optimizer->SetNumberOfIterations(200);
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver(itk::IterationEvent(), observer);
constexpr unsigned int numberOfLevels = 1;
RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
shrinkFactorsPerLevel.SetSize(1);
shrinkFactorsPerLevel[0] = 1;
RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
smoothingSigmasPerLevel.SetSize(1);
smoothingSigmasPerLevel[0] = 0;
registration->SetNumberOfLevels(numberOfLevels);
registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);
try
{
registration->Update();
std::cout << "Optimizer stop condition: "
<< registration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch (const itk::ExceptionObject & err)
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
TransformType::ParametersType finalParameters = transform->GetParameters();
const double finalAngle = finalParameters[0];
const double finalTranslationX = finalParameters[1];
const double finalTranslationY = finalParameters[2];
const double rotationCenterX =
registration->GetOutput()->Get()->GetFixedParameters()[0];
const double rotationCenterY =
registration->GetOutput()->Get()->GetFixedParameters()[1];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
const double finalAngleInDegrees = finalAngle * 180.0 / itk::Math::pi;
std::cout << "Result = " << std::endl;
std::cout << " Angle (radians) " << finalAngle << std::endl;
std::cout << " Angle (degrees) " << finalAngleInDegrees << std::endl;
std::cout << " Translation X = " << finalTranslationX << std::endl;
std::cout << " Translation Y = " << finalTranslationY << std::endl;
std::cout << " Fixed Center X = " << rotationCenterX << std::endl;
std::cout << " Fixed Center Y = " << rotationCenterY << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
TransformType::MatrixType matrix = transform->GetMatrix();
TransformType::OffsetType offset = transform->GetOffset();
std::cout << "Matrix = " << std::endl << matrix << std::endl;
std::cout << "Offset = " << std::endl << offset << std::endl;
using ResampleFilterType =
itk::ResampleImageFilter<MovingImageType, FixedImageType>;
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(transform);
resample->SetInput(movingImageReader->GetOutput());
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resample->SetOutputOrigin(fixedImage->GetOrigin());
resample->SetOutputSpacing(fixedImage->GetSpacing());
resample->SetOutputDirection(fixedImage->GetDirection());
resample->SetDefaultPixelValue(100);
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 DifferenceImageType = itk::Image<float, Dimension>;
using DifferenceFilterType = itk::
SubtractImageFilter<FixedImageType, FixedImageType, DifferenceImageType>;
DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
using OutputPixelType = unsigned char;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
using RescalerType =
itk::RescaleIntensityImageFilter<DifferenceImageType, OutputImageType>;
RescalerType::Pointer intensityRescaler = RescalerType::New();
intensityRescaler->SetOutputMinimum(0);
intensityRescaler->SetOutputMaximum(255);
difference->SetInput1(fixedImageReader->GetOutput());
difference->SetInput2(resample->GetOutput());
resample->SetDefaultPixelValue(1);
intensityRescaler->SetInput(difference->GetOutput());
using WriterType = itk::ImageFileWriter<OutputImageType>;
WriterType::Pointer writer2 = WriterType::New();
writer2->SetInput(intensityRescaler->GetOutput());
try
{
// Compute the difference image between the
// fixed and moving image after registration.
if (argc > 5)
{
writer2->SetFileName(argv[5]);
writer2->Update();
}
// Compute the difference image between the
// fixed and resampled moving image after registration.
TransformType::Pointer identityTransform = TransformType::New();
identityTransform->SetIdentity();
resample->SetTransform(identityTransform);
if (argc > 4)
{
writer2->SetFileName(argv[4]);
writer2->Update();
}
}
catch (const itk::ExceptionObject & excp)
{
std::cerr << "Error while writing difference images" << std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}3.6.3 Similarity Transform in 2D
相似变换可以看作是旋转、平动和均匀(各向同性)缩放的组合。它保留角度并将线性映射到线性映射。这个transform是在toolkit中实现的，它来源于一个刚性的2D转换，并添加了一个缩放参数。
在使用这种转换时，应该注意缩放和转换不是独立的这一事实。就像旋转可以被局部看作平移一样，缩放也会导致局部位移。缩放一般是根据坐标的原点来进行的。但是，我们已经看到在旋转的情况下它是多么的模糊。因此，此转换还允许用户设置特定的中心。这个中心用于旋转和缩放。
#include "itkImageRegistrationMethodv4.h"
#include "itkMeanSquaresImageToImageMetricv4.h"
#include "itkRegularStepGradientDescentOptimizerv4.h"
#include "itkCenteredTransformInitializer.h"
#include "itkSimilarity2DTransform.h"//头文件
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkIdentityTransform.h"
#include "itkCommand.h"
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::RegularStepGradientDescentOptimizerv4<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() << std::endl;
}
};
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 [differenceBeforeRegistration] ";
std::cerr << " [differenceAfterRegistration] ";
std::cerr << " [steplength] ";
std::cerr << " [initialScaling] [initialAngle] ";
std::cerr << 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 TransformType = itk::Similarity2DTransform<double>;//唯一模板参数是空间坐标的表示类型。
using OptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
using MetricType =
itk::MeanSquaresImageToImageMetricv4<FixedImageType, MovingImageType>;
using RegistrationType = itk::
ImageRegistrationMethodv4<FixedImageType, MovingImageType, TransformType>;
MetricType::Pointer metric = MetricType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetMetric(metric);
registration->SetOptimizer(optimizer);
TransformType::Pointer transform = TransformType::New();
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]);
registration->SetFixedImage(fixedImageReader->GetOutput());
registration->SetMovingImage(movingImageReader->GetOutput());
using TransformInitializerType =
itk::CenteredTransformInitializer<TransformType,
FixedImageType,
MovingImageType>;//初始旋转中心和缩放比例以及初始平移的合理值。
TransformInitializerType::Pointer initializer =
TransformInitializerType::New();
initializer->SetTransform(transform);
initializer->SetFixedImage(fixedImageReader->GetOutput());
initializer->SetMovingImage(movingImageReader->GetOutput());
initializer->MomentsOn();
initializer->InitializeTransform();
double initialScale = 1.0;
if (argc > 7)
{
initialScale = std::stod(argv[7]);
}
double initialAngle = 0.0;
if (argc > 8)
{
initialAngle = std::stod(argv[8]);
}
transform->SetScale(initialScale);//初始化transform剩余参数
transform->SetAngle(initialAngle);
registration->SetInitialTransform(transform);//将transform对象设置到registration方法，初始化蚕食用于初始化registration
registration->InPlaceOn();//initialized transform设置到output transform object或者“grafted”到registration process
using OptimizerScalesType = OptimizerType::ScalesType;//scaling、rotation、translation的单位尺度是不同的
OptimizerScalesType optimizerScales(transform->GetNumberOfParameters());
const double translationScale = 1.0 / 100.0;
optimizerScales[0] = 10.0;//放缩
optimizerScales[1] = 1.0;//旋转
optimizerScales[2] = translationScale;//平移，平移的尺度比较小
optimizerScales[3] = translationScale;
optimizer->SetScales(optimizerScales);
double steplength = 1.0;
if (argc > 6)
{
steplength = std::stod(argv[6]);
}
optimizer->SetLearningRate(steplength);//itk::RegularStepGradientDescentOptimizerv4。下面我们定义优化参数，即初始学习率(步长)、最小步长和迭代次数。最后两个作为优化的停止标准。
optimizer->SetMinimumStepLength(0.0001);
optimizer->SetNumberOfIterations(500);
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver(itk::IterationEvent(), observer);
constexpr unsigned int numberOfLevels = 1;
RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
shrinkFactorsPerLevel.SetSize(1);
shrinkFactorsPerLevel[0] = 1;
RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
smoothingSigmasPerLevel.SetSize(1);
smoothingSigmasPerLevel[0] = 0;
registration->SetNumberOfLevels(numberOfLevels);
registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);
try
{
registration->Update();
std::cout << "Optimizer stop condition: "
<< registration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch (const itk::ExceptionObject & err)
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
TransformType::ParametersType finalParameters = transform->GetParameters();
const double finalScale = finalParameters[0];
const double finalAngle = finalParameters[1];
const double finalTranslationX = finalParameters[2];
const double finalTranslationY = finalParameters[3];
const double rotationCenterX =
registration->GetOutput()->Get()->GetFixedParameters()[0];
const double rotationCenterY =
registration->GetOutput()->Get()->GetFixedParameters()[1];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
const double finalAngleInDegrees = finalAngle * 180.0 / itk::Math::pi;
std::cout << std::endl;
std::cout << "Result = " << std::endl;
std::cout << " Scale = " << finalScale << std::endl;
std::cout << " Angle (radians) = " << finalAngle << std::endl;
std::cout << " Angle (degrees) = " << finalAngleInDegrees << std::endl;
std::cout << " Translation X = " << finalTranslationX << std::endl;
std::cout << " Translation Y = " << finalTranslationY << std::endl;
std::cout << " Fixed Center X = " << rotationCenterX << std::endl;
std::cout << " Fixed Center Y = " << rotationCenterY << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
using ResampleFilterType =
itk::ResampleImageFilter<MovingImageType, FixedImageType>;
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
resampler->SetTransform(transform);
resampler->SetInput(movingImageReader->GetOutput());
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resampler->SetOutputOrigin(fixedImage->GetOrigin());
resampler->SetOutputSpacing(fixedImage->GetSpacing());
resampler->SetOutputDirection(fixedImage->GetDirection());
resampler->SetDefaultPixelValue(100);
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(resampler->GetOutput());
writer->SetInput(caster->GetOutput());
writer->Update();
using DifferenceFilterType =
itk::SubtractImageFilter<FixedImageType, FixedImageType, FixedImageType>;
DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
using RescalerType =
itk::RescaleIntensityImageFilter<FixedImageType, OutputImageType>;
RescalerType::Pointer intensityRescaler = RescalerType::New();
intensityRescaler->SetInput(difference->GetOutput());
intensityRescaler->SetOutputMinimum(0);
intensityRescaler->SetOutputMaximum(255);
difference->SetInput1(fixedImageReader->GetOutput());
difference->SetInput2(resampler->GetOutput());
resampler->SetDefaultPixelValue(1);
WriterType::Pointer writer2 = WriterType::New();
writer2->SetInput(intensityRescaler->GetOutput());
if (argc > 5)
{
writer2->SetFileName(argv[5]);
writer2->Update();
}
using IdentityTransformType = itk::IdentityTransform<double, Dimension>;
IdentityTransformType::Pointer identity = IdentityTransformType::New();
if (argc > 4)
{
resampler->SetTransform(identity);
writer2->SetFileName(argv[4]);
writer2->Update();
}
return EXIT_SUCCESS;
}VersorRigid3DTransform的参数空间不是一个向量空间，因为加法在versor分量空间中不是一个封闭的操作。因此，我们需要使用Versor组合操作来更新参数数组的前三个组件(旋转参数)，使用向量加法来更新参数数组的后三个组件(平移参数)。
在ITK的前一个版本中，需要一个特殊的优化器ITK::VersorRigid3DTransformOptimizer来进行配准，以处理versor计算。幸运的是，在ITKv4中，itk::RegularStepGradientDescentOptimizerv4可以用于vector和versor transform优化，因为在新的配准框架中，更新参数的任务被委托给了moving transform本身。UpdateTransformParameters方法是在itk::Transform类中作为一个虚函数实现的，所有派生的Transform类都可以有自己的该函数的实现。由于这一事实，更新函数重新实现了versor转换，因此它可以处理旋转参数的versor组合。
#include "itkImageRegistrationMethodv4.h"
#include "itkMeanSquaresImageToImageMetricv4.h"
#include "itkVersorRigid3DTransform.h"//头文件
#include "itkCenteredTransformInitializer.h"
#include "itkRegularStepGradientDescentOptimizerv4.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkExtractImageFilter.h"
#include "itkCommand.h"
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::RegularStepGradientDescentOptimizerv4<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() << std::endl;
}
};
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 [differenceBeforeRegistration] ";
std::cerr << " [differenceAfterRegistration] ";
std::cerr << " [sliceBeforeRegistration] ";
std::cerr << " [sliceDifferenceBeforeRegistration] ";
std::cerr << " [sliceDifferenceAfterRegistration] ";
std::cerr << " [sliceAfterRegistration] " << std::endl;
return EXIT_FAILURE;
}
constexpr unsigned int Dimension = 3;
using PixelType = float;
using FixedImageType = itk::Image<PixelType, Dimension>;
using MovingImageType = itk::Image<PixelType, Dimension>;
using TransformType = itk::VersorRigid3DTransform<double>;//唯一模板参数是空间坐标的表示类型。
using OptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
using MetricType =
itk::MeanSquaresImageToImageMetricv4<FixedImageType, MovingImageType>;
using RegistrationType = itk::
ImageRegistrationMethodv4<FixedImageType, MovingImageType, TransformType>;
MetricType::Pointer metric = MetricType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetMetric(metric);
registration->SetOptimizer(optimizer);
TransformType::Pointer initialTransform = TransformType::New();//此转换将被初始化，并在配准过程开始时使用其初始参数。
using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
//input images取自readers。这里没有必要显式地在readers上调用Update()，因为itk::CenteredTransformInitializer将把它作为计算的一部分来执行。
fixedImageReader->SetFileName(argv[1]);
movingImageReader->SetFileName(argv[2]);
registration->SetFixedImage(fixedImageReader->GetOutput());
registration->SetMovingImage(movingImageReader->GetOutput());
using TransformInitializerType =
itk::CenteredTransformInitializer<TransformType,
FixedImageType,
MovingImageType>;//下面的代码实例化初始化器的类型。该类在fixed和moving图像类型以及transform类型上模板化。
TransformInitializerType::Pointer initializer =
TransformInitializerType::New();//通过调用New()方法并将结果赋值给智能指针来构造初始化器。
initializer->SetTransform(initialTransform);//initializer连接到transform、fixed、moving图像。
initializer->SetFixedImage(fixedImageReader->GetOutput());
initializer->SetMovingImage(movingImageReader->GetOutput());
initializer->MomentsOn();//使用质心配准，GeometryOn()方法为几何中心
initializer->InitializeTransform();//最后，中心和平移的计算由InitializeTransform()方法触发。结果值将直接传递给transform。
//transform的rotation部分使用itk::Versor进行初始化，它只是一个单位四元数。从transform traits中可以得到VersorType。versor本身定义了用于表示旋转轴的vector的类型。可以将此特征提取为VectorType。下面的代码行创建了一个versor对象，并通过传递一个旋转轴和一个角度来初始化其参数。
using VersorType = TransformType::VersorType;
using VectorType = VersorType::VectorType;
VersorType rotation;
VectorType axis;
axis[0] = 0.0;
axis[1] = 0.0;
axis[2] = 1.0;
constexpr double angle = 0;
rotation.Set(axis, angle);
initialTransform->SetRotation(rotation);
registration->SetInitialTransform(initialTransform);//将initialized transform设置到registration，因此可以使用其初始参数来初始化配准过程。
using OptimizerScalesType = OptimizerType::ScalesType;
OptimizerScalesType optimizerScales(
initialTransform->GetNumberOfParameters());
const double translationScale = 1.0 / 1000.0;
optimizerScales[0] = 1.0;
optimizerScales[1] = 1.0;
optimizerScales[2] = 1.0;
optimizerScales[3] = translationScale;
optimizerScales[4] = translationScale;
optimizerScales[5] = translationScale;
optimizer->SetScales(optimizerScales);
optimizer->SetNumberOfIterations(200);
optimizer->SetLearningRate(0.2);
optimizer->SetMinimumStepLength(0.001);
optimizer->SetReturnBestParametersAndValue(true);
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver(itk::IterationEvent(), observer);
constexpr unsigned int numberOfLevels = 1;
RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
shrinkFactorsPerLevel.SetSize(1);
shrinkFactorsPerLevel[0] = 1;
RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
smoothingSigmasPerLevel.SetSize(1);
smoothingSigmasPerLevel[0] = 0;
registration->SetNumberOfLevels(numberOfLevels);
registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);
try
{
registration->Update();
std::cout << "Optimizer stop condition: "
<< registration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch (const itk::ExceptionObject & err)
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
const TransformType::ParametersType finalParameters =
registration->GetOutput()->Get()->GetParameters();
const double versorX = finalParameters[0];
const double versorY = finalParameters[1];
const double versorZ = finalParameters[2];
const double finalTranslationX = finalParameters[3];
const double finalTranslationY = finalParameters[4];
const double finalTranslationZ = finalParameters[5];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
std::cout << std::endl << std::endl;
std::cout << "Result = " << std::endl;
std::cout << " versor X = " << versorX << std::endl;
std::cout << " versor Y = " << versorY << std::endl;
std::cout << " versor Z = " << versorZ << std::endl;
std::cout << " Translation X = " << finalTranslationX << std::endl;
std::cout << " Translation Y = " << finalTranslationY << std::endl;
std::cout << " Translation Z = " << finalTranslationZ << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
TransformType::Pointer finalTransform = TransformType::New();
finalTransform->SetFixedParameters(
registration->GetOutput()->Get()->GetFixedParameters());
finalTransform->SetParameters(finalParameters);
TransformType::MatrixType matrix = finalTransform->GetMatrix();
TransformType::OffsetType offset = finalTransform->GetOffset();
std::cout << "Matrix = " << std::endl << matrix << std::endl;
std::cout << "Offset = " << std::endl << offset << std::endl;
using ResampleFilterType =
itk::ResampleImageFilter<MovingImageType, FixedImageType>;
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
resampler->SetTransform(finalTransform);
resampler->SetInput(movingImageReader->GetOutput());
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resampler->SetOutputOrigin(fixedImage->GetOrigin());
resampler->SetOutputSpacing(fixedImage->GetSpacing());
resampler->SetOutputDirection(fixedImage->GetDirection());
resampler->SetDefaultPixelValue(100);
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(resampler->GetOutput());
writer->SetInput(caster->GetOutput());
writer->Update();
using DifferenceFilterType =
itk::SubtractImageFilter<FixedImageType, FixedImageType, FixedImageType>;
DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
using RescalerType =
itk::RescaleIntensityImageFilter<FixedImageType, OutputImageType>;
RescalerType::Pointer intensityRescaler = RescalerType::New();
intensityRescaler->SetInput(difference->GetOutput());
intensityRescaler->SetOutputMinimum(0);
intensityRescaler->SetOutputMaximum(255);
difference->SetInput1(fixedImageReader->GetOutput());
difference->SetInput2(resampler->GetOutput());
resampler->SetDefaultPixelValue(1);
WriterType::Pointer writer2 = WriterType::New();
writer2->SetInput(intensityRescaler->GetOutput());
if (argc > 5)
{
writer2->SetFileName(argv[5]);
writer2->Update();
}
using IdentityTransformType = itk::IdentityTransform<double, Dimension>;
IdentityTransformType::Pointer identity = IdentityTransformType::New();
if (argc > 4)
{
resampler->SetTransform(identity);
writer2->SetFileName(argv[4]);
writer2->Update();
}
using OutputSliceType = itk::Image<OutputPixelType, 2>;
using ExtractFilterType =
itk::ExtractImageFilter<OutputImageType, OutputSliceType>;
ExtractFilterType::Pointer extractor = ExtractFilterType::New();
extractor->SetDirectionCollapseToSubmatrix();
extractor->InPlaceOn();
FixedImageType::RegionType inputRegion =
fixedImage->GetLargestPossibleRegion();
FixedImageType::SizeType size = inputRegion.GetSize();
FixedImageType::IndexType start = inputRegion.GetIndex();
size[2] = 0;
start[2] = 90;
FixedImageType::RegionType desiredRegion;
desiredRegion.SetSize(size);
desiredRegion.SetIndex(start);
extractor->SetExtractionRegion(desiredRegion);
using SliceWriterType = itk::ImageFileWriter<OutputSliceType>;
SliceWriterType::Pointer sliceWriter = SliceWriterType::New();
sliceWriter->SetInput(extractor->GetOutput());
if (argc > 6)
{
extractor->SetInput(caster->GetOutput());
resampler->SetTransform(identity);
sliceWriter->SetFileName(argv[6]);
sliceWriter->Update();
}
if (argc > 7)
{
extractor->SetInput(intensityRescaler->GetOutput());
resampler->SetTransform(identity);
sliceWriter->SetFileName(argv[7]);
sliceWriter->Update();
}
if (argc > 8)
{
resampler->SetTransform(finalTransform);
sliceWriter->SetFileName(argv[8]);
sliceWriter->Update();
}
if (argc > 9)
{
extractor->SetInput(caster->GetOutput());
resampler->SetTransform(finalTransform);
sliceWriter->SetFileName(argv[9]);
sliceWriter->Update();
}
return EXIT_SUCCESS;
}You are strongly encouraged to run the example code, since only in this way can you gain first-hand experience with the behavior of the registration process. Once again, this is a simple reflection of the philosophy that we put forward in this book:
If you can not replicate it, then it does not exist!
We have seen enough published papers with pretty pictures, presenting results that in practice are impossible to replicate. That is vanity, not science.
3.6.5 Centered Initialized Affine Transform
#include "itkImageRegistrationMethodv4.h"
#include "itkMeanSquaresImageToImageMetricv4.h"
#include "itkRegularStepGradientDescentOptimizerv4.h"
#include "itkCenteredTransformInitializer.h"
#include "itkAffineTransform.h"//仿射变换的头文件
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkCommand.h"
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::RegularStepGradientDescentOptimizerv4<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();
vnl_matrix<double> p(2, 2);
p[0][0] = (double)optimizer->GetCurrentPosition()[0];
p[0][1] = (double)optimizer->GetCurrentPosition()[1];
p[1][0] = (double)optimizer->GetCurrentPosition()[2];
p[1][1] = (double)optimizer->GetCurrentPosition()[3];
vnl_svd<double> svd(p);
vnl_matrix<double> r(2, 2);
r = svd.U() * vnl_transpose(svd.V());
double angle = std::asin(r[1][0]);
std::cout << " AffineAngle: " << angle * 180.0 / itk::Math::pi
<< std::endl;
}
};
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::endl;
std::cerr << " outputImagefile [differenceBeforeRegistration] "
<< std::endl;
std::cerr << " [differenceAfterRegistration] " << std::endl;
std::cerr << " [stepLength] [maxNumberOfIterations] " << 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 TransformType = itk::AffineTransform<double, Dimension>;//此类的模板参数是空间坐标和空间尺寸的表示类型。
using OptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
using MetricType =
itk::MeanSquaresImageToImageMetricv4<FixedImageType, MovingImageType>;
using RegistrationType = itk::
ImageRegistrationMethodv4<FixedImageType, MovingImageType, TransformType>;
MetricType::Pointer metric = MetricType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetMetric(metric);
registration->SetOptimizer(optimizer);
TransformType::Pointer transform = TransformType::New();//构建transform对象，在配准开始前被初始化
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]);
registration->SetFixedImage(fixedImageReader->GetOutput());
registration->SetMovingImage(movingImageReader->GetOutput());
using TransformInitializerType =
itk::CenteredTransformInitializer<TransformType,
FixedImageType,
MovingImageType>;//使用itk::CenteredTransformInitializer助手类来计算初始旋转中心和平动的合理值。初始化器被设置为使用每个图像的质心作为初始对应校正。
TransformInitializerType::Pointer initializer =
TransformInitializerType::New();
initializer->SetTransform(transform);
initializer->SetFixedImage(fixedImageReader->GetOutput());
initializer->SetMovingImage(movingImageReader->GetOutput());
initializer->MomentsOn();
initializer->InitializeTransform();
registration->SetInitialTransform(transform);//将transform对象传递给registration，通过在配准过程中更新其参数，将其嫁接到配准过滤器的输出变换。
registration->InPlaceOn();
double translationScale = 1.0 / 1000.0;
if (argc > 8)
{
translationScale = std::stod(argv[8]);
}
using OptimizerScalesType = OptimizerType::ScalesType;
OptimizerScalesType optimizerScales(transform->GetNumberOfParameters());
//请记住，缩放、旋转和平移中的单位比例是非常不同的，我们利用了优化器提供的比例功能。我们知道，参数数组的前N×N个元素对应于旋转矩阵因子，后N为与矩阵相乘后要应用的平移分量。
optimizerScales[0] = 1.0;
optimizerScales[1] = 1.0;
optimizerScales[2] = 1.0;
optimizerScales[3] = 1.0;
optimizerScales[4] = translationScale;
optimizerScales[5] = translationScale;
optimizer->SetScales(optimizerScales);
double steplength = 1.0;
if (argc > 6)
{
steplength = std::stod(argv[6]);
}
unsigned int maxNumberOfIterations = 300;
if (argc > 7)
{
maxNumberOfIterations = std::stoi(argv[7]);
}
//并对优化方法的常用参数进行了设置。在本例中，我们像以前一样使用itk::RegularStepGradientDescentOptimizerv4。下面，我们定义学习率(初始步长)、最小步长和迭代次数等优化参数。最后两个作为优化的停止标准。
optimizer->SetLearningRate(steplength);
optimizer->SetMinimumStepLength(0.0001);
optimizer->SetNumberOfIterations(maxNumberOfIterations);
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver(itk::IterationEvent(), observer);
constexpr unsigned int numberOfLevels = 1;
RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
shrinkFactorsPerLevel.SetSize(1);
shrinkFactorsPerLevel[0] = 1;
RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
smoothingSigmasPerLevel.SetSize(1);
smoothingSigmasPerLevel[0] = 0;
registration->SetNumberOfLevels(numberOfLevels);
registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);
try
{
//最后，我们通过调用Update（）方法来触发配准方法的执行。如果抛出任何异常，则将调用放置在try / catch块中。
registration->Update();
std::cout << "Optimizer stop condition: "
<< registration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch (const itk::ExceptionObject & err)
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
//优化收敛后，我们从配准方法中恢复参数。我们还可以使用GetValue（）方法恢复度量的最终值，并使用GetCurrentIteration（）方法恢复最终的迭代次数。
const TransformType::ParametersType finalParameters =
registration->GetOutput()->Get()->GetParameters();
const double finalRotationCenterX = transform->GetCenter()[0];
const double finalRotationCenterY = transform->GetCenter()[1];
const double finalTranslationX = finalParameters[4];
const double finalTranslationY = finalParameters[5];
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
std::cout << "Result = " << std::endl;
std::cout << " Center X = " << finalRotationCenterX << std::endl;
std::cout << " Center Y = " << finalRotationCenterY << std::endl;
std::cout << " Translation X = " << finalTranslationX << std::endl;
std::cout << " Translation Y = " << finalTranslationY << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
vnl_matrix<double> p(2, 2);
p[0][0] = (double)finalParameters[0];
p[0][1] = (double)finalParameters[1];
p[1][0] = (double)finalParameters[2];
p[1][1] = (double)finalParameters[3];
vnl_svd<double> svd(p);
vnl_matrix<double> r(2, 2);
r = svd.U() * vnl_transpose(svd.V());
double angle = std::asin(r[1][0]);
const double angleInDegrees = angle * 180.0 / itk::Math::pi;
std::cout << " Scale 1 = " << svd.W(0) << std::endl;
std::cout << " Scale 2 = " << svd.W(1) << std::endl;
std::cout << " Angle (degrees) = " << angleInDegrees << std::endl;
using ResampleFilterType =
itk::ResampleImageFilter<MovingImageType, FixedImageType>;
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
resampler->SetTransform(transform);
resampler->SetInput(movingImageReader->GetOutput());
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
resampler->SetOutputOrigin(fixedImage->GetOrigin());
resampler->SetOutputSpacing(fixedImage->GetSpacing());
resampler->SetOutputDirection(fixedImage->GetDirection());
resampler->SetDefaultPixelValue(100);
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(resampler->GetOutput());
writer->SetInput(caster->GetOutput());
writer->Update();
using DifferenceFilterType =
itk::SubtractImageFilter<FixedImageType, FixedImageType, FixedImageType>;
DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
difference->SetInput1(fixedImageReader->GetOutput());
difference->SetInput2(resampler->GetOutput());
WriterType::Pointer writer2 = WriterType::New();
using RescalerType =
itk::RescaleIntensityImageFilter<FixedImageType, OutputImageType>;
RescalerType::Pointer intensityRescaler = RescalerType::New();
intensityRescaler->SetInput(difference->GetOutput());
intensityRescaler->SetOutputMinimum(0);
intensityRescaler->SetOutputMaximum(255);
writer2->SetInput(intensityRescaler->GetOutput());
resampler->SetDefaultPixelValue(1);
if (argc > 5)
{
writer2->SetFileName(argv[5]);
writer2->Update();
}
using IdentityTransformType = itk::IdentityTransform<double, Dimension>;
IdentityTransformType::Pointer identity = IdentityTransformType::New();
if (argc > 4)
{
resampler->SetTransform(identity);
writer2->SetFileName(argv[4]);
writer2->Update();
}
return EXIT_SUCCESS;
}
271
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
