Apply GradientRecursiveGaussianImageFilter on Image With Vector type#

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 images (gradient along the X and Y directions). In this example, the input is an image having two values per pixel (vector type). One value is coming from the inputed image (first image); the second is the inverse of the inputed image (second image).

The output of GradientRecursiveGaussianImageFilter will then have four values per pixel. These are in the following order: gradient along X axis (first image), gradient along Y axis (first image), gradient along X axis (second image), gradient along Y axis (second image).

Note that the Python example will only work with ITK 4.8.0 or greater. The ComoposeImagefilter could be used to cast from unsigned char to double type without explicitely using the CastImageFilter; but in Python this possibilty is not allowed. This reduces the number of combinations that need to be wrapped, which helps reducing the size of the ITK binaries.

Results#

Input image

Input image (first image)#

Gradient along X direction (first image)

Gradient along X direction (first image)#

Gradient along Y direction (first image)

Gradient along Y direction (first image)#

Gradient along X direction (second image)

Gradient along X direction (second image)#

Gradient along Y direction (second image)

Gradient along Y direction (second image)#

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.8.0"):
    print("ITK 4.8.0 is required (see example documentation).")
    sys.exit(1)

parser = argparse.ArgumentParser(
    description="Apply Gradient Recursive Gaussian With Vector Input."
)
parser.add_argument("input_image")
parser.add_argument("output_image_1x")
parser.add_argument("output_image_1y")
parser.add_argument("output_image_2x")
parser.add_argument("output_image_2y")
args = parser.parse_args()

filenames = [
    args.output_image_1x,
    args.output_image_1y,
    args.output_image_2x,
    args.output_image_2y,
]

ImageDimension = 2
VectorDimension = 2
CovDimension = 4

PixelType = itk.UC
ImageType = itk.Image[PixelType, ImageDimension]
FloatPixelType = itk.F
FloatImageType = itk.Image[FloatPixelType, ImageDimension]
VecPixelType = itk.Vector[FloatPixelType, VectorDimension]
VecImageType = itk.Image[VecPixelType, ImageDimension]
CovPixelType = itk.CovariantVector[FloatPixelType, CovDimension]
CovImageType = itk.Image[CovPixelType, ImageDimension]

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

# Invert the input image
InvertType = itk.InvertIntensityImageFilter[ImageType, ImageType]
inverter = InvertType.New()
inverter.SetInput(reader.GetOutput())

# Cast the image to double type.
CasterType = itk.CastImageFilter[ImageType, FloatImageType]
caster = CasterType.New()
caster2 = CasterType.New()

# Create an image, were each pixel has 2 values: first value is the value
# coming from the input image, second value is coming from the inverted
# image
ComposeType = itk.ComposeImageFilter[FloatImageType, VecImageType]
composer = ComposeType.New()
caster.SetInput(reader.GetOutput())
composer.SetInput(0, caster.GetOutput())
caster2.SetInput(inverter.GetOutput())
composer.SetInput(1, caster2.GetOutput())

# Apply the gradient filter on the two images, this will return and image
# with 4 values per pixel: two X and Y gradients
FilterType = itk.GradientRecursiveGaussianImageFilter[VecImageType, CovImageType]
gradientfilter = FilterType.New()
gradientfilter.SetInput(composer.GetOutput())

# Set up the filter to select each gradient
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())

for i in range(4):
    indexSelectionFilter.SetIndex(i)
    writer.SetFileName(filenames[i])
    writer.Update()

C++#

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

int
main(int argc, char * argv[])
{
  if (argc != 6)
  {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0];
    std::cerr << " <InputFileName> <OutputFileName1X> <OutputFileName1Y> <OutputFileName2X> <OutputFileName2Y>";
    std::cerr << std::endl;
    return EXIT_FAILURE;
  }

  const char * inputFileName = argv[1];
  const char * outputFileName1X = argv[2];
  const char * outputFileName1Y = argv[3];
  const char * outputFileName2X = argv[4];
  const char * outputFileName2Y = argv[5];

  const char * filenames[4];
  filenames[0] = outputFileName1X;
  filenames[1] = outputFileName1Y;
  filenames[2] = outputFileName2X;
  filenames[3] = outputFileName2Y;

  constexpr unsigned int ImageDimension = 2;
  constexpr unsigned int VectorDimension = 2;
  constexpr unsigned int CovDimension = 4;

  using PixelType = unsigned char;
  using ImageType = itk::Image<PixelType, ImageDimension>;
  using DoublePixelType = double;
  using DoubleImageType = itk::Image<DoublePixelType, ImageDimension>;
  using VecPixelType = itk::Vector<DoublePixelType, VectorDimension>;
  using VecImageType = itk::Image<VecPixelType, ImageDimension>;
  using CovPixelType = itk::CovariantVector<DoublePixelType, CovDimension>;
  using CovImageType = itk::Image<CovPixelType, ImageDimension>;

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

  // Invert the input image
  using InvertType = itk::InvertIntensityImageFilter<ImageType, ImageType>;
  auto inverter = InvertType::New();
  inverter->SetInput(input);

  // Cast the image to double type.
  using CasterType = itk::CastImageFilter<ImageType, DoubleImageType>;
  auto caster = CasterType::New();
  auto caster2 = CasterType::New();

  // Create an image, were each pixel has 2 values: first value is the value
  // coming from the input image, second value is coming from the inverted
  // image
  using ComposeType = itk::ComposeImageFilter<DoubleImageType, VecImageType>;
  auto composer = ComposeType::New();
  caster->SetInput(input);
  composer->SetInput(0, caster->GetOutput());
  caster2->SetInput(inverter->GetOutput());
  composer->SetInput(1, caster2->GetOutput());

  // Apply the gradient filter on the two images, this will return and image
  // with 4 values per pixel: two X and Y gradients
  using FilterType = itk::GradientRecursiveGaussianImageFilter<VecImageType, CovImageType>;
  auto filter = FilterType::New();
  filter->SetInput(composer->GetOutput());

  // Set up the filter to select each gradient
  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 < 4; ++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;
    }
  }

  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.