Segment With Geodesic Active Contour Level Set

Synopsis

Segments structures in images based on a user supplied edge potential map.

This example is the same as the one in the ITK software guide. It was re-organized and does not contain any comments about the executed filters. Please refer to the guide for more details.

Results

Input image

Input image

Evolution of the segmentation depending on the number of iterations (10 to 500).

Used parameters:

  • Left ventricle: 81 114 5.0 1.0 -0.5 3.0 2.0

  • Right ventricle: 99 114 5.0 1.0 -0.5 3.0 2.0

  • White matter: 56 92 5.0 1.0 -0.3 2.0 10.0

  • Gray matter: 40 90 5.0 .5 -0.3 2.0 10.0

Code

C++

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkGeodesicActiveContourLevelSetImageFilter.h"
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkFastMarchingImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"

int
main(int argc, char * argv[])
{
  if (argc != 11)
  {
    std::cerr << "Usage: " << argv[0];
    std::cerr << " <InputFileName>  <OutputFileName>";
    std::cerr << " <seedX> <seedY> <InitialDistance>";
    std::cerr << " <Sigma> <SigmoidAlpha> <SigmoidBeta>";
    std::cerr << " <PropagationScaling> <NumberOfIterations>" << std::endl;
    return EXIT_FAILURE;
  }

  const char * inputFileName = argv[1];
  const char * outputFileName = argv[2];
  const int    seedPosX = std::stoi(argv[3]);
  const int    seedPosY = std::stoi(argv[4]);

  const double initialDistance = std::stod(argv[5]);
  const double sigma = std::stod(argv[6]);
  const double alpha = std::stod(argv[7]);
  const double beta = std::stod(argv[8]);
  const double propagationScaling = std::stod(argv[9]);
  const double numberOfIterations = std::stoi(argv[10]);
  const double seedValue = -initialDistance;

  constexpr unsigned int Dimension = 2;

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

  using ReaderType = itk::ImageFileReader<InputImageType>;
  using WriterType = itk::ImageFileWriter<OutputImageType>;

  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(inputFileName);

  using SmoothingFilterType = itk::CurvatureAnisotropicDiffusionImageFilter<InputImageType, InputImageType>;
  SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
  smoothing->SetTimeStep(0.125);
  smoothing->SetNumberOfIterations(5);
  smoothing->SetConductanceParameter(9.0);
  smoothing->SetInput(reader->GetOutput());

  using GradientFilterType = itk::GradientMagnitudeRecursiveGaussianImageFilter<InputImageType, InputImageType>;
  GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();
  gradientMagnitude->SetSigma(sigma);
  gradientMagnitude->SetInput(smoothing->GetOutput());

  using SigmoidFilterType = itk::SigmoidImageFilter<InputImageType, InputImageType>;
  SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
  sigmoid->SetOutputMinimum(0.0);
  sigmoid->SetOutputMaximum(1.0);
  sigmoid->SetAlpha(alpha);
  sigmoid->SetBeta(beta);
  sigmoid->SetInput(gradientMagnitude->GetOutput());

  using FastMarchingFilterType = itk::FastMarchingImageFilter<InputImageType, InputImageType>;
  FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();

  using GeodesicActiveContourFilterType = itk::GeodesicActiveContourLevelSetImageFilter<InputImageType, InputImageType>;
  GeodesicActiveContourFilterType::Pointer geodesicActiveContour = GeodesicActiveContourFilterType::New();
  geodesicActiveContour->SetPropagationScaling(propagationScaling);
  geodesicActiveContour->SetCurvatureScaling(1.0);
  geodesicActiveContour->SetAdvectionScaling(1.0);
  geodesicActiveContour->SetMaximumRMSError(0.02);
  geodesicActiveContour->SetNumberOfIterations(numberOfIterations);
  geodesicActiveContour->SetInput(fastMarching->GetOutput());
  geodesicActiveContour->SetFeatureImage(sigmoid->GetOutput());

  using ThresholdingFilterType = itk::BinaryThresholdImageFilter<InputImageType, OutputImageType>;
  ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
  thresholder->SetLowerThreshold(-1000.0);
  thresholder->SetUpperThreshold(0.0);
  thresholder->SetOutsideValue(itk::NumericTraits<OutputPixelType>::min());
  thresholder->SetInsideValue(itk::NumericTraits<OutputPixelType>::max());
  thresholder->SetInput(geodesicActiveContour->GetOutput());

  using NodeContainer = FastMarchingFilterType::NodeContainer;
  using NodeType = FastMarchingFilterType::NodeType;

  InputImageType::IndexType seedPosition;
  seedPosition[0] = seedPosX;
  seedPosition[1] = seedPosY;

  NodeContainer::Pointer seeds = NodeContainer::New();
  NodeType               node;
  node.SetValue(seedValue);
  node.SetIndex(seedPosition);

  seeds->Initialize();
  seeds->InsertElement(0, node);

  fastMarching->SetTrialPoints(seeds);
  fastMarching->SetSpeedConstant(1.0);

  using CastFilterType = itk::RescaleIntensityImageFilter<InputImageType, OutputImageType>;

  CastFilterType::Pointer caster1 = CastFilterType::New();
  CastFilterType::Pointer caster2 = CastFilterType::New();
  CastFilterType::Pointer caster3 = CastFilterType::New();
  CastFilterType::Pointer caster4 = CastFilterType::New();

  WriterType::Pointer writer1 = WriterType::New();
  WriterType::Pointer writer2 = WriterType::New();
  WriterType::Pointer writer3 = WriterType::New();
  WriterType::Pointer writer4 = WriterType::New();

  caster1->SetInput(smoothing->GetOutput());
  writer1->SetInput(caster1->GetOutput());
  writer1->SetFileName("GeodesicActiveContourImageFilterOutput1.png");
  caster1->SetOutputMinimum(itk::NumericTraits<OutputPixelType>::min());
  caster1->SetOutputMaximum(itk::NumericTraits<OutputPixelType>::max());
  writer1->Update();

  caster2->SetInput(gradientMagnitude->GetOutput());
  writer2->SetInput(caster2->GetOutput());
  writer2->SetFileName("GeodesicActiveContourImageFilterOutput2.png");
  caster2->SetOutputMinimum(itk::NumericTraits<OutputPixelType>::min());
  caster2->SetOutputMaximum(itk::NumericTraits<OutputPixelType>::max());
  writer2->Update();

  caster3->SetInput(sigmoid->GetOutput());
  writer3->SetInput(caster3->GetOutput());
  writer3->SetFileName("GeodesicActiveContourImageFilterOutput3.png");
  caster3->SetOutputMinimum(itk::NumericTraits<OutputPixelType>::min());
  caster3->SetOutputMaximum(itk::NumericTraits<OutputPixelType>::max());
  writer3->Update();

  caster4->SetInput(fastMarching->GetOutput());
  writer4->SetInput(caster4->GetOutput());
  writer4->SetFileName("GeodesicActiveContourImageFilterOutput4.png");
  caster4->SetOutputMinimum(itk::NumericTraits<OutputPixelType>::min());
  caster4->SetOutputMaximum(itk::NumericTraits<OutputPixelType>::max());

  InputImageType::Pointer inImage = reader->GetOutput();
  fastMarching->SetOutputDirection(inImage->GetDirection());
  fastMarching->SetOutputOrigin(inImage->GetOrigin());
  fastMarching->SetOutputRegion(inImage->GetBufferedRegion());
  fastMarching->SetOutputSpacing(inImage->GetSpacing());

  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName(outputFileName);
  writer->SetInput(thresholder->GetOutput());
  try
  {
    writer->Update();
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  std::cout << std::endl;
  std::cout << "Max. no. iterations: " << geodesicActiveContour->GetNumberOfIterations() << std::endl;
  std::cout << "Max. RMS error: " << geodesicActiveContour->GetMaximumRMSError() << std::endl;
  std::cout << std::endl;
  std::cout << "No. elpased iterations: " << geodesicActiveContour->GetElapsedIterations() << std::endl;
  std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() << std::endl;

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

  using InternalWriterType = itk::ImageFileWriter<InputImageType>;

  InternalWriterType::Pointer mapWriter = InternalWriterType::New();
  mapWriter->SetInput(fastMarching->GetOutput());
  mapWriter->SetFileName("GeodesicActiveContourImageFilterOutput4.mha");
  try
  {
    mapWriter->Update();
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  InternalWriterType::Pointer speedWriter = InternalWriterType::New();
  speedWriter->SetInput(sigmoid->GetOutput());
  speedWriter->SetFileName("GeodesicActiveContourImageFilterOutput3.mha");
  try
  {
    speedWriter->Update();
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  InternalWriterType::Pointer gradientWriter = InternalWriterType::New();
  gradientWriter->SetInput(gradientMagnitude->GetOutput());
  gradientWriter->SetFileName("GeodesicActiveContourImageFilterOutput2.mha");
  try
  {
    gradientWriter->Update();
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

Python

#!/usr/bin/env python
import itk

if len(sys.argv) != 11:
    print(
        "Usage: " + sys.argv[0] + "<InputFileName>  <OutputFileName>"
        " <seedX> <seedY> <InitialDistance> <Sigma> <SigmoidAlpha> "
        "<SigmoidBeta> <PropagationScaling> <NumberOfIterations>")
    sys.exit(1)

inputFileName = sys.argv[1]
outputFileName = sys.argv[2]
seedPosX = int(sys.argv[3])
seedPosY = int(sys.argv[4])

initialDistance = float(sys.argv[5])
sigma = float(sys.argv[6])
alpha = float(sys.argv[7])
beta = float(sys.argv[8])
propagationScaling = float(sys.argv[9])
numberOfIterations = int(sys.argv[10])
seedValue = -initialDistance

Dimension = 2

InputPixelType = itk.F
OutputPixelType = itk.UC

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

ReaderType = itk.ImageFileReader[InputImageType]
WriterType = itk.ImageFileWriter[OutputImageType]

reader = ReaderType.New()
reader.SetFileName(inputFileName)

SmoothingFilterType = itk.CurvatureAnisotropicDiffusionImageFilter[
    InputImageType, InputImageType]
smoothing = SmoothingFilterType.New()
smoothing.SetTimeStep(0.125)
smoothing.SetNumberOfIterations(5)
smoothing.SetConductanceParameter(9.0)
smoothing.SetInput(reader.GetOutput())

GradientFilterType = itk.GradientMagnitudeRecursiveGaussianImageFilter[
    InputImageType, InputImageType]
gradientMagnitude = GradientFilterType.New()
gradientMagnitude.SetSigma(sigma)
gradientMagnitude.SetInput(smoothing.GetOutput())

SigmoidFilterType = itk.SigmoidImageFilter[InputImageType, InputImageType]
sigmoid = SigmoidFilterType.New()
sigmoid.SetOutputMinimum(0.0)
sigmoid.SetOutputMaximum(1.0)
sigmoid.SetAlpha(alpha)
sigmoid.SetBeta(beta)
sigmoid.SetInput(gradientMagnitude.GetOutput())

FastMarchingFilterType = itk.FastMarchingImageFilter[
    InputImageType, InputImageType]
fastMarching = FastMarchingFilterType.New()

GeoActiveContourFilterType = itk.GeodesicActiveContourLevelSetImageFilter[
    InputImageType, InputImageType, InputPixelType]
geodesicActiveContour = GeoActiveContourFilterType.New()
geodesicActiveContour.SetPropagationScaling(propagationScaling)
geodesicActiveContour.SetCurvatureScaling(1.0)
geodesicActiveContour.SetAdvectionScaling(1.0)
geodesicActiveContour.SetMaximumRMSError(0.02)
geodesicActiveContour.SetNumberOfIterations(numberOfIterations)
geodesicActiveContour.SetInput(fastMarching.GetOutput())
geodesicActiveContour.SetFeatureImage(sigmoid.GetOutput())

ThresholdingFilterType = itk.BinaryThresholdImageFilter[
    InputImageType, OutputImageType]
thresholder = ThresholdingFilterType.New()
thresholder.SetLowerThreshold(-1000.0)
thresholder.SetUpperThreshold(0.0)
thresholder.SetOutsideValue(itk.NumericTraits[OutputPixelType].min())
thresholder.SetInsideValue(itk.NumericTraits[OutputPixelType].max())
thresholder.SetInput(geodesicActiveContour.GetOutput())

seedPosition = itk.Index[Dimension]()
seedPosition[0] = seedPosX
seedPosition[1] = seedPosY

node = itk.LevelSetNode[InputPixelType, Dimension]()
node.SetValue(seedValue)
node.SetIndex(seedPosition)

seeds = itk.VectorContainer[
    itk.UI, itk.LevelSetNode[InputPixelType, Dimension]].New()
seeds.Initialize()
seeds.InsertElement(0, node)

fastMarching.SetTrialPoints(seeds)
fastMarching.SetSpeedConstant(1.0)

CastFilterType = itk.RescaleIntensityImageFilter[
    InputImageType, OutputImageType]

caster1 = CastFilterType.New()
caster2 = CastFilterType.New()
caster3 = CastFilterType.New()
caster4 = CastFilterType.New()

writer1 = WriterType.New()
writer2 = WriterType.New()
writer3 = WriterType.New()
writer4 = WriterType.New()

caster1.SetInput(smoothing.GetOutput())
writer1.SetInput(caster1.GetOutput())
writer1.SetFileName("GeodesicActiveContourImageFilterOutput1.png")
caster1.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
caster1.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
writer1.Update()

caster2.SetInput(gradientMagnitude.GetOutput())
writer2.SetInput(caster2.GetOutput())
writer2.SetFileName("GeodesicActiveContourImageFilterOutput2.png")
caster2.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
caster2.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
writer2.Update()

caster3.SetInput(sigmoid.GetOutput())
writer3.SetInput(caster3.GetOutput())
writer3.SetFileName("GeodesicActiveContourImageFilterOutput3.png")
caster3.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
caster3.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())
writer3.Update()

caster4.SetInput(fastMarching.GetOutput())
writer4.SetInput(caster4.GetOutput())
writer4.SetFileName("GeodesicActiveContourImageFilterOutput4.png")
caster4.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min())
caster4.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max())

fastMarching.SetOutputSize(
    reader.GetOutput().GetBufferedRegion().GetSize())

writer = WriterType.New()
writer.SetFileName(outputFileName)
writer.SetInput(thresholder.GetOutput())
writer.Update()

print(
    "Max. no. iterations: " +
    str(geodesicActiveContour.GetNumberOfIterations()) + "\n")
print(
    "Max. RMS error: " +
    str(geodesicActiveContour.GetMaximumRMSError()) + "\n")
print(
    "No. elpased iterations: " +
    str(geodesicActiveContour.GetElapsedIterations()) + "\n")
print("RMS change: " + str(geodesicActiveContour.GetRMSChange()) + "\n")

writer4.Update()

InternalWriterType = itk.ImageFileWriter[InputImageType]

mapWriter = InternalWriterType.New()
mapWriter.SetInput(fastMarching.GetOutput())
mapWriter.SetFileName("GeodesicActiveContourImageFilterOutput4.mha")
mapWriter.Update()

speedWriter = InternalWriterType.New()
speedWriter.SetInput(sigmoid.GetOutput())
speedWriter.SetFileName("GeodesicActiveContourImageFilterOutput3.mha")
speedWriter.Update()

gradientWriter = InternalWriterType.New()
gradientWriter.SetInput(gradientMagnitude.GetOutput())
gradientWriter.SetFileName("GeodesicActiveContourImageFilterOutput2.mha")
gradientWriter.Update()

Classes demonstrated