# Why are empty spaces observed after performing Deformable registration ?

**URL:** <https://discourse.itk.org/t/why-are-empty-spaces-observed-after-performing-deformable-registration/1810>\
**Category:** Beginner Questions\
**Tags:** simpleitk, registration, python\
**Created:** [April 24, 2019, 9:49am UTC](https://discourse.itk.org/t/why-are-empty-spaces-observed-after-performing-deformable-registration/1810 "2019-04-24T09:49:02Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![EE18D504\_SINDHURA\_C](https://discourse.itk.org/letter_avatar_proxy/v4/letter/e/b3f665/32.png) [@EE18D504\_SINDHURA\_C](https://discourse.itk.org/u/EE18D504_SINDHURA_C)\
**Post date:** [April 24, 2019, 9:49am UTC](https://discourse.itk.org/t/why-are-empty-spaces-observed-after-performing-deformable-registration/1810/1 "2019-04-24T09:49:02Z")

</div>

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](https://discourse.itk.org/uploads/default/original/1X/075019288979fce21b534a0905084c8d2b776054.png)

```
  #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();

```

}

---

<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:** [April 29, 2019, 4:23pm UTC](https://discourse.itk.org/t/why-are-empty-spaces-observed-after-performing-deformable-registration/1810/2 "2019-04-29T16:23:26Z")

</div>

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](https://itk.org/Doxygen/html/classitk_1_1ResampleImageFilter.html#a227397f3786a0ec4d3107cb73033b4be). 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](https://discourse.itk.org/t/inverse-of-bspline-transform/496) for possible alternatives.