Warp an Image Using a Deformation Field#

Synopsis#

Use WarpImageFilter and a deformation field for resampling an image.

Results#

Code#

Python#

#!/usr/bin/env python

import itk
import argparse

parser = argparse.ArgumentParser(description="Wrap An Image Using A Deformation Field.")
parser.add_argument("input_image")
parser.add_argument("displacement_field")
parser.add_argument("output_image")
args = parser.parse_args()

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(args.input_image)

fieldReader = itk.ImageFileReader[DisplacementFieldType].New()
fieldReader.SetFileName(args.displacement_field)
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(args.output_image)

writer.Update()

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

  const auto input = itk::ReadImage<ImageType>(inputFileName);

  const auto deformationField = itk::ReadImage<DisplacementFieldType>(displacementFieldFileName);

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

  auto filter = FilterType::New();

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

  auto interpolator = InterpolatorType::New();

  filter->SetInterpolator(interpolator);

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

  filter->SetDisplacementField(deformationField);

  filter->SetInput(input);

  try
  {
    itk::WriteImage(filter->GetOutput(), outputFileName);
  }
  catch (const itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

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: