图像配准的基本过程如下：
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 中还提供了使用多分辨策略的配准框架，弹性配准框架等等，现在刚刚开始研究，其它的慢慢研究。