ShapeDetectionLevelSetFilter.cxx

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

//  Software Guide : BeginCommandLineArgs
//    INPUTS:  {BrainProtonDensitySlice.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput1.png}
//    ARGUMENTS:    81 114 5 1.0  -0.5  3.0   .05 1
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput1Smoothing.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput1GradientMagnitude.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput1Sigmoid.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput1FastMarching.png}
//  Software Guide : EndCommandLineArgs
//  Software Guide : BeginCommandLineArgs
//    INPUTS:  {BrainProtonDensitySlice.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput2.png}
//    ARGUMENTS:    99 114 5 1.0  -0.5  3.0   .05 1
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput2Smoothing.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput2GradientMagnitude.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput2Sigmoid.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput2FastMarching.png}
//  Software Guide : EndCommandLineArgs
//  Software Guide : BeginCommandLineArgs
//    INPUTS:  {BrainProtonDensitySlice.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput3.png}
//    ARGUMENTS:    56 92 5 1.0  -0.3  2.0   .05 1
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput3Smoothing.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput3GradientMagnitude.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput3Sigmoid.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput3FastMarching.png}
//  Software Guide : EndCommandLineArgs
//  Software Guide : BeginCommandLineArgs
//    INPUTS:  {BrainProtonDensitySlice.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput4.png}
//    ARGUMENTS:    40 90 5 0.5  -0.3  2.0   .05 1
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput4Smoothing.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput4GradientMagnitude.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput4Sigmoid.png}
//    OUTPUTS: {ShapeDetectionLevelSetFilterOutput4FastMarching.png}
//  Software Guide : EndCommandLineArgs

// Software Guide : BeginLatex
//
// The use of the \doxygen{ShapeDetectionLevelSetImageFilter} is illustrated
// in the following example.  The implementation of this filter in ITK is
// based on the paper by Malladi et al \cite{Malladi1995}.  In this
// implementation, the governing differential equation has an additional
// curvature-based term. This term acts as a smoothing term where areas of
// high curvature, assumed to be due to noise, are smoothed out. Scaling
// parameters are used to control the tradeoff between the expansion term and
// the smoothing term.  One consequence of this additional curvature term is
// that the fast marching algorithm is no longer applicable, because the
// contour is no longer guaranteed to always be expanding. Instead, the level
// set function is updated iteratively.
//
// The ShapeDetectionLevelSetImageFilter expects two inputs,
// the first being an initial Level Set in the form of an
// \doxygen{Image}, and the second being a feature image. For this algorithm,
// the feature image is an edge potential image that basically
// follows the same rules applicable to the speed image used for the
// FastMarchingImageFilter discussed in
// Section~\ref{sec:FastMarchingImageFilter}.
//
// In this example we use an FastMarchingImageFilter to produce the initial
// level set as the distance function to a set of user-provided seeds. The
// FastMarchingImageFilter is run with a constant speed value which enables
// us to employ this filter as a distance map calculator.
//
// \begin{figure} \center
// \includegraphics[width=\textwidth]{ShapeDetectionCollaborationDiagram1}
// \itkcaption[ShapeDetectionLevelSetImageFilter collaboration
// diagram]{Collaboration diagram for the ShapeDetectionLevelSetImageFilter
// applied to a segmentation task.}
// \label{fig:ShapeDetectionCollaborationDiagram}
// \end{figure}
//
// Figure~\ref{fig:ShapeDetectionCollaborationDiagram} shows the major
// components involved in the application of the
// ShapeDetectionLevelSetImageFilter to a segmentation task. The first stage
// involves smoothing using the
// \doxygen{CurvatureAnisotropicDiffusionImageFilter}. The smoothed image is
// passed as the input for the
// \doxygen{GradientMagnitudeRecursiveGaussianImageFilter} and then to the
// \doxygen{SigmoidImageFilter} in order to produce the edge potential image.
// A set of user-provided seeds is passed to an FastMarchingImageFilter in
// order to compute the distance map. A constant value is subtracted from
// this map in order to obtain a level set in which the \emph{zero set}
// represents the initial contour. This level set is also passed as input to
// the ShapeDetectionLevelSetImageFilter.
//
// Finally, the level set at the output of the
// ShapeDetectionLevelSetImageFilter is passed to an
// BinaryThresholdImageFilter in order to produce a binary mask representing
// the segmented object.
//
// Let's start by including the headers of the main filters involved in the
// preprocessing.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
// Software Guide : EndCodeSnippet

