Why are empty spaces observed after performing Deformable registration ?

Beginner Questions
, ,
1

I am trying to register two CT scans of the brain. I have used B-spline transformation code. The output I get after the registration is having some spaces as shown below. How to deal with it? How to apply the inverse transformation?
Please find the code and output image.

image
image.png662Ă—650 32.3 KB 

  #include "itkImageRegistrationMethod.h"
  #include "itkMattesMutualInformationImageToImageMetric.h"

   #include "itkBSplineTransform.h"
   #include "itkLBFGSBOptimizer.h"

   #include "itkImageFileReader.h"
   #include "itkImageFileWriter.h"

    #include "itkResampleImageFilter.h"
    #include "itkCastImageFilter.h"
    #include "itkSquaredDifferenceImageFilter.h"

int main()
{
const    unsigned int    ImageDimension = 3;
typedef  unsigned char     PixelType;
typedef itk::Image< PixelType, ImageDimension >  FixedImageType;
typedef itk::Image< PixelType, ImageDimension >  MovingImageType;

const unsigned int SpaceDimension = ImageDimension;
const unsigned int SplineOrder = 3;
typedef double CoordinateRepType;
typedef itk::BSplineTransform<
	CoordinateRepType,
	SpaceDimension,
	SplineOrder >     TransformType;

typedef itk::LBFGSBOptimizer       OptimizerType;
/*typedef itk::MeanSquaresImageToImageMetric<
	FixedImageType,
	MovingImageType >    MetricType;*/
typedef itk::MattesMutualInformationImageToImageMetric<
	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();
registration->SetMetric(metric);
registration->SetOptimizer(optimizer);
registration->SetInterpolator(interpolator);

TransformType::Pointer  transform = TransformType::New();
registration->SetTransform(transform);

typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
FixedImageReaderType::Pointer  fixedImageReader = FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
fixedImageReader->SetFileName("C:\\VISUAL STUDIO\\projects\\bspline3dregistration1\\bin\\Debug\\images\\subsampled\\Normal_sampled1.mhd");
movingImageReader->SetFileName("C:\\VISUAL STUDIO\\projects\\bspline3dregistration1\\bin\\Debug\\images\\subsampled\\Zone3_sampled.mhd");
FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
registration->SetFixedImage(fixedImage);
registration->SetMovingImage(movingImageReader->GetOutput());
fixedImageReader->Update();
FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
registration->SetFixedImageRegion(fixedRegion);

typedef TransformType::RegionType RegionType;
unsigned int numberOfGridNodes = 4;
TransformType::PhysicalDimensionsType   fixedPhysicalDimensions;
TransformType::MeshSizeType             meshSize;
TransformType::OriginType               fixedOrigin;
for (unsigned int i = 0; i < SpaceDimension; i++)
{
	fixedOrigin[i] = fixedImage->GetOrigin()[i];
	fixedPhysicalDimensions[i] = fixedImage->GetSpacing()[i] *
		static_cast<double>(
			fixedImage->GetLargestPossibleRegion().GetSize()[i] - 1);
}
meshSize.Fill(numberOfGridNodes - SplineOrder);
transform->SetTransformDomainOrigin(fixedOrigin);
transform->SetTransformDomainPhysicalDimensions(
	fixedPhysicalDimensions);
transform->SetTransformDomainMeshSize(meshSize);
transform->SetTransformDomainDirection(fixedImage->GetDirection());
typedef TransformType::ParametersType     ParametersType;
const unsigned int numberOfParameters =
	transform->GetNumberOfParameters();
ParametersType parameters(numberOfParameters);
parameters.Fill(0.0);
transform->SetParameters(parameters);

registration->SetInitialTransformParameters(transform->GetParameters());

std::cout << "Intial Parameters = " << std::endl;
std::cout << transform->GetParameters() << std::endl;

OptimizerType::BoundSelectionType boundSelect(transform->GetNumberOfParameters());
OptimizerType::BoundValueType upperBound(transform->GetNumberOfParameters());
OptimizerType::BoundValueType lowerBound(transform->GetNumberOfParameters());
boundSelect.Fill(0);
upperBound.Fill(0.0);
lowerBound.Fill(0.0);
optimizer->SetBoundSelection(boundSelect);
optimizer->SetUpperBound(upperBound);
optimizer->SetLowerBound(lowerBound);
optimizer->SetCostFunctionConvergenceFactor(1e+12);
optimizer->SetProjectedGradientTolerance(1.0);
optimizer->SetMaximumNumberOfIterations(200);
optimizer->SetMaximumNumberOfEvaluations(100);
optimizer->SetMaximumNumberOfCorrections(5);

metric->SetNumberOfHistogramBins(50);
//FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
const unsigned int numberOfPixels = fixedRegion.GetNumberOfPixels();
metric->ReinitializeSeed(76926294);
metric->SetUseExplicitPDFDerivatives(atoi("False"));
metric->SetUseCachingOfBSplineWeights(atoi("True"));

/*metric->SetNumberOfHistogramBins(50);
const unsigned int numberOfSamples =
	static_cast<unsigned int>(fixedRegion.GetNumberOfPixels() * 20.0 / 100.0);
metric->SetNumberOfSpatialSamples(numberOfSamples);
metric->ReinitializeSeed(76926294);*/

//double samplingPercentage = 0.05;
//registration->SetMetricSamplingPercentage(samplingPercentage);
metric->SetNumberOfSpatialSamples(10000L);

std::cout << std::endl << "Starting Registration" << std::endl;
try
{
	registration->Update();
	std::cout << "Optimizer stop condition = "
		<< registration->GetOptimizer()->GetStopConditionDescription()
		<< std::endl;
}
catch (itk::ExceptionObject & err)
{
	std::cerr << "ExceptionObject caught !" << std::endl;
	std::cerr << err << std::endl;
	return EXIT_FAILURE;
}
OptimizerType::ParametersType finalParameters =
	registration->GetLastTransformParameters();
std::cout << "Last Transform Parameters" << std::endl;
std::cout << finalParameters << std::endl;

transform->SetParameters(finalParameters);

typedef itk::ResampleImageFilter<
	MovingImageType,
	FixedImageType >    ResampleFilterType;
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform(transform);
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  unsigned char  OutputPixelType;
typedef itk::Image< OutputPixelType, ImageDimension > OutputImageType;
typedef itk::CastImageFilter<
	FixedImageType,
	OutputImageType > CastFilterType;
typedef itk::ImageFileWriter< OutputImageType >  WriterType;
WriterType::Pointer      writer = WriterType::New();
CastFilterType::Pointer  caster = CastFilterType::New();
writer->SetFileName("C:\\VISUAL STUDIO\\projects\\bspline3dregistration1\\bin\\Debug\\images\\RegOutput\\Normal_Zone3s.mhd");
caster->SetInput(resample->GetOutput());
writer->SetInput(caster->GetOutput());
writer->Update();

}

2

The “spaces” appear because the images you are registering do not cover the same spatial extent (at least after registration). The missing space has to filled by some value, see SetDefaultPixelValue. The default value is 0 by default, but for CT images more appropriate default is -1000 or -1024.

BSpline transform is not necessarily invertible. See this discussion for possible alternatives.

1 Like