ITK配准：基于DeformationFieldTransform的配准

本文采用DeformationFieldTransform进行图像的配准。

注意在ITKv4中 DeformationFieldTransform被命名为 DisplacementFieldTransform。

1、DeformationFieldTransform提供一个进行局部、稠密和高纬度的变换场，也称为位移场。该位移场存储着向量的 移位，向量的维度记为N。在给定 的输入点上面增加一个位移以实现变换！（Transformation is performed at a given point by adding the displacement at that point to the input point.）

T(x, p), x为所在位置, p 为在x点上面的局部参数（ the local parameter at position x）。 For a 2D example:

x = (x0, x1), p = (p0, p1)

then T(x, p) is defined as:

T(x, p) = (T0, T1) = (x0+p0, x1+p1)

变换期间, 超出边界的输入点返回的是0位移（zero displacement）

2、位移场是通过一张itkImage来定义的,所以需要在使用之前进行设置，通过 SetDisplacementField来实现设置。 该图像和输入图像、输出图像的维度是一致的。此外，需要对位移场图像设置一个类型VectorInterpolateImageFunction 的插值器。默认情况下为VectorLinearInterpolateImageFunction ，也可以通过 SetInterpolator进行重载。

实现过程与之前所撰写的两篇文章为互补。在刚性配准该文章中，并不是就有位移场进行配准的，而基于手动设置配准点的图像配准则是与本文相同的思路；但是在设置滤波器与插值器进行关联时，本文所采用的是刚性配准文章中所用的ResampleFilterType，而不是“基于手动设置配准点的配准”中的WarpImageFilterType所定义的warpImageFilter。

代码如下：

