Process a 2D Slice of a 3D Image

Synopsis

This example illustrates the common task of extracting a 2D slice from a 3D volume. Perform some processing on that slice and then paste it on an output volume of the same size as the volume from the input.

Results

Code

Python

#!/usr/bin/env python

import sys
import itk

if len(sys.argv) != 4:
    print("Usage: " + sys.argv[0] + " input3DImageFile  output3DImageFile  sliceNumber")
    sys.exit(1)
inputFilename = sys.argv[1]
outputFilename = sys.argv[2]

Dimension = 3
PixelType = itk.ctype("short")
ImageType = itk.Image[PixelType, Dimension]

reader = itk.ImageFileReader[ImageType].New()
reader.SetFileName(inputFilename)
reader.Update()
inputImage = reader.GetOutput()

extractFilter = itk.ExtractImageFilter.New(inputImage)
extractFilter.SetDirectionCollapseToSubmatrix()

# set up the extraction region [one slice]
inputRegion = inputImage.GetBufferedRegion()
size = inputRegion.GetSize()
size[2] = 1  # we extract along z direction
start = inputRegion.GetIndex()
sliceNumber = int(sys.argv[3])
start[2] = sliceNumber
desiredRegion = inputRegion
desiredRegion.SetSize(size)
desiredRegion.SetIndex(start)

extractFilter.SetExtractionRegion(desiredRegion)
pasteFilter = itk.PasteImageFilter.New(inputImage)
medianFilter = itk.MedianImageFilter.New(extractFilter)
pasteFilter.SetSourceImage(medianFilter.GetOutput())
pasteFilter.SetDestinationImage(inputImage)
pasteFilter.SetDestinationIndex(start)

indexRadius = size
indexRadius[0] = 1  # radius along x
indexRadius[1] = 1  # radius along y
indexRadius[2] = 0  # radius along z
medianFilter.SetRadius(indexRadius)
medianFilter.UpdateLargestPossibleRegion()
medianImage = medianFilter.GetOutput()
pasteFilter.SetSourceRegion(medianImage.GetBufferedRegion())

itk.imwrite(pasteFilter.GetOutput(), outputFilename)

C++

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkExtractImageFilter.h"
#include "itkPasteImageFilter.h"
#include "itkMedianImageFilter.h"

