Warp An Image Using A Deformation Field

Synopsis

Use WarpImageFilter and a deformation field for resampling an image.

Results

Input image

Input image

Output image

Output image

Code

C++

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkVector.h"
#include "itkWarpImageFilter.h"

int
main(int argc, char * argv[])
{
  if (argc < 4)
  {
    std::cerr << "Usage: \n" << argv[0] << "<InputFileName> <DisplacementFieldFileName> <OutputFileName>" << std::endl;
    return EXIT_FAILURE;
  }

  const char * inputFileName = argv[1];
  const char * displacementFieldFileName = argv[2];
  const char * outputFileName = argv[3];

  constexpr unsigned int Dimension = 2;

  using VectorComponentType = float;
  using VectorPixelType = itk::Vector<VectorComponentType, Dimension>;

  using DisplacementFieldType = itk::Image<VectorPixelType, Dimension>;

  using PixelType = unsigned char;
  using ImageType = itk::Image<PixelType, Dimension>;

  using ReaderType = itk::ImageFileReader<ImageType>;
  using WriterType = itk::ImageFileWriter<ImageType>;

  using FieldReaderType = itk::ImageFileReader<DisplacementFieldType>;

  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(inputFileName);

  FieldReaderType::Pointer fieldReader = FieldReaderType::New();
  fieldReader->SetFileName(displacementFieldFileName);
  fieldReader->Update();

  DisplacementFieldType::ConstPointer deformationField = fieldReader->GetOutput();

  using FilterType = itk::WarpImageFilter<ImageType, ImageType, DisplacementFieldType>;

  FilterType::Pointer filter = FilterType::New();

  using InterpolatorType = itk::LinearInterpolateImageFunction<ImageType, double>;

  InterpolatorType::Pointer interpolator = InterpolatorType::New();

  filter->SetInterpolator(interpolator);

  filter->SetOutputSpacing(deformationField->GetSpacing());
  filter->SetOutputOrigin(deformationField->GetOrigin());
  filter->SetOutputDirection(deformationField->GetDirection());

  filter->SetDisplacementField(deformationField);

  filter->SetInput(reader->GetOutput());

  WriterType::Pointer writer = WriterType::New();
  writer->SetInput(filter->GetOutput());
  writer->SetFileName(outputFileName);

  try
  {
    writer->Update();
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

Python

#!/usr/bin/env python

import sys
import itk

if len(sys.argv) != 4:
    print('Usage: ' + sys.argv[0] +
          ' <InputFileName> <DisplacementFieldFileName> <OutputFileName>')
    sys.exit(1)

inputFileName = sys.argv[1]
displacementFieldFileName = sys.argv[2]
outputFileName = sys.argv[3]

Dimension = 2

VectorComponentType = itk.F
VectorPixelType = itk.Vector[VectorComponentType, Dimension]

DisplacementFieldType = itk.Image[VectorPixelType, Dimension]

PixelType = itk.UC
ImageType = itk.Image[PixelType, Dimension]

reader = itk.ImageFileReader[ImageType].New()
reader.SetFileName(inputFileName)

fieldReader = itk.ImageFileReader[DisplacementFieldType].New()
fieldReader.SetFileName(displacementFieldFileName)
fieldReader.Update()

deformationField = fieldReader.GetOutput()

warpFilter = \
        itk.WarpImageFilter[ImageType, ImageType, DisplacementFieldType].New()

interpolator = itk.LinearInterpolateImageFunction[ImageType, itk.D].New()

warpFilter.SetInterpolator(interpolator)

warpFilter.SetOutputSpacing(deformationField.GetSpacing())
warpFilter.SetOutputOrigin(deformationField.GetOrigin())
warpFilter.SetOutputDirection(deformationField.GetDirection())

warpFilter.SetDisplacementField(deformationField)

warpFilter.SetInput(reader.GetOutput())

writer = itk.ImageFileWriter[ImageType].New()
writer.SetInput(warpFilter.GetOutput())
writer.SetFileName(outputFileName)

writer.Update()

Classes demonstrated

template<typename TInputImage, typename TOutputImage, typename TDisplacementField>
class WarpImageFilter : public itk::ImageToImageFilter<TInputImage, TOutputImage>

Warps an image using an input displacement field.

WarpImageFilter warps an existing image with respect to a given displacement field.

A displacement field is represented as a image whose pixel type is some vector type with at least N elements, where N is the dimension of the input image. The vector type must support element access via operator [].

The output image is produced by inverse mapping: the output pixels are mapped back onto the input image. This scheme avoids the creation of any holes and overlaps in the output image.

Each vector in the displacement field represent the distance between a geometric point in the input space and a point in the output space such that:

p_{in} = p_{out} + d

Typically the mapped position does not correspond to an integer pixel position in the input image. Interpolation via an image function is used to compute values at non-integer positions. The default interpolation typed used is the LinearInterpolateImageFunction. The user can specify a particular interpolation function via SetInterpolator(). Note that the input interpolator must derive from base class InterpolateImageFunction.

Position mapped to outside of the input image buffer are assigned a edge padding value.

The LargetPossibleRegion for the output is inherited from the input displacement field. The output image spacing, origin and orientation may be set via SetOutputSpacing, SetOutputOrigin and SetOutputDirection. The default are respectively a vector of 1’s, a vector of 0’s and an identity matrix.

This class is templated over the type of the input image, the type of the output image and the type of the displacement field.

The input image is set via SetInput. The input displacement field is set via SetDisplacementField.

This filter is implemented as a multithreaded filter.

Warning

This filter assumes that the input type, output type and displacement field type all have the same number of dimensions.

ITK Sphinx Examples:

See itk::WarpImageFilter for additional documentation.