ITK学习笔记（一）ITK的使用

ITK官方文档

第一个ITK程序

1、CMakeLists.txt

# This is the root ITK CMakeLists file.
cmake_minimum_required(VERSION 3.10)
 
# This project is designed to be built outside the Insight source tree.
project(ITK_demo)
 
# Find VTK
set(ITK_DIR D:/ProgramFiles/ITK-5.1.1/lib/cmake/ITK-5.1) # add wmz
find_package(ITK REQUIRED)
include_directories(${ITK_INCLUDE_DIRS})
message("ITK dir = ${ITK_INCLUDE_DIRS}")
message("ITK lib = ${ITK_LIBRARIES}") 

include(${ITK_USE_FILE})
#aux_source_directory(src SRC_LIST) 
set(SRC_LIST
    ./itk_demo.cpp)

add_executable(itk_demo ${SRC_LIST} ) 
target_link_libraries(itk_demo ${ITK_LIBRARIES})

关于 include(${ITK_USE_FILE}) 的说明可以在 UseITK.cmake 中找到：

#  -------------
#
# This file is not part of the ITK API.  It exists purely as an
# implementation detail.  This CMake module may change from version to
# version without notice, or even be removed.
#
# We mean it.
#

# This file sets up include directories, link directories, IO settings and
# compiler settings for a project to use ITK.  It should not be
# included directly, but rather through the ITK_USE_FILE setting
# obtained from ITKConfig.cmake.

2、 测试数据

测试数据下载路径：https://github.com/InsightSoftwareConsortium/ITK/tree/master/Examples/Data
其实编译ITK时的目录下就有需要的测试数据，比如我的ITK-5.1.1目录下。

\ITK-5.1.1\Examples\Data

3、代码

作为第一个示例程序本来应该写一个很简单的像HelloWorld的程序，但是一些比较简单的官网的程序 要么依赖VTK，要么版本高于ITK5.1.1.
所以就找了一个比较长的程序，是一个配准的程序。
代码来自：https://github.com/InsightSoftwareConsortium/ITK/blob/master/Examples/RegistrationITKv4/MultiResImageRegistration1.cxx
我找了一个其他人做过的中文注释版

#include "itkImageRegistrationMethodv4.h"
#include "itkTranslationTransform.h"
#include "itkMeanSquaresImageToImageMetricv4.h"
#include "itkRegularStepGradientDescentOptimizerv4.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkPNGImageIOFactory.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkSubtractImageFilter.h"
/*****************************************************************************************************************
* 本例子是一个图像配准的Demo
*   0、创建了一个Command对象，用于监控配准的过程，被后面的对象调用
*   1、首先要定义像素的维度以及像素类型：进进而链接参考图像以及浮动图像
*   2、定义框架的基本组件：
*               确定变换种类：TransformType：二维变换
*               确定优化方法：OptimizerType：梯度下降
*               确定相似度度量：MetricType：链接两个图像：浮动图像以及参考图像
*   3、创建图像组件，并且通过创建上述框架，进而进行设置（链接）
*   4、设置插值方法：LinearInterpolateImageFunction并且链接在一起
*   5、6：通过ImageFileReader方法进行读取，链接到 registration并更新
*   7、针对前面的TransformType进行实例化：平移变换用于配准SetInitialTransformParameters：用于设置初始值
*   8、针对优化方法的设置：OptimizerType：前面在创建的时候已经设置了其梯度下降方法，此步骤用于对其微调：初始步长，收敛公差，最大迭代次数
*   9、通过RegistrationParameterScalesFromPhysicalShift：将每一个配准要素链接到配准方法中执行，
*   10、实例化Common对象，监控配准过程的执行，触发配准过程--迭代
*   11、通过update函数触发配准的执行
*   12、配准结果定义空间变换的参数序列：其结果由GetLastTransformParameters( )获得并且输出
*           X、Y的变换：TranslationAlongX；TranslationAlongY
*           迭代次数：numberOfIterations
*           最后的结果：bestValue
*       通过CompositeTransform：AddTransform将转换添加到堆栈的背面，并且拥有可优化的参数。
*           也就是说：添加堆栈，副本？？
*   13、14、15、16、ResampleFilterType方法：
*           用变换参数将两幅图像进行叠加比较，并设置重采样滤波器：输入两幅图像
*           输出的是一个变换
*           对滤波器进行相关参数的设置：大小、原点、间距、位置
*           并通过CastFilterType：setInput：weiter进行相关的输出
*           此时：这个图象就是配准结束后的图像
*   17、通过itk::SubtractImageFilter对两幅图像进行比较：
*           fixedImageReader；resampler
*   18、对图像进行处理 itk::RescaleIntensityImageFilter：调节一下亮度；并进行输出
*   19、一致性转发计算参考图像与正在移动图像之间的不同，输出图片5
******************************************************************************************************************/
/*CommandIterationUpdate 类：
    继承Command，监视配准过程的执行。每调用一次，输出相应参数
        object类指向事件的观察者
        Execute方法，类似cellbake，回转
        observer方法：

  */
