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 < 4
	typedef   float           VectorComponentType;
#else
	typedef   double          VectorComponentType;
#endif
	typedef   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 < 4
	typedef itk::DeformationFieldSource<DeformationFieldType>  DeformationFieldSourceType;
#else
	typedef itk::LandmarkDisplacementFieldSource<DeformationFieldType>  DeformationFieldSourceType;
#endif
	DeformationFieldSourceType::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 < 4
	typedef itk::DeformationFieldTransform<VectorComponentType, 2>  DeformationFieldTransformType;
#else
	typedef itk::DisplacementFieldTransform<VectorComponentType, 2>  DeformationFieldTransformType;
#endif
	DeformationFieldTransformType::Pointer deformationFieldTransform = DeformationFieldTransformType::New();
	//设置位移场的变换
#if ITK_VERSION_MAJOR < 4
	deformationFieldTransform->SetDeformationField( deformationFieldSource->GetOutput() );
#else
	deformationFieldTransform->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 output
	typedef 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 square
	ImageType::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 field
	typedef 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 square
	ImageType::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 field
	typedef itk::ImageFileWriter<  ImageType  > WriterType;
	WriterType::Pointer writer = WriterType::New();
	writer->SetInput (  image );
	writer->SetFileName( "moving.png" );
	writer->Update();
}

结果图如下所示：