Resample Segmented Image

Synopsis

Resample (stretch or compress) a label image resulting form segmentation.

For more on why label image interpolation is necessary and how it works, see the associated Insight Journal article.

Results

Input label image.

Input label image

QuickView output.

Resample with label image gaussian interpolation.

QuickView output.

Resample with nearest neighbor interpolation.

Code

Python

#!/usr/bin/env python

import itk
import sys

if len(sys.argv) != 6:
    print(
        f"Usage: {sys.argv[0]} input_image_file spacing_fraction sigma_fraction output_image_file_label_image_interpolator output_image_file_nearest_neighbor_interpolator"
    )
    sys.exit(1)

input_image_file = sys.argv[1]
spacing_fraction = float(sys.argv[2])
sigma_fraction = float(sys.argv[3])
output_image_file_label_image_interpolator = sys.argv[4]
output_image_file_nearest_neighbor_interpolator = sys.argv[5]

input_image = itk.imread(input_image_file)

resize_filter = itk.ResampleImageFilter.New(input_image)

input_spacing = itk.spacing(input_image)
output_spacing = [s * spacing_fraction for s in input_spacing]
resize_filter.SetOutputSpacing(output_spacing)

input_size = itk.size(input_image)
output_size = [
    int(s * input_spacing[dim] / spacing_fraction) for dim, s in enumerate(input_size)
]
resize_filter.SetSize(output_size)

gaussian_interpolator = itk.LabelImageGaussianInterpolateImageFunction.New(input_image)
sigma = [s * sigma_fraction for s in output_spacing]
gaussian_interpolator.SetSigma(sigma)
gaussian_interpolator.SetAlpha(3.0)
resize_filter.SetInterpolator(gaussian_interpolator)

itk.imwrite(resize_filter, output_image_file_label_image_interpolator)

nearest_neighbor_interpolator = itk.NearestNeighborInterpolateImageFunction.New(
    input_image
)
resize_filter.SetInterpolator(nearest_neighbor_interpolator)

itk.imwrite(resize_filter, output_image_file_nearest_neighbor_interpolator)

C++

#include <iostream>
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkLabelImageGaussianInterpolateImageFunction.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkResampleImageFilter.h"

int
main(int argc, char * argv[])
{
  if (argc != 6)
  {
    std::cerr << "Usage: " << argv[0]
              << " inputImageFile spacingFraction sigmaFraction outputImageFileLabelImageInterpolator "
                 "outputImageFileNearestNeighborInterpolator"
              << std::endl;

    return EXIT_FAILURE;
  }
  const char * const inputImageFile = argv[1];
  const double       spacingFraction = std::stod(argv[2]);
  const double       sigmaFraction = std::stod(argv[3]);
  const char * const outputImageFileLabelImageInterpolator = argv[4];
  const char * const outputImageFileNearestNeighborInterpolator = argv[5];

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

  using ReaderType = itk::ImageFileReader<ImageType>;
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(inputImageFile);
  reader->Update();

  using ResampleFilterType = itk::ResampleImageFilter<ImageType, ImageType>;
  ResampleFilterType::Pointer resizeFilter = ResampleFilterType::New();
  resizeFilter->SetInput(reader->GetOutput());

  //     Compute and set the output size
  //
  //     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
  //     or we specify new spacing and compute new height and width
  //     and computations that follows need to be modified a little (as it is

  const ImageType::SpacingType inputSpacing{ reader->GetOutput()->GetSpacing() };
  ImageType::SpacingType       outputSpacing;
  for (unsigned int dim = 0; dim < Dimension; ++dim)
  {
    outputSpacing[dim] = inputSpacing[dim] * spacingFraction;
  }
  resizeFilter->SetOutputSpacing(outputSpacing);

  const ImageType::RegionType inputRegion{ reader->GetOutput()->GetLargestPossibleRegion() };
  const ImageType::SizeType   inputSize{ inputRegion.GetSize() };
  ImageType::SizeType         outputSize;
  for (unsigned int dim = 0; dim < Dimension; ++dim)
  {
    outputSize[dim] = inputSize[dim] * inputSpacing[dim] / spacingFraction;
  }
  resizeFilter->SetSize(outputSize);

  using GaussianInterpolatorType = itk::LabelImageGaussianInterpolateImageFunction<ImageType, double>;
  GaussianInterpolatorType::Pointer   gaussianInterpolator = GaussianInterpolatorType::New();
  GaussianInterpolatorType::ArrayType sigma;
  for (unsigned int dim = 0; dim < Dimension; ++dim)
  {
    sigma[dim] = outputSpacing[dim] * sigmaFraction;
  }
  gaussianInterpolator->SetSigma(sigma);
  gaussianInterpolator->SetAlpha(3.0);
  resizeFilter->SetInterpolator(gaussianInterpolator);

  using WriterType = itk::ImageFileWriter<ImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetInput(resizeFilter->GetOutput());
  writer->SetFileName(outputImageFileLabelImageInterpolator);
  try
  {
    writer->Update();
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  using NearestNeighborInterpolatorType = itk::NearestNeighborInterpolateImageFunction<ImageType, double>;
  NearestNeighborInterpolatorType::Pointer nearestNeighborInterpolator = NearestNeighborInterpolatorType::New();
  resizeFilter->SetInterpolator(nearestNeighborInterpolator);

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

  return EXIT_SUCCESS;
}

Classes demonstrated

template<typename TInputImage, typename TCoordRep = double, typename TPixelCompare = std::less<typename itk::NumericTraits<typename TInputImage::PixelType>::RealType>>
class LabelImageGaussianInterpolateImageFunction : public itk::GaussianInterpolateImageFunction<TInputImage, TCoordRep>

Interpolation function for multi-label images that implicitly smooths each unique value in the image corresponding to each label set element and returns the corresponding label set element with the largest weight.

This filter is an alternative to nearest neighbor interpolation for multi-label images. Given a multi-label image I with label set L, this function returns a label at the non-voxel position I(x), based on the following rule

I(x) = \arg\max_{l \in L} (G_\sigma * I_l)(x)

Where I_l is the l-th binary component of the multilabel image. In other words, each label in the multi-label image is convolved with a Gaussian, and the label for which the response is largest is returned. For sigma=0, this is just nearest neighbor interpolation.

This class defines an N-dimensional Gaussian interpolation function for label using the vnl error function. The two parameters associated with this function are:

  • Sigma - a scalar array of size ImageDimension determining the width of the interpolation function.

  • Alpha - a scalar specifying the cutoff distance over which the function is calculated.

Note

The input image can be of any type, but the number of unique intensity values in the image will determine the amount of memory needed to complete each interpolation.

Author

Paul Yushkevich

Author

Nick Tustison

ITK Sphinx Examples:

See itk::LabelImageGaussianInterpolateImageFunction for additional documentation.