Image Read Extract Filter Insert Write

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

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;
    }

  typedef short PixelType;

  typedef itk::Image< PixelType, 3 > ImageType;

  typedef itk::ImageFileReader< ImageType > ReaderType;
  typedef itk::ImageFileWriter< ImageType > WriterType;

  // 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 );

  typedef itk::ExtractImageFilter< ImageType, ImageType >
    ExtractFilterType;
  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 = atoi( argv[3] );
  start[2] = sliceNumber;
  ImageType::RegionType desiredRegion;
  desiredRegion.SetSize(  size  );
  desiredRegion.SetIndex( start );

  extractFilter->SetExtractionRegion( desiredRegion );
  typedef itk::PasteImageFilter< ImageType > PasteFilterType;
  PasteFilterType::Pointer pasteFilter = PasteFilterType::New();
  typedef itk::MedianImageFilter< ImageType,
                                  ImageType > MedianFilterType;
  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;
}

Python

#!/usr/bin/env python
import itk

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

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

# Here we recover the file names from the command line arguments
inputFilename = sys.argv[1]
outputFilename = sys.argv[2]
reader = itk.ImageFileReader[ImageType].New()
reader.SetFileName(inputFilename)
reader.Update()

extractFilter = itk.ExtractImageFilter[ImageType, ImageType].New()
extractFilter.SetInput(reader.GetOutput())
extractFilter.SetDirectionCollapseToSubmatrix()

# set up the extraction region [one slice]
inputImage = reader.GetOutput()
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[ImageType].New()
medianFilter = itk.MedianImageFilter[ImageType, ImageType].New()
extractFilter.SetInput(inputImage)
medianFilter.SetInput(extractFilter.GetOutput())
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())

writer = itk.ImageFileWriter[ImageType].New()
writer.SetFileName(outputFilename)
writer.SetInput(pasteFilter.GetOutput())
writer.Update()

Classes demonstrated

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

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 ThreadedGenerateData() method for its implementation.

Note
This filter is derived from InPlaceImageFilter. When the input to this filter matched the output requirested 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
Wiki Examples:

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

Paste an image into another image.

PasteImageFilter allows you to take a section of one image and paste into another image. The SetDestinationIndex() method prescribes where in the first input to start pasting data from the second input. The SetSourceRegion method prescribes the section of the second image to paste into the first. 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.

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

Wiki Examples:

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

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

Neighborhood

NeighborhoodOperator

NeighborhoodIterator

Wiki Examples:

See itk::MedianImageFilter for additional documentation.