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#

QuickView output. — Resample with label image gaussian interpolation.#

Code#

Python#

#!/usr/bin/env python

import itk
import argparse

parser = argparse.ArgumentParser(description="Resample Segmented Image.")
parser.add_argument("input_image")
parser.add_argument("spacing_fraction", type=float)
parser.add_argument("sigma_fraction", type=float)
parser.add_argument("output_image_file_label_image_interpolator")
parser.add_argument("output_image_file_nearest_neighbor_interpolator")
args = parser.parse_args()

input_image = itk.imread(args.input_image)

resize_filter = itk.ResampleImageFilter.New(input_image)

input_spacing = itk.spacing(input_image)
output_spacing = [s * args.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] / args.spacing_fraction)
    for dim, s in enumerate(input_size)
]
resize_filter.SetSize(output_size)

gaussian_interpolator = itk.LabelImageGaussianInterpolateImageFunction.New(input_image)
sigma = [s * args.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, args.output_image_file_label_image_interpolator)

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

itk.imwrite(resize_filter, args.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>;

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

  using ResampleFilterType = itk::ResampleImageFilter<ImageType, ImageType>;
  auto resizeFilter = ResampleFilterType::New();
  resizeFilter->SetInput(input);

  //     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{ input->GetSpacing() };
  ImageType::SpacingType       outputSpacing;
  for (unsigned int dim = 0; dim < Dimension; ++dim)
  {
    outputSpacing[dim] = inputSpacing[dim] * spacingFraction;
  }
  resizeFilter->SetOutputSpacing(outputSpacing);

  const ImageType::RegionType inputRegion{ input->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>;
  auto                                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);

  try
  {
    itk::WriteImage(resizeFilter->GetOutput(), outputImageFileLabelImageInterpolator, true); // compress
  }
  catch (const itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

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

  try
  {
    itk::WriteImage(resizeFilter->GetOutput(), outputImageFileNearestNeighborInterpolator, true); // compress
  }
  catch (const 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:

All ITK Sphinx Examples

Resample Segmented Image

See itk::LabelImageGaussianInterpolateImageFunction for additional documentation.