ITK-基于相似性度量的二维平移配准算法-EW帮帮网

作者：翟天保Steven
版权声明：著作权归作者所有，商业转载请联系作者获得授权，非商业转载请注明出处

前言

‌‌图像配准是医学图像处理中一项基础且重要的任务，它的目标是将不同时间、不同模态或不同视角下获取的图像对齐到同一坐标系下。ITK (Insight Segmentation and Registration Toolkit) 作为开源的医学图像处理工具包，提供了强大的图像配准功能。

本文将通过一段完整的ITK代码（参考ITK官方例程），详细介绍二维图像平移配准的实现过程。

环境准备

参见：Windows下用CMake编译ITK及配置测试_itk配置-CSDN博客

功能说明

图像配准包含几个核心要素：：

变换模型：浮动图像与固定图像之间的变换，常见的有平移、旋转、仿射和非线性变换等。
相似性度量：量化固定图像与浮动图像之间的匹配程度，如均方误差、互信息等。
优化器：寻找最佳变换参数使相似性度量达到极值。
插值方法：处理变换过程中的像素重采样问题。

本文案例使用的是平移变换模型和均方误差 (Mean Squares) 相似性度量。

代码模块解析

1. 配准框架搭建

代码首先定义了图像维度 (2D) 和像素类型 (float)，然后创建了配准所需的核心组件：

变换模型：使用itk::TranslationTransform定义二维平移变换
相似性度量：使用itk::MeanSquaresImageToImageMetricv4计算均方误差
优化器：使用itk::RegularStepGradientDescentOptimizerv4实现梯度下降优化
配准方法：使用itk::ImageRegistrationMethodv4整合所有组件

2. 图像读取与预处理

通过itk::ImageFileReader读取固定图像和浮动图像，其中固定图像是作为参考的标准图像，浮动图像是需要变换的图像。代码中还注册了 PNG 图像 IO 工厂，确保能正确处理 PNG 格式图像。

3. 优化过程配置

优化器参数配置是配准的关键部分：

学习率 (Learning Rate)：控制每次迭代的步长大小，本例设为 4
最小步长 (Minimum Step Length)：当步长小于该值时停止迭代，本例设为 0.001
松弛因子 (Relaxation Factor)：用于调整步长，本例设为 0.5
最大迭代次数 (Number of Iterations)：设置为 200 次

此外，代码中还创建了一个迭代观察者 (CommandIterationUpdate)，用于在每次迭代时输出当前的迭代次数、度量值和变换参数，便于监控配准过程。

4. 多分辨率策略

虽然本例中配置了多分辨率策略，但将层级数设为 1，实际上并未启用多分辨率配准。多分辨率策略是一种常用的配准技巧，它从低分辨率图像开始配准，逐步过渡到高分辨率，有助于避免局部最优解。

5. 配准结果处理

配准完成后，代码获取最终的变换参数，并使用itk::ResampleImageFilter将浮动图像重采样到固定图像的空间坐标系中。然后通过类型转换和强度重映射，将结果保存为 PNG 图像。

6. 差异图像计算

为了直观展示配准效果，代码计算了配准前后的差异图像：

配准前差异：使用恒等变换重采样浮动图像，与固定图像相减
配准后差异：使用配准得到的变换重采样浮动图像，与固定图像相减

差异图像中，像素值越接近 0 表示匹配越好，通过对比这两张差异图像可以明显看出配准效果。

使用说明

需安装ITK库及相关依赖。
将ITK官方项目中Examples\Data中的图像文件（文件名字与例程中名字一致），复制到你项目的路径下，并更改代码。
编译时需链接 ITK的库文件。

完整代码

#include <itkImageRegistrationMethodv4.h>
#include <itkTranslationTransform.h>
#include <itkMeanSquaresImageToImageMetricv4.h>
#include <itkRegularStepGradientDescentOptimizerv4.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>
#include <itkCastImageFilter.h>
#include <itkRescaleIntensityImageFilter.h>
#include <itkSubtractImageFilter.h>
#include <itkImageIOFactory.h>
#include <itkPNGImageIOFactory.h>
#include <itkJPEGImageIOFactory.h>

