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:
#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:
- I have printed the size of displacement field by:
std::cout << "size of displacement field is: " << cpuDisplacementField->GetLargestPossibleRegion().GetSize() << std::endl;
and the information is:
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.
- 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:
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:
#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”.