ITK-hello word配准

#include "itkImageRegistrationMethod.h"
#include "itkTranslationTransform.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkSubtractImageFilter.h"
//
#include <iostream>
using namespace std;

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[]) {
//    if(argc < 4) {
//        std::cerr << "Missing Parameters " << std::endl;
//        std::cerr << "Usage: " << argv[0];
//        std::cerr << " fixedImageFile  movingImageFile ";
//        std::cerr << "outputImagefile [differenceImageAfter]";
//        std::cerr << "[differenceImageBefore]" << std::endl;
//        return EXIT_FAILURE;
//    }
//实例化图像类型
    const    unsigned int    Dimension = 2;
    typedef  float           PixelType;

    typedef itk::Image< PixelType, Dimension >  FixedImageType;
    typedef itk::Image< PixelType, Dimension >  MovingImageType;
//实例化变换类型
    typedef itk::TranslationTransform< double, Dimension > 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();
    TransformType::Pointer      transform     = TransformType::New();
    OptimizerType::Pointer      optimizer     = OptimizerType::New();
    InterpolatorType::Pointer   interpolator  = InterpolatorType::New();
    RegistrationType::Pointer   registration  = RegistrationType::New();
//初始化配准对象参数
    registration->SetMetric(metric);
    registration->SetOptimizer(optimizer);
    registration->SetTransform(transform);
    registration->SetInterpolator(interpolator);
//创建数据源对象用于读取图像
    typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
    typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
    FixedImageReaderType::Pointer  fixedImageReader  = FixedImageReaderType::New();
    MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();

    fixedImageReader->SetFileName("../data/BrainProtonDensitySliceBorder20.png");//金标准
    movingImageReader->SetFileName("../data/BrainProtonDensitySliceR10X13Y17.png");//待配准的
//初始化配准中的参考图像及待配准图像
    registration->SetFixedImage(fixedImageReader->GetOutput());
    registration->SetMovingImage(movingImageReader->GetOutput());
//触发读操作
    fixedImageReader->Update();
    registration->SetFixedImageRegion(
        fixedImageReader->GetOutput()->GetBufferedRegion());

    typedef RegistrationType::ParametersType ParametersType;
    ParametersType initialParameters(transform->GetNumberOfParameters());

    initialParameters[0] = 0.0;  // Initial offset in mm along X
    initialParameters[1] = 0.0;  // Initial offset in mm along Y

    registration->SetInitialTransformParameters(initialParameters);

    optimizer->SetMaximumStepLength(4.00);
    optimizer->SetMinimumStepLength(0.01);

    optimizer->SetNumberOfIterations(200);

    CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
    optimizer->AddObserver(itk::IterationEvent(), observer);
    try {
        registration->Update();
    } catch(itk::ExceptionObject &err) {
        std::cerr << "ExceptionObject caught !" << std::endl;
        std::cerr << err << std::endl;
        return EXIT_FAILURE;
    }

    ParametersType finalParameters = registration->GetLastTransformParameters();

    const double TranslationAlongX = finalParameters[0];
    const double TranslationAlongY = finalParameters[1];
    const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
    const double bestValue = optimizer->GetValue();
    //
    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;

    typedef itk::ResampleImageFilter <
    MovingImageType,
    FixedImageType >    ResampleFilterType;

    ResampleFilterType::Pointer resampler = ResampleFilterType::New();
    resampler->SetInput(movingImageReader->GetOutput());

    resampler->SetTransform(registration->GetOutput()->Get());

    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);
//将参考图像类型映射为输出图像类型，便于对图像进行写操作
    typedef unsigned char                            OutputPixelType;
    typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
    typedef itk::CastImageFilter <
    FixedImageType,
    OutputImageType >          CastFilterType;
    //实例化数据槽，用于将重采样结果写入
    typedef itk::ImageFileWriter< OutputImageType >  WriterType;

    WriterType::Pointer      writer =  WriterType::New();
    CastFilterType::Pointer  caster =  CastFilterType::New();
    writer->SetFileName("../save/resample.png");
    caster->SetInput(resampler->GetOutput());
    writer->SetInput(caster->GetOutput());
    writer->Update();

    typedef itk::SubtractImageFilter <
    FixedImageType,
    FixedImageType,
    FixedImageType > DifferenceFilterType;

    DifferenceFilterType::Pointer difference = DifferenceFilterType::New();

    difference->SetInput1(fixedImageReader->GetOutput());
    difference->SetInput2(resampler->GetOutput());

    typedef itk::RescaleIntensityImageFilter <
    FixedImageType,
    OutputImageType >   RescalerType;

    RescalerType::Pointer intensityRescaler = RescalerType::New();

    intensityRescaler->SetInput(difference->GetOutput());
    intensityRescaler->SetOutputMinimum(0);
    intensityRescaler->SetOutputMaximum(255);

    resampler->SetDefaultPixelValue(1);
    
//创建数据源，将相减后的图像写入图像文件
    WriterType::Pointer writer2 = WriterType::New();
    writer2->SetInput(intensityRescaler->GetOutput());
    writer2->SetFileName("../save/xiangjian.png");
    writer2->Update();
//创建变换函数，将变换后的图像写入
    TransformType::Pointer identityTransform = TransformType::New();
    identityTransform->SetIdentity();
    resampler->SetTransform(identityTransform);
    writer2->SetFileName("../save/bianhuan.png");
    writer2->Update();

    return EXIT_SUCCESS;
}