class CommandIterationUpdate : public itk::Command
{
public:
    typedef CommandIterationUpdate   Self;
    typedef itk::Command             Superclass;
    typedef itk::SmartPointer<Self>  Pointer;
    itkNewMacro(Self);//宏，包装了所有的new（）所有代码
protected:
    CommandIterationUpdate() {};

public:

    typedef itk::RegularStepGradientDescentOptimizerv4<double> OptimizerType;
    typedef const OptimizerType* OptimizerPointer;

    void Execute(itk::Object* caller, const itk::EventObject& event) ITK_OVERRIDE
    {
        Execute((const itk::Object*)caller, event);
    }
    //Object表示激活事件的对象，event表示被激活的事件
    void Execute(const itk::Object* object, const itk::EventObject& event) ITK_OVERRIDE
    {
        OptimizerPointer optimizer = static_cast<OptimizerPointer>(object);
        //checkEvent表示是否观察的对象
        if (!itk::IterationEvent().CheckEvent(&event))
        {
            return;
        }

        std::cout << optimizer->GetCurrentIteration() << " = ";
        std::cout << optimizer->GetValue() << " : ";
        std::cout << optimizer->GetCurrentPosition() << std::endl;
    }

};


int main()
{


    //1、定义图像的维度以及像素执行
    const    unsigned int    Dimension = 2;//定义维度
    typedef  float           PixelType;//图像像素类型
    typedef itk::Image< PixelType, Dimension >  FixedImageType;//输入数据的类型：参考图像 
    typedef itk::Image< PixelType, Dimension >  MovingImageType;//浮动图像



    //2、定义配准框架的基本组件：变换、优化、测度配准组件
    //用于2D图像的一个刚性配准，变换的唯一参数是：空间坐标类型
   //配准
    typedef itk::TranslationTransform< double, Dimension > TransformType;//把参考图像的空间映射到待配准图像的映射
   //优化
    typedef itk::RegularStepGradientDescentOptimizerv4<double> OptimizerType;//优化算法：牛顿梯度下降法
   //度量
    typedef itk::MeanSquaresImageToImageMetricv4<//相似度测量：均方根
        FixedImageType,
        MovingImageType >    MetricType;



    //3、该组件用用于连接其他组件
    typedef itk::ImageRegistrationMethodv4<
        FixedImageType,
        MovingImageType,
        TransformType   >    RegistrationType;

    MetricType::Pointer         metric = MetricType::New();
    OptimizerType::Pointer      optimizer = OptimizerType::New();
    RegistrationType::Pointer   registration = RegistrationType::New();
    //连接组件：变换、优化组件
    registration->SetMetric(metric);
    registration->SetOptimizer(optimizer);

    //4、插值方法
    typedef itk::LinearInterpolateImageFunction<//选择校对机类型，校对机会对配准图像在非网格位置的程度进行评估
        FixedImageType,
        double > FixedLinearInterpolatorType;
    typedef itk::LinearInterpolateImageFunction<
        MovingImageType,
        double > MovingLinearInterpolatorType;
    FixedLinearInterpolatorType::Pointer fixedInterpolator =//每一个配准要素需要其new创建
        FixedLinearInterpolatorType::New();
    MovingLinearInterpolatorType::Pointer movingInterpolator =
        MovingLinearInterpolatorType::New();
    metric->SetFixedInterpolator(fixedInterpolator);
    metric->SetMovingInterpolator(movingInterpolator);

    //5、设置待配准图像以及变换区域
    typedef itk::ImageFileReader< FixedImageType  >   FixedImageReaderType;
    typedef itk::ImageFileReader< MovingImageType >   MovingImageReaderType;
    FixedImageReaderType::Pointer   fixedImageReader = FixedImageReaderType::New();
    MovingImageReaderType::Pointer  movingImageReader = MovingImageReaderType::New();


    //6、读图像！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！
    itk::PNGImageIOFactory::RegisterOneFactory();
    fixedImageReader->SetFileName("E:\\documents\\vs2019\\itk_demo\\build\\RelWithDebInfo\\data\\BrainProtonDensitySliceBorder20.png");//输入图像文件
    movingImageReader->SetFileName("E:\\documents\\vs2019\\itk_demo\\build\\RelWithDebInfo\\data\\BrainProtonDensitySliceShifted13x17y.png");
    //因为图像是从文件读取的，所以下面方法是用于获取图像数据的
    //需要 itk::ImageRegistrationMethod 从 file readers 的输出获取输入
    registration->SetFixedImage(fixedImageReader->GetOutput());
    registration->SetMovingImage(movingImageReader->GetOutput());
    //更新reader，确保其有效
   // fixedImageReader->Update();
   // movingImageReader->Update();


    //7、平移变换用于配准
    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);
    registration->SetMovingInitialTransform(movingInitialTransform);

    //8、准备执行配准方法：对优化器参数进行微调
    TransformType::Pointer   identityTransform = TransformType::New();
    identityTransform->SetIdentity();

    registration->SetFixedInitialTransform(identityTransform);
    //初始振幅的长度用SetMaximumStepLength( ) 定义
    //建立迭代的次数需要谨慎。最大数用SetNumberOfIterations()定义：
    optimizer->SetLearningRate(4);
    optimizer->SetMinimumStepLength(0.001);//优化器的收敛公差
    optimizer->SetRelaxationFactor(0.5);

    //9、将每一个配准要素连接到配准方法执行中
    bool useEstimator = false;

    //useEstimator = atoi(argv[6]) != 0;


    if (useEstimator)
    {


        typedef itk::RegistrationParameterScalesFromPhysicalShift<MetricType> ScalesEstimatorType;
        ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();
        scalesEstimator->SetMetric(metric);
        scalesEstimator->SetTransformForward(true);
        optimizer->SetScalesEstimator(scalesEstimator);
        optimizer->SetDoEstimateLearningRateOnce(true);
    }
    optimizer->SetNumberOfIterations(200);//最大迭代次数


    //10、实例化commend对象，监视配准过程的执行，并处触发配准过程
    CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
    optimizer->AddObserver(itk::IterationEvent(), observer);


    const 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);



    //11、通过调用Update函数触发配准执行
    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;
    }

    //12、配准结果是一系列定义空间变换的参数序列，结果由get获得
    TransformType::ConstPointer transform = registration->GetTransform();
    TransformType::ParametersType finalParameters = transform->GetParameters();
    const double TranslationAlongX = finalParameters[0];//队列中每个元素对应着沿着一个空间维度的平移
    const double TranslationAlongY = finalParameters[1];

    //优化器能够询问抵达收敛的迭代的实际次数并通过GetCurrentIteration()返回出来
    const unsigned int numberOfIterations = optimizer->GetCurrentIteration();//迭代次数
    //最终参数集合的图像量规值通过优化器的GetValue();
    const double bestValue = optimizer->GetValue();//最优化的度量

    //将上述输出
    std::cout << "Result = " << std::endl;
    std::cout << " Translation X = " << TranslationAlongX << std::endl;//输出移动X的值
    std::cout << " Translation Y = " << TranslationAlongY << std::endl;//输出移动Y的值
    std::cout << " Iterations    = " << numberOfIterations << std::endl;//输出迭代次数
    std::cout << " Metric value  = " << bestValue << std::endl;//输出优化的度量

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


    //13、用变换结果将待配准图映射到参考图像中
    typedef itk::ResampleImageFilter<
        MovingImageType,
        FixedImageType >    ResampleFilterType;

    //14、创建一个重采样滤波器，输入待配准图像
    ResampleFilterType::Pointer resampler = ResampleFilterType::New();
    resampler->SetInput(movingImageReader->GetOutput());
    //配准函数生成的变换也作为重采样滤波器的输入被传递
    resampler->SetTransform(outputCompositeTransform);

    //15、ResampleImageFilter要求指定额外的参数，特别是输出图像的间 距、原点和大小
    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);

    //16、滤波器的输出被传递给一个在文件中存储图像的writer
    typedef unsigned char OutputPixelType;

    typedef itk::Image< OutputPixelType, Dimension > OutputImageType;

    typedef itk::CastImageFilter<//转化重采样的像素类型到最终的writer类型
        FixedImageType,
        OutputImageType >          CastFilterType;
    typedef itk::ImageFileWriter< OutputImageType >  WriterType;
    //调用new函数创建新的滤波器
    WriterType::Pointer      writer = WriterType::New();
    CastFilterType::Pointer  caster = CastFilterType::New();
    writer->SetFileName("E:\\documents\\vs2019\\itk_demo\\build\\RelWithDebInfo\\output\\RegistrationITKv4Moving13x17yInputType.png");//写到文件夹位置
    caster->SetInput(resampler->GetOutput());
    writer->SetInput(caster->GetOutput());
    writer->Update();//触发更新

    //17、参照图像和被变换的待配准图像很容易用itk::SubtractImageFilter比较
    //pixel-wise滤波器 计算两幅输入的同源像素的不同：
    typedef itk::SubtractImageFilter<
        FixedImageType,
        FixedImageType,
        FixedImageType > DifferenceFilterType;

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

    difference->SetInput1(fixedImageReader->GetOutput());//不同
    difference->SetInput2(resampler->GetOutput());

    //18、两幅图像的不同也许比较暗，我们用下面方法对其进行调节亮度，使之更加的明显
    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("E:\\documents\\vs2019\\itk_demo\\build\\RelWithDebInfo\\output\\Moving13x17yInputType.png");
    writer2->Update();

    //设置了一致性转换，计算参考图像的不同

    resampler->SetTransform(identityTransform);

    writer2->SetFileName("E:\\documents\\vs2019\\itk_demo\\build\\RelWithDebInfo\\output\\DifferenceBeforeRegistration.png");
    writer2->Update();

    return EXIT_SUCCESS;

}

4、结果

输入图像

输出图像