Compute Curvature Anisotropic Diffusion

Synopsis

Perform anisotropic diffusion on an image.

Results

Input image

Output image

Code

Python

#!/usr/bin/env python

import itk
import argparse

parser = argparse.ArgumentParser(description="Compute Curvature Anisotropic Diffusion.")
parser.add_argument("input_image")
parser.add_argument("output_image")
parser.add_argument("number_of_iterations", type=int)
parser.add_argument("time_step", type=float)
parser.add_argument("conductance", type=float)
args = parser.parse_args()

InputPixelType = itk.F
OutputPixelType = itk.UC
Dimension = 2

InputImageType = itk.Image[InputPixelType, Dimension]
OutputImageType = itk.Image[OutputPixelType, Dimension]

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

FilterType = itk.CurvatureAnisotropicDiffusionImageFilter[
    InputImageType, InputImageType
]
curvatureAnisotropicDiffusionFilter = FilterType.New()

curvatureAnisotropicDiffusionFilter.SetInput(reader.GetOutput())
curvatureAnisotropicDiffusionFilter.SetNumberOfIterations(args.number_of_iterations)
curvatureAnisotropicDiffusionFilter.SetTimeStep(args.time_step)
curvatureAnisotropicDiffusionFilter.SetConductanceParameter(args.conductance)

RescaleFilterType = itk.RescaleIntensityImageFilter[InputImageType, OutputImageType]
rescaler = RescaleFilterType.New()
rescaler.SetInput(curvatureAnisotropicDiffusionFilter.GetOutput())

outputPixelTypeMinimum = itk.NumericTraits[OutputPixelType].min()
outputPixelTypeMaximum = itk.NumericTraits[OutputPixelType].max()

rescaler.SetOutputMinimum(outputPixelTypeMinimum)
rescaler.SetOutputMaximum(outputPixelTypeMaximum)

WriterType = itk.ImageFileWriter[OutputImageType]
writer = WriterType.New()
writer.SetFileName(args.output_image)
writer.SetInput(rescaler.GetOutput())

writer.Update()

C++

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"

int
main(int argc, char * argv[])
{
  if (argc != 6)
  {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0];
    std::cerr << " <InputFileName> <OutputFileName>";
    std::cerr << " <numberOfIterations> <timeStep> <conductance>";
    std::cerr << std::endl;
    return EXIT_FAILURE;
  }

  const char * inputFileName = argv[1];
  const char * outputFileName = argv[2];

  constexpr unsigned int Dimension = 2;

  using InputPixelType = float;
  using InputImageType = itk::Image<InputPixelType, Dimension>;
  using OutputPixelType = unsigned char;
  using OutputImageType = itk::Image<OutputPixelType, Dimension>;

  const int            numberOfIterations = std::stoi(argv[3]);
  const InputPixelType timeStep = std::stod(argv[4]);
  const InputPixelType conductance = std::stod(argv[5]);

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

  using FilterType = itk::CurvatureAnisotropicDiffusionImageFilter<InputImageType, InputImageType>;
  FilterType::Pointer filter = FilterType::New();
  filter->SetInput(reader->GetOutput());
  filter->SetNumberOfIterations(numberOfIterations);
  filter->SetTimeStep(timeStep);
  filter->SetConductanceParameter(conductance);

  using RescaleType = itk::RescaleIntensityImageFilter<InputImageType, OutputImageType>;
  RescaleType::Pointer rescaler = RescaleType::New();
  rescaler->SetInput(filter->GetOutput());
  rescaler->SetOutputMinimum(itk::NumericTraits<OutputPixelType>::min());
  rescaler->SetOutputMaximum(itk::NumericTraits<OutputPixelType>::max());

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

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

  return EXIT_SUCCESS;
}

Classes demonstrated

template<typename TInputImage, typename TOutputImage>
class CurvatureAnisotropicDiffusionImageFilter : public itk::AnisotropicDiffusionImageFilter<TInputImage, TOutputImage>

This filter performs anisotropic diffusion on a scalar itk::Image using the modified curvature diffusion equation (MCDE).

For detailed information on anisotropic diffusion and the MCDE see itkAnisotropicDiffusionFunction and itkCurvatureNDAnisotropicDiffusionFunction.

The default time step for this filter is set to the maximum theoretically stable value: 0.5 / 2^N, where N is the dimensionality of the image. For a 2D image, this means valid time steps are below 0.1250. For a 3D image, valid time steps are below 0.0625.

Inputs and Outputs

The input and output to this filter must be a scalar itk::Image with numerical pixel types (float or double). A user defined type which correctly defines arithmetic operations with floating point accuracy should also give correct results.

Parameters

Please first read all the documentation found in AnisotropicDiffusionImageFilter and AnisotropicDiffusionFunction. Also see CurvatureNDAnisotropicDiffusionFunction.

See

AnisotropicDiffusionImageFilter

See

AnisotropicDiffusionFunction

See

CurvatureNDAnisotropicDiffusionFunction

See itk::CurvatureAnisotropicDiffusionImageFilter for additional documentation.