int
main(int argc, char ** argv)
{
  // Verify the number of parameters in the command line
  if (argc <= 3)
  {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0] << " input3DImageFile  output3DImageFile " << std::endl;
    std::cerr << " sliceNumber " << std::endl;
    return EXIT_FAILURE;
  }

  using PixelType = short;

  using ImageType = itk::Image<PixelType, 3>;

  using ReaderType = itk::ImageFileReader<ImageType>;
  using WriterType = itk::ImageFileWriter<ImageType>;

  // Here we recover the file names from the command line arguments
  const char * inputFilename = argv[1];
  const char * outputFilename = argv[2];

  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(inputFilename);
  try
  {
    reader->Update();
  }
  catch (itk::ExceptionObject & err)
  {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
  }

  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName(outputFilename);

  using ExtractFilterType = itk::ExtractImageFilter<ImageType, ImageType>;
  ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
  extractFilter->SetDirectionCollapseToSubmatrix();

  // set up the extraction region [one slice]
  const ImageType *     inputImage = reader->GetOutput();
  ImageType::RegionType inputRegion = inputImage->GetBufferedRegion();
  ImageType::SizeType   size = inputRegion.GetSize();
  size[2] = 1; // we extract along z direction
  ImageType::IndexType start = inputRegion.GetIndex();
  const unsigned int   sliceNumber = std::stoi(argv[3]);
  start[2] = sliceNumber;
  ImageType::RegionType desiredRegion;
  desiredRegion.SetSize(size);
  desiredRegion.SetIndex(start);

  extractFilter->SetExtractionRegion(desiredRegion);
  using PasteFilterType = itk::PasteImageFilter<ImageType>;
  PasteFilterType::Pointer pasteFilter = PasteFilterType::New();
  using MedianFilterType = itk::MedianImageFilter<ImageType, ImageType>;
  MedianFilterType::Pointer medianFilter = MedianFilterType::New();
  extractFilter->SetInput(inputImage);
  medianFilter->SetInput(extractFilter->GetOutput());
  pasteFilter->SetSourceImage(medianFilter->GetOutput());
  pasteFilter->SetDestinationImage(inputImage);
  pasteFilter->SetDestinationIndex(start);

  ImageType::SizeType indexRadius;
  indexRadius[0] = 1; // radius along x
  indexRadius[1] = 1; // radius along y
  indexRadius[2] = 0; // radius along z
  medianFilter->SetRadius(indexRadius);
  medianFilter->UpdateLargestPossibleRegion();
  const ImageType * medianImage = medianFilter->GetOutput();
  pasteFilter->SetSourceRegion(medianImage->GetBufferedRegion());
  writer->SetInput(pasteFilter->GetOutput());

  try
  {
    writer->Update();
  }
  catch (itk::ExceptionObject & err)
  {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

Classes demonstrated

template<typename TInputImage, typename TOutputImage>
class ExtractImageFilter : public itk::InPlaceImageFilter<TInputImage, TOutputImage>

Decrease the image size by cropping the image to the selected region bounds.

ExtractImageFilter changes the image boundary of an image by removing pixels outside the target region. The target region must be specified.

ExtractImageFilter also collapses dimensions so that the input image may have more dimensions than the output image (i.e. 4-D input image to a 3-D output image). To specify what dimensions to collapse, the ExtractionRegion must be specified. For any dimension dim where ExtractionRegion.Size[dim] = 0, that dimension is collapsed. The index to collapse on is specified by ExtractionRegion.Index[dim]. For example, we have a image 4D = a 4x4x4x4 image, and we want to get a 3D image, 3D = a 4x4x4 image, specified as [x,y,z,2] from 4D (i.e. the 3rd “time” slice from 4D). The ExtractionRegion.Size = [4,4,4,0] and ExtractionRegion.Index = [0,0,0,2].

The number of dimension in ExtractionRegion.Size and Index must = InputImageDimension. The number of non-zero dimensions in ExtractionRegion.Size must = OutputImageDimension.

The output image produced by this filter will have the same origin as the input image, while the ImageRegion of the output image will start at the starting index value provided in the ExtractRegion parameter. If you are looking for a filter that will re-compute the origin of the output image, and provide an output image region whose index is set to zeros, then you may want to use the RegionOfInterestImageFilter. The output spacing is is simply the collapsed version of the input spacing.

Determining the direction of the collapsed output image from an larger dimensional input space is an ill defined problem in general. It is required that the application developer select the desired transformation strategy for collapsing direction cosines. It is REQUIRED that a strategy be explicitly requested (i.e. there is no working default). Direction Collapsing Strategies: 1) DirectionCollapseToUnknown(); This is the default and the filter can not run when this is set. The reason is to explicitly force the application developer to define their desired behavior. 1) DirectionCollapseToIdentity(); Output has identity direction no matter what 2) DirectionCollapseToSubmatrix(); Output direction is the sub-matrix if it is positive definite, else throw an exception.

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

Note

This filter is derived from InPlaceImageFilter. When the input to this filter matched the output requested region, like with streaming filter for input, then setting this filter to run in-place will result in no copying of the bulk pixel data.

See

CropImageFilter

ITK Sphinx Examples:

Subclassed by itk::CropImageFilter< TInputImage, TOutputImage >

See itk::ExtractImageFilter for additional documentation.
template<typename TInputImage, typename TSourceImage = TInputImage, typename TOutputImage = TInputImage>
class PasteImageFilter : public itk::InPlaceImageFilter<TInputImage, TOutputImage>

Paste an image (or a constant value) into another image.

PasteImageFilter allows a region in a destination image to be filled with a source image or a constant pixel value. The SetDestinationIndex() method prescribes where in the destination input to start pasting data from the source input. The SetSourceRegion method prescribes the section of the second image to paste into the first. When a constant pixel value is set, the SourceRegion describes the size of the region filled. If the output requested region does not include the SourceRegion after it has been repositioned to DestinationIndex, then the output will just be a copy of the input.

This filter supports running “InPlace” to efficiently reuse the destination image buffer for the output, removing the need to copy the destination pixels to the output.

When the source image has a lower dimension than the destination image then the DestinationSkipAxes parameter specifies which axes in the destination image are set to 1 when copying the region or filling with a constant.

The two inputs and output image will have the same pixel type.

ITK Sphinx Examples:

See itk::PasteImageFilter for additional documentation.
template<typename TInputImage, typename TOutputImage>
class MedianImageFilter : public itk::BoxImageFilter<TInputImage, TOutputImage>

Applies a median filter to an image.

Computes an image where a given pixel is the median value of the the pixels in a neighborhood about the corresponding input pixel.

A median filter is one of the family of nonlinear filters. It is used to smooth an image without being biased by outliers or shot noise.

This filter requires that the input pixel type provides an operator<() (LessThan Comparable).

See

Image

See

Neighborhood

See

NeighborhoodOperator

See

NeighborhoodIterator

ITK Sphinx Examples:

See itk::MedianImageFilter for additional documentation.