ITKv4框架 配准2D图像
本例子为itk官方示例ImageRegistration1.cxx。
关键字： 2D图像配准，单模，刚性配准（本处只用了平移）
1、整体配准流程
如上图，
- Fixed image 和 Moving image都从本地读入。
- 配准器包含4个部分，
1）Interpolator： 插值器。
ITKv4新引入了virtual images domains 的概念，该空间为定义了统一resolution和物理坐标系的物理空间。Fixed image 和 Moving image 都变换到该空间，然后进行相似度度量计算，变换过程中需要进行插值操作。本处只有Moving image需要进行插值，是因为通常将Fixed image domains直接作为virtual image domains。
2） Metric ：相似度度量器。
3）Optimizer: 优化器，本例使用RegularStepGradientDescentOptimizerv4。
4） Transform：变换。 - 配准完成之后，只得到了变换关系。resample filter根据变换关系，将Moving image 变换到配准后的位置。
- Subtract Filter同于显示两幅图像的不同，上面的Subtract filter显示配准后的Moving image 与Fixed image的不同，下面的Subtract filter显示的是原始Moving image 与 Fixed image的不同。
2、配准结果
A 官方图像示例
B 自测图像示例
3.源代码及部分注释
#include <iostream>
#include "itkImageRegistrationMethodv4.h" //配准框架
#include "itkTranslationTransform.h" //平移变换
#include "itkMeanSquaresImageToImageMetricv4.h" //度量标准：选择为MeanSquares
#include "itkRegularStepGradientDescentOptimizerv4.h" //优化器:选择为梯度下降
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h" //重采样器
#include "itkCastImageFilter.h" //图像类型转换：例如将float类型转为char
#include "itkRescaleIntensityImageFilter.h" //图像增强显示 more visible
#include "itkSubtractImageFilter.h" //图像差异比较：选择为 直接像素级相减
#include "itkSmartPointer.h" //智能指针声明
#include "itkPNGImageIOFactory.h" //图像工厂注册
#include "itkBMPImageIOFactory.h"
#include "itkTIFFImageIOFactory.h"
#include "itkJPEGImageIOFactory.h"
#include "importLib_ITK.h" //ITK库导入
//定义一个观察器，输出配准的中间迭代信息
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[])
{
itk::PNGImageIOFactory::RegisterOneFactory();
itk::TIFFImageIOFactory::RegisterOneFactory();
itk::BMPImageIOFactory::RegisterOneFactory();
itk::BMPImageIOFactory::RegisterOneFactory();
const char * fixedImageFile = R"(testPicture/BrainProtonDensitySliceBorder20.png)";
const char * movingImageFile = R"(testPicture/BrainProtonDensitySliceShifted13x17y.png)";
const char * outputImageFile = R"(testPicture/outputImageFile.png)";
const char * differenceImageAfterFile = R"(testPicture/differenceImageAfter.png)";
const char * differenceImageBeforeFile = R"(testPicture/differenceImageBefore.png)";
const char * differenceUseEstimator = R"(yes)";
constexpr unsigned int Dimension = 2;
using PixelType = float;
using FixedImageType = itk::Image<PixelType, Dimension>;
using MovingImageType = itk::Image<PixelType, Dimension>;
using TransformType = itk::TranslationTransform<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);
//线性插值器，使得在非网格位置也能进行比较
using FixedLinearInterpolatorType =
itk::LinearInterpolateImageFunction<FixedImageType, double>;
using MovingLinearInterpolatorType =
itk::LinearInterpolateImageFunction<MovingImageType, double>;
FixedLinearInterpolatorType::Pointer fixedInterpolator =
FixedLinearInterpolatorType::New();
MovingLinearInterpolatorType::Pointer movingInterpolator =
MovingLinearInterpolatorType::New();
metric->SetFixedInterpolator(fixedInterpolator);
metric->SetMovingInterpolator(movingInterpolator);
//文件读写器
using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();
fixedImageReader->SetFileName(fixedImageFile);
movingImageReader->SetFileName(movingImageFile);
registration->SetFixedImage(fixedImageReader->GetOutput());
registration->SetMovingImage(movingImageReader->GetOutput());
/**获取参数个数，二维平移只有两个参数*/
TransformType::Pointer movingInitialTransform = TransformType::New();
TransformType::ParametersType initialParameters(
movingInitialTransform->GetNumberOfParameters());
initialParameters[0] = 0.0; // Initial offset in mm along X
initialParameters[1] = 0.0; // Initial offset in mm along Y
movingInitialTransform->SetParameters(initialParameters);
/**设定 fixed image 和 moving image 从virtual image domain 取出时的初始变换
* 默认设定virtual image 与 fixed image 相同，因此FixedInitialTransform设为单位矩阵
*/
registration->SetMovingInitialTransform(movingInitialTransform);
TransformType::Pointer identityTransform = TransformType::New();
identityTransform->SetIdentity();
registration->SetFixedInitialTransform(identityTransform);
optimizer->SetLearningRate(4); //学习率（步长）
optimizer->SetMinimumStepLength(0.001);// 停止优化条件，步长到达该值时认为已经收敛
optimizer->SetRelaxationFactor(0.5); //方向导数激变时，优化器认为达到一个极值点，以该值减少步长
//估算器，使用后更精准
bool useEstimator = false;
useEstimator = (differenceUseEstimator) != 0;
if (useEstimator)
{
using ScalesEstimatorType =
itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;
ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();
scalesEstimator->SetMetric(metric);
scalesEstimator->SetTransformForward(true);
optimizer->SetScalesEstimator(scalesEstimator);
optimizer->SetDoEstimateLearningRateOnce(true);
}
optimizer->SetNumberOfIterations(200); //最多迭代次数
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver(itk::IterationEvent(), observer);
/**ITKv4配准框架可以配准多层（空间收缩，平滑等）图像，
* 本例子只需要配准平移，因此只需要一层
*/
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 std::exception& err)
{
std::cout << "ExceptionObject caught !" <<err.what()<< std::endl;
return EXIT_FAILURE;
}
/************************************************************************/
/* 获得配准结果 */
/************************************************************************/
TransformType::ConstPointer transform = registration->GetTransform();
TransformType::ParametersType finalParameters = transform->GetParameters();
const double TranslationAlongX = finalParameters[0];
const double TranslationAlongY = finalParameters[1];
const unsigned int numberOfIterations
= optimizer->GetCurrentIteration(); //获取停止迭代时的迭代次数
const double bestValue = optimizer->GetValue();//停止时的Metric值
// Print out results
std::cout << "Result = " << std::endl;
std::cout << " Translation X = " << TranslationAlongX << std::endl;
std::cout << " Translation Y = " << TranslationAlongY << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
/************************************************************************/
/* 使用重采样器将moving image 变换到配准后的位置 */
/************************************************************************/
/** 在图像上进行映射，将moving image 映射到 fixed image
* 在此示例中，我们没有使用输出变换的直接初始化，因此，移动初始变换的
* 参数未反映在配准滤波器的输出参数中。 因此，
* 需要复合转换来将初始转换和输出转换连接在一起。
*/
using CompositeTransformType = itk::CompositeTransform<double, Dimension>;
CompositeTransformType::Pointer outputCompositeTransform =
CompositeTransformType::New();
outputCompositeTransform->AddTransform(
movingInitialTransform); //使outputCompositeTransform初始值与movingimage初始值一致？
outputCompositeTransform->AddTransform(
registration->GetModifiableTransform());
using ResampleFilterType =
itk::ResampleImageFilter<MovingImageType, FixedImageType>;
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
resampler->SetInput(movingImageReader->GetOutput());
resampler->SetTransform(outputCompositeTransform);
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(outputImageFile);
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());
using RescalerType =
itk::RescaleIntensityImageFilter<FixedImageType, OutputImageType>;
RescalerType::Pointer intensityRescaler = RescalerType::New();
//对图像放缩到0-255以便于观察
intensityRescaler->SetInput(difference->GetOutput());
intensityRescaler->SetOutputMinimum(0);
intensityRescaler->SetOutputMaximum(255);
resampler->SetDefaultPixelValue(1);
WriterType::Pointer writer2 = WriterType::New();
writer2->SetInput(intensityRescaler->GetOutput());
writer2->SetFileName(differenceImageAfterFile);
writer2->Update();
resampler->SetTransform(identityTransform);
writer2->SetFileName(differenceImageBeforeFile);
writer2->Update();
return EXIT_SUCCESS;
}
337
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
