Segment Blood Vessels With Multi-Scale Hessian-Based Measure#

Synopsis#

Segment blood vessels with multi-scale Hessian-based measure.

Results#

Input image

Input image#

Output image

Output image#

Code#

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 with multi-scale Hessian-based measure."
)
parser.add_argument("input_image")
parser.add_argument("output_image")
parser.add_argument("--sigma_minimum", type=float, default=1.0)
parser.add_argument("--sigma_maximum", type=float, default=10.0)
parser.add_argument("--number_of_sigma_steps", type=int, default=10)
args = parser.parse_args()

input_image = itk.imread(args.input_image, itk.F)

ImageType = type(input_image)
Dimension = input_image.GetImageDimension()
HessianPixelType = itk.SymmetricSecondRankTensor[itk.D, Dimension]
HessianImageType = itk.Image[HessianPixelType, Dimension]

objectness_filter = itk.HessianToObjectnessMeasureImageFilter[
    HessianImageType, ImageType
].New()
objectness_filter.SetBrightObject(False)
objectness_filter.SetScaleObjectnessMeasure(False)
objectness_filter.SetAlpha(0.5)
objectness_filter.SetBeta(1.0)
objectness_filter.SetGamma(5.0)

multi_scale_filter = itk.MultiScaleHessianBasedMeasureImageFilter[
    ImageType, HessianImageType, ImageType
].New()
multi_scale_filter.SetInput(input_image)
multi_scale_filter.SetHessianToMeasureFilter(objectness_filter)
multi_scale_filter.SetSigmaStepMethodToLogarithmic()
multi_scale_filter.SetSigmaMinimum(args.sigma_minimum)
multi_scale_filter.SetSigmaMaximum(args.sigma_maximum)
multi_scale_filter.SetNumberOfSigmaSteps(args.number_of_sigma_steps)

OutputPixelType = itk.UC
OutputImageType = itk.Image[OutputPixelType, Dimension]

rescale_filter = itk.RescaleIntensityImageFilter[ImageType, OutputImageType].New()
rescale_filter.SetInput(multi_scale_filter)

itk.imwrite(rescale_filter.GetOutput(), args.output_image)

C++#

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkHessianToObjectnessMeasureImageFilter.h"
#include "itkMultiScaleHessianBasedMeasureImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"

int
main(int argc, char * argv[])
{
  if (argc < 3)
  {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0];
    std::cerr << " <InputFileName> <OutputFileName>";
    std::cerr << " [SigmaMinimum] [SigmaMaximum]";
    std::cerr << " [NumberOfSigmaSteps]";
    std::cerr << std::endl;
    return EXIT_FAILURE;
  }

  const char * inputFileName = argv[1];
  const char * outputFileName = argv[2];
  double       sigmaMinimum = 1.0;
  if (argc > 3)
  {
    sigmaMinimum = std::stod(argv[3]);
  }
  double sigmaMaximum = 10.0;
  if (argc > 4)
  {
    sigmaMaximum = std::stod(argv[4]);
  }
  unsigned int numberOfSigmaSteps = 10;
  if (argc > 5)
  {
    numberOfSigmaSteps = std::stoi(argv[5]);
  }

  constexpr unsigned int Dimension = 2;

  using PixelType = float;
  using ImageType = itk::Image<PixelType, Dimension>;

  const auto input = itk::ReadImage<ImageType>(inputFileName);

  using HessianPixelType = itk::SymmetricSecondRankTensor<double, Dimension>;
  using HessianImageType = itk::Image<HessianPixelType, Dimension>;
  using ObjectnessFilterType = itk::HessianToObjectnessMeasureImageFilter<HessianImageType, ImageType>;
  auto objectnessFilter = ObjectnessFilterType::New();
  objectnessFilter->SetBrightObject(false);
  objectnessFilter->SetScaleObjectnessMeasure(false);
  objectnessFilter->SetAlpha(0.5);
  objectnessFilter->SetBeta(1.0);
  objectnessFilter->SetGamma(5.0);

  using MultiScaleEnhancementFilterType =
    itk::MultiScaleHessianBasedMeasureImageFilter<ImageType, HessianImageType, ImageType>;
  auto multiScaleEnhancementFilter = MultiScaleEnhancementFilterType::New();
  multiScaleEnhancementFilter->SetInput(input);
  multiScaleEnhancementFilter->SetHessianToMeasureFilter(objectnessFilter);
  multiScaleEnhancementFilter->SetSigmaStepMethodToLogarithmic();
  multiScaleEnhancementFilter->SetSigmaMinimum(sigmaMinimum);
  multiScaleEnhancementFilter->SetSigmaMaximum(sigmaMaximum);
  multiScaleEnhancementFilter->SetNumberOfSigmaSteps(numberOfSigmaSteps);

  using OutputImageType = itk::Image<unsigned char, Dimension>;
  using RescaleFilterType = itk::RescaleIntensityImageFilter<ImageType, OutputImageType>;
  auto rescaleFilter = RescaleFilterType::New();
  rescaleFilter->SetInput(multiScaleEnhancementFilter->GetOutput());

  try
  {
    itk::WriteImage(rescaleFilter->GetOutput(), outputFileName);
  }
  catch (const itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

Classes demonstrated#

template<typename TInputImage, typename THessianImage, typename TOutputImage = TInputImage>
class MultiScaleHessianBasedMeasureImageFilter : public itk::ImageToImageFilter<TInputImage, TOutputImage>

A filter to enhance structures using Hessian eigensystem-based measures in a multiscale framework.

The filter evaluates a Hessian-based enhancement measure, such as vesselness or objectness, at different scale levels. The Hessian-based measure is computed from the Hessian image at each scale level and the best response is selected.

Minimum and maximum sigma value can be set using SetMinSigma and SetMaxSigma methods respectively. The number of scale levels is set using SetNumberOfSigmaSteps method. Exponentially distributed scale levels are computed within the bound set by the minimum and maximum sigma values

The filter computes a second output image (accessed by the GetScalesOutput method) containing the scales at which each pixel gave the best response.

This code was contributed in the Insight Journal paper: “Generalizing vesselness with respect to dimensionality and shape” by Antiga L. https://www.insight-journal.org/browse/publication/175

Author

Luca Antiga Ph.D. Medical Imaging Unit, Bioengineering Department, Mario Negri Institute, Italy.

See

HessianToObjectnessMeasureImageFilter

See

Hessian3DToVesselnessMeasureImageFilter

See

HessianSmoothed3DToVesselnessMeasureImageFilter

See

HessianRecursiveGaussianImageFilter

See

SymmetricEigenAnalysisImageFilter

See

SymmetricSecondRankTensor

See itk::MultiScaleHessianBasedMeasureImageFilter for additional documentation.