Apply GradientRecursiveGaussianImageFilter

Synopsis

Computes the gradient of an image by convolution with the first derivative of a Gaussian.

The output of the GradientRecursiveGaussianImageFilter is composed of two imagescomponents (gradients along the X and Y directions). In this example, we store these components using an image with CovariantVector pixel type. Covariant vectors are the natural representation for gradients since they are the equivalent of normals to iso-values manifolds.

This example shows also how to extract the X and Y gradients as images, and how to compute the magnitude of the CovariantVector image.

Note that the covariant vector types were only added to the VectorIndexSelectionCastImageFilter Python wrapping in ITK 4.7.

Results

Input image

Input image

Gradient along X direction

Gradient along X direction

Gradient along Y direction

Gradient along Y direction

Output magnitude

The image of the magnitude of the vector containing the X and Y values.

Code

C++

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImage.h"
#include "itkGradientRecursiveGaussianImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkVectorIndexSelectionCastImageFilter.h"
#include "itkVectorMagnitudeImageFilter.h"

int
main(int argc, char * argv[])
{
  if (argc != 5)
  {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0];
    std::cerr << " <InputFileName> <OutputFileNameX> <OutputFileNameY> <OutputFileNameMagnitude>";
    std::cerr << std::endl;
    return EXIT_FAILURE;
  }

  const char * inputFileName = argv[1];
  const char * outputFileNameX = argv[2];
  const char * outputFileNameY = argv[3];
  const char * outputFileNameMagnitude = argv[4];

  const char * filenames[2];
  filenames[0] = outputFileNameX;
  filenames[1] = outputFileNameY;

  constexpr unsigned int Dimension = 2;

  // Input and output are png files, use unsigned char
  using PixelType = unsigned char;
  using ImageType = itk::Image<PixelType, Dimension>;
  // Double type for GradientRecursiveGaussianImageFilter
  using DoublePixelType = double;
  using DoubleImageType = itk::Image<DoublePixelType, Dimension>;
  // The output of GradientRecursiveGaussianImageFilter
  // are images of the gradient along X and Y, so the type of
  // the output is a covariant vector of dimension 2 (X, Y)
  using CovPixelType = itk::CovariantVector<DoublePixelType, Dimension>;
  using CovImageType = itk::Image<CovPixelType, Dimension>;

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

  using FilterType = itk::GradientRecursiveGaussianImageFilter<ImageType, CovImageType>;
  FilterType::Pointer filter = FilterType::New();
  filter->SetInput(reader->GetOutput());

  // Allows to select the X or Y output images
  using IndexSelectionType = itk::VectorIndexSelectionCastImageFilter<CovImageType, DoubleImageType>;
  IndexSelectionType::Pointer indexSelectionFilter = IndexSelectionType::New();
  indexSelectionFilter->SetInput(filter->GetOutput());

  // Rescale for png output
  using RescalerType = itk::RescaleIntensityImageFilter<DoubleImageType, ImageType>;
  RescalerType::Pointer rescaler = RescalerType::New();
  rescaler->SetOutputMinimum(itk::NumericTraits<PixelType>::min());
  rescaler->SetOutputMaximum(itk::NumericTraits<PixelType>::max());
  rescaler->SetInput(indexSelectionFilter->GetOutput());

  using WriterType = itk::ImageFileWriter<ImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetInput(rescaler->GetOutput());

  // Write the X and Y images
  for (int i = 0; i < 2; ++i)
  {

    writer->SetFileName(filenames[i]);
    indexSelectionFilter->SetIndex(i);

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

  // Compute the magnitude of the vector and output the image
  using MagnitudeType = itk::VectorMagnitudeImageFilter<CovImageType, DoubleImageType>;
  MagnitudeType::Pointer magnitudeFilter = MagnitudeType::New();
  magnitudeFilter->SetInput(filter->GetOutput());

  // Rescale for png output
  rescaler->SetInput(magnitudeFilter->GetOutput());

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

  return EXIT_SUCCESS;
}

Python

#!/usr/bin/env python
import itk

from distutils.version import StrictVersion as VS
if VS(itk.Version.GetITKVersion()) < VS("4.7.0"):
    print("ITK 4.7.0 is required (see example documentation).")
    sys.exit(1)

if len(sys.argv) != 5:
    print("Usage: " + sys.argv[0] +
          " [InputFileName] [OutputFileNameX] [OutputFileNameY]" +
          " [OutputFileNameMagnitude]")
    sys.exit(1)

inputFileName = sys.argv[1]
outputFileNameX = sys.argv[2]
outputFileNameY = sys.argv[3]
outputFileNameMagnitude = sys.argv[4]

Dimension = 2

filenames = []
filenames.append(outputFileNameX)
filenames.append(outputFileNameY)

# Input and output are png files, use unsigned char
PixelType = itk.UC
ImageType = itk.Image[PixelType, Dimension]
# Float type for GradientRecursiveGaussianImageFilter
FloatPixelType = itk.F
FloatImageType = itk.Image[FloatPixelType, Dimension]
# The output of GradientRecursiveGaussianImageFilter
# are images of the gradient along X and Y, so the type of
# the output is a covariant vector of dimension 2 (X, Y)
CovPixelType = itk.CovariantVector[FloatPixelType, Dimension]
CovImageType = itk.Image[CovPixelType, Dimension]

ReaderType = itk.ImageFileReader[ImageType]
reader = ReaderType.New()
reader.SetFileName(inputFileName)

FilterType = itk.GradientRecursiveGaussianImageFilter[ImageType, CovImageType]
gradientFilter = FilterType.New()
gradientFilter.SetInput(reader.GetOutput())

# Allows to select the X or Y output images
IndexSelectionType = itk.VectorIndexSelectionCastImageFilter[
    CovImageType, FloatImageType]
indexSelectionFilter = IndexSelectionType.New()
indexSelectionFilter.SetInput(gradientFilter.GetOutput())

# Rescale for png output
RescalerType = itk.RescaleIntensityImageFilter[FloatImageType, ImageType]
rescaler = RescalerType.New()
rescaler.SetOutputMinimum(itk.NumericTraits[PixelType].min())
rescaler.SetOutputMaximum(itk.NumericTraits[PixelType].max())
rescaler.SetInput(indexSelectionFilter.GetOutput())

WriterType = itk.ImageFileWriter[ImageType]
writer = WriterType.New()
writer.SetInput(rescaler.GetOutput())

# Write the X and Y images
for i in range(2):
    writer.SetFileName(filenames[i])
    indexSelectionFilter.SetIndex(i)
    writer.Update()

# Compute the magnitude of the vector and output the image
MagnitudeType = itk.VectorMagnitudeImageFilter[
    CovImageType, FloatImageType]
magnitudeFilter = MagnitudeType.New()
magnitudeFilter.SetInput(gradientFilter.GetOutput())

# Rescale for png output
rescaler.SetInput(magnitudeFilter.GetOutput())

writer.SetFileName(outputFileNameMagnitude)
writer.SetInput(rescaler.GetOutput())
writer.Update()

Classes demonstrated