Singlephase Chan and Vese Sparse Field Level Set Segmentation

Warning

Fix Problem Contains problems not fixed from original wiki.

Synopsis

Single-phase Chan And Vese Sparse Field Level Set Segmentation.

Results

Note

Help Wanted Implementation of Results for sphinx examples containing this message. Reconfiguration of CMakeList.txt may be necessary. Write An Example <https://itk.org/ITKExamples/Documentation/Contribute/WriteANewExample.html>

Code

C++

// The use of the ScalarChanAndVeseSparseLevelSetImageFilter is
// illustrated in the following example. The implementation of this filter in
// ITK is based on the paper by Chan And Vese.  This
// implementation extends the functionality of the
// level-set filters in ITK by using region-based variational techniques. These methods
// do not rely on the presence of edges in the images.
//
// ScalarChanAndVeseSparseLevelSetImageFilter expects two inputs.  The first is
// an initial level set in the form of an \doxygen{Image}. The second input
// is a feature image. For this algorithm, the feature image is the original
// raw or preprocessed image. Several parameters are required by the algorithm
// for regularization and weights of different energy terms. The user is encouraged to
// change different parameter settings to optimize the code example on their images.
//
// Let's start by including the headers of the main filters involved in the
// preprocessing.
//


#include "itkScalarChanAndVeseSparseLevelSetImageFilter.h"
#include "itkScalarChanAndVeseLevelSetFunctionData.h"
#include "itkConstrainedRegionBasedLevelSetFunctionSharedData.h"
#include "itkFastMarchingImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImage.h"
#include "itkAtanRegularizedHeavisideStepFunction.h"

