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 argparse

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

template<typename TPixel>
class Hessian3DToVesselnessMeasureImageFilter : public itk::ImageToImageFilter<Image<SymmetricSecondRankTensor<double, 3>, 3>, Image<TPixel, 3>>

Line filter to provide a vesselness measure for tubular objects from the hessian matrix.

The filter takes as input an image of hessian pixels (SymmetricSecondRankTensor pixels) and preserves pixels that have eigen values \lambda_3 close to 0 and \lambda_2 and \lambda_1 as large negative values (for bright tubular structures).

| \lambda_1 | < | \lambda_2 | < | \lambda_3 |

  • Bright tubular structures will have low \lambda_1 and large negative values of \lambda_2 and \lambda_3.

  • Conversely dark tubular structures will have a low value of \lambda_1 and large positive values of \lambda_2 and \lambda_3.

  • Bright plate like structures have low values of \lambda_1 and \lambda_2 and large negative values of \lambda_3

  • Dark plate like structures have low values of \lambda_1 and \lambda_2 and large positive values of \lambda_3

  • Bright spherical (blob) like structures have all three eigen values as large negative numbers

  • Dark spherical (blob) like structures have all three eigen values as large positive numbers

This filter is used to discriminate the Bright tubular structures.

Notes:

The filter takes into account that the eigen values play a crucial role in discriminating shape and orientation of structures.

http://www.image.med.osaka-u.ac.jp/member/yoshi/paper/linefilter.pdf

References:

”3D Multi-scale line filter for segmentation and visualization of

curvilinear structures in medical images”, Yoshinobu Sato, Shin Nakajima, Hideki Atsumi, Thomas Koller, Guido Gerig, Shigeyuki Yoshida, Ron Kikinis.

See

HessianRecursiveGaussianImageFilter

See

SymmetricEigenAnalysisImageFilter

See

SymmetricSecondRankTensor

See itk::Hessian3DToVesselnessMeasureImageFilter for additional documentation.
template<typename TInputImage, typename TOutputImage = Image<SymmetricSecondRankTensor<typename NumericTraits<typename TInputImage::PixelType>::RealType, TInputImage::ImageDimension>, TInputImage::ImageDimension>>
class HessianRecursiveGaussianImageFilter : public itk::ImageToImageFilter<TInputImage, TOutputImage>

Computes the Hessian matrix of an image by convolution with the Second and Cross derivatives of a Gaussian.

This filter is implemented using the recursive gaussian filters

See itk::HessianRecursiveGaussianImageFilter for additional documentation.