ITK 配准框架示例
图像配准的基本过程如下：
1.指定用于评估配准效果的相似度或误差测度；
2.指定一个变换模型，如刚体变换、仿射变换、弹性变换(elastic)、流体变换或B-样条等；
3.指定插值策略，如最邻近插值(nearst neighbour)、三线性插值(trilinear)、sinc插值等；
4.寻找变换参数以最大化相似性测度。
如下图所示：
配准框架的基本流程如下：
1．输入待配准的两幅图像，参考图 Fixed Image,浮动图 Moving Image.
2．对参考图的指定区域进行几何坐标变换(Transform) 得到新的区域 ;
3．通过插值方法(Image Interpolator)得到浮动图在上一步新区域的坐标；
4．相似性测度模块计算参考图和插值图之间的相似度，是一个关于几何变换参数的函数；
5．相似度函数输入到优化模块中进行最优化计算得到最终变换参数，这个过程一般通过
迭代来实现，即重复2~4步直到取得最大值；
6．整个配准算法模块输出浮动图在最优变换下的插值图像。
配准过程是一个优化问题，配准过程每进行一次迭代，得到一测度值，将该测度值与
我们所设定的值进行比较，如果达到预期的效果则停止迭代，得到最终配准结果。当然，
迭代可能无限制进行，所以我们还需要设置一迭代上限。
下面以一个简单例子对 ITK 中的配准框架进行说明：
ITK 的配准框架由如下几个组件组成：变换组件、相似性测度组件、插值组件和优化组件；各个组件通过一个称为“配准方法”的组件连接到一起，形成一个一个管道结构。
//代码：
- #if defined(_MSC_VER)
- #pragma warning ( disable : 4786 )
- #endif
- //本例演示了使用刚体变换进行 2D 图像配准, 使用的变换为 itk::CenteredRigid2DTransform
- #include "itkImageRegistrationMethod.h"
- #include "itkMeanSquaresImageToImageMetric.h"
- #include "itkLinearInterpolateImageFunction.h"
- #include "itkRegularStepGradientDescentOptimizer.h"
- #include "itkImage.h"
- #include "itkCenteredRigid2DTransform.h"
- #include "itkImageFileReader.h"
- #include "itkImageFileWriter.h"
- #include "itkResampleImageFilter.h"
- #include "itkSubtractImageFilter.h"
- #include "itkRescaleIntensityImageFilter.h"
- #include "itkCommand.h"
- /********************************************************************
- *类 CommandIterationUpdate：创建 Command/Observer 设计模式的简单方式
- *本例中它的主要功能为监视配准过程的执行，每迭代一次，调用一次Execute() 方法，输出相应参数
- *此处涉及三个类: itk::Object, itk::Command，itk::EventObject。
- *Object 是大多数类的基类，该类维护一个着一个指针链表,指针指向事件观察者(event observer).
- *Observer 的角色由 Command 类扮演, Observers 通过一个 Object 注册自己,声明当一特定事件发
- *生时. 配准过程通过一 itk::Optimizer 执行一个迭代过程来控制. 大多数 optimizer 在每次迭代
- *后都会调用一 itk::IterationEvent. 当一个事件发生时, 该对象遍历 registred observers(Command)
- *链表,并检查哪一个对该事件感兴趣, 当找到一个对该事件感兴趣的 observer 时,则调用与其关联的
- *Execute() 方法. 在这方面, Execute() 应该看做是一个 callback,所以它应该 遵循一般 callback
- *的原则,如不应该执行计算机量大的任务,应该只完成一些快速且简单的任务.
- ********************************************************************/
- class CommandIterationUpdate : public itk::Command
- {
- public:
- typedef CommandIterationUpdate Self;
- typedef itk::Command Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- //itkNewMacro: 该宏包装了 New() 方法的所有代码
- itkNewMacro( Self );
- protected:
- CommandIterationUpdate() {};
- public:
- typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
- typedef const OptimizerType * OptimizerPointer;
- void Execute(itk::Object *caller, const itk::EventObject & event)
- {
- Execute( (const itk::Object *)caller, event);
- }
- //object: 指向激活事件的对象；event: 需要被激活的事件
- void Execute(const itk::Object * object, const itk::EventObject & event)
- {
- OptimizerPointer optimizer = dynamic_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( 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 [differenceAfterRegistration] ";
- std::cerr << " [differenceBeforeRegistration] ";
- std::cerr << " [initialStepLength] "<< std::endl;
- return EXIT_FAILURE;
- }
- //1.定义待配准图像类型: 维数, 像素类型
- const unsigned int Dimension = 2;
- typedef unsigned char PixelType;
- typedef itk::Image< PixelType, Dimension > FixedImageType;
- typedef itk::Image< PixelType, Dimension > MovingImageType;
- //2.定义配准框架的基本组件: 变换, 优化, 测度, 插值, 配准组件
- //用于 2D 图像的一个刚体变换,该变换的唯一参数为: 空间坐标的类型
- typedef itk::CenteredRigid2DTransform< double > 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();
- OptimizerType::Pointer optimizer = OptimizerType::New();
- InterpolatorType::Pointer interpolator = InterpolatorType::New();
- RegistrationType::Pointer registration = RegistrationType::New();
- //3.使用配准组件将：变换, 优化, 测度, 插值 四个基本组件连接至一起
- registration->SetMetric( metric );
- registration->SetOptimizer( optimizer );
- registration->SetInterpolator( interpolator );
- TransformType::Pointer transform = TransformType::New();
- registration->SetTransform( transform );
- //4.设置待配准图像及变换区域
- //定义待配准两幅输入图像的 reader
- typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
- typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
- FixedImageReaderType::Pointer fixedImageReader
- = FixedImageReaderType::New();
- MovingImageReaderType::Pointer movingImageReader
- = MovingImageReaderType::New();
- fixedImageReader->SetFileName( argv[1] );
- movingImageReader->SetFileName( argv[2] );
- //本例, 固定图像与浮动图像都来自文件
- //需要 itk::ImageRegistrationMethod 从 file readers 的输出获取输入
- registration->SetFixedImage( fixedImageReader->GetOutput() );
- registration->SetMovingImage( movingImageReader->GetOutput() );
- //更新 readers,确保在初始化变换时图像参数(size,origin,spacing)有效.
- fixedImageReader->Update();
- //可以将配准限制在固定图像的一特定区域,通过将其作为测度(metric)计算的输入.
- //该区域通过方法SetFixedImageRegion() 定义. 通过使用这个特征, 可以避免图像
- //中某些不需要的的对象影响配准输出,或者减少计算时间,本例中使用整幅图像.
- //此区域通过 BufferedRegion 标识.
- registration->SetFixedImageRegion(
- fixedImageReader->GetOutput()->GetBufferedRegion() );
- //更新固定,浮动图像 reader, 使其 origin, size, spacing 有效,
- fixedImageReader->Update();
- movingImageReader->Update();
- //5.设置各组件的参数，主要是变换、优化组件的参数
- //本例使用的为变换为: 先进行一个旋转变换,再进行一平移变换.
- typedef FixedImageType::SpacingType SpacingType;
- typedef FixedImageType::PointType OriginType;
- typedef FixedImageType::RegionType RegionType;
- typedef FixedImageType::SizeType SizeType;
- FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
- const SpacingType fixedSpacing = fixedImage->GetSpacing();
- const OriginType fixedOrigin = fixedImage->GetOrigin();
- const RegionType fixedRegion = fixedImage->GetLargestPossibleRegion();
- const SizeType fixedSize = fixedRegion.GetSize();
- //本例使用固定图像的中心点做为旋转中心
- //使用固定图与浮动图中心之间距离向量作为平移量
- //首先计算固定图像中心点
- TransformType::InputPointType centerFixed;
- centerFixed[0] = fixedOrigin[0] + fixedSpacing[0] * fixedSize[0] / 2.0;
- centerFixed[1] = fixedOrigin[1] + fixedSpacing[1] * fixedSize[1] / 2.0;
- MovingImageType::Pointer movingImage = movingImageReader->GetOutput();
- const SpacingType movingSpacing = movingImage->GetSpacing();
- const OriginType movingOrigin = movingImage->GetOrigin();
- const RegionType movingRegion = movingImage->GetLargestPossibleRegion();
- const SizeType movingSize = movingRegion.GetSize();
- //然后计算浮动图像中心点
- TransformType::InputPointType centerMoving;
- centerMoving[0] = movingOrigin[0] + movingSpacing[0] * movingSize[0] / 2.0;
- centerMoving[1] = movingOrigin[1] + movingSpacing[1] * movingSize[1] / 2.0;
- //设置旋转中心
- transform->SetCenter( centerFixed );
- //设置平移量
- transform->SetTranslation( centerMoving - centerFixed );
- //设置初始旋转角度
- transform->SetAngle( 0.0 );
- //然后将当前的变换参数作为初始参数, 传递给配准过程
- registration->SetInitialTransformParameters( transform->GetParameters() );
- //设置优化组件参数.
- //注意: 旋转和变换的"单位刻度" 是不同的,旋转用弧度度量, 平移以毫米度量
- typedef OptimizerType::ScalesType OptimizerScalesType;
- OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
- const double translationScale = 1.0 / 1000.0;
- optimizerScales[0] = 1.0;
- optimizerScales[1] = translationScale;
- optimizerScales[2] = translationScale;
- optimizerScales[3] = translationScale;
- optimizerScales[4] = translationScale;
- optimizer->SetScales( optimizerScales );
- //本例使用的优化器是 itk::RegularStepGradientDescentOptimizer,
- //梯度下降法的一种,下面定义优化参数:
- //relaxation fator, initial step length, minimal step length, number of iterations
- double initialStepLength = 0.1; //初始步长
- if( argc > 6 )
- {
- initialStepLength = atof( argv[6] );
- }
- optimizer->SetRelaxationFactor( 0.6 ); //松驰因子
- optimizer->SetMaximumStepLength( initialStepLength ); //最大步长
- //最小步长和迭代最大次数, 用来限制优化过程迭代的执行
- optimizer->SetMinimumStepLength( 0.001 ); //最小步长
- optimizer->SetNumberOfIterations( 200 ); //最大迭代次数,以防止无限次迭代
- /*
- *应该注意: 配准过程是一个优化的过程, 优化过程每迭代一次,
- *就与相似性测度进行一次比较
- *当相似性测度值达到我们设定的预期值,或者达到设定的迭代上限时
- *就停止迭代,得到最终结果.
- */
- //6.实例化一 Command/Observer 对象, 监视配准过程的执行, 并触发配准过程的执行.
- //实例化一个 Command/Observer 对象, 将其注册为 optimizer 的观察者
- CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
- optimizer->AddObserver( itk::IterationEvent(), observer );
- try
- {
- registration->StartRegistration(); //触发配准过程的执行
- }
- catch( itk::ExceptionObject & err )
- {
- std::cerr << "ExceptionObject caught !" << std::endl;
- std::cerr << err << std::endl;
- return EXIT_FAILURE;
- }
- //取得配准执行完毕的最终参数
- OptimizerType::ParametersType finalParameters =
- registration->GetLastTransformParameters();
- const double finalAngle = finalParameters[0];
- const double finalRotationCenterX = finalParameters[1];
- const double finalRotationCenterY = finalParameters[2];
- const double finalTranslationX = finalParameters[3];
- const double finalTranslationY = finalParameters[4];
- const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
- const double bestValue = optimizer->GetValue();
- const double finalAngleInDegrees = finalAngle * 180.0 / vnl_math::pi;
- std::cout << "Result = " << std::endl;
- std::cout << " Angle (radians) = " << finalAngle << std::endl;
- std::cout << " Angle (degrees) = " << finalAngleInDegrees << std::endl;
- std::cout << " Center X = " << finalRotationCenterX << std::endl;
- std::cout << " Center Y = " << finalRotationCenterY << std::endl;
- std::cout << " Translation X = " << finalTranslationX << std::endl;
- std::cout << " Translation Y = " << finalTranslationY << std::endl;
- std::cout << " Iterations = " << numberOfIterations << std::endl;
- std::cout << " Metric value = " << bestValue << std::endl;
- //7.重采样, 将变换后的浮动图像映射到固定图像空间中,保存配准结果
- //一般情况,最后一步为: 利用最终的变换结果将浮动图像映射到固定图像空间
- //可通过 itk::ResampleImageFilter 完成
- typedef itk::ResampleImageFilter< MovingImageType,
- FixedImageType > ResampleFilterType;
- TransformType::Pointer finalTransform = TransformType::New();
- finalTransform->SetParameters( finalParameters );
- finalTransform->SetFixedParameters( transform->GetFixedParameters() );
- //设置重采样过滤器的相应参数
- ResampleFilterType::Pointer resample = ResampleFilterType::New();
- resample->SetTransform( finalTransform );
- //将浮动图连接至重采样过滤器的输入端
- resample->SetInput( movingImageReader->GetOutput() );
- resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
- resample->SetOutputOrigin( fixedImage->GetOrigin() );
- resample->SetOutputSpacing( fixedImage->GetSpacing() );
- resample->SetOutputDirection( fixedImage->GetDirection() );
- //设置默认灰度值,用来 "突出" 显示映射后在浮动图像之外的区域
- resample->SetDefaultPixelValue( 100 );
- //保存配准后图像
- typedef itk::ImageFileWriter< FixedImageType > WriterFixedType;
- WriterFixedType::Pointer writer = WriterFixedType::New();
- writer->SetFileName( argv[3] );
- writer->SetInput( resample->GetOutput() ); //重采样过滤器的输出
- try
- {
- writer->Update();
- }
- catch( itk::ExceptionObject & excp )
- {
- std::cerr << "ExceptionObject while writing the resampled image !" << std::endl;
- std::cerr << excp << std::endl;
- return EXIT_FAILURE;
- }
- //通过 itk::SubtractImageFilter 比较"经过变换的浮动图像"与"固定图像"之间的差异
- typedef itk::Image< float, Dimension > DifferenceImageType;
- typedef itk::SubtractImageFilter< FixedImageType,
- FixedImageType, DifferenceImageType > DifferenceFilterType;
- DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
- typedef unsigned char OutputPixelType;
- typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
- //由于两幅图像之间的差异可能对应非常小的亮度值
- //使用 itk::RescaleIntensityImageFilter 重新调节亮度值,使其更加明显
- typedef itk::RescaleIntensityImageFilter< DifferenceImageType,
- OutputImageType > RescalerType;
- RescalerType::Pointer intensityRescaler = RescalerType::New();
- intensityRescaler->SetOutputMinimum( 0 );
- intensityRescaler->SetOutputMaximum( 255 );
- difference->SetInput1( fixedImageReader->GetOutput() );
- difference->SetInput2( resample->GetOutput() );
- resample->SetDefaultPixelValue( 1 );
- intensityRescaler->SetInput( difference->GetOutput() );
- typedef itk::ImageFileWriter< OutputImageType > WriterType;
- WriterType::Pointer writer2 = WriterType::New();
- writer2->SetInput( intensityRescaler->GetOutput() );
- try
- {
- if( argc > 4 )
- {
- //保存配准后, 浮动图与固定图之间差异的对比
- writer2->SetFileName( argv[4] );
- writer2->Update();
- }
- //保存配准前, 浮动图与固定图之间差异的对比, 这里仍然需要重采样
- //因为浮动图与固定图不一定有相同的 origin, size, spacing.
- //可以过一个恒等变换(Identity Transform)完成
- TransformType::Pointer identityTransform = TransformType::New();
- identityTransform->SetIdentity();
- resample->SetTransform( identityTransform );
- if( argc > 5 )
- {
- writer2->SetFileName( argv[5] );
- writer2->Update();
- }
- }
- catch( itk::ExceptionObject & excp )
- {
- std::cerr << "Error while writing difference images" << std::endl;
- std::cerr << excp << std::endl;
- return EXIT_FAILURE;
- }
- return EXIT_SUCCESS;
- }
除去细节，可以看出，使用 ITK 的配准框架进行医学图像配准主要遵循以下步骤：
1.定义待配准图像类型: 维数, 像素类型
2.定义配准框架的基本组件: 变换, 优化, 测度, 插值, 配准组件
3.使用配准组件将：变换, 优化, 测度, 插值 四个基本组件连接至一起
4.设置待配准图像及变换区域，有需要时进行适当处理
5.设置各组件的参数，变换、优化、测度，不同方法设置不同
6.实例化一个 Command/Observer 对象, 监视配准过程的执行, 并触发配准过程的执行.
7.重采样, 将变换后的浮动图像映射到固定图像空间中，保存配准结果
各个配准组件组成的管道结构如下所示：
Command/Observer 观察者模式的执行过程如下图所示：
本例的 CMakeLists.txt 内容为：
CMAKE_MINIMUM_REQUIRED(VERSION 2.4)
PROJECT(ImageRegistration5)
FIND_PACKAGE(ITK)
IF(ITK_FOUND)
INCLUDE(${ITK_USE_FILE})
ELSE(ITK_FOUND)
MESSAGE(FATAL_ERROR
"ITK not found. Please set ITK_DIR.")
ENDIF(ITK_FOUND)
ADD_EXECUTABLE(ImageRegistration5 ImageRegistration5.cxx )
TARGET_LINK_LIBRARIES(ImageRegistration5 ITKIO ITKNumerics)
程序所使用的两幅示例图像，固定图像，浮动图像分别如下：
程序执行后的过程及结果为：
Debug>ImageRegistration5.exe BrainProtonDensitySliceBorder20.png BrainProtonDensitySliceRotated10.png output.png after.png befor.png
0 2098.33 [0.0252835, 110.5, 128.5, 0.0912944, -0.0320326]
1 1766.73 [0.0618533, 110.501, 128.502, 0.177712, -0.0665167]
2 1380.97 [0.104328, 110.501, 128.5, 0.135417, -0.0669017]
3 877.601 [0.122006, 110.501, 128.494, 0.0784302, -0.0690638]
4 647.837 [0.140307, 110.5, 128.487, 0.0219416, -0.0639346]
5 410.321 [0.162961, 110.497, 128.479, -0.0313332, -0.0501843]
6 153.093 [0.189839, 110.49, 128.483, -0.0120674, -0.00086769]
7 115.912 [0.177587, 110.487, 128.494, 0.0427617, 0.0171322]
8 61.0641 [0.177564, 110.489, 128.488, 0.0095898, 0.00463817]
9 60.3813 [0.17645, 110.495, 128.489, 0.0217426, -0.0286424]
10 61.411 [0.177373, 110.491, 128.488, 0.0108792, -0.0103758]
11 60.4302 [0.177628, 110.488, 128.488, 0.00999985, 0.0108734]
12 60.5575 [0.177367, 110.49, 128.488, 0.00930347, -0.00186576]
13 60.3963 [0.177692, 110.489, 128.488, 0.011102, 0.00556995]
14 60.3989 [0.177429, 110.489, 128.488, 0.00980777, 0.0011698]
15 60.3822 [0.177691, 110.489, 128.488, 0.0109154, 0.00368068]
16 60.3942 [0.177487, 110.489, 128.488, 0.0103087, 0.00215544]
17 60.3781 [0.177817, 110.489, 128.488, 0.00869253, 0.00202342]
18 60.4265 [0.17767, 110.489, 128.488, 0.00967173, 0.00209358]
19 60.3867 [0.177458, 110.489, 128.488, 0.0106296, 0.00194103]
//配准完毕后的最终参数
Result =
Angle (radians) = 0.177458
Angle (degrees) = 10.1676
Center X = 110.489
Center Y = 128.488
Translation X = 0.0106296
Translation Y = 0.00194103
Iterations = 21
Metric value = 60.3809
配准后的图像，配准前两幅图像的差异，配准后的差异：
可以看到，配准后的两幅图像仍然有稍许没有对齐的部分。可以通过更复杂的测度，优化方法进行配准。
ITK 提供了大量不同的几何变换、测度方法、插值算法以及优化算法，主要的不同在于不同组件的参数设置。比较复杂的方法参数设置是件困难的事情，通常需要大量实验才能确定合适的参数，而且不同的实际问题参数设置也不尽相同，比较重要的参数通常需要在其它参数设置好以后再进行微调，才能得到合适的值。ITK 中还提供了使用多分辨策略的配准框架，弹性配准框架等等，现在刚刚开始研究，其它的慢慢研究。
转自zhangcunli的博文
1099
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