int
main(int argc, char ** argv)
{
  if (argc < 6)
  {
    std::cerr << "Missing arguments" << std::endl;
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0] << " featureImage outputImage";
    std::cerr << " startx starty seedValue" << std::endl;
    return EXIT_FAILURE;
  }

  unsigned int nb_iteration = 500;
  double       rms = 0.;
  double       epsilon = 1.;
  double       curvature_weight = 0.;
  double       area_weight = 0.;
  double       reinitialization_weight = 0.;
  double       volume_weight = 0.;
  double       volume = 0.;
  double       l1 = 1.;
  double       l2 = 1.;

  //
  //  We now define the image type using a particular pixel type and
  //  dimension. In this case the \code{float} type is used for the pixels
  //  due to the requirements of the smoothing filter.
  //
  constexpr unsigned int Dimension = 2;
  using ScalarPixelType = float;
  using InternalImageType = itk::Image<ScalarPixelType, Dimension>;


  using DataHelperType = itk::ScalarChanAndVeseLevelSetFunctionData<InternalImageType, InternalImageType>;

  using SharedDataHelperType =
    itk::ConstrainedRegionBasedLevelSetFunctionSharedData<InternalImageType, InternalImageType, DataHelperType>;

  using LevelSetFunctionType =
    itk::ScalarChanAndVeseLevelSetFunction<InternalImageType, InternalImageType, SharedDataHelperType>;


  //  We declare now the type of the numerically discretized Step and Delta functions that
  //  will be used in the level-set computations for foreground and background regions
  //
  using DomainFunctionType = itk::AtanRegularizedHeavisideStepFunction<ScalarPixelType, ScalarPixelType>;

  DomainFunctionType::Pointer domainFunction = DomainFunctionType::New();
  domainFunction->SetEpsilon(epsilon);

  // We instantiate reader and writer types in the following lines.
  //
  using ReaderType = itk::ImageFileReader<InternalImageType>;
  using WriterType = itk::ImageFileWriter<InternalImageType>;

  ReaderType::Pointer reader = ReaderType::New();
  WriterType::Pointer writer = WriterType::New();

  reader->SetFileName(argv[1]);
  reader->Update();

  writer->SetFileName(argv[2]);

  InternalImageType::Pointer featureImage = reader->GetOutput();

  //  We declare now the type of the FastMarchingImageFilter that
  //  will be used to generate the initial level set in the form of a distance
  //  map.
  //
  using FastMarchingFilterType = itk::FastMarchingImageFilter<InternalImageType, InternalImageType>;

  FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();

  //  The FastMarchingImageFilter requires the user to provide a seed
  //  point from which the level set will be generated. The user can actually
  //  pass not only one seed point but a set of them. Note the the
  //  FastMarchingImageFilter is used here only as a helper in the
  //  determination of an initial level set. We could have used the
  //  \doxygen{DanielssonDistanceMapImageFilter} in the same way.
  //
  //  The seeds are passed stored in a container. The type of this
  //  container is defined as \code{NodeContainer} among the
  //  FastMarchingImageFilter traits.
  //
  using NodeContainer = FastMarchingFilterType::NodeContainer;
  using NodeType = FastMarchingFilterType::NodeType;

  NodeContainer::Pointer seeds = NodeContainer::New();

  InternalImageType::IndexType seedPosition;

  seedPosition[0] = std::stoi(argv[3]);
  seedPosition[1] = std::stoi(argv[4]);

  const double initialDistance = std::stod(argv[5]);

  NodeType node;

  const double seedValue = -initialDistance;

  node.SetValue(seedValue);
  node.SetIndex(seedPosition);

  //  The list of nodes is initialized and then every node is inserted using
  //  the \code{InsertElement()}.
  //
  seeds->Initialize();
  seeds->InsertElement(0, node);


  //  The set of seed nodes is passed now to the
  //  FastMarchingImageFilter with the method
  //  \code{SetTrialPoints()}.
  //
  fastMarching->SetTrialPoints(seeds);


  //  Since the FastMarchingImageFilter is used here just as a
  //  Distance Map generator. It does not require a speed image as input.
  //  Instead the constant value $1.0$ is passed using the
  //  \code{SetSpeedConstant()} method.
  //
  fastMarching->SetSpeedConstant(1.0);

  //  The FastMarchingImageFilter requires the user to specify the
  //  size of the image to be produced as output. This is done using the
  //  \code{SetOutputSize()}. Note that the size is obtained here from the
  //  output image of the smoothing filter. The size of this image is valid
  //  only after the \code{Update()} methods of this filter has been called
  //  directly or indirectly.
  //
  fastMarching->SetOutputSize(featureImage->GetBufferedRegion().GetSize());
  fastMarching->Update();

  //  We declare now the type of the ScalarChanAndVeseSparseLevelSetImageFilter that
  //  will be used to generate a segmentation.
  //

  using MultiLevelSetType = itk::ScalarChanAndVeseSparseLevelSetImageFilter<InternalImageType,
                                                                            InternalImageType,
                                                                            InternalImageType,
                                                                            LevelSetFunctionType,
                                                                            SharedDataHelperType>;

  MultiLevelSetType::Pointer levelSetFilter = MultiLevelSetType::New();

  //  We set the function count to 1 since a single level-set is being evolved.
  //
  levelSetFilter->SetFunctionCount(1);

  //  Set the feature image and initial level-set image as output of the
  //  fast marching image filter.
  //
  levelSetFilter->SetFeatureImage(featureImage);
  levelSetFilter->SetLevelSet(0, fastMarching->GetOutput());

  //  Once activiated the level set evolution will stop if the convergence
  //  criteria or if the maximum number of iterations is reached.  The
  //  convergence criteria is defined in terms of the root mean squared (RMS)
  //  change in the level set function. The evolution is said to have
  //  converged if the RMS change is below a user specified threshold.  In a
  //  real application is desirable to couple the evolution of the zero set
  //  to a visualization module allowing the user to follow the evolution of
  //  the zero set. With this feedback, the user may decide when to stop the
  //  algorithm before the zero set leaks through the regions of low gradient
  //  in the contour of the anatomical structure to be segmented.
  //
  levelSetFilter->SetNumberOfIterations(nb_iteration);
  levelSetFilter->SetMaximumRMSError(rms);

  //  Often, in real applications, images have different pixel resolutions. In such
  //  cases, it is best to use the native spacings to compute derivatives etc rather
  //  than sampling the images.
  //
  levelSetFilter->SetUseImageSpacing(1);

  //  For large images, we may want to compute the level-set over the initial supplied
  //  level-set image. This saves a lot of memory.
  //
  levelSetFilter->SetInPlace(false);

  //  For the level set with phase 0, set different parameters and weights. This may
  //  to be set in a loop for the case of multiple level-sets evolving simultaneously.
  //
  levelSetFilter->GetDifferenceFunction(0)->SetDomainFunction(domainFunction);
  levelSetFilter->GetDifferenceFunction(0)->SetCurvatureWeight(curvature_weight);
  levelSetFilter->GetDifferenceFunction(0)->SetAreaWeight(area_weight);
  levelSetFilter->GetDifferenceFunction(0)->SetReinitializationSmoothingWeight(reinitialization_weight);
  levelSetFilter->GetDifferenceFunction(0)->SetVolumeMatchingWeight(volume_weight);
  levelSetFilter->GetDifferenceFunction(0)->SetVolume(volume);
  levelSetFilter->GetDifferenceFunction(0)->SetLambda1(l1);
  levelSetFilter->GetDifferenceFunction(0)->SetLambda2(l2);

  levelSetFilter->Update();

  writer->SetInput(levelSetFilter->GetOutput());

  try
  {
    writer->Update();
  }
  catch (itk::ExceptionObject & excep)
  {
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
    return -1;
  }

  return EXIT_SUCCESS;
}

Classes demonstrated

template<typename TInputImage, typename TFeatureImage, typename TOutputImage, typename TFunction = ScalarChanAndVeseLevelSetFunction<TInputImage, TFeatureImage>, class TSharedData = typename TFunction::SharedDataType, typename TIdCell = unsigned int>
class ScalarChanAndVeseSparseLevelSetImageFilter : public itk::MultiphaseSparseFiniteDifferenceImageFilter<TInputImage, TFeatureImage, TOutputImage, TFunction, TIdCell>

Sparse implementation of the Chan and Vese multiphase level set image filter.

This code was adapted from the paper: 

   "An active contour model without edges"
    T. Chan and L. Vese.
    In Scale-Space Theories in Computer Vision, pages 141-151, 1999.

This code was taken from the Insight Journal paper: 

"Cell Tracking using Coupled Active Surfaces for Nuclei and Membranes"
http://www.insight-journal.org/browse/publication/642
https://hdl.handle.net/10380/3055
Author

Mosaliganti K., Smith B., Gelas A., Gouaillard A., Megason S.

That is based on the papers: 

"Level Set Segmentation: Active Contours without edge"
http://www.insight-journal.org/browse/publication/322
https://hdl.handle.net/1926/1532

and

"Level set segmentation using coupled active surfaces"
http://www.insight-journal.org/browse/publication/323
https://hdl.handle.net/1926/1533

ITK Sphinx Examples:

See itk::ScalarChanAndVeseSparseLevelSetImageFilter for additional documentation.