Problem in getting displacement field

Hi, I used landmark warping code for deforming two sets of landmarks. The output is desirable and I wanted to get displacement field from the code. So I added the following part to the code but the displacement field is zero when I open the file in paraview. Does anyone know the problem? I really appreciate it if you could help me with that!

	// Write the deformation field
	using FieldVectorType = itk::Vector< float, Dimension >;
	using DeformationFieldType = itk::Image< FieldVectorType, Dimension >;
	{
		typedef itk::ImageFileWriter<  DeformationFieldType  > WriterType;
		WriterType::Pointer writer = WriterType::New();
		writer->SetInput(deformer->GetOutput());
		writer->SetFileName(argv[5]);
		writer->Update();
	}
//  Software Guide : BeginLatex
//
//  This example illustrates how to deform an image using a KernelBase spline
//  and two sets of landmarks.
//
//  \index{WarpImageFilter}
//  \index{LandmarkDisplacementFieldSource}
//
//  Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkVector.h"
#include "itkImage.h"
#include "itkLandmarkDisplacementFieldSource.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkWarpImageFilter.h"
// Software Guide : EndCodeSnippet
#include <fstream>
int main(int argc, char * argv[])
{
	if (argc < 4)
	{
		std::cerr << "Missing Parameters " << std::endl;
		std::cerr << "Usage: " << argv[0];
		std::cerr << " landmarksFile fixedImage ";
		std::cerr << "movingImage deformedMovingImage" << std::endl;
		return EXIT_FAILURE;
	}
	constexpr unsigned int Dimension = 3;
	using VectorComponentType = float;
	using VectorType = itk::Vector< VectorComponentType, Dimension >;
	using DisplacementFieldType = itk::Image< VectorType, Dimension >;
	using PixelType = signed short;
	using FixedImageType = itk::Image< PixelType, Dimension >;
	using MovingImageType = itk::Image< PixelType, Dimension >;
	using FixedReaderType = itk::ImageFileReader< FixedImageType  >;
	using MovingReaderType = itk::ImageFileReader< MovingImageType >;
	using MovingWriterType = itk::ImageFileWriter< MovingImageType >;
	FixedReaderType::Pointer fixedReader = FixedReaderType::New();
	fixedReader->SetFileName(argv[2]);
	try
	{
		fixedReader->Update();
	}
	catch (itk::ExceptionObject & excp)
	{
		std::cerr << "Exception thrown " << std::endl;
		std::cerr << excp << std::endl;
		return EXIT_FAILURE;
	}
	MovingReaderType::Pointer movingReader = MovingReaderType::New();
	MovingWriterType::Pointer movingWriter = MovingWriterType::New();
	movingReader->SetFileName(argv[3]);
	movingWriter->SetFileName(argv[4]);
	FixedImageType::ConstPointer fixedImage = fixedReader->GetOutput();
	using DisplacementSourceType = itk::LandmarkDisplacementFieldSource<
		DisplacementFieldType >;
	DisplacementSourceType::Pointer deformer = DisplacementSourceType::New();
	deformer->SetOutputSpacing(fixedImage->GetSpacing());
	deformer->SetOutputOrigin(fixedImage->GetOrigin());
	deformer->SetOutputRegion(fixedImage->GetLargestPossibleRegion());
	deformer->SetOutputDirection(fixedImage->GetDirection());
	//  Create source and target landmarks.
	//
	using LandmarkContainerType = DisplacementSourceType::LandmarkContainer;
	using LandmarkPointType = DisplacementSourceType::LandmarkPointType;
	LandmarkContainerType::Pointer sourceLandmarks = LandmarkContainerType::New();
	LandmarkContainerType::Pointer targetLandmarks = LandmarkContainerType::New();
	LandmarkPointType sourcePoint;
	LandmarkPointType targetPoint;
	std::ifstream pointsFile;
	pointsFile.open(argv[1]);
	unsigned int pointId = 0;
	pointsFile >> sourcePoint;
	pointsFile >> targetPoint;
	while (!pointsFile.fail())
	{
		sourceLandmarks->InsertElement(pointId, sourcePoint);
		targetLandmarks->InsertElement(pointId, targetPoint);
		pointId++;
		pointsFile >> sourcePoint;
		pointsFile >> targetPoint;
	}
	pointsFile.close();
	deformer->SetSourceLandmarks(sourceLandmarks);
	deformer->SetTargetLandmarks(targetLandmarks);
	try
	{
		deformer->UpdateLargestPossibleRegion();
	}
	catch (itk::ExceptionObject & excp)
	{
		std::cerr << "Exception thrown " << std::endl;
		std::cerr << excp << std::endl;
		return EXIT_FAILURE;
	}

	//////////////////////////////////////Added By User////////////////////////////////////////////////////////
	// Write the deformation field
	using FieldVectorType = itk::Vector< float, Dimension >;
	using DeformationFieldType = itk::Image< FieldVectorType, Dimension >;
	{
		typedef itk::ImageFileWriter<  DeformationFieldType  > WriterType;
		WriterType::Pointer writer = WriterType::New();
		writer->SetInput(deformer->GetOutput());
		writer->SetFileName(argv[5]);
		writer->Update();
	}
//////////////////////////////////////////////////

	DisplacementFieldType::ConstPointer displacementField = deformer->GetOutput();
	using FilterType = itk::WarpImageFilter< MovingImageType,
		MovingImageType,
		DisplacementFieldType  >;
	FilterType::Pointer warper = FilterType::New();
	using InterpolatorType = itk::LinearInterpolateImageFunction<
		MovingImageType, double >;
	InterpolatorType::Pointer interpolator = InterpolatorType::New();
	warper->SetInterpolator(interpolator);
	warper->SetOutputSpacing(displacementField->GetSpacing());
	warper->SetOutputOrigin(displacementField->GetOrigin());
	warper->SetDisplacementField(displacementField);
	warper->SetInput(movingReader->GetOutput());
	movingWriter->SetInput(warper->GetOutput());
	try
	{
		movingWriter->Update();
	}
	catch (itk::ExceptionObject & excp)
	{
		std::cerr << "Exception thrown " << std::endl;
		std::cerr << excp << std::endl;
		return EXIT_FAILURE;
	}
	return EXIT_SUCCESS;
}

Hi @basal,

Nothing jumps out as incorrect in the provided code.

A few items to check:

  • Are the landmarks different?
  • If the displacement field file is opened in a text/hex editor, are all the values 0.0?