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#

Python#

#!/usr/bin/env python

import sys
import itk
import argparse

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)

parser = argparse.ArgumentParser(description="Apply Gradient Recursive Gaussian.")
parser.add_argument("input_image")
parser.add_argument("output_image_x")
parser.add_argument("output_image_y")
parser.add_argument("output_image_magnitude")
args = parser.parse_args()

Dimension = 2

filenames = []
filenames.append(args.output_image_x)
filenames.append(args.output_image_y)

# 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(args.input_image)

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(args.output_image_magnitude)
writer.SetInput(rescaler.GetOutput())
writer.Update()

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

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

  using FilterType = itk::GradientRecursiveGaussianImageFilter<ImageType, CovImageType>;
  auto filter = FilterType::New();
  filter->SetInput(input);

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

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

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

    try
    {
      itk::WriteImage(rescaler->GetOutput(), filenames[i]);
    }
    catch (const 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>;
  auto magnitudeFilter = MagnitudeType::New();
  magnitudeFilter->SetInput(filter->GetOutput());

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

  try
  {
    itk::WriteImage(rescaler->GetOutput(), outputFileNameMagnitude);
  }
  catch (const itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

Classes demonstrated#

template<typename TInputImage, typename TOutputImage = Image<CovariantVector<typename NumericTraits<typename TInputImage::PixelType>::RealType, TInputImage::ImageDimension>, TInputImage::ImageDimension>>
class GradientRecursiveGaussianImageFilter : public itk::ImageToImageFilter<TInputImage, TOutputImage>

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

This filter is implemented using the recursive gaussian filters.

This filter supports both scalar and vector pixel types within the input image, including VectorImage type.

ITK Sphinx Examples:

See itk::GradientRecursiveGaussianImageFilter for additional documentation.