//  Software Guide : BeginLatex
//
// The edge potential map is generated using these filters as in the previous
// example.
//
//  Software Guide : EndLatex


//  Software Guide : BeginLatex
//
//  We will need the Image class, the FastMarchingImageFilter class and the
//  ShapeDetectionLevelSetImageFilter class. Hence we include their headers
//  here.
//
//  Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
#include "itkFastMarchingImageFilter.h"
#include "itkShapeDetectionLevelSetImageFilter.h"
// Software Guide : EndCodeSnippet


//  Software Guide : BeginLatex
//
//  The level set resulting from the ShapeDetectionLevelSetImageFilter will
//  be thresholded at the zero level in order to get a binary image
//  representing the segmented object. The BinaryThresholdImageFilter is used
//  for this purpose.
//
//  Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
#include "itkBinaryThresholdImageFilter.h"
// Software Guide : EndCodeSnippet


//  Reading and writing images will be done with the \doxygen{ImageFileReader}
//  and \doxygen{ImageFileWriter}.

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"

#include "itksys/SystemTools.hxx"


//  The RescaleIntensityImageFilter is used to renormailize the output
//  of filters before sending them to files.
#include "itkRescaleIntensityImageFilter.h"


int
main(int argc, char * argv[])
{
  if (argc < 11)
  {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " inputImage  outputImage";
    std::cerr << " seedX seedY InitialDistance";
    std::cerr << " Sigma SigmoidAlpha SigmoidBeta ";
    std::cerr << " curvatureScaling propagationScaling" << std::endl;
    return EXIT_FAILURE;
  }
  const std::string inputImageFile(argv[1]);
  const std::string outputImageFile(argv[2]);


  //  Software Guide : BeginLatex
  //
  //  We now define the image type using a particular pixel type and a
  //  dimension. In this case the \code{float} type is used for the pixels
  //  due to the requirements of the smoothing filter.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  using InternalPixelType = float;
  constexpr unsigned int Dimension = 2;
  using InternalImageType = itk::Image<InternalPixelType, Dimension>;
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The output image, on the other hand, is declared to be binary.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  using OutputPixelType = unsigned char;
  using OutputImageType = itk::Image<OutputPixelType, Dimension>;
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The type of the BinaryThresholdImageFilter filter is instantiated below
  //  using the internal image type and the output image type.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  using ThresholdingFilterType =
    itk::BinaryThresholdImageFilter<InternalImageType, OutputImageType>;
  auto thresholder = ThresholdingFilterType::New();
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The upper threshold of the BinaryThresholdImageFilter is set
  //  to $0.0$ in order to display the zero set of the resulting level
  //  set. The lower threshold is set to a large negative number in order to
  //  ensure that the interior of the segmented object will appear
  //  inside the binary region.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  thresholder->SetLowerThreshold(-1000.0);
  thresholder->SetUpperThreshold(0.0);

  thresholder->SetOutsideValue(0);
  thresholder->SetInsideValue(255);
  // Software Guide : EndCodeSnippet


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

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

  reader->SetFileName(inputImageFile);
  writer->SetFileName(outputImageFile);


  //  The RescaleIntensityImageFilter type is declared below. This filter will
  //  renormalize image before sending them to writers.
  //
  using CastFilterType =
    itk::RescaleIntensityImageFilter<InternalImageType, OutputImageType>;

  //  Software Guide : BeginLatex
  //
  //  The CurvatureAnisotropicDiffusionImageFilter type is instantiated using
  //  the internal image type.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  using SmoothingFilterType =
    itk::CurvatureAnisotropicDiffusionImageFilter<InternalImageType,
                                                  InternalImageType>;
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The filter is instantiated by invoking the \code{New()} method and
  //  assigning the result to a \doxygen{SmartPointer}.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  auto smoothing = SmoothingFilterType::New();
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The types of the GradientMagnitudeRecursiveGaussianImageFilter and
  //  SigmoidImageFilter are instantiated using the internal
  //  image type.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  using GradientFilterType =
    itk::GradientMagnitudeRecursiveGaussianImageFilter<InternalImageType,
                                                       InternalImageType>;

  using SigmoidFilterType =
    itk::SigmoidImageFilter<InternalImageType, InternalImageType>;
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The corresponding filter objects are created with the method
  //  \code{New()}.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  auto gradientMagnitude = GradientFilterType::New();
  auto sigmoid = SigmoidFilterType::New();
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The minimum and maximum values of the SigmoidImageFilter output are
  //  defined with the methods \code{SetOutputMinimum()} and
  //  \code{SetOutputMaximum()}. In our case, we want these two values to be
  //  $0.0$ and $1.0$ respectively in order to get a nice speed image to feed
  //  to the FastMarchingImageFilter. Additional details on the use of the
  //  SigmoidImageFilter are presented in
  //  Section~\ref{sec:IntensityNonLinearMapping}.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  sigmoid->SetOutputMinimum(0.0);
  sigmoid->SetOutputMaximum(1.0);
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  We now declare the type of the FastMarchingImageFilter that
  //  will be used to generate the initial level set in the form of a distance
  //  map.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  using FastMarchingFilterType =
    itk::FastMarchingImageFilter<InternalImageType, InternalImageType>;
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  Next we construct one filter of this class using the \code{New()}
  //  method.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  auto fastMarching = FastMarchingFilterType::New();
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  In the following lines we instantiate the type of the
  //  ShapeDetectionLevelSetImageFilter and create an object of this type
  //  using the \code{New()} method.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  using ShapeDetectionFilterType =
    itk::ShapeDetectionLevelSetImageFilter<InternalImageType,
                                           InternalImageType>;
  auto shapeDetection = ShapeDetectionFilterType::New();
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The filters are now connected in a pipeline indicated in
  //  Figure~\ref{fig:ShapeDetectionCollaborationDiagram} with the following
  //  code.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  smoothing->SetInput(reader->GetOutput());
  gradientMagnitude->SetInput(smoothing->GetOutput());
  sigmoid->SetInput(gradientMagnitude->GetOutput());

  shapeDetection->SetInput(fastMarching->GetOutput());
  shapeDetection->SetFeatureImage(sigmoid->GetOutput());

  thresholder->SetInput(shapeDetection->GetOutput());

  writer->SetInput(thresholder->GetOutput());
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The CurvatureAnisotropicDiffusionImageFilter requires a couple
  //  of parameters to be defined. The following are typical values for $2D$
  //  images. However they may have to be adjusted depending on the amount of
  //  noise present in the input image. This filter has been discussed in
  //  Section~\ref{sec:GradientAnisotropicDiffusionImageFilter}.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  smoothing->SetTimeStep(0.125);
  smoothing->SetNumberOfIterations(5);
  smoothing->SetConductanceParameter(9.0);
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The GradientMagnitudeRecursiveGaussianImageFilter performs the
  //  equivalent of a convolution with a Gaussian kernel followed by a
  //  derivative operator. The sigma of this Gaussian can be used to control
  //  the range of influence of the image edges. This filter has been
  //  discussed in
  //  Section~\ref{sec:GradientMagnitudeRecursiveGaussianImageFilter}.
  //
  //  \index{itk::Gradient\-Magnitude\-Recursive\-Gaussian\-Image\-Filter!SetSigma()}
  //
  //  Software Guide : EndLatex

  const double sigma = std::stod(argv[6]);

  // Software Guide : BeginCodeSnippet
  gradientMagnitude->SetSigma(sigma);
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The SigmoidImageFilter requires two parameters that define the linear
  //  transformation to be applied to the sigmoid argument. These parameters
  //  have been discussed in Sections~\ref{sec:IntensityNonLinearMapping} and
  //  \ref{sec:FastMarchingImageFilter}.
  //
  //  Software Guide : EndLatex

  const double alpha = std::stod(argv[7]);
  const double beta = std::stod(argv[8]);

  // Software Guide : BeginCodeSnippet
  sigmoid->SetAlpha(alpha);
  sigmoid->SetBeta(beta);
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  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
  //  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.
  //
  //  \index{itk::FastMarchingImageFilter!Multiple seeds}
  //
  //  The seeds are stored in a container. The type of this
  //  container is defined as \code{NodeContainer} among the
  //  FastMarchingImageFilter traits.
  //
  //  \index{itk::FastMarchingImageFilter!Nodes}
  //  \index{itk::FastMarchingImageFilter!NodeContainer}
  //  \index{itk::FastMarchingImageFilter!NodeType}
  //
  //  Software Guide : EndLatex

  //  Software Guide : BeginCodeSnippet
  using NodeContainer = FastMarchingFilterType::NodeContainer;
  using NodeType = FastMarchingFilterType::NodeType;
  auto seeds = NodeContainer::New();
  //  Software Guide : EndCodeSnippet


  InternalImageType::IndexType seedPosition;

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


  //  Software Guide : BeginLatex
  //
  //  Nodes are created as stack variables and initialized with a value and
  //  an \doxygen{Index} position. Note that we assign the negative of the
  //  value of the user-provided distance to the unique node of the seeds
  //  passed to the FastMarchingImageFilter. In this way, the value
  //  will increment as the front is propagated, until it reaches the zero
  //  value corresponding to the contour. After this, the front will continue
  //  propagating until it fills up the entire image. The initial distance is
  //  taken from the command line arguments. The rule of thumb for the user
  //  is to select this value as the distance from the seed points at which
  //  the initial contour should be.
  //
  //  \index{itk::FastMarchingImageFilter!Seed initialization}
  //
  //  Software Guide : EndLatex

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

  // Software Guide : BeginCodeSnippet
  NodeType     node;
  const double seedValue = -initialDistance;

  node.SetValue(seedValue);
  node.SetIndex(seedPosition);
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The list of nodes is initialized and then every node is inserted using
  //  \code{InsertElement()}.
  //
  //  Software Guide : EndLatex

  //  Software Guide : BeginCodeSnippet
  seeds->Initialize();
  seeds->InsertElement(0, node);
  //  Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The set of seed nodes is now passed to the FastMarchingImageFilter with
  //  the method \code{SetTrialPoints()}.
  //
  //  \index{itk::FastMarchingImageFilter!SetTrialPoints()}
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  fastMarching->SetTrialPoints(seeds);
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  Since the FastMarchingImageFilter is used here only 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.
  //
  //  Software Guide : EndLatex

  //  Software Guide : BeginCodeSnippet
  fastMarching->SetSpeedConstant(1.0);
  //  Software Guide : EndCodeSnippet


  //  Here we configure all the writers required to see the intermediate
  //  outputs of the pipeline. This is added here only for
  //  pedagogical/debugging purposes. These intermediate output are normaly
  //  not required. Only the output of the final thresholding filter should
  //  be relevant.  Observing intermediate output is helpful in the process
  //  of fine tuning the parameters of filters in the pipeline.
  //
  auto caster1 = CastFilterType::New();
  auto caster2 = CastFilterType::New();
  auto caster3 = CastFilterType::New();
  auto caster4 = CastFilterType::New();

  auto writer1 = WriterType::New();
  auto writer2 = WriterType::New();
  auto writer3 = WriterType::New();
  auto writer4 = WriterType::New();

  const std::string outputImageFilePrefix =
    itksys::SystemTools::GetFilenameWithoutExtension(outputImageFile);
  caster1->SetInput(smoothing->GetOutput());
  writer1->SetInput(caster1->GetOutput());
  writer1->SetFileName(outputImageFilePrefix + "Smoothing.png");
  caster1->SetOutputMinimum(0);
  caster1->SetOutputMaximum(255);
  writer1->Update();

  caster2->SetInput(gradientMagnitude->GetOutput());
  writer2->SetInput(caster2->GetOutput());
  writer2->SetFileName(outputImageFilePrefix + "GradientMagnitude.png");
  caster2->SetOutputMinimum(0);
  caster2->SetOutputMaximum(255);
  writer2->Update();

  caster3->SetInput(sigmoid->GetOutput());
  writer3->SetInput(caster3->GetOutput());
  writer3->SetFileName(outputImageFilePrefix + "Sigmoid.png");
  caster3->SetOutputMinimum(0);
  caster3->SetOutputMaximum(255);
  writer3->Update();

  caster4->SetInput(fastMarching->GetOutput());
  writer4->SetInput(caster4->GetOutput());
  writer4->SetFileName(outputImageFilePrefix + "FastMarching.png");
  caster4->SetOutputMinimum(0);
  caster4->SetOutputMaximum(255);


  //  Software Guide : BeginLatex
  //
  //  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 have been called
  //  directly or indirectly.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  fastMarching->SetOutputSize(
    reader->GetOutput()->GetBufferedRegion().GetSize());
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  ShapeDetectionLevelSetImageFilter provides two parameters to control
  //  the competition between the propagation or expansion term and the
  //  curvature smoothing term. The methods \code{SetPropagationScaling()}
  //  and \code{SetCurvatureScaling()} defines the relative weighting between
  //  the two terms. In this example, we will set the propagation scaling to
  //  one and let the curvature scaling be an input argument. The larger the
  //  the curvature scaling parameter the smoother the resulting
  //  segmentation.  However, the curvature scaling parameter should not be
  //  set too large, as it will draw the contour away from the shape
  //  boundaries.
  //
  //  \index{itk::Shape\-Detection\-LevelSet\-Image\-Filter!SetPropagationScaling()}
  //  \index{itk::Segmentation\-Level\-Set\-Image\-Filter!SetPropagationScaling()}
  //  \index{itk::Shape\-Detection\-Level\-Set\-Image\-Filter!SetCurvatureScaling()}
  //  \index{itk::Segmentation\-Level\-Set\-Image\-Filter!SetCurvatureScaling()}
  //
  //  Software Guide : EndLatex

  const double curvatureScaling = std::stod(argv[9]);
  const double propagationScaling = std::stod(argv[10]);

  //  Software Guide : BeginCodeSnippet
  shapeDetection->SetPropagationScaling(propagationScaling);
  shapeDetection->SetCurvatureScaling(curvatureScaling);
  //  Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  Once activated, the level set evolution will stop if the convergence
  //  criteria or the maximum number of iterations is reached.  The
  //  convergence criteria are 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, it 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.
  //
  //  \index{itk::Shape\-Detection\-Level\-Set\-Image\-Filter!SetMaximumRMSError()}
  //  \index{itk::Shape\-Detection\-Level\-Set\-Image\-Filter!SetNumberOfIterations()}
  //  \index{itk::Segmentation\-Level\-Set\-Image\-Filter!SetMaximumRMSError()}
  //  \index{itk::Segmentation\-Level\-Set\-Image\-Filter!SetNumberOfIterations()}
  //
  //  Software Guide : EndLatex

  //  Software Guide : BeginCodeSnippet
  shapeDetection->SetMaximumRMSError(0.02);
  shapeDetection->SetNumberOfIterations(800);
  //  Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  The invocation of the \code{Update()} method on the writer triggers the
  //  execution of the pipeline.  As usual, the call is placed in a
  //  \code{try/catch} block should any errors occur or exceptions be thrown.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  try
  {
    writer->Update();
  }
  catch (const itk::ExceptionObject & excep)
  {
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
    return EXIT_FAILURE;
  }
  // Software Guide : EndCodeSnippet

  // Print out some useful information
  std::cout << std::endl;
  std::cout << "Max. no. iterations: "
            << shapeDetection->GetNumberOfIterations() << std::endl;
  std::cout << "Max. RMS error: " << shapeDetection->GetMaximumRMSError()
            << std::endl;
  std::cout << std::endl;
  std::cout << "No. elpased iterations: "
            << shapeDetection->GetElapsedIterations() << std::endl;
  std::cout << "RMS change: " << shapeDetection->GetRMSChange() << std::endl;

  writer4->Update();


  // The following writer type is used to save the output of the time-crossing
  // map in a file with apropiate pixel representation. The advantage of
  // saving this image in native format is that it can be used with a viewer
  // to help determine an appropriate threshold to be used on the output of
  // the fastmarching filter.
  //
  using InternalWriterType = itk::ImageFileWriter<InternalImageType>;

  auto mapWriter = InternalWriterType::New();
  mapWriter->SetInput(fastMarching->GetOutput());
  mapWriter->SetFileName("ShapeDetectionLevelSetFilterOutput4.mha");
  mapWriter->Update();

  auto speedWriter = InternalWriterType::New();
  speedWriter->SetInput(sigmoid->GetOutput());
  speedWriter->SetFileName("ShapeDetectionLevelSetFilterOutput3.mha");
  speedWriter->Update();

  auto gradientWriter = InternalWriterType::New();
  gradientWriter->SetInput(gradientMagnitude->GetOutput());
  gradientWriter->SetFileName("ShapeDetectionLevelSetFilterOutput2.mha");
  gradientWriter->Update();


  //  Software Guide : BeginLatex
  //
  //  Let's now run this example using as input the image
  //  \code{BrainProtonDensitySlice.png} provided in the directory
  //  \code{Examples/Data}. We can easily segment the major anatomical
  //  structures by providing seeds in the appropriate locations.
  //  Table~\ref{tab:ShapeDetectionLevelSetFilterOutput} presents the
  //  parameters used for some structures.  For all of the examples
  //  illustrated in this table, the propagation scaling was set to $1.0$, and
  //  the curvature scaling set to 0.05.
  //
  //  \begin{table}
  //  \begin{center}
  //  \begin{tabular}{|l|c|c|c|c|c|c|}
  //  \hline
  //  Structure    & Seed Index & Distance & $\sigma$ & $\alpha$ & $\beta$ &
  //  Output Image \\ \hline Left Ventricle  & $(81,114)$ & 5.0 & 1.0 & -0.5
  //  & 3.0  & First  in Figure \ref{fig:ShapeDetectionLevelSetFilterOutput2}
  //  \\  \hline Right Ventricle &
  //  $(99,114)$ & 5.0 & 1.0 & -0.5 & 3.0  & Second in Figure
  //  \ref{fig:ShapeDetectionLevelSetFilterOutput2} \\  \hline White matter &
  //  $(56, 92)$ & 5.0 & 1.0 & -0.3 & 2.0  & Third  in Figure
  //  \ref{fig:ShapeDetectionLevelSetFilterOutput2} \\  \hline Gray matter &
  //  $(40, 90)$ & 5.0 & 0.5 & -0.3 & 2.0  & Fourth in Figure
  //  \ref{fig:ShapeDetectionLevelSetFilterOutput2} \\  \hline \end{tabular}
  //  \end{center}
  //  \itkcaption[ShapeDetection example parameters]{Parameters used for
  //  segmenting some brain structures shown in
  //  Figure~\ref{fig:ShapeDetectionLevelSetFilterOutput} using the filter
  //  ShapeDetectionLevelSetFilter. All of them used a propagation
  //  scaling of $1.0$ and curvature scaling of
  //  $0.05$.\label{tab:ShapeDetectionLevelSetFilterOutput}}
  //  \end{table}
  //
  //  Figure~\ref{fig:ShapeDetectionLevelSetFilterOutput} presents the
  //  intermediate outputs of the pipeline illustrated in
  //  Figure~\ref{fig:ShapeDetectionCollaborationDiagram}. They are from left
  //  to right: the output of the anisotropic diffusion filter, the gradient
  //  magnitude of the smoothed image and the sigmoid of the gradient
  //  magnitude which is finally used as the edge potential for the
  //  ShapeDetectionLevelSetImageFilter.
  //
  //  Notice that in Figure~\ref{fig:ShapeDetectionLevelSetFilterOutput2} the
  //  segmented shapes are rounder than in
  //  Figure~\ref{fig:FastMarchingImageFilterOutput2} due to the effects of
  //  the curvature term in the driving equation. As with the previous
  //  example, segmentation of the gray matter is still problematic.
  //
  // \begin{figure} \center
  // \includegraphics[height=0.40\textheight]{BrainProtonDensitySlice}
  // \includegraphics[height=0.40\textheight]{ShapeDetectionLevelSetFilterOutput1Smoothing}
  // \includegraphics[height=0.40\textheight]{ShapeDetectionLevelSetFilterOutput1GradientMagnitude}
  // \includegraphics[height=0.40\textheight]{ShapeDetectionLevelSetFilterOutput1Sigmoid}
  // \itkcaption[ShapeDetectionLevelSetImageFilter intermediate output]{Images
  // generated by the segmentation process based on the
  // ShapeDetectionLevelSetImageFilter. From left to right and top to bottom:
  // input image to be segmented, image smoothed with an edge-preserving
  // smoothing filter, gradient magnitude of the smoothed image, sigmoid of
  // the gradient magnitude. This last image, the sigmoid, is used to compute
  // the speed term for the front propagation.}
  // \label{fig:ShapeDetectionLevelSetFilterOutput}
  // \end{figure}
  //
  //  A larger number of iterations is reguired for segmenting large
  //  structures since it takes longer for the front to propagate and cover
  //  the structure. This drawback can be easily mitigated by setting many
  //  seed points in the initialization of the
  //  FastMarchingImageFilter. This will generate an initial level
  //  set much closer in shape to the object to be segmented and hence
  //  require fewer iterations to fill and reach the edges of the anatomical
  //  structure.
  //
  // \begin{figure} \center
  // \includegraphics[width=0.24\textwidth]{ShapeDetectionLevelSetFilterOutput1}
  // \includegraphics[width=0.24\textwidth]{ShapeDetectionLevelSetFilterOutput2}
  // \includegraphics[width=0.24\textwidth]{ShapeDetectionLevelSetFilterOutput3}
  // \includegraphics[width=0.24\textwidth]{ShapeDetectionLevelSetFilterOutput4}
  // \itkcaption[ShapeDetectionLevelSetImageFilter segmentations]{Images
  // generated by the segmentation process based on the
  // ShapeDetectionLevelSetImageFilter. From left to right: segmentation of
  // the left ventricle, segmentation of the right ventricle, segmentation of
  // the white matter, attempt of segmentation of the gray matter.}
  // \label{fig:ShapeDetectionLevelSetFilterOutput2}
  // \end{figure}
  //
  //  Software Guide : EndLatex

  return EXIT_SUCCESS;
}