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

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

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

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

  // Cast the image to double type.
  using CasterType = itk::CastImageFilter<ImageType, DoubleImageType>;
  CasterType::Pointer caster = CasterType::New();
  CasterType::Pointer 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>;
  ComposeType::Pointer 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
  using FilterType = itk::GradientRecursiveGaussianImageFilter<VecImageType, CovImageType>;
  FilterType::Pointer filter = FilterType::New();
  filter->SetInput(composer->GetOutput());

  // Set up the filter to select each gradient
  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 < 4; ++i)
  {
    indexSelectionFilter->SetIndex(i);
    writer->SetFileName(filenames[i]);

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

if len(sys.argv) != 6:
    print(
        "Usage: " + sys.argv[0] +
        " [InputFileName] [OutputFileName1X] [OutputFileName1Y]" +
        " [OutputFileName2X] [OutputFileName2Y]")
    sys.exit(1)

inputFileName = sys.argv[1]
outputFileName1X = sys.argv[2]
outputFileName1Y = sys.argv[3]
outputFileName2X = sys.argv[4]
outputFileName2Y = sys.argv[5]

filenames = [
    outputFileName1X,
    outputFileName1Y,
    outputFileName2X,
    outputFileName2Y]

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(inputFileName)

# 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()

Classes demonstrated