#include "itkImageFileReader.h"#include "itkImageFileWriter.h"#include "itkImage.h"#include "itkVector.h"#if ITK_VERSION_MAJOR < 4#include "itkDeformationFieldTransform.h"#include "itkDeformationFieldSource.h"#else#include "itkVectorLinearInterpolateImageFunction.h"#include "itkDisplacementFieldTransform.h"#include "itkLandmarkDisplacementFieldSource.h"#endif#include "itkResampleImageFilter.h"const     unsigned int   Dimension = 2;typedef   unsigned char  PixelType;typedef   itk::Image< PixelType, Dimension > ImageType;static void CreateFixedImage(ImageType::Pointer image);static void CreateMovingImage(ImageType::Pointer image);int main(int argc, char * argv[]){#if ITK_VERSION_MAJOR < 4typedef   float           VectorComponentType;#elsetypedef   double          VectorComponentType;#endiftypedef   itk::Vector< VectorComponentType, Dimension >    VectorType;typedef   itk::Image< VectorType,  Dimension >   DeformationFieldType;//创建图像ImageType::Pointer fixedImage = ImageType::New();CreateFixedImage(fixedImage);ImageType::Pointer movingImage = ImageType::New();CreateMovingImage(movingImage);//定义位移场的场源类型#if ITK_VERSION_MAJOR < 4typedef itk::DeformationFieldSource<DeformationFieldType>  DeformationFieldSourceType;#elsetypedef itk::LandmarkDisplacementFieldSource<DeformationFieldType>  DeformationFieldSourceType;#endifDeformationFieldSourceType::Pointer deformationFieldSource = DeformationFieldSourceType::New();deformationFieldSource->SetOutputSpacing( fixedImage->GetSpacing() );deformationFieldSource->SetOutputOrigin(  fixedImage->GetOrigin() );deformationFieldSource->SetOutputRegion(  fixedImage->GetLargestPossibleRegion() );deformationFieldSource->SetOutputDirection( fixedImage->GetDirection() );//  Create source and target landmarks.typedef DeformationFieldSourceType::LandmarkContainerPointer   LandmarkContainerPointer;typedef DeformationFieldSourceType::LandmarkContainer          LandmarkContainerType;typedef DeformationFieldSourceType::LandmarkPointType          LandmarkPointType;LandmarkContainerType::Pointer sourceLandmarks = LandmarkContainerType::New();LandmarkContainerType::Pointer targetLandmarks = LandmarkContainerType::New();LandmarkPointType sourcePoint;LandmarkPointType targetPoint;//添加地标点sourcePoint[0] = 40;sourcePoint[1] = 40;targetPoint[0] = 20;targetPoint[1] = 20;sourceLandmarks->InsertElement( 0, sourcePoint );targetLandmarks->InsertElement( 0, targetPoint );sourcePoint[0] = 40;sourcePoint[1] = 60;targetPoint[0] = 20;targetPoint[1] = 80;sourceLandmarks->InsertElement( 1, sourcePoint );targetLandmarks->InsertElement( 1, targetPoint );sourcePoint[0] = 60;sourcePoint[1] = 40;targetPoint[0] = 80;targetPoint[1] = 20;sourceLandmarks->InsertElement( 2, sourcePoint );targetLandmarks->InsertElement( 2, targetPoint );sourcePoint[0] = 60;sourcePoint[1] = 60;targetPoint[0] = 80;targetPoint[1] = 80;sourceLandmarks->InsertElement( 3, sourcePoint );targetLandmarks->InsertElement( 3, targetPoint );//为位移场场源设置地标点deformationFieldSource->SetSourceLandmarks( sourceLandmarks.GetPointer() );deformationFieldSource->SetTargetLandmarks( targetLandmarks.GetPointer() );deformationFieldSource->UpdateLargestPossibleRegion();//将设置好的位移场场源进行保存// Write the deformation field{typedef itk::ImageFileWriter<  DeformationFieldType  > WriterType;WriterType::Pointer writer = WriterType::New();writer->SetInput (  deformationFieldSource->GetOutput() );writer->SetFileName( "deformation.mhd" );writer->Update();}//设置位移场的变换类型#if ITK_VERSION_MAJOR < 4typedef itk::DeformationFieldTransform<VectorComponentType, 2>  DeformationFieldTransformType;#elsetypedef itk::DisplacementFieldTransform<VectorComponentType, 2>  DeformationFieldTransformType;#endifDeformationFieldTransformType::Pointer deformationFieldTransform = DeformationFieldTransformType::New();//设置位移场的变换#if ITK_VERSION_MAJOR < 4deformationFieldTransform->SetDeformationField( deformationFieldSource->GetOutput() );#elsedeformationFieldTransform->SetDisplacementField( deformationFieldSource->GetOutput() );//设置位移场，将场源与其关联#endif  typedef itk::ResampleImageFilter<ImageType, ImageType, VectorComponentType >    ResampleFilterType;ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();resampleFilter->SetInput( movingImage );resampleFilter->SetTransform( deformationFieldTransform );resampleFilter->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );resampleFilter->SetOutputOrigin(  fixedImage->GetOrigin() );resampleFilter->SetOutputSpacing( fixedImage->GetSpacing() );resampleFilter->SetOutputDirection( fixedImage->GetDirection() );resampleFilter->SetDefaultPixelValue( 200 );resampleFilter->GetOutput();// Write the outputtypedef itk::ImageFileWriter<  ImageType  > WriterType;WriterType::Pointer writer = WriterType::New();writer->SetInput (  resampleFilter->GetOutput() );writer->SetFileName( "output.png" );writer->Update();return EXIT_SUCCESS;}void CreateFixedImage(ImageType::Pointer image){// Create a black image with a white squareImageType::IndexType start;start.Fill(0);ImageType::SizeType size;size.Fill(100);ImageType::RegionType region;region.SetSize(size);region.SetIndex(start);image->SetRegions(region);image->Allocate();image->FillBuffer(0);itk::ImageRegionIterator<ImageType> imageIterator(image,region);while(!imageIterator.IsAtEnd()){if(imageIterator.GetIndex()[0] > 40 && imageIterator.GetIndex()[0] < 60 &&imageIterator.GetIndex()[1] > 40 && imageIterator.GetIndex()[1] < 60){imageIterator.Set(255);}++imageIterator;}// Write the deformation fieldtypedef itk::ImageFileWriter<  ImageType  > WriterType;WriterType::Pointer writer = WriterType::New();writer->SetInput (  image );writer->SetFileName( "fixed.png" );writer->Update();}void CreateMovingImage(ImageType::Pointer image){// Create a black image with a white squareImageType::IndexType start;start.Fill(0);ImageType::SizeType size;size.Fill(100);ImageType::RegionType region;region.SetSize(size);region.SetIndex(start);image->SetRegions(region);image->Allocate();image->FillBuffer(0);itk::ImageRegionIterator<ImageType> imageIterator(image,region);while(!imageIterator.IsAtEnd()){if(imageIterator.GetIndex()[0] > 20 && imageIterator.GetIndex()[0] < 80 &&imageIterator.GetIndex()[1] > 20 && imageIterator.GetIndex()[1] < 80){imageIterator.Set(100);}++imageIterator;}// Write the deformation fieldtypedef itk::ImageFileWriter<  ImageType  > WriterType;WriterType::Pointer writer = WriterType::New();writer->SetInput (  image );writer->SetFileName( "moving.png" );writer->Update();}

结果图如下所示：