Segment Blood Vessels

Synopsis

The example takes an image (say MRA image), computes the vesselness measure of the image using the HessianRecursiveGaussianImageFilter and the Hessian3DToVesselnessMeasureImageFilter. The goal is to detect bright, tubular structures in the image.

Note that since the algorithm is based on the Hessian, it will also identify black tubular structures.

An important parameter is the Sigma that determines the amount of smoothing applied during Hessian estimation. A larger Sigma will decrease the identification of noise or small structures as vessels. In this example, the Sigma is large enough only vessels comprising the Circle of Willis and other large vessels are segmented.

Results

Input image

Input head MRA.

Output image

Output vessel segmentation.

Code

C++

#include "itkHessian3DToVesselnessMeasureImageFilter.h"
#include "itkHessianRecursiveGaussianImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

int
main(int argc, char * argv[])
{
  if (argc < 3)
  {
    std::cerr << "Usage: <InputImage> <OutputImage> [Sigma] [Alpha1] [Alpha2]" << std::endl;
    return EXIT_FAILURE;
  }
  const char * inputImage = argv[1];
  const char * outputImage = argv[2];
  double       sigma = 1.0;
  if (argc > 3)
  {
    sigma = std::atof(argv[3]);
  }
  double alpha1 = 0.5;
  if (argc > 4)
  {
    alpha1 = std::atof(argv[4]);
  }
  double alpha2 = 2.0;
  if (argc > 5)
  {
    alpha2 = std::atof(argv[5]);
  }

  constexpr unsigned int Dimension = 3;
  using PixelType = float;
  using ImageType = itk::Image<PixelType, Dimension>;

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

  using HessianFilterType = itk::HessianRecursiveGaussianImageFilter<ImageType>;
  HessianFilterType::Pointer hessianFilter = HessianFilterType::New();
  hessianFilter->SetInput(reader->GetOutput());
  hessianFilter->SetSigma(sigma);

  using VesselnessMeasureFilterType = itk::Hessian3DToVesselnessMeasureImageFilter<PixelType>;
  VesselnessMeasureFilterType::Pointer vesselnessFilter = VesselnessMeasureFilterType::New();
  vesselnessFilter->SetInput(hessianFilter->GetOutput());
  vesselnessFilter->SetAlpha1(alpha1);
  vesselnessFilter->SetAlpha2(alpha2);

  using WriterType = itk::ImageFileWriter<ImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetInput(vesselnessFilter->GetOutput());
  writer->SetFileName(outputImage);

  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("5.0.0"):
    print("ITK 5.0.0 or newer is required.")
    sys.exit(1)

parser = argparse.ArgumentParser(description='Segment blood vessels.')
parser.add_argument('input_image')
parser.add_argument('output_image')
parser.add_argument('--sigma', type=float, default=1.0)
parser.add_argument('--alpha1', type=float, default=0.5)
parser.add_argument('--alpha2', type=float, default=2.0)
args = parser.parse_args()

input_image = itk.imread(args.input_image, itk.ctype('float'))

hessian_image = itk.hessian_recursive_gaussian_image_filter(input_image, sigma=args.sigma)

vesselness_filter = itk.Hessian3DToVesselnessMeasureImageFilter[itk.ctype('float')].New()
vesselness_filter.SetInput(hessian_image)
vesselness_filter.SetAlpha1(args.alpha1)
vesselness_filter.SetAlpha2(args.alpha2)

itk.imwrite(vesselness_filter, args.output_image)

Classes demonstrated