“hello world“ Registration

#include<iostream>

#include<itkImageRegistrationMethodv4.h>
#include<itkTranslationTransform.h>
#include<itkMeanSquaresImageToImageMetricv4.h>
#include<itkRegularStepGradientDescentOptimizerv4.h>

#include <itkImage.h>
#include<itkNiftiImageIO.h>
#include<itkImageFileReader.h>
#include<itkImageFileWriter.h>

#include<itkResampleImageFilter.h>
#include<itkCastImageFilter.h>
#include<itkRescaleIntensityImageFilter.h>
#include<itkSubtractImageFilter.h>
#include<itkSmartPointer.h>

#include<itkPNGImageIOFactory.h>
#include<itkBMPImageIOFactory.h>
#include<itkTIFFImageIOFactory.h>
#include<itkJPEGImageIOFactory.h>


using ImageType = itk::Image<short, 3>;


// 读取nii数据
template<typename image_type, typename image_pointer>
void readData(const std::string& file_path, image_pointer& out_image) {
    using imageIOType = itk::NiftiImageIO;
    using readerType = itk::ImageFileReader<image_type>;

    auto reader = readerType::New();
    auto nifitiIO = imageIOType::New();

    reader->SetImageIO(nifitiIO);
    reader->SetFileName(file_path);
    try
    {
        reader->Update();
    }
    catch (const itk::ExceptionObject& e)
    {
        std::cout << e.what() << std::endl;
    }

    out_image = reader->GetOutput();
}

// 写入nii数据
template<typename image_type, typename image_pointer>
void writeData(const image_pointer& write_image, const std::string& write_path) {
    using writerType = itk::ImageFileWriter<image_type>;
    auto writer = writerType::New();

    using niiIOType = itk::NiftiImageIO;
    auto niiIO = niiIOType::New();
    writer->SetImageIO(niiIO);

    writer->SetInput(write_image);
    writer->SetFileName(write_path);
    try
    {
        writer->Update();
    }
    catch (const itk::ExceptionObject& e)
    {
        std::cout << e.what() << std::endl;
    }
}


#include<itkCommand.h>
#include<iomanip>

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 << std::setiosflags(std::ios::fixed)<<std::setw(4) << "Iteration "<< optimizer->GetCurrentIteration() << ": opt val = ";
        std::cout <<std::setprecision(2)<<optimizer->GetValue() << ",\tresult = ";
        std::cout << std::setprecision(2)<<optimizer->GetCurrentPosition() << std::endl;
    }
};



int main()
{
    itk::PNGImageIOFactory::RegisterOneFactory();
    itk::TIFFImageIOFactory::RegisterOneFactory();
    itk::BMPImageIOFactory::RegisterOneFactory();
    itk::JPEGImageIOFactory::RegisterOneFactory();

    std::string fixedImagePath
    ("D:\\lib\\itk53_20230526_vs2019\\source\\Examples\\Data\\BrainProtonDensitySliceBorder20.png");
    std::string movingImagePath("D:\\lib\\itk53_20230526_vs2019\\source\\Examples\\Data\\BrainProtonDensitySliceShifted13x17y.png");
    std::string outputImagePath("./outputImage.png");
    std::string differenceImageAfterPath("./differenceImageAfter.png");
    std::string differenceImageBeforePath("./differenceImageBefore.png");

    bool useEstimator = true;

    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>;
    // Software Guide : EndCodeSnippet

    using OptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;
    using MetricType =
        itk::MeanSquaresImageToImageMetricv4<FixedImageType, MovingImageType>;

    using RegistrationType = itk::
        ImageRegistrationMethodv4<FixedImageType, MovingImageType, TransformType>;

    auto metric = MetricType::New();
    auto optimizer = OptimizerType::New();
    auto registration = RegistrationType::New();

    registration->SetMetric(metric);
    registration->SetOptimizer(optimizer);

    using fixedLinearInterpolatorType = itk::LinearInterpolateImageFunction<FixedImageType, double>;
    using movingLinearInterpolatorType = itk::LinearInterpolateImageFunction<MovingImageType, double>;
    auto fixedInterpolator = fixedLinearInterpolatorType::New();
    auto movingInterpolator = movingLinearInterpolatorType::New();

    metric->SetFixedInterpolator(fixedInterpolator);
    metric->SetMovingInterpolator(movingInterpolator);

    using fixedImageReadType = itk::ImageFileReader<FixedImageType>;
    using movingImageReadType = itk::ImageFileReader<MovingImageType>;

    auto fixedImageReader = fixedImageReadType::New();
    auto movingImageReader = movingImageReadType::New();

    fixedImageReader->SetFileName(fixedImagePath);
    movingImageReader->SetFileName(movingImagePath);

    registration->SetFixedImage(fixedImageReader->GetOutput());
    registration->SetMovingImage(movingImageReader->GetOutput());

    TransformType::Pointer movingInitialTransform = TransformType::New();
    TransformType::ParametersType initialParameters(movingInitialTransform->GetNumberOfParameters());
    initialParameters[0] = 0.0;
    initialParameters[1] = 0.0;
    
    movingInitialTransform->SetParameters(initialParameters);

    registration->SetMovingInitialTransform(movingInitialTransform);

    TransformType::Pointer identityTransform = TransformType::New();
    identityTransform->SetIdentity();
    registration->SetFixedInitialTransform(identityTransform);

    optimizer->SetLearningRate(4);
    optimizer->SetMinimumStepLength(1e-3);
    optimizer->SetRelaxationFactor(0.5);

    if (useEstimator)
    {
        using ScaleEstimatorType = itk::RegistrationParameterScalesFromPhysicalShift<MetricType>;
        ScaleEstimatorType::Pointer scaleEstimator = ScaleEstimatorType::New();

        scaleEstimator->SetMetric(metric);
        scaleEstimator->SetTransformForward(true);
        optimizer->SetScalesEstimator(scaleEstimator);
        optimizer->SetDoEstimateLearningRateOnce(true);
    }

    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::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();
    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;


    using CompositeTransformType = itk::CompositeTransform<double, Dimension>;
    CompositeTransformType::Pointer outputCompositeTransform = CompositeTransformType::New();
    outputCompositeTransform->AddTransform(movingInitialTransform);
    outputCompositeTransform->AddTransform(registration->GetModifiableTransform());

    //std::cout << outputCompositeTransform->GetParameters() << std::endl;

    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(outputImagePath);
    caster->SetInput(resampler->GetOutput());

    writer->SetInput(caster->GetOutput());
    try 
    {
        writer->Update();
    }
    catch (const itk::ExceptionObject& e)
    {
        std::cout << e.what() << std::endl;
    }

    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();
    //
    intensityRescaler->SetInput(difference->GetOutput());
    intensityRescaler->SetOutputMinimum(0);
    intensityRescaler->SetOutputMaximum(255);
    resampler->SetDefaultPixelValue(1);

    WriterType::Pointer writer2 = WriterType::New();
    writer2->SetInput(intensityRescaler->GetOutput());
    writer2->SetFileName(differenceImageAfterPath);
    writer2->Update();

    resampler->SetTransform(identityTransform);
    writer2->SetFileName(differenceImageBeforePath);
    writer2->Update();


    return 0;
}