# Question of the pipeline of MultiResolutionPDEDeformableRegistration

**URL:** https://discourse.itk.org/t/question-of-the-pipeline-of-multiresolutionpdedeformableregistration/2749
**Category:** Beginner Questions
**Created:** [February 20, 2020, 8:50am UTC](https://discourse.itk.org/t/question-of-the-pipeline-of-multiresolutionpdedeformableregistration/2749 "2020-02-20T08:50:00Z")
**Posts on this page:** 3
**Page:** 1

<div class="post-metadata">

### Author: ![zhang-qiang-github](https://discourse.itk.org/user_avatar/discourse.itk.org/zhang-qiang-github/32/1183_2.png) [@zhang-qiang-github](https://discourse.itk.org/u/zhang-qiang-github)
#### Post date: [February 20, 2020, 8:50am UTC](https://discourse.itk.org/t/question-of-the-pipeline-of-multiresolutionpdedeformableregistration/2749/1 "2020-02-20T08:50:00Z")

</div>

Currently, I need to use GPU based Demons on my project. Unfortunately, the GPU based Demons don’t work for multi-resolution framework. Thus, I implement myself multi-resolution code for GPU based Demons. My code is:

```auto
#include <itkImage.h>
#include <itkGPUImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkCastImageFilter.h>
#include <itkHistogramMatchingImageFilter.h>
#include <itkDemonsRegistrationFilter.h>
#include <itkGPUDemonsRegistrationFilter.h>
#include <itkMultiResolutionPDEDeformableRegistration.h>
#include <itkDisplacementFieldTransform.h>
#include <itkBSplineInterpolateImageFunction.h>
#include <ctime>
#include <cstdlib>
#include <itkMultiResolutionPyramidImageFilter.h>
#include <itkResampleImageFilter.h>

clock_t start, end;

//
// The following section of code implements a Command observer
// that will monitor the evolution of the registration process.
// This observer has a layer of intelligence, for deciding what
// MaximumRMS convergence criteria to use at every resolution level.
//
class CommandIterationUpdate : public itk::Command
{
public:
	using Self = CommandIterationUpdate;
	using Superclass = itk::Command;
	using Pointer = itk::SmartPointer<CommandIterationUpdate>;
	itkNewMacro(CommandIterationUpdate);
protected:
	CommandIterationUpdate() = default;

	using InternalImageType = itk::Image< float, 2 >;
	using VectorPixelType = itk::Vector< float, 2 >;
	using DisplacementFieldType = itk::Image< VectorPixelType, 2 >;

	using RegistrationFilterType = itk::DemonsRegistrationFilter<
		InternalImageType,
		InternalImageType,
		DisplacementFieldType>;
	int iter = 0;

public:

	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
	{
		const auto * filter = static_cast<const RegistrationFilterType *>(object);
		if (!(itk::IterationEvent().CheckEvent(&event)))
		{
			return;
		}
		this->iter = this->iter + 1;
		std::cout << "iteration: " << iter << " || metric:" << filter->GetMetric() << std::endl;
	}
};

void main()
{
	start = clock();
	constexpr unsigned int Dimension = 3;
	using PixelType = float;
	//using PixelType = double;
	using InternalPixelType = float;

	using ImageType = itk::Image<PixelType, Dimension >;
	using InternalImageType = itk::Image< InternalPixelType, Dimension >;
	using ImageCasterType = itk::CastImageFilter< ImageType, InternalImageType >;

	using ReaderType = itk::ImageFileReader<ImageType>;
	ReaderType::Pointer targetReader = ReaderType::New();
	targetReader->SetFileName("E:\\COPD\\code\\Registration\\ITK_gpu\\fixImg.mha");
	targetReader->Update();
	ReaderType::Pointer sourceReader = ReaderType::New();
	sourceReader->SetFileName("E:\\COPD\\code\\Registration\\ITK_gpu\\movingImg.mha");
	sourceReader->Update();

	ImageCasterType::Pointer targetImageCaster = ImageCasterType::New();
	ImageCasterType::Pointer sourceImageCaster = ImageCasterType::New();
	targetImageCaster->SetInput(targetReader->GetOutput());
	sourceImageCaster->SetInput(sourceReader->GetOutput());

	using MatchingFilterType = itk::HistogramMatchingImageFilter<
		InternalImageType, InternalImageType >;

	MatchingFilterType::Pointer matcher = MatchingFilterType::New();

	matcher->SetInput(sourceImageCaster->GetOutput());
	matcher->SetReferenceImage(targetImageCaster->GetOutput());
	matcher->SetNumberOfHistogramLevels(1024);
	matcher->SetNumberOfMatchPoints(7);
	matcher->ThresholdAtMeanIntensityOn();

	// setup the deformation field and filter
	using VectorPixelType = itk::Vector< float, Dimension >;
	using CPUDisplacementFieldType = itk::Image< VectorPixelType, Dimension >;
	CPUDisplacementFieldType::Pointer cpuDisplacementField = CPUDisplacementFieldType::New();

	using GPUImageType = itk::GPUImage<PixelType, Dimension>;
	using GPUDisplacementFieldType = itk::GPUImage< VectorPixelType, Dimension >;
	GPUDisplacementFieldType::Pointer gpuDisplacementField = GPUDisplacementFieldType::New();

	using CPUTOGPUImageFilter = itk::CastImageFilter<InternalImageType, GPUImageType>;
	CPUTOGPUImageFilter::Pointer movingImageCPUToGPUFilter = CPUTOGPUImageFilter::New();
	CPUTOGPUImageFilter::Pointer fixImageCPUToGPUFilter = CPUTOGPUImageFilter::New();

	using CPUToGPUDisplacementFieldFilter = itk::CastImageFilter<CPUDisplacementFieldType, GPUDisplacementFieldType>;
	using GPUToCPUDisplacementFieldFilter = itk::CastImageFilter<GPUDisplacementFieldType, CPUDisplacementFieldType>;
	CPUToGPUDisplacementFieldFilter::Pointer cpuToGpuDisplacementFieldFilter = CPUToGPUDisplacementFieldFilter::New();
	GPUToCPUDisplacementFieldFilter::Pointer gpuToCpuDisplacementFieldFilter = GPUToCPUDisplacementFieldFilter::New();

	using FieldExpanderType = itk::ResampleImageFilter<CPUDisplacementFieldType, CPUDisplacementFieldType>;
	FieldExpanderType::Pointer fieldExpander = FieldExpanderType::New();

	using RegistrationFilterType = itk::GPUDemonsRegistrationFilter<
		GPUImageType,
		GPUImageType,
		GPUDisplacementFieldType>;

	using ImagePyramidType = itk::MultiResolutionPyramidImageFilter<InternalImageType, InternalImageType>;
	ImagePyramidType::Pointer movingImagePyramid = ImagePyramidType::New();
	ImagePyramidType::Pointer fixImagePyramid = ImagePyramidType::New();

	unsigned int nIterations[4] = {200, 100, 100, 50};
	int numberOfLevel = 4;
	movingImagePyramid->SetNumberOfLevels(numberOfLevel);
	fixImagePyramid->SetNumberOfLevels(numberOfLevel);

	InternalImageType::Pointer fixImage = targetImageCaster->GetOutput();
	InternalImageType::Pointer movingImage = matcher->GetOutput();

	movingImage->SetRequestedRegionToLargestPossibleRegion();
	fixImage->SetRequestedRegionToLargestPossibleRegion();

	movingImagePyramid->SetInput(movingImage);
	movingImagePyramid->UpdateLargestPossibleRegion();

	fixImagePyramid->SetInput(fixImage);
	fixImagePyramid->UpdateLargestPossibleRegion();

	bool lastShrinkFactorsAllOnes = false;

	for (int i = 0; i < numberOfLevel; i++)
	{
		//std::cout << "-------------------------level: " << i << "------------------------" << std::endl;
		if (i == 0)
		{
			//demonsFilter->SetInitialDisplacementField(nullptr);
		}
		else
		{
			// Resample the field to be the same size as the fixed image
			// at the current level
			fieldExpander->SetInput(cpuDisplacementField);
			auto fi = fixImagePyramid->GetOutput(i);
			auto sizex = fi->GetLargestPossibleRegion().GetSize();
			fieldExpander->SetSize(fi->GetLargestPossibleRegion().GetSize());
			fieldExpander->SetOutputStartIndex(fi->GetLargestPossibleRegion().GetIndex());
			fieldExpander->SetOutputOrigin(fi->GetOrigin());
			fieldExpander->SetOutputSpacing(fi->GetSpacing());
			fieldExpander->SetOutputDirection(fi->GetDirection());
			fieldExpander->UpdateLargestPossibleRegion(); 
			fieldExpander->SetInput(nullptr);
			cpuDisplacementField = fieldExpander->GetOutput();
			cpuDisplacementField->DisconnectPipeline();

			cpuToGpuDisplacementFieldFilter->SetInput(cpuDisplacementField);
			cpuToGpuDisplacementFieldFilter->UpdateLargestPossibleRegion();
			gpuDisplacementField = cpuToGpuDisplacementFieldFilter->GetOutput();
			gpuDisplacementField->DisconnectPipeline();
		}
		movingImageCPUToGPUFilter->SetInput(movingImagePyramid->GetOutput(i));
		movingImageCPUToGPUFilter->GetOutput()->SetRequestedRegionToLargestPossibleRegion();
		movingImageCPUToGPUFilter->UpdateLargestPossibleRegion();
		GPUImageType::Pointer movingGPUImage = GPUImageType::New();
		movingGPUImage = movingImageCPUToGPUFilter->GetOutput();
		movingGPUImage->SetRequestedRegionToLargestPossibleRegion();
		movingGPUImage->DisconnectPipeline();

		fixImageCPUToGPUFilter->SetInput(fixImagePyramid->GetOutput(i));
		fixImageCPUToGPUFilter->GetOutput()->SetRequestedRegionToLargestPossibleRegion();
		fixImageCPUToGPUFilter->UpdateLargestPossibleRegion();
		GPUImageType::Pointer fixGPUImage = GPUImageType::New();
		fixGPUImage = fixImageCPUToGPUFilter->GetOutput();
		fixGPUImage->SetRequestedRegionToLargestPossibleRegion();
		fixGPUImage->DisconnectPipeline();

		RegistrationFilterType::Pointer DemonsFilter = RegistrationFilterType::New();
		CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
		DemonsFilter->AddObserver(itk::IterationEvent(), observer);
		DemonsFilter->SetFixedImage(fixGPUImage);
		DemonsFilter->SetMovingImage(movingGPUImage);
		if (i != 0)
		{
			gpuDisplacementField->DisconnectPipeline();
			DemonsFilter->SetInitialDisplacementField(gpuDisplacementField);
		}
		else
		{
			DemonsFilter->SetInitialDisplacementField(nullptr);
		}

		DemonsFilter->SetNumberOfIterations(nIterations[i]);
		DemonsFilter->SetStandardDeviations(1.0);
		DemonsFilter->UpdateLargestPossibleRegion();

		gpuDisplacementField = DemonsFilter->GetOutput();
		gpuDisplacementField->DisconnectPipeline();
		gpuToCpuDisplacementFieldFilter->SetInput(gpuDisplacementField);
		gpuToCpuDisplacementFieldFilter->Update();
		cpuDisplacementField = gpuToCpuDisplacementFieldFilter->GetOutput();
		cpuDisplacementField->SetRequestedRegionToLargestPossibleRegion();

		movingImagePyramid->GetOutput(i)->ReleaseData();
		fixImagePyramid->GetOutput(i)->ReleaseData();

	}
	fieldExpander->SetInput(nullptr);
	fieldExpander->GetOutput()->ReleaseData();

	end = clock();
	double endtime = (double)(end - start) / CLOCKS_PER_SEC;
	std::cout << "Totol time:" << endtime << std::endl;

	using DisplacementWrite = itk::ImageFileWriter<CPUDisplacementFieldType>;
	DisplacementWrite::Pointer displacementWrite = DisplacementWrite::New();
	displacementWrite->SetInput(cpuDisplacementField);
	displacementWrite->SetFileName("Wrap\\displacement.mha");
	displacementWrite->Update();
}

```

I can run the above code without bug, and the code refers to itkMultiResolutionPDEDeformableRegistration. I have saved the displacement field of 0-, 1-, 2-th level pyramid, and they look like normal. Also, the result of [200, 100, 100, 50] is better than that of 1 level pyramid (50). However, I find the result is worse than the provided exampe (DeformableRegistration16.cxx).

Does the pipeline of my code is wrong?

In addition, I have two questions:

1. I have printed the size of displacement field by:

```auto
std::cout << "size of displacement field is: " << cpuDisplacementField->GetLargestPossibleRegion().GetSize() << std::endl;

```

and the information is:

```auto
size of displacement field is: [348, 512, 512]

```

But, I think the displacement field should be [348, 512, 512, 3]. The last 3 means the displacement of X, Y, Z.

1. Let’s say, the displacement field is [100, 100, 100] in the 1-th level pyramid, and it should be [200, 200, 200] in the 2-th level pyramid. In the itkMultiResolutionPDEDeformableRegistration, the displacement is interploted to the new size by:

```auto
      typename FloatImageType::Pointer fi =
        m_FixedImagePyramid->GetOutput(fixedLevel);
      m_FieldExpander->SetSize(
        fi->GetLargestPossibleRegion().GetSize() );
      m_FieldExpander->SetOutputStartIndex(
        fi->GetLargestPossibleRegion().GetIndex() );
      m_FieldExpander->SetOutputOrigin( fi->GetOrigin() );
      m_FieldExpander->SetOutputSpacing( fi->GetSpacing() );
      m_FieldExpander->SetOutputDirection( fi->GetDirection() );

      m_FieldExpander->UpdateLargestPossibleRegion();
      m_FieldExpander->SetInput(nullptr);
      tempField = m_FieldExpander->GetOutput();
      tempField->DisconnectPipeline();

      m_RegistrationFilter->SetInitialDisplacementField(tempField);

```

However, I think the value of displacement field should be multiplied by 2 for the displacement when up-resampling. The displacement is [x, y, z] when the size of image is [100, 100, 100]. Thus, the displacement should be [2x, 2y, 2z] when the size of image is [200, 200, 200]. Am I right?

* * *

Update:

Finally, I make the performance of myself code similar with DeformableRegistration16.cxx.  
As I said for the second point, the displacement field should be multipied by a factor (for example: 2) when up-resampling. I add part of code as when the cpuDisplacementField is up-resampled:

```
		using DisplacementIterator = itk::ImageRegionIterator<CPUDisplacementFieldType>;
		DisplacementIterator it(cpuDisplacementField, cpuDisplacementField->GetRequestedRegion());
		int k = 0;
		while (!it.IsAtEnd())
		{
			it.Set(it.Get() * 2);
			++it;
			k++;
		}

```

And the performance of above code is more or less better than that of DeformableRegistration16.cxx.

Also, I still feel confuzed that why the displacement field don’t multiply a factor in the MultiResolutionPDEDeformableRegistration. In my code, if I do not multiply a factor for displacement field when up-sampleing, the result would be very worse.

The whole code is updated as:

```auto
#include <itkImage.h>
#include <itkGPUImage.h>
#include <itkImageFileReader.h>
#include <itkImageFileWriter.h>
#include <itkCastImageFilter.h>
#include <itkHistogramMatchingImageFilter.h>
#include <itkDemonsRegistrationFilter.h>
#include <itkGPUDemonsRegistrationFilter.h>
#include <itkMultiResolutionPDEDeformableRegistration.h>
#include <itkDisplacementFieldTransform.h>
#include <itkBSplineInterpolateImageFunction.h>
#include <ctime>
#include <cstdlib>
#include <itkMultiResolutionPyramidImageFilter.h>
#include <itkResampleImageFilter.h>
#include "itkImageIterator.h"
#include "itkImageRegionIterator.h"
clock_t start, end;

unsigned int RmsCounter = 0;

//
// The following section of code implements a Command observer
// that will monitor the evolution of the registration process.
// This observer has a layer of intelligence, for deciding what
// MaximumRMS convergence criteria to use at every resolution level.
//
class CommandIterationUpdate : public itk::Command
{
public:
	using Self = CommandIterationUpdate;
	using Superclass = itk::Command;
	using Pointer = itk::SmartPointer<CommandIterationUpdate>;
	itkNewMacro(CommandIterationUpdate);
protected:
	CommandIterationUpdate() = default;

	using InternalImageType = itk::Image< float, 2 >;
	using VectorPixelType = itk::Vector< float, 2 >;
	using DisplacementFieldType = itk::Image< VectorPixelType, 2 >;

	using RegistrationFilterType = itk::DemonsRegistrationFilter<
		InternalImageType,
		InternalImageType,
		DisplacementFieldType>;
	int iter = 0;

public:

	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
	{
		const auto * filter = static_cast<const RegistrationFilterType *>(object);
		if (!(itk::IterationEvent().CheckEvent(&event)))
		{
			return;
		}
		this->iter = this->iter + 1;
		std::cout << "iteration: " << iter << " || metric:" << filter->GetMetric() << std::endl;
	}
};

//
// The following command observer reports the progress of the registration
// inside a given resolution level.
//
class CommandResolutionLevelUpdate : public itk::Command
{
public:
	using Self = CommandResolutionLevelUpdate;
	using Superclass = itk::Command;
	using Pointer = itk::SmartPointer<Self>;
	itkNewMacro(Self);

protected:
	CommandResolutionLevelUpdate() = default;

public:
	void Execute(itk::Object *caller, const itk::EventObject & event) override
	{
		Execute((const itk::Object *)caller, event);
	}
	void Execute(const itk::Object *, const itk::EventObject &) override
	{
		std::cout << "----------------------------------" << std::endl;
		RmsCounter = RmsCounter + 1;
		std::cout << "----------------------------------" << std::endl;
	}
};

char path[] = "E:\\COPD\\code\\Registration\\ITK_gpu\\";

char * fileCat(char * name1, char * name2)
{
	char name[1000];
	strcpy(name, name1);
	strcat(name, name2);
	return name;
}

void main()
{
	start = clock();
	constexpr unsigned int Dimension = 3;
	using PixelType = float;
	//using PixelType = double;
	using InternalPixelType = float;

	using ImageType = itk::Image<PixelType, Dimension >;
	using InternalImageType = itk::Image< InternalPixelType, Dimension >;
	using ImageCasterType = itk::CastImageFilter< ImageType, InternalImageType >;

	using ReaderType = itk::ImageFileReader<ImageType>;
	ReaderType::Pointer targetReader = ReaderType::New();
	targetReader->SetFileName("E:\\COPD\\code\\Registration\\ITK_gpu\\fixImg.mha");
	targetReader->UpdateLargestPossibleRegion();
	ReaderType::Pointer sourceReader = ReaderType::New();
	sourceReader->SetFileName("E:\\COPD\\code\\Registration\\ITK_gpu\\movingImg.mha");
	sourceReader->UpdateLargestPossibleRegion();

	ImageCasterType::Pointer targetImageCaster = ImageCasterType::New();
	ImageCasterType::Pointer sourceImageCaster = ImageCasterType::New();
	targetImageCaster->SetInput(targetReader->GetOutput());
	targetImageCaster->UpdateLargestPossibleRegion();
	sourceImageCaster->SetInput(sourceReader->GetOutput());
	sourceImageCaster->UpdateLargestPossibleRegion();

	using MatchingFilterType = itk::HistogramMatchingImageFilter<
		InternalImageType, InternalImageType >;

	MatchingFilterType::Pointer matcher = MatchingFilterType::New();

	matcher->SetInput(sourceImageCaster->GetOutput());
	matcher->SetReferenceImage(targetImageCaster->GetOutput());
	matcher->SetNumberOfHistogramLevels(1024);
	matcher->SetNumberOfMatchPoints(7);
	matcher->ThresholdAtMeanIntensityOn();
	matcher->UpdateLargestPossibleRegion();

	// setup the deformation field and filter
	using VectorPixelType = itk::Vector< float, Dimension >;
	using CPUDisplacementFieldType = itk::Image< VectorPixelType, Dimension >;
	CPUDisplacementFieldType::Pointer cpuDisplacementField = CPUDisplacementFieldType::New();

	using GPUImageType = itk::GPUImage<PixelType, Dimension>;
	using GPUDisplacementFieldType = itk::GPUImage< VectorPixelType, Dimension >;
	GPUDisplacementFieldType::Pointer gpuDisplacementField = GPUDisplacementFieldType::New();

	using CPUTOGPUImageFilter = itk::CastImageFilter<InternalImageType, GPUImageType>;
	CPUTOGPUImageFilter::Pointer movingImageCPUToGPUFilter = CPUTOGPUImageFilter::New();
	CPUTOGPUImageFilter::Pointer fixImageCPUToGPUFilter = CPUTOGPUImageFilter::New();

	using CPUToGPUDisplacementFieldFilter = itk::CastImageFilter<CPUDisplacementFieldType, GPUDisplacementFieldType>;
	using GPUToCPUDisplacementFieldFilter = itk::CastImageFilter<GPUDisplacementFieldType, CPUDisplacementFieldType>;
	CPUToGPUDisplacementFieldFilter::Pointer cpuToGpuDisplacementFieldFilter = CPUToGPUDisplacementFieldFilter::New();
	GPUToCPUDisplacementFieldFilter::Pointer gpuToCpuDisplacementFieldFilter = GPUToCPUDisplacementFieldFilter::New();
	
	using DisplacementFieldWrite = itk::ImageFileWriter<CPUDisplacementFieldType>;

	using FieldExpanderType = itk::ResampleImageFilter<CPUDisplacementFieldType, CPUDisplacementFieldType>;
	FieldExpanderType::Pointer fieldExpander = FieldExpanderType::New();

	using RegistrationFilterType = itk::GPUDemonsRegistrationFilter<
		GPUImageType,
		GPUImageType,
		GPUDisplacementFieldType>;
	//RegistrationFilterType::Pointer demonsFilter = RegistrationFilterType::New();
	RegistrationFilterType::Pointer DemonsFilter = RegistrationFilterType::New();
	CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
	DemonsFilter->AddObserver(itk::IterationEvent(), observer);

	using ImagePyramidType = itk::MultiResolutionPyramidImageFilter<InternalImageType, InternalImageType>;
	ImagePyramidType::Pointer movingImagePyramid = ImagePyramidType::New();
	ImagePyramidType::Pointer fixImagePyramid = ImagePyramidType::New();

	unsigned int nIterations[4] = {256, 256, 128, 30};
	int numberOfLevel = 4;
	movingImagePyramid->SetNumberOfLevels(numberOfLevel);
	fixImagePyramid->SetNumberOfLevels(numberOfLevel);

	InternalImageType::Pointer fixImage = targetImageCaster->GetOutput();
	InternalImageType::Pointer movingImage = matcher->GetOutput();

	movingImage->SetRequestedRegionToLargestPossibleRegion();
	fixImage->SetRequestedRegionToLargestPossibleRegion();

	movingImagePyramid->SetInput(movingImage);
	movingImagePyramid->UpdateLargestPossibleRegion();
	auto movingSize = movingImage->GetLargestPossibleRegion().GetSize();
// std::cout << "moving image size" << movingSize << std::endl;
// std::cout << "image size of level 0" << movingImagePyramid->GetOutput(0)->GetLargestPossibleRegion().GetSize() << std::endl;
// std::cout << "image size of level 1" << movingImagePyramid->GetOutput(1)->GetLargestPossibleRegion().GetSize() << std::endl;
// std::cout << "image size of level 2" << movingImagePyramid->GetOutput(2)->GetLargestPossibleRegion().GetSize() << std::endl;
// std::cout << "image size of level 3" << movingImagePyramid->GetOutput(3)->GetLargestPossibleRegion().GetSize() << std::endl;

	fixImagePyramid->SetInput(fixImage);
	fixImagePyramid->UpdateLargestPossibleRegion();

	bool lastShrinkFactorsAllOnes = false;

	for (int i = 0; i < numberOfLevel; i++)
	{
		//std::cout << "-------------------------level: " << i << "------------------------" << std::endl;
		if (i == 0)
		{
			//demonsFilter->SetInitialDisplacementField(nullptr);
		}
		else
		{
// char a[10];
// std::cout << "please input and continue" << std::endl;
// scanf("%s", a);
			// Resample the field to be the same size as the fixed image
			// at the current level
			DisplacementFieldWrite::Pointer displacementBeforResampleWrite = DisplacementFieldWrite::New();
			displacementBeforResampleWrite->SetInput(cpuDisplacementField);
			displacementBeforResampleWrite->SetFileName("E:\\COPD\\code\\Registration\\ITK_gpu\\Wrap\\displacement_beforeResample.mha");
			displacementBeforResampleWrite->Update();

			fieldExpander->SetInput(cpuDisplacementField);
			auto fi = fixImagePyramid->GetOutput(i);
			auto sizex = fi->GetLargestPossibleRegion().GetSize();
			fieldExpander->SetSize(fi->GetLargestPossibleRegion().GetSize());
			fieldExpander->SetOutputStartIndex(fi->GetLargestPossibleRegion().GetIndex());
			fieldExpander->SetOutputOrigin(fi->GetOrigin());
			fieldExpander->SetOutputSpacing(fi->GetSpacing());
			fieldExpander->SetOutputDirection(fi->GetDirection());
			fieldExpander->UpdateLargestPossibleRegion(); // 仅插值，不需要乘以一个系数吗？
			fieldExpander->SetInput(nullptr);
			cpuDisplacementField = fieldExpander->GetOutput();
			cpuDisplacementField->DisconnectPipeline();

			std::cout << "multiply 2 start ..." << std::endl;
			std::cout << "image size is: " << cpuDisplacementField->GetLargestPossibleRegion().GetSize() << std::endl;
			//using DisplacementIterator = itk::ImageIterator<CPUDisplacementFieldType>;
			using DisplacementIterator = itk::ImageRegionIterator<CPUDisplacementFieldType>;
			DisplacementIterator it(cpuDisplacementField, cpuDisplacementField->GetRequestedRegion());
			int k = 0;
			while (!it.IsAtEnd())
			{
				it.Set(it.Get() * 2);
				++it;
				k++;
			}
			std::cout << "multiply 2 finish, element number: "<< k << std::endl;
			cpuToGpuDisplacementFieldFilter->SetInput(cpuDisplacementField);
			cpuToGpuDisplacementFieldFilter->UpdateLargestPossibleRegion();
			gpuDisplacementField = cpuToGpuDisplacementFieldFilter->GetOutput();
			gpuDisplacementField->DisconnectPipeline();
			//demonsFilter->SetInitialDisplacementField(gpuDisplacementField);

			DisplacementFieldWrite::Pointer displacementAfterResampleWrite = DisplacementFieldWrite::New();
			displacementAfterResampleWrite->SetInput(cpuDisplacementField);
			displacementAfterResampleWrite->SetFileName("E:\\COPD\\code\\Registration\\ITK_gpu\\Wrap\\displacement_afterResample.mha");
			displacementAfterResampleWrite->Update();
		}
		movingImageCPUToGPUFilter->SetInput(movingImagePyramid->GetOutput(i));
		movingImageCPUToGPUFilter->GetOutput()->SetRequestedRegionToLargestPossibleRegion();
		movingImageCPUToGPUFilter->UpdateLargestPossibleRegion();
		GPUImageType::Pointer movingGPUImage = GPUImageType::New();
		movingGPUImage = movingImageCPUToGPUFilter->GetOutput();
		movingGPUImage->SetRequestedRegionToLargestPossibleRegion();
		movingGPUImage->DisconnectPipeline();
		auto size1 = movingGPUImage->GetLargestPossibleRegion().GetSize();
		//demonsFilter->SetMovingImage(movingImageCPUToGPUFilter->GetOutput());

		fixImageCPUToGPUFilter->SetInput(fixImagePyramid->GetOutput(i));
		fixImageCPUToGPUFilter->GetOutput()->SetRequestedRegionToLargestPossibleRegion();
		fixImageCPUToGPUFilter->UpdateLargestPossibleRegion();
		GPUImageType::Pointer fixGPUImage = GPUImageType::New();
		fixGPUImage = fixImageCPUToGPUFilter->GetOutput();
		fixGPUImage->SetRequestedRegionToLargestPossibleRegion();
		fixGPUImage->DisconnectPipeline();
		auto size2 = fixGPUImage->GetLargestPossibleRegion().GetSize();
		//demonsFilter->SetFixedImage(fixImageCPUToGPUFilter->GetOutput());

		//demonsFilter->SetNumberOfIterations(nIterations[i]);
		//demonsFilter->GetOutput()->SetRequestedRegionToLargestPossibleRegion();
// using FixImageWrite = itk::ImageFileWriter<InternalImageType>;
// FixImageWrite::Pointer fixImageWrite = FixImageWrite::New();
// auto fixImg = fixImageCPUToGPUFilter->GetOutput();
// fixImageWrite->SetInput(fixGPUImage);
// //fixImageWrite->SetFileName(fileCat(path, "Wrap\\fixImage.mha"));
// fixImageWrite->SetFileName("E:\\COPD\\code\\Registration\\ITK_gpu\\Wrap\\fixImage.mha");
// fixImageWrite->Update();

// using MovingImageWrite = itk::ImageFileWriter<InternalImageType>;
// MovingImageWrite::Pointer movingImageWrite = MovingImageWrite::New();
// movingImageWrite->SetInput(movingGPUImage);
// //movingImageWrite->SetFileName(fileCat(path, "Wrap\\movingImage.mha"));
// movingImageWrite->SetFileName("E:\\COPD\\code\\Registration\\ITK_gpu\\Wrap\\movingImage.mha");
// movingImageWrite->Update();

		DemonsFilter->SetFixedImage(fixGPUImage);
		DemonsFilter->SetMovingImage(movingGPUImage);
		//DemonsFilter->GetOutput()->SetRequestedRegionToLargestPossibleRegion();
		auto sizexxx = gpuDisplacementField->GetLargestPossibleRegion().GetSize();
		//DemonsFilter->SetInitialDisplacementField(nullptr);
		if (i != 0)
		{
			gpuDisplacementField->DisconnectPipeline();
			DemonsFilter->SetInitialDisplacementField(gpuDisplacementField);
		}
		else
		{
			DemonsFilter->SetInitialDisplacementField(nullptr);
		}

		DemonsFilter->SetNumberOfIterations(nIterations[i]);
		DemonsFilter->SetStandardDeviations(1.0);
		DemonsFilter->UpdateLargestPossibleRegion();
		//auto rms = DemonsFilter->GetMaximumRMSError();
		//std::cout << "Maximum RMS Error: " << rms << std::endl;
// auto fixImageSize = DemonsFilter->GetFixedImage()->GetLargestPossibleRegion().GetSize();
// std::cout << "fix image size: " << fixImageSize << std::endl;
// auto movingImageSize = DemonsFilter->GetMovingImage()->GetLargestPossibleRegion().GetSize();
// std::cout << "moving image size: " << movingImageSize << std::endl;
// auto displacementSize = DemonsFilter->GetOutput()->GetLargestPossibleRegion().GetSize();
// std::cout << "final displacement size: " << displacementSize << std::endl;
// auto initialDisplacement = DemonsFilter->GetInitialDisplacementField();
// if (initialDisplacement == nullptr)
// {
// }
// else
// {
// auto initDisplacementSize = DemonsFilter->GetInitialDisplacementField()->GetLargestPossibleRegion().GetSize();
// std::cout << "init displacement size: " << initDisplacementSize << std::endl;
// }
		

		gpuDisplacementField = DemonsFilter->GetOutput();
		auto size11 = gpuDisplacementField->GetLargestPossibleRegion().GetSize();
		//std::cout << "size of displacement: " << size11[0] << " " << size11[1] << " " << size11[2] << std::endl;

		lastShrinkFactorsAllOnes = true;
		for (int idim = 0; idim < Dimension; idim++)
		{
			auto a = fixImagePyramid->GetSchedule()[i];
			if (fixImagePyramid->GetSchedule()[i][idim] > 1)
			{
				lastShrinkFactorsAllOnes = false;
				break;
			}
		}
		//demonsFilter->UpdateLargestPossibleRegion();
		//gpuDisplacementField = demonsFilter->GetOutput();
		//auto size10 = cpuDisplacementField->GetLargestPossibleRegion().GetSize();
		gpuDisplacementField->DisconnectPipeline();
		gpuToCpuDisplacementFieldFilter->SetInput(gpuDisplacementField);
		gpuToCpuDisplacementFieldFilter->Update();
		cpuDisplacementField = gpuToCpuDisplacementFieldFilter->GetOutput();
		cpuDisplacementField->SetRequestedRegionToLargestPossibleRegion();
		auto size9 = cpuDisplacementField->GetLargestPossibleRegion().GetSize();

		DisplacementFieldWrite::Pointer displacementModifyWrite = DisplacementFieldWrite::New();
		displacementModifyWrite->SetInput(cpuDisplacementField);
		displacementModifyWrite->SetFileName("E:\\COPD\\code\\Registration\\ITK_gpu\\Wrap\\displacement_modify.mha");
		displacementModifyWrite->Update();

		movingImagePyramid->GetOutput(i)->ReleaseData();
		fixImagePyramid->GetOutput(i)->ReleaseData();

	}
	if (!lastShrinkFactorsAllOnes)
	{
		// Some of the last shrink factors are not one
		// graft the output of the expander filter to
		// to output of this filter

		// resample the field to the same size as the fixed image
	}
	else
	{

	}
	fieldExpander->SetInput(nullptr);
	fieldExpander->GetOutput()->ReleaseData();
	//demonsFilter->SetInput(nullptr);
	auto o1 = cpuDisplacementField;
	//auto o1 = demonsFilter->GetOutput();
	//demonsFilter->GetOutput()->ReleaseData();
	//auto output = demonsFilter->GetOutput();

	end = clock();
	double endtime = (double)(end - start) / CLOCKS_PER_SEC;
	std::cout << "Totol time:" << endtime << std::endl;

// using DisplacementWrite = itk::ImageFileWriter<CPUDisplacementFieldType>;
// DisplacementWrite::Pointer displacementWrite = DisplacementWrite::New();
// displacementWrite->SetInput(cpuDisplacementField);
// displacementWrite->SetFileName(fileCat(path, "Wrap\\displacement.mha"));
// displacementWrite->Update();
// auto size1 = cpuDisplacementField->GetLargestPossibleRegion().GetSize();
// std::cout <<"size of displacement: "<< size1[0] << " " << size1[1] << " " << size1[2] << std::endl;

	using DisplacementFieldTransformType = itk::DisplacementFieldTransform<InternalPixelType, Dimension>;
	auto displacementTransform = DisplacementFieldTransformType::New();
	displacementTransform->SetDisplacementField(cpuDisplacementField);

	// compute the output (warped) image
	using InterpolatorPrecisionType = double;
	using WarperType = itk::ResampleImageFilter< ImageType, ImageType, InterpolatorPrecisionType, InternalPixelType >;
	using InterpolatorType = itk::BSplineInterpolateImageFunction< ImageType, InterpolatorPrecisionType >;

	WarperType::Pointer warper = WarperType::New();

	InterpolatorType::Pointer interpolator = InterpolatorType::New();

	warper->SetInterpolator(interpolator);
	warper->SetTransform(displacementTransform);
	warper->SetSize(targetReader->GetOutput()->GetLargestPossibleRegion().GetSize());

	//std::cout << "start warp image ..." << std::endl;
	ReaderType::Pointer warperImgReader = ReaderType::New();
	using WriterType = itk::ImageFileWriter<ImageType>;
	WriterType::Pointer warperImgWriter = WriterType::New();

	warperImgReader->SetFileName(fileCat(path, "Wrap\\0.mha"));
	warperImgReader->Update();
	warper->SetInput(warperImgReader->GetOutput());
	warper->Update();
	warperImgWriter->SetFileName(fileCat(path, "Wrap\\wraped_0.mha"));
	warperImgWriter->SetInput(warper->GetOutput());
	warperImgWriter->Update();

	warperImgReader->SetFileName(fileCat(path, "Wrap\\1.mha"));
	warperImgReader->Update();
	warper->SetInput(warperImgReader->GetOutput());
	warper->Update();
	warperImgWriter->SetFileName(fileCat(path, "Wrap\\wraped_1.mha"));
	warperImgWriter->SetInput(warper->GetOutput());
	warperImgWriter->Update();
}

```

Hope ITK can develop the “GPUMultiResolutionPDEDeformableRegistration”.

---

<div class="post-metadata">

### Author: ![dzenanz](https://discourse.itk.org/user_avatar/discourse.itk.org/dzenanz/32/1093_2.png) [@dzenanz](https://discourse.itk.org/u/dzenanz)
#### Post date: [February 20, 2020, 2:26pm UTC](https://discourse.itk.org/t/question-of-the-pipeline-of-multiresolutionpdedeformableregistration/2749/2 "2020-02-20T14:26:48Z")

</div>

> [@zhang-qiang-github](#):
>
> But, I think the displacement field should be [348, 512, 512, 3]. The last 3 means the displacement of X, Y, Z.

If the pixel type is `itk::Vector<3>` (or some other 3-component type), then the image dimensions are [348, 512, 512]. But the memory layout is the same as that of [3, 348, 512, 512] image of scalar type.

> [@zhang-qiang-github](#):
>
> Hope ITK can develop the “GPUMultiResolutionPDEDeformableRegistration”.

ITK community consists of all of us who use ITK, which includes you. If `GPUMultiResolutionPDEDeformableRegistration` is important for you, you can propose its addition to the library by following the [contribution instructions](https://github.com/InsightSoftwareConsortium/ITK/blob/master/CONTRIBUTING.md). If you decided to do so, please add a unit test which demonstrates the use of the class and does regression testing of the result. Possibly be inspired by [this test](https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Registration/PDEDeformable/test/itkMultiResolutionPDEDeformableRegistrationTest.cxx).

---

<div class="post-metadata">

### Author: ![zhang-qiang-github](https://discourse.itk.org/user_avatar/discourse.itk.org/zhang-qiang-github/32/1183_2.png) [@zhang-qiang-github](https://discourse.itk.org/u/zhang-qiang-github)
#### Post date: [February 20, 2020, 3:00pm UTC](https://discourse.itk.org/t/question-of-the-pipeline-of-multiresolutionpdedeformableregistration/2749/3 "2020-02-20T15:00:31Z")

</div>

Why the displacement field don’t multiply a factor when up-resasmpling?

I am happy to implement a GPUMultiResolutionPDEDeformableRegistration. But, to be honest, I am a newbie for ITK. I fear that I cannot implement it. But, I would like to have a try. Hope I can do it.