Upsample an Image#

Synopsis#

Upsample an image.

Results#

Input image

Input image#

Output image

Output image#

Code#

C++#

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkIdentityTransform.h"
#include "itkBSplineInterpolateImageFunction.h"
#include "itkResampleImageFilter.h"

#include <iostream>

int
main(int argc, char * argv[])
{
  if (argc != 5)
  {
    std::cerr << "Usage: " << std::endl
              << argv[0] << " inputImageFile outputImageFile nNewWidth nNewHeight" << std::endl;

    return EXIT_FAILURE;
  }

  // Typedef's for pixel, image, reader and writer types
  using T_InputPixel = unsigned char;
  using T_Image = itk::Image<T_InputPixel, 2>;

  const auto input = itk::ReadImage<T_Image>(argv[1]);

  // Doesn't work for RGB pixels
  // using T_OutputPixel = unsigned char;
  // using T_InputPixel = itk::CovariantVector<unsigned char, 3>;
  // using T_OutputPixel = itk::CovariantVector<unsigned char, 3>;

  // Typedefs for the different (numerous!) elements of the "resampling"

  // Identity transform.
  // We don't want any transform on our image except rescaling which is not
  // specified by a transform but by the input/output spacing as we will see
  // later.
  // So no transform will be specified.
  using T_Transform = itk::IdentityTransform<double, 2>;

  // If ITK resampler determines there is something to interpolate which is
  // usually the case when upscaling (!) then we must specify the interpolation
  // algorithm. In our case, we want bicubic interpolation. One way to implement
  // it is with a third order b-spline. So the type is specified here and the
  // order will be specified with a method call later on.
  using T_Interpolator = itk::BSplineInterpolateImageFunction<T_Image, double, double>;

  // The resampler type itself.
  using T_ResampleFilter = itk::ResampleImageFilter<T_Image, T_Image>;

  // Prepare the resampler.

  // Instantiate the transform and specify it should be the id transform.
  auto _pTransform = T_Transform::New();
  _pTransform->SetIdentity();

  // Instantiate the b-spline interpolator and set it as the third order
  // for bicubic.
  auto _pInterpolator = T_Interpolator::New();
  _pInterpolator->SetSplineOrder(3);

  // Instantiate the resampler. Wire in the transform and the interpolator.
  auto _pResizeFilter = T_ResampleFilter::New();
  _pResizeFilter->SetTransform(_pTransform);
  _pResizeFilter->SetInterpolator(_pInterpolator);

  // Set the output origin. You may shift the original image "inside" the
  // new image size by specifying something else than 0.0, 0.0 here.

  const double vfOutputOrigin[2] = { 0.0, 0.0 };
  _pResizeFilter->SetOutputOrigin(vfOutputOrigin);

  //     Compute and set the output spacing
  //     Compute the output spacing from input spacing and old and new sizes.
  //
  //     The computation must be so that the following holds:
  //
  //     new width         old x spacing
  //     ----------   =   ---------------
  //     old width         new x spacing
  //
  //
  //     new height         old y spacing
  //    ------------  =   ---------------
  //     old height         new y spacing
  //
  //     So either we specify new height and width and compute new spacings (as
  //     we do here) or we specify new spacing and compute new height and width
  //     and computations that follows need to be modified a little (as it is
  //     done at step 2 there:
  //       https://itk.org/Wiki/ITK/Examples/DICOM/ResampleDICOM)
  //
  unsigned int nNewWidth = std::stoi(argv[3]);
  unsigned int nNewHeight = std::stoi(argv[4]);

  // Fetch original image size.
  const T_Image::RegionType & inputRegion = input->GetLargestPossibleRegion();
  const T_Image::SizeType &   vnInputSize = inputRegion.GetSize();
  unsigned int                nOldWidth = vnInputSize[0];
  unsigned int                nOldHeight = vnInputSize[1];

  // Fetch original image spacing.
  const T_Image::SpacingType & vfInputSpacing = input->GetSpacing();
  // Will be {1.0, 1.0} in the usual
  // case.

  double vfOutputSpacing[2];
  vfOutputSpacing[0] = vfInputSpacing[0] * (double)nOldWidth / (double)nNewWidth;
  vfOutputSpacing[1] = vfInputSpacing[1] * (double)nOldHeight / (double)nNewHeight;

  // Set the output spacing. If you comment out the following line, the original
  // image will be simply put in the upper left corner of the new image without
  // any scaling.
  _pResizeFilter->SetOutputSpacing(vfOutputSpacing);

  // Set the output size as specified on the command line.

  itk::Size<2> vnOutputSize = { { nNewWidth, nNewHeight } };
  _pResizeFilter->SetSize(vnOutputSize);

  // Specify the input.

  _pResizeFilter->SetInput(input);

  // Write the result
  itk::WriteImage(_pResizeFilter->GetOutput(), argv[2]);

  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:

See itk::ResampleImageFilter for additional documentation.
template<typename TImageType, typename TCoordRep = double, typename TCoefficientType = double>
class BSplineInterpolateImageFunction : public itk::InterpolateImageFunction<TImageType, TCoordRep>

Evaluates the B-Spline interpolation of an image. Spline order may be from 0 to 5.

This class defines N-Dimension B-Spline transformation. It is based on:

[1] M. Unser,
    "Splines: A Perfect Fit for Signal and Image Processing,"
    IEEE Signal Processing Magazine, vol. 16, no. 6, pp. 22-38,
    November 1999.
[2] M. Unser, A. Aldroubi and M. Eden,
    "B-Spline Signal Processing: Part I--Theory,"
    IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 821-832,
    February 1993.
[3] M. Unser, A. Aldroubi and M. Eden,
    "B-Spline Signal Processing: Part II--Efficient Design and Applications,"
     IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 834-848,
     February 1993.
And code obtained from bigwww.epfl.ch by Philippe Thevenaz

The B spline coefficients are calculated through the BSplineDecompositionImageFilter

Limitations: Spline order must be between 0 and 5. Spline order must be set before setting the image. Uses mirror boundary conditions. Requires the same order of Spline for each dimension. Spline is determined in all dimensions, cannot selectively pick dimension for calculating spline.

See

BSplineDecompositionImageFilter

See itk::BSplineInterpolateImageFunction for additional documentation.