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;
}