// 迭代过程更新
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()
{
    // 指定输入输出文件名
    const char* fixedImageFile = "ExampleData\\BrainProtonDensitySliceBorder20.png";
    const char* movingImageFile = "ExampleData\\BrainProtonDensitySliceShifted13x17y.png";
    const char* outputImageFile = "Output\\ImageRegistration1Output.png";
    const char* differenceAfterFile = "Output\\ImageRegistration1DifferenceAfter.png";
    const char* differenceBeforeFile = "Output\\ImageRegistration1DifferenceBefore.png";

    // 创建PNG工厂
    itk::PNGImageIOFactory::RegisterOneFactory();

    // 定义图像类型和维度
    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>;

    // 创建配准组件
    auto metric = MetricType::New();
    auto optimizer = OptimizerType::New();
    auto registration = RegistrationType::New();

    // 配置插值器
    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 FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
    using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
    auto fixedImageReader = FixedImageReaderType::New();
    auto movingImageReader = MovingImageReaderType::New();
    fixedImageReader->SetFileName(fixedImageFile);
    movingImageReader->SetFileName(movingImageFile);

    // 连接配准组件
    registration->SetMetric(metric);
    registration->SetOptimizer(optimizer);
    registration->SetFixedImage(fixedImageReader->GetOutput());
    registration->SetMovingImage(movingImageReader->GetOutput());

    // 设置初始变换
    auto movingInitialTransform = TransformType::New();
    TransformType::ParametersType initialParameters(movingInitialTransform->GetNumberOfParameters());
    initialParameters[0] = 0.0; // X轴平移
    initialParameters[1] = 0.0; // Y轴平移
    movingInitialTransform->SetParameters(initialParameters);
    registration->SetMovingInitialTransform(movingInitialTransform);

    // 设置固定图像初始变换（恒等变换）
    auto identityTransform = TransformType::New();
    identityTransform->SetIdentity();
    registration->SetFixedInitialTransform(identityTransform);

    // 配置优化器参数
    optimizer->SetLearningRate(4);
    optimizer->SetMinimumStepLength(0.001);
    optimizer->SetRelaxationFactor(0.5);
    optimizer->SetNumberOfIterations(200);

    // 连接迭代观察者
    auto 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 << "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;

    // 构建输出变换（初始变换+优化后的变换）
    using CompositeTransformType = itk::CompositeTransform<double, Dimension>;
    auto outputCompositeTransform = CompositeTransformType::New();
    outputCompositeTransform->AddTransform(movingInitialTransform);
    outputCompositeTransform->AddTransform(registration->GetModifiableTransform());

    // 重采样浮动图像到固定图像空间
    using ResampleFilterType = itk::ResampleImageFilter<MovingImageType, FixedImageType>;
    auto 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>;
    auto caster = CastFilterType::New();
    auto writer = WriterType::New();
    caster->SetInput(resampler->GetOutput());
    writer->SetInput(caster->GetOutput());
    writer->SetFileName(outputImageFile);
    writer->Update();

    // 计算配准后的差异图像
    using DifferenceFilterType = itk::SubtractImageFilter<FixedImageType, FixedImageType, FixedImageType>;
    auto difference = DifferenceFilterType::New();
    difference->SetInput1(fixedImageReader->GetOutput());
    difference->SetInput2(resampler->GetOutput());
    using RescalerType = itk::RescaleIntensityImageFilter<FixedImageType, OutputImageType>;
    auto intensityRescaler = RescalerType::New();
    intensityRescaler->SetInput(difference->GetOutput());
    intensityRescaler->SetOutputMinimum(0);
    intensityRescaler->SetOutputMaximum(255);
    resampler->SetDefaultPixelValue(1);
    auto writer2 = WriterType::New();
    writer2->SetInput(intensityRescaler->GetOutput());
    writer2->SetFileName(differenceAfterFile);
    writer2->Update();

    // 计算配准前的差异图像（使用恒等变换重采样）
    resampler->SetTransform(identityTransform);
    difference->SetInput2(resampler->GetOutput());
    writer2->SetFileName(differenceBeforeFile);
    writer2->Update();

    return EXIT_SUCCESS;
}