ITK-居中变换（配准）

#include "itkImageRegistrationMethod.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkCenteredRigid2DTransform.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:
    typedef  CommandIterationUpdate   Self;
    typedef  itk::Command             Superclass;
    typedef itk::SmartPointer<Self>   Pointer;
    itkNewMacro(Self);

  protected:
    CommandIterationUpdate() {};

  public:
    typedef itk::RegularStepGradientDescentOptimizer  OptimizerType;
    typedef   const OptimizerType                    *OptimizerPointer;

    void Execute(itk::Object *caller, const itk::EventObject &event) ITK_OVERRIDE {
        Execute((const itk::Object *)caller, event);
    }

    void Execute(const itk::Object *object, const itk::EventObject &event) ITK_OVERRIDE {
        OptimizerPointer 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[]) {
    //实例化图像类
    const    unsigned int    Dimension = 2;
    typedef  unsigned char   PixelType;

    typedef itk::Image< PixelType, Dimension >  FixedImageType;
    typedef itk::Image< PixelType, Dimension >  MovingImageType;

    typedef itk::CenteredRigid2DTransform< double > TransformType;
    typedef itk::RegularStepGradientDescentOptimizer       OptimizerType;
    typedef itk::MeanSquaresImageToImageMetric <
    FixedImageType,
    MovingImageType >    MetricType;
    typedef itk:: LinearInterpolateImageFunction <
    MovingImageType,
    double          >    InterpolatorType;
    typedef itk::ImageRegistrationMethod <
    FixedImageType,
    MovingImageType >    RegistrationType;
//创建配准所需的对象
    MetricType::Pointer         metric        = MetricType::New();
    OptimizerType::Pointer      optimizer     = OptimizerType::New();
    InterpolatorType::Pointer   interpolator  = InterpolatorType::New();
    RegistrationType::Pointer   registration  = RegistrationType::New();
//实例化配准对象
    registration->SetMetric(metric);
    registration->SetOptimizer(optimizer);
    registration->SetInterpolator(interpolator);
    TransformType::Pointer  transform = TransformType::New();
    registration->SetTransform(transform);
//实例化数据源，两个输入图像
    typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
    typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;

    FixedImageReaderType::Pointer  fixedImageReader  = FixedImageReaderType::New();
    MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();

    fixedImageReader->SetFileName("../data/BrainT1SliceBorder20.png");
    movingImageReader->SetFileName("../data/BrainProtonDensitySliceR10X13Y17.png");

//设置待配准图像和参考图像
    registration->SetFixedImage(fixedImageReader->GetOutput());
    registration->SetMovingImage(movingImageReader->GetOutput());
    fixedImageReader->Update();
    registration->SetFixedImageRegion(
        fixedImageReader->GetOutput()->GetBufferedRegion());
    fixedImageReader->Update();
    movingImageReader->Update();

//计算参考图像的中心，以及得到最大有效区域
    typedef FixedImageType::SpacingType    SpacingType;
    typedef FixedImageType::PointType      OriginType;
    typedef FixedImageType::RegionType     RegionType;
    typedef FixedImageType::SizeType       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;
    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;
    centerMoving[0] = movingOrigin[0] + movingSpacing[0] * movingSize[0] / 2.0;
    centerMoving[1] = movingOrigin[1] + movingSpacing[1] * movingSize[1] / 2.0;

    transform->SetCenter(centerFixed);//以参考图像中心为旋转中心
    //以两幅图像重心相减的结果作为变换参数
    transform->SetTranslation(centerMoving - centerFixed);
    transform->SetAngle(0.0);
//初始化配准的变换参数
    registration->SetInitialTransformParameters(transform->GetParameters());
//设置优化器在旋转过程中的缩放比例以及迭代过程中的初始步长，收敛不长，迭代次数等
    typedef OptimizerType::ScalesType       OptimizerScalesType;
    OptimizerScalesType optimizerScales(transform->GetNumberOfParameters());
    const double translationScale = 1.0 / 1000.0;
    optimizerScales[0] = 1.0;
    optimizerScales[1] = translationScale;
    optimizerScales[2] = translationScale;
    optimizerScales[3] = translationScale;
    optimizerScales[4] = translationScale;
    optimizer->SetScales(optimizerScales);
    double initialStepLength = 0.1;
    optimizer->SetRelaxationFactor(0.6);
    optimizer->SetMaximumStepLength(initialStepLength);
    optimizer->SetMinimumStepLength(0.001);
    optimizer->SetNumberOfIterations(200);

    //创建优化器的监听器，并跟踪迭代过程
    CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
    optimizer->AddObserver(itk::IterationEvent(), observer);

    try {
        registration->Update();//触发配准操作
        std::cout << "Optimizer stop condition: "
                  << registration->GetOptimizer()->GetStopConditionDescription()
                  << std::endl;
    } catch(itk::ExceptionObject &err) {
        std::cerr << "ExceptionObject caught !" << std::endl;
        std::cerr << err << std::endl;
        return EXIT_FAILURE;
    }
//获取并打印配准结果参数
    OptimizerType::ParametersType finalParameters =
        registration->GetLastTransformParameters();
    const double finalAngle           = finalParameters[0];
    const double finalRotationCenterX = finalParameters[1];
    const double finalRotationCenterY = finalParameters[2];
    const double finalTranslationX    = finalParameters[3];
    const double finalTranslationY    = finalParameters[4];
    const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
    const double bestValue = optimizer->GetValue();
    // Print out results
    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 << " 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;
//对待配准图像进行重采样操作。以参考图像的参数为输出图像的参数
    typedef itk::ResampleImageFilter <
    MovingImageType,
    FixedImageType >    ResampleFilterType;

    TransformType::Pointer finalTransform = TransformType::New();

    finalTransform->SetParameters(finalParameters);
    finalTransform->SetFixedParameters(transform->GetFixedParameters());

    ResampleFilterType::Pointer resample = ResampleFilterType::New();

    resample->SetTransform(finalTransform);
    resample->SetInput(movingImageReader->GetOutput());

    resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
    resample->SetOutputOrigin(fixedImage->GetOrigin());
    resample->SetOutputSpacing(fixedImage->GetSpacing());
    resample->SetOutputDirection(fixedImage->GetDirection());
    resample->SetDefaultPixelValue(100);
//创建数据槽，将重采样结果保存
    typedef itk::ImageFileWriter< FixedImageType >  WriterFixedType;
    WriterFixedType::Pointer      writer =  WriterFixedType::New();
    writer->SetFileName("../save/resample.png");
    writer->SetInput(resample->GetOutput());
    try {
        writer->Update();
    } catch(itk::ExceptionObject &excp) {
        std::cerr << "ExceptionObject while writing the resampled image !" << std::endl;
        std::cerr << excp << std::endl;
        return EXIT_FAILURE;
    }
//利用图像相减滤波器获取参数图像和重采样变换后的待配准原图像见得差异
    typedef itk::Image< float, Dimension > DifferenceImageType;
    typedef itk::SubtractImageFilter <
    FixedImageType,
    FixedImageType,
    DifferenceImageType > DifferenceFilterType;
    DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
    typedef  unsigned char  OutputPixelType;
    typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
    //利用像素强度映射滤波器，将输入图像像素映射到输出图像范围
    typedef itk::RescaleIntensityImageFilter <
    DifferenceImageType,
    OutputImageType >   RescalerType;
    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());
//创建数据草
    typedef itk::ImageFileWriter< OutputImageType >  WriterType;
    WriterType::Pointer      writer2 =  WriterType::New();
    writer2->SetInput(intensityRescaler->GetOutput());
    try {
        //配准后的参考图像和待配准图间的差异写入
        writer2->SetFileName("../save/after.png");
        writer2->Update();
       //配准前的参考图像和待配准图间的差异写入
        TransformType::Pointer identityTransform = TransformType::New();
        identityTransform->SetIdentity();
        resample->SetTransform(identityTransform);
        writer2->SetFileName("../save/before.png");
        writer2->Update();

    } catch(itk::ExceptionObject &excp) {
        std::cerr << "Error while writing difference images" << std::endl;
        std::cerr << excp << std::endl;
        return EXIT_FAILURE;
    }
//
    return EXIT_SUCCESS;
}