Resample an Image#

Synopsis#

Resample an image.

Results#

Code#

Python#

#!/usr/bin/env python

import sys
import itk

if len(sys.argv) != 4:
    print("Usage: " + sys.argv[0] + " <input_image> <output_image> <scale>")
    sys.exit(1)

input_file_name = sys.argv[1]
output_file_name = sys.argv[2]
scale = float(sys.argv[3])

input_image = itk.imread(input_file_name)
input_size = itk.size(input_image)
input_spacing = itk.spacing(input_image)
input_origin = itk.origin(input_image)
Dimension = input_image.GetImageDimension()

# We will scale the objects in the image by the factor `scale`; that is they
# will be shrunk (scale < 1.0) or enlarged (scale > 1.0).  However, the number
# of pixels for each dimension of the output image will equal the corresponding
# number of pixels in the input image, with cropping or padding as necessary.
# Furthermore, the physical distance between adjacent pixels will be the same
# in the input and the output images.  In contrast, if you want to change the
# resolution of the image without changing the represented physical size of the
# objects in the image, omit the transform and instead supply:
#
# output_size = [int(input_size[d] * scale) for d in range(Dimension)]
# output_spacing = [input_spacing[d] / scale for d in range(Dimension)]
# output_origin = [input_origin[d] + 0.5 * (output_spacing[d] - input_spacing[d])
#                  for d in range(Dimension)]

output_size = input_size
output_spacing = input_spacing
output_origin = input_origin
scale_transform = itk.ScaleTransform[itk.D, Dimension].New()
scale_transform_parameters = scale_transform.GetParameters()
for i in range(len(scale_transform_parameters)):
    scale_transform_parameters[i] = scale
scale_transform_center = [float(int(s / 2)) for s in input_size]
scale_transform.SetParameters(scale_transform_parameters)
scale_transform.SetCenter(scale_transform_center)

interpolator = itk.LinearInterpolateImageFunction.New(input_image)

resampled = itk.resample_image_filter(
    input_image,
    transform=scale_transform,
    interpolator=interpolator,
    size=output_size,
    output_spacing=output_spacing,
    output_origin=output_origin,
)

itk.imwrite(resampled, output_file_name)

C++#

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkResampleImageFilter.h"
#include "itkScaleTransform.h"

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

  const char * inputFileName = argv[1];
  const char * outputFileName = argv[2];
  const float  scale = std::stod(argv[3]);

  constexpr unsigned int Dimension = 2;

  using PixelType = unsigned char;
  using ImageType = itk::Image<PixelType, Dimension>;
  using ScalarType = double;
  using IndexValueType = typename itk::Index<Dimension>::IndexValueType;

  const typename ImageType::Pointer     inputImage = itk::ReadImage<ImageType>(inputFileName);
  const typename ImageType::RegionType  inputRegion = inputImage->GetLargestPossibleRegion();
  const typename ImageType::SizeType    inputSize = inputRegion.GetSize();
  const typename ImageType::SpacingType inputSpacing = inputImage->GetSpacing();
  const typename ImageType::PointType   inputOrigin = inputImage->GetOrigin();

  /*
   * We will scale the objects in the image by the factor `scale`; that is they
   * will be shrunk (scale < 1.0) or enlarged (scale > 1.0).  However, the
   * number of pixels for each dimension of the output image will equal the
   * corresponding number of pixels in the input image, with padding (if
   * shrunk) or cropping (if enlarged) as necessary.  Furthermore, the physical
   * distance between adjacent pixels will be the same in the input and the
   * output images.  In contrast, if you want to change the resolution of the
   * image without changing the represented physical size of the objects in the
   * image, omit the transform and instead use:
   *
   *   outputSize[d] = inputSize[d] * scale;
   *   outputSpacing[d] = inputSpacing[d] / scale;
   *   outputOrigin[d] = inputOrigin[d] + 0.5 * (outputSpacing[d] - inputSpacing[d]);
   *
   * in the loop over dimensions.
   */

  typename ImageType::SizeType    outputSize = inputSize;
  typename ImageType::SpacingType outputSpacing = inputSpacing;
  typename ImageType::PointType   outputOrigin = inputOrigin;

  using ScaleTransformType = itk::ScaleTransform<ScalarType, Dimension>;
  auto scaleTransform = ScaleTransformType::New();

  typename ScaleTransformType::ParametersType scaleTransformParameters = scaleTransform->GetParameters();
  itk::Point<ScalarType, Dimension>           scaleTransformCenter;
  for (unsigned int d = 0; d < Dimension; ++d)
  {
    scaleTransformParameters[d] = scale;
    scaleTransformCenter[d] = static_cast<ScalarType>(static_cast<IndexValueType>(inputSize[d] / 2));
  }
  scaleTransform->SetParameters(scaleTransformParameters);
  scaleTransform->SetCenter(scaleTransformCenter);

  using LinearInterpolatorType = itk::LinearInterpolateImageFunction<ImageType, ScalarType>;
  auto interpolator = LinearInterpolatorType::New();

  using ResampleFilterType = itk::ResampleImageFilter<ImageType, ImageType>;
  auto resampleFilter = ResampleFilterType::New();

  resampleFilter->SetInput(inputImage);
  resampleFilter->SetTransform(scaleTransform);
  resampleFilter->SetInterpolator(interpolator);
  resampleFilter->SetSize(outputSize);
  resampleFilter->SetOutputSpacing(outputSpacing);
  resampleFilter->SetOutputOrigin(outputOrigin);

  try
  {
    itk::WriteImage(resampleFilter->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 TInterpolatorPrecisionType = double, typename TTransformPrecisionType = TInterpolatorPrecisionType> class ResampleImageFilter : public itk::ImageToImageFilter<TInputImage, TOutputImage>

Resample an image via a coordinate transform.

ResampleImageFilter resamples an existing image through some coordinate transform, interpolating via some image function. The class is templated over the types of the input and output images.

Note that the choice of interpolator function can be important. This function is set via SetInterpolator(). The default is LinearInterpolateImageFunction<InputImageType, TInterpolatorPrecisionType>, which is reasonable for ordinary medical images. However, some synthetic images have pixels drawn from a finite prescribed set. An example would be a mask indicating the segmentation of a brain into a small number of tissue types. For such an image, one does not want to interpolate between different pixel values, and so NearestNeighborInterpolateImageFunction< InputImageType, TCoordRep > would be a better choice.

If an sample is taken from outside the image domain, the default behavior is to use a default pixel value. If different behavior is desired, an extrapolator function can be set with SetExtrapolator().

Output information (spacing, size and direction) for the output image should be set. This information has the normal defaults of unit spacing, zero origin and identity direction. Optionally, the output information can be obtained from a reference image. If the reference image is provided and UseReferenceImage is On, then the spacing, origin and direction of the reference image will be used.

Since this filter produces an image which is a different size than its input, it needs to override several of the methods defined in ProcessObject in order to properly manage the pipeline execution model. In particular, this filter overrides ProcessObject::GenerateInputRequestedRegion() and ProcessObject::GenerateOutputInformation().

This filter is implemented as a multithreaded filter. It provides a DynamicThreadedGenerateData() method for its implementation.

Warning

For multithreading, the TransformPoint method of the user-designated coordinate transform must be threadsafe.

ITK Sphinx Examples: