ITK-基于Mattes互信息的二维多模态配准算法

翟天保Steven

于 2025-10-10 15:41:59 发布

阅读量713

点赞数 10

CC 4.0 BY-SA版权

分类专栏： ITK 文章标签：算法

本文链接：https://blog.csdn.net/zhaitianbao/article/details/152798253

ITK 专栏收录该内容

25 篇文章

订阅专栏

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

前言

‌‌在医学图像处理中，多模态图像配准是实现不同成像技术（如 MRI 的 T1 序列与质子密度序列、CT 与 PET）图像融合的核心步骤。由于不同模态图像的灰度分布差异显著（例如 T1 加权 MRI 侧重解剖结构，质子密度 MRI 侧重组织水分含量），传统基于灰度相似性的度量（如均方误差 MSE）不再适用。

‌‌ITK（Insight Segmentation and Registration Toolkit）作为开源医学图像处理标杆工具包，提供了专为多模态场景设计的Mattes 互信息相似度度量，可有效捕捉不同模态图像间的统计相关性。本文将基于 ITK 官方例程，详细讲解采用 Mattes 互信息度量与平移变换的二维多模态配准算法实现，重点解决多模态图像因患者轻微移动导致的 “平移类错位” 问题。

环境准备

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

功能说明

图像配准的核心是寻找最优空间变换，使浮动图像（待对齐图像）与固定图像（参考图像）的相似度最高。本方案针对多模态场景设计，核心组件及原理如下：

1. 空间变换：2D 平移变换（TranslationTransform）

平移变换是多模态配准中最基础的刚性变换（不改变图像尺度与旋转角度），仅包含2 个自由度，专门解决 “纯平移” 类错位：

平移：X 轴、Y 轴方向的位移调整（单位：像素 / 物理毫米，取决于图像间距）。
适用场景：多模态医学影像中因患者扫描时轻微平移（如呼吸、肢体移动）导致的图像错位（如 MRI T1 切片与质子密度切片的位置偏差）。

2. 相似度度量：Mattes 互信息（Mattes Mutual Information）

Mattes 互信息是多模态配准的 “核心利器”，其本质是通过计算两图像灰度分布的统计依赖关系衡量相似度，而非直接比较灰度值：

原理：通过构建两图像的联合灰度直方图，估算互信息值 —— 互信息值越大，说明两图像的空间对应关系越紧密（即对齐效果越好）。
优势：不依赖灰度分布一致性，完美适配多模态场景（如 T1 MRI 与质子密度 MRI、CT 与 MRI 的配准）。
关键参数：numberOfBins（直方图箱数量）—— 控制联合直方图的分辨率，通常取 20~30（本文设为 24，平衡计算精度与效率）。

3. 优化器：正则步长梯度下降（RegularStepGradientDescentv4）

优化器的目标是 “最大化 Mattes 互信息值”（ITK 内部通过负互信息转换为最小化问题），找到最优平移参数。本方案配置的优化器核心特性：

自适应步长：初始步长较大（快速接近最优解），迭代中逐步减小（避免在最优解附近震荡）；
松弛因子：设为 0.8，进一步平滑步长调整，提升收敛稳定性；
早停机制：当步长小于MinimumStepLength（本文 0.001）或达到NumberOfIterations（本文 200）时停止；
结果回溯：启用ReturnBestParametersAndValueOn()，确保返回迭代过程中最优的参数（而非最后一次迭代参数）。

4. 采样策略：随机采样（Random Sampling）

多模态配准中全像素计算互信息效率较低，本方案采用20% 随机采样：

原理：从图像中随机选取部分像素点计算互信息，在保证精度的前提下降低计算量；
可重现性：设置固定随机种子（121213），确保每次运行结果一致，便于调试与验证。

5. 可视化工具：棋盘图（CheckerBoard Image）

为直观对比配准效果，代码生成配准前后的棋盘图：

原理：将固定图像与浮动图像按棋盘格模式交替拼接（如 2x2 像素块）；
效果：配准前棋盘格明显错位，配准后棋盘格无缝对齐，可直接通过视觉判断配准精度。

代码模块解析

本节按 “配准流水线” 拆解代码，清晰说明各模块的功能与逻辑关联。

1. 配准框架搭建：定义核心组件并关联

首先定义各组件的类型（适配 2D 浮点型图像），并实例化关联，构建完整配准链路：

// 图像维度与像素类型
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::MattesMutualInformationImageToImageMetricv4<FixedImageType, MovingImageType>; // 互信息度量
using RegistrationType = itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType, TransformType>; // 配准器

// 实例化并关联
auto optimizer = OptimizerType::New();
auto metric = MetricType::New();
auto registration = RegistrationType::New();
registration->SetOptimizer(optimizer); // 配准器绑定优化器
registration->SetMetric(metric);       // 配准器绑定度量

2. 图像读取与 IO 配置：支持多模态图像加载

// 注册PNG IO工厂，支持PNG格式读写
itk::PNGImageIOFactory::RegisterOneFactory();

// 图像读取器
using FixedImageReaderType = itk::ImageFileReader<FixedImageType>;
using MovingImageReaderType = itk::ImageFileReader<MovingImageType>;
auto fixedImageReader = FixedImageReaderType::New();
auto movingImageReader = MovingImageReaderType::New();

// 设置多模态图像路径（固定：T1序列；浮动：质子密度序列）
fixedImageReader->SetFileName("ExampleData\\BrainT1SliceBorder20.png");
movingImageReader->SetFileName("ExampleData\\BrainProtonDensitySliceShifted13x17y.png");

// 配准器绑定输入图像
registration->SetFixedImage(fixedImageReader->GetOutput());
registration->SetMovingImage(movingImageReader->GetOutput());

3. 度量与优化器参数配置：适配多模态场景

（1）Mattes 互信息度量配置

关闭梯度滤波器（多模态图像梯度特征差异大，梯度计算易引入噪声），设置直方图箱数量：

metric->SetNumberOfHistogramBins(24);  // 直方图箱数量（核心参数）
metric->SetUseMovingImageGradientFilter(false); // 关闭浮动图像梯度滤波
metric->SetUseFixedImageGradientFilter(false);  // 关闭固定图像梯度滤波

（2）优化器参数配置

设置学习率、迭代次数等关键参数，确保稳定收敛：

optimizer->SetLearningRate(8.00);          // 初始学习率（步长）
optimizer->SetMinimumStepLength(0.001);    // 最小步长（早停条件1）
optimizer->SetNumberOfIterations(200);     // 最大迭代次数（早停条件2）
optimizer->ReturnBestParametersAndValueOn();// 启用最优参数回溯
optimizer->SetRelaxationFactor(0.8);       // 松弛因子（平滑步长调整）

（3）迭代观察者：监控配准过程

自定义CommandIterationUpdate类，实时输出迭代次数、当前度量值与平移参数，便于调试：

class CommandIterationUpdate : public itk::Command
{
public:
    itkNewMacro(Self); // ITK智能指针构造
    using OptimizerType = itk::RegularStepGradientDescentOptimizerv4<double>;

    void Execute(const itk::Object* object, const itk::EventObject& event) override
    {
        auto optimizer = static_cast<const OptimizerType*>(object);
        if (itk::IterationEvent().CheckEvent(&event))
        {
            // 输出：迭代次数 + 互信息值（负的，因ITK转为最小化问题） + 平移参数
            std::cout << optimizer->GetCurrentIteration() << "   ";
            std::cout << optimizer->GetValue() << "   ";
            std::cout << optimizer->GetCurrentPosition() << std::endl;
        }
    }
};

// 注册观察者到优化器
auto observer = CommandIterationUpdate::New();
optimizer->AddObserver(itk::IterationEvent(), observer);

4. 采样策略配置：平衡效率与精度

设置随机采样比例与种子，降低计算量并确保结果可重现：

// 采样策略：随机采样
registration->SetMetricSamplingStrategy(RegistrationType::MetricSamplingStrategyEnum::RANDOM);
registration->SetMetricSamplingPercentage(0.20); // 采样20%像素
registration->MetricSamplingReinitializeSeed(121213); // 固定随机种子

5. 配准执行与异常处理

调用Update()执行配准，通过try-catch捕获 ITK 异常（如文件路径错误、参数非法）：

try
{
    registration->Update();
    std::cout << "优化器停止条件: " << registration->GetOptimizer()->GetStopConditionDescription() << std::endl;
}
catch (const itk::ExceptionObject& err)
{
    std::cerr << "捕获到异常: " << err << std::endl;
    return EXIT_FAILURE;
}

6. 结果提取与打印

获取最终平移参数、迭代次数与最优度量值，验证配准效果：

// 提取最终平移参数（X、Y方向）
TransformType::ParametersType finalParameters = registration->GetOutput()->Get()->GetParameters();
double TranslationAlongX = finalParameters[0];
double TranslationAlongY = finalParameters[1];

// 提取迭代信息
unsigned int numberOfIterations = optimizer->GetCurrentIteration();
double bestValue = optimizer->GetValue();

// 打印结果
std::cout << "结果 = " << std::endl;
std::cout << " X方向平移 = " << TranslationAlongX << std::endl;
std::cout << " Y方向平移 = " << TranslationAlongY << std::endl;
std::cout << " 迭代次数 = " << numberOfIterations << std::endl;
std::cout << " 最优度量值 = " << bestValue << std::endl;

7. 重采样与结果保存

将浮动图像通过最优变换重采样到固定图像坐标系，转换为 8 位无符号 char 类型并保存为 PNG：

// 重采样滤波器（将浮动图像对齐到固定图像空间）
using ResampleFilterType = itk::ResampleImageFilter<MovingImageType, FixedImageType>;
auto resample = ResampleFilterType::New();
resample->SetTransform(registration->GetTransform()); // 应用最优变换
resample->SetInput(movingImageReader->GetOutput());
resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize()); // 匹配固定图像尺寸
resample->SetOutputOrigin(fixedImage->GetOrigin()); // 匹配固定图像原点
resample->SetOutputSpacing(fixedImage->GetSpacing()); // 匹配固定图像间距
resample->SetDefaultPixelValue(100); // 边界填充像素值

// 类型转换（float -> unsigned char）与保存
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(resample->GetOutput());
writer->SetInput(caster->GetOutput());
writer->SetFileName("Output\\ImageRegistration4Output.png");
writer->Update();

8. 棋盘图生成：可视化配准效果

生成配准前后的棋盘图，直观对比对齐效果：

// 棋盘图滤波器
using CheckerBoardFilterType = itk::CheckerBoardImageFilter<FixedImageType>;
auto checker = CheckerBoardFilterType::New();
checker->SetInput1(fixedImage); // 输入1：固定图像
checker->SetInput2(resample->GetOutput()); // 输入2：浮动图像（重采样后）

// 配准前棋盘图（用恒等变换，即未对齐状态）
auto identityTransform = TransformType::New();
identityTransform->SetIdentity();
resample->SetTransform(identityTransform);
writer->SetFileName("Output\\ImageRegistration4CheckerboardBefore.png");
writer->SetInput(caster->GetInput()); // caster输入已绑定棋盘图
writer->Update();

// 配准后棋盘图（用最优变换）
resample->SetTransform(registration->GetTransform());
writer->SetFileName("Output\\ImageRegistration4CheckerboardAfter.png");
writer->Update();

使用说明

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

完整代码

#include <itkImageRegistrationMethodv4.h>
#include <itkTranslationTransform.h>
#include <itkMeanSquaresImageToImageMetricv4.h>
#include <itkRegularStepGradientDescentOptimizerv4.h>
#include <itkMattesMutualInformationImageToImageMetricv4.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkResampleImageFilter.h>
#include <itkCastImageFilter.h>
#include <itkImageIOFactory.h>
#include <itkPNGImageIOFactory.h>
#include <itkCommand.h>
#include <itkCheckerBoardImageFilter.h>
#include <itkMersenneTwisterRandomVariateGenerator.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(int argc, char* argv[])
{
    // 固定输入输出文件名
    const char* fixedImageFile = "ExampleData\\BrainT1SliceBorder20.png";
    const char* movingImageFile = "ExampleData\\BrainProtonDensitySliceShifted13x17y.png";
    const char* outputImageFile = "Output\\ImageRegistration4Output.png";
    const int defaultPixelValue = 100;  // 默认像素值参数
    const char* checkerboardBeforeFile = "Output\\ImageRegistration4CheckerboardBefore.png";
    const char* checkerboardAfterFile = "Output\\ImageRegistration4CheckerboardAfter.png";
    const unsigned int numberOfBins = 24;  // 直方图箱数量参数

    // 创建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 RegistrationType = itk::ImageRegistrationMethodv4<FixedImageType, MovingImageType, TransformType>; // 配准器类型
    using MetricType = itk::MattesMutualInformationImageToImageMetricv4<FixedImageType, MovingImageType>;  // 互信息度量类型

    auto optimizer = OptimizerType::New();    // 创建优化器实例
    auto registration = RegistrationType::New();  // 创建配准器实例
    registration->SetOptimizer(optimizer);   // 为配准器设置优化器

    auto metric = MetricType::New();  // 创建度量实例
    registration->SetMetric(metric);  // 为配准器设置度量
    metric->SetNumberOfHistogramBins(numberOfBins);  // 设置直方图箱数量
    metric->SetUseMovingImageGradientFilter(false);  // 不使用移动图像梯度滤波器
    metric->SetUseFixedImageGradientFilter(false);   // 不使用固定图像梯度滤波器

    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->SetFixedImage(fixedImageReader->GetOutput());  // 为配准器设置固定图像
    registration->SetMovingImage(movingImageReader->GetOutput());// 为配准器设置移动图像

    optimizer->SetLearningRate(8.00);          // 设置学习率
    optimizer->SetMinimumStepLength(0.001);    // 设置最小步长
    optimizer->SetNumberOfIterations(200);     // 设置最大迭代次数
    optimizer->ReturnBestParametersAndValueOn();// 启用返回最佳参数和值
    optimizer->SetRelaxationFactor(0.8);  // 设置松弛因子

    // 创建命令观察者并将其注册到优化器
    auto observer = CommandIterationUpdate::New();
    optimizer->AddObserver(itk::IterationEvent(), observer);

    // 无收缩和平滑的单级配准过程
    constexpr unsigned int numberOfLevels = 1;

    RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
    shrinkFactorsPerLevel.SetSize(1);
    shrinkFactorsPerLevel[0] = 1;  // 收缩因子（1表示不收缩）

    RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
    smoothingSigmasPerLevel.SetSize(1);
    smoothingSigmasPerLevel[0] = 0;  // 平滑 sigma（0表示不平滑）

    registration->SetNumberOfLevels(numberOfLevels);
    registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
    registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);

    RegistrationType::MetricSamplingStrategyEnum samplingStrategy = RegistrationType::MetricSamplingStrategyEnum::RANDOM;  // 采样策略：随机
    double samplingPercentage = 0.20;  // 采样百分比（20%）
    registration->SetMetricSamplingStrategy(samplingStrategy);  // 设置采样策略
    registration->SetMetricSamplingPercentage(samplingPercentage);  // 设置采样百分比
    registration->MetricSamplingReinitializeSeed(121213);  // 设置随机数种子（保证结果可重现）

    try
    {
        registration->Update();  // 执行配准
        std::cout << "优化器停止条件: "
            << registration->GetOptimizer()->GetStopConditionDescription()
            << std::endl;
    }
    catch (const itk::ExceptionObject& err)
    {
        std::cerr << "捕获到异常!" << std::endl;
        std::cerr << err << std::endl;
        return EXIT_FAILURE;
    }

    // 输出结果
    TransformType::ParametersType finalParameters = registration->GetOutput()->Get()->GetParameters();  // 获取最终变换参数
    double TranslationAlongX = finalParameters[0];  // X方向平移
    double TranslationAlongY = finalParameters[1];  // Y方向平移
    unsigned int numberOfIterations = optimizer->GetCurrentIteration();  // 实际迭代次数
    double bestValue = optimizer->GetValue();  // 最佳度量值 
    std::cout << std::endl;
    std::cout << "结果 = " << std::endl;
    std::cout << " X方向平移 = " << TranslationAlongX << std::endl;
    std::cout << " Y方向平移 = " << TranslationAlongY << std::endl;
    std::cout << " 迭代次数 = " << numberOfIterations << std::endl;
    std::cout << " 度量值 = " << bestValue << std::endl;
    std::cout << " 停止条件 = " << optimizer->GetStopConditionDescription() << std::endl;

    using ResampleFilterType = itk::ResampleImageFilter<MovingImageType, FixedImageType>;  // 重采样滤波器类型
    auto resample = ResampleFilterType::New();  // 创建重采样滤波器
    resample->SetTransform(registration->GetTransform());  // 设置变换
    resample->SetInput(movingImageReader->GetOutput());    // 设置输入图像
    FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();  // 获取固定图像
    resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());  // 设置输出大小
    resample->SetOutputOrigin(fixedImage->GetOrigin());  // 设置输出原点
    resample->SetOutputSpacing(fixedImage->GetSpacing());  // 设置输出间距
    resample->SetOutputDirection(fixedImage->GetDirection());  // 设置输出方向
    resample->SetDefaultPixelValue(defaultPixelValue);  // 设置默认像素值

    using OutputPixelType = unsigned char;  // 输出像素类型
    using OutputImageType = itk::Image<OutputPixelType, Dimension>;  // 输出图像类型
    using CastFilterType = itk::CastImageFilter<FixedImageType, OutputImageType>;  // 类型转换滤波器类型
    using WriterType = itk::ImageFileWriter<OutputImageType>;  // 写入器类型

    auto writer = WriterType::New();  // 创建写入器
    auto caster = CastFilterType::New();  // 创建类型转换滤波器

    writer->SetFileName(outputImageFile);  // 设置输出文件名
    caster->SetInput(resample->GetOutput());  // 设置转换输入
    writer->SetInput(caster->GetOutput());    // 设置写入输入
    writer->Update();  // 执行写入

    // 生成配准前后的棋盘图
    using CheckerBoardFilterType = itk::CheckerBoardImageFilter<FixedImageType>;  // 棋盘图滤波器类型
    auto checker = CheckerBoardFilterType::New();  // 创建棋盘图滤波器
    checker->SetInput1(fixedImage);  // 设置输入1：固定图像
    checker->SetInput2(resample->GetOutput());  // 设置输入2：重采样后的移动图像
    caster->SetInput(checker->GetOutput());  // 设置转换输入为棋盘图
    writer->SetInput(caster->GetOutput());    // 设置写入输入

    // 配准前的棋盘图（使用恒等变换）
    auto identityTransform = TransformType::New();
    identityTransform->SetIdentity();  // 恒等变换
    resample->SetTransform(identityTransform);  // 使用恒等变换
    writer->SetFileName(checkerboardBeforeFile);  // 设置配准前棋盘图文件名
    writer->Update();  // 写入配准前棋盘图

    // 配准后的棋盘图（使用配准得到的变换）
    resample->SetTransform(registration->GetTransform());  // 使用配准变换
    writer->SetFileName(checkerboardAfterFile);  // 设置配准后棋盘图文件名
    writer->Update();  // 写入配准后棋盘图

    return EXIT_SUCCESS;
}