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];
  char * sigma = NULL;
  if( argc > 3 )
    {
    sigma = argv[3];
    }
  char * alpha1 = NULL;
  if( argc > 4 )
    {
    alpha1 = argv[4];
    }
  char * alpha2 = NULL;
  if( argc > 5 )
    {
    alpha2 = argv[5];
    }

  const   unsigned int        Dimension = 3;
  typedef double              InputPixelType;
  typedef float               OutputPixelType;
  typedef itk::Image< InputPixelType, Dimension >   InputImageType;
  typedef itk::Image< OutputPixelType, Dimension >  OutputImageType;

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

  typedef itk::HessianRecursiveGaussianImageFilter< InputImageType >
    HessianFilterType;
  HessianFilterType::Pointer hessianFilter = HessianFilterType::New();
  hessianFilter->SetInput( reader->GetOutput() );
  if( sigma )
    {
    hessianFilter->SetSigma( atof( sigma ) );
    }

  typedef itk::Hessian3DToVesselnessMeasureImageFilter< OutputPixelType >
    VesselnessMeasureFilterType;
  VesselnessMeasureFilterType::Pointer vesselnessFilter = VesselnessMeasureFilterType::New();
  vesselnessFilter->SetInput( hessianFilter->GetOutput() );
  if( alpha1 )
    {
    vesselnessFilter->SetAlpha1( atof( alpha1 ) );
    }
  if( alpha2 )
    {
    vesselnessFilter->SetAlpha2( atof( alpha2 ) );
    }

  typedef itk::ImageFileWriter< OutputImageType > WriterType;
  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;
}

Classes demonstrated

template <typename TPixel>
class itk::Hessian3DToVesselnessMeasureImageFilter

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 discrimintaitng shape and orientation of structures.

http://www.spl.harvard.edu/archive/spl-pre2007/pages/papers/yoshi/

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

SymmetricEigenAnalysisImageFilter

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 itk::HessianRecursiveGaussianImageFilter

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.