ITK  5.1.0
Insight Toolkit
Examples/Segmentation/CurvesLevelSetImageFilter.cxx
/*=========================================================================
*
* Copyright Insight Software Consortium
*
* 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 : BeginLatex
//
// WORK IN PROGRESS: THIS WAS TAKEN FROM THE GEODESIC ACTIVE CONTOURS.
// IT NEED TO BE REWORKED TO MATCH THE CURVESLEVELSET FILTER.
//
// The use of the \doxygen{CurvesLevelSetImageFilter} is
// illustrated in the following example. The implementation of this filter in
// ITK is based on the paper by Caselles \cite{Caselles1997}. This
// implementation extends the functionality of the
// \doxygen{ShapeDetectionLevelSetImageFilter} by the addition of a third
// avection term which attracts the level set to the object boundaries.
//
// CurvesLevelSetImageFilter 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 an edge
// potential image that basically follows the same rules used for the
// ShapeDetectionLevelSetImageFilter discussed in
// Section~\ref{sec:ShapeDetectionLevelSetFilter}. The configuration of this
// example is quite similar to the example on the use of the
// ShapeDetectionLevelSetImageFilter. We omit most of the redundant
// description. A look at the code will reveal the great degree of similarity
// between both examples.
//
// \begin{figure} \center
// \includegraphics[width=\textwidth]{CurvessCollaborationDiagram1}
// \itkcaption[CurvesLevelSetImageFilter collaboration
// diagram]{Collaboration diagram for the CurvesLevelSetImageFilter
// applied to a segmentation task.}
// \label{fig:CurvessCollaborationDiagram}
// \end{figure}
//
// Figure~\ref{fig:CurvessCollaborationDiagram} shows the major
// components involved in the application of the
// CurvesLevelSetImageFilter to a segmentation task.
// This pipeline is quite similar to the one used by the
// ShapeDetectionLevelSetImageFilter in
// section~\ref{sec:ShapeDetectionLevelSetFilter}.
//
// The pipeline involves a first stage of smoothing using the
// \doxygen{CurvatureAnisotropicDiffusionImageFilter}. The smoothed image is
// passed as the input to 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 a
// \doxygen{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
// CurvesLevelSetImageFilter.
//
// Finally, the level set generated by the
// CurvesLevelSetImageFilter is passed to a
// \doxygen{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
// Software Guide : EndCodeSnippet
int
main(int argc, char * argv[])
{
if (argc < 10)
{
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 << " PropagationScaling" << std::endl;
return EXIT_FAILURE;
}
// Software Guide : BeginLatex
//
// 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.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using InternalPixelType = float;
constexpr unsigned int Dimension = 2;
using InternalImageType = itk::Image<InternalPixelType, Dimension>;
// Software Guide : EndCodeSnippet
// The following lines instantiate the thresholding filter that will
// process the final level set at the output of the
// CurvesLevelSetImageFilter.
//
using OutputPixelType = unsigned char;
using OutputImageType = itk::Image<OutputPixelType, Dimension>;
using ThresholdingFilterType =
ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
thresholder->SetLowerThreshold(-1000.0);
thresholder->SetUpperThreshold(0.0);
thresholder->SetOutsideValue(0);
thresholder->SetInsideValue(255);
// We instantiate reader and writer types in the following lines.
//
ReaderType::Pointer reader = ReaderType::New();
WriterType::Pointer writer = WriterType::New();
reader->SetFileName(argv[1]);
writer->SetFileName(argv[2]);
// The RescaleIntensityImageFilter type is declared below. This filter will
// renormalize image before sending them to writers.
//
using CastFilterType =
// The \doxygen{CurvatureAnisotropicDiffusionImageFilter} type is
// instantiated using the internal image type.
//
using SmoothingFilterType =
SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
// The types of the
// GradientMagnitudeRecursiveGaussianImageFilter and
// SigmoidImageFilter are instantiated using the internal image
// type.
//
using GradientFilterType =
InternalImageType>;
using SigmoidFilterType =
GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();
SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
// 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
// the \code{FastMarchingImageFilter}. Additional details on the user of the
// \doxygen{SigmoidImageFilter} are presented in
// section~\ref{sec:IntensityNonLinearMapping}.
sigmoid->SetOutputMinimum(0.0);
sigmoid->SetOutputMaximum(1.0);
// 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 =
// Next we construct one filter of this class using the \code{New()}
// method.
//
FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();
// Software Guide : BeginLatex
//
// In the following lines we instantiate the type of the
// CurvesLevelSetImageFilter and create an object of this
// type using the \code{New()} method.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using CurvesFilterType =
CurvesFilterType::Pointer geodesicActiveContour = CurvesFilterType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// For the CurvesLevelSetImageFilter, scaling parameters
// are used to trade off between the propagation (inflation), the
// curvature (smoothing) and the advection terms. These parameters are set
// using methods \code{SetPropagationScaling()},
// \code{SetCurvatureScaling()} and \code{SetAdvectionScaling()}. In this
// example, we will set the curvature and advection scales to one and let
// the propagation scale be a command-line argument.
//
// \index{itk::Geodesic\-Active\-Contour\-LevelSet\-Image\-Filter!SetPropagationScaling()}
// \index{itk::Segmentation\-Level\-Set\-Image\-Filter!SetPropagationScaling()}
// \index{itk::Geodesic\-Active\-Contour\-LevelSet\-Image\-Filter!SetCurvatureScaling()}
// \index{itk::Segmentation\-Level\-Set\-Image\-Filter!SetCurvatureScaling()}
// \index{itk::Geodesic\-Active\-Contour\-LevelSet\-Image\-Filter!SetAdvectionScaling()}
// \index{itk::Segmentation\-Level\-Set\-Image\-Filter!SetAdvectionScaling()}
//
// Software Guide : EndLatex
const double propagationScaling = std::stod(argv[9]);
// Software Guide : BeginCodeSnippet
geodesicActiveContour->SetPropagationScaling(propagationScaling);
geodesicActiveContour->SetCurvatureScaling(1.0);
geodesicActiveContour->SetAdvectionScaling(1.0);
// Software Guide : EndCodeSnippet
// 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.
geodesicActiveContour->SetMaximumRMSError(0.02);
geodesicActiveContour->SetNumberOfIterations(800);
// Software Guide : BeginLatex
//
// The filters are now connected in a pipeline indicated in
// Figure~\ref{fig:CurvessCollaborationDiagram} using the
// following lines:
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
smoothing->SetInput(reader->GetOutput());
gradientMagnitude->SetInput(smoothing->GetOutput());
sigmoid->SetInput(gradientMagnitude->GetOutput());
geodesicActiveContour->SetInput(fastMarching->GetOutput());
geodesicActiveContour->SetFeatureImage(sigmoid->GetOutput());
thresholder->SetInput(geodesicActiveContour->GetOutput());
writer->SetInput(thresholder->GetOutput());
// Software Guide : EndCodeSnippet
// The CurvatureAnisotropicDiffusionImageFilter requires a couple of
// parameter 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}.
smoothing->SetTimeStep(0.125);
smoothing->SetNumberOfIterations(5);
smoothing->SetConductanceParameter(9.0);
// 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}.
const double sigma = std::stod(argv[6]);
gradientMagnitude->SetSigma(sigma);
// The SigmoidImageFilter requires two parameters that define the linear
// transformation to be applied to the sigmoid argument. This parameters
// have been discussed in Sections~\ref{sec:IntensityNonLinearMapping} and
// \ref{sec:FastMarchingImageFilter}.
const double alpha = std::stod(argv[7]);
const double beta = std::stod(argv[8]);
sigmoid->SetAlpha(alpha);
sigmoid->SetBeta(beta);
// 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();
seedPosition[0] = std::stoi(argv[3]);
seedPosition[1] = std::stoi(argv[4]);
// Nodes are created as stack variables and initialized with a value and an
// \doxygen{Index} position. Note that here we assign the value of minus 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 here 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 she want the initial
// contour to be.
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);
// 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.
//
CastFilterType::Pointer caster1 = CastFilterType::New();
CastFilterType::Pointer caster2 = CastFilterType::New();
CastFilterType::Pointer caster3 = CastFilterType::New();
CastFilterType::Pointer caster4 = CastFilterType::New();
WriterType::Pointer writer1 = WriterType::New();
WriterType::Pointer writer2 = WriterType::New();
WriterType::Pointer writer3 = WriterType::New();
WriterType::Pointer writer4 = WriterType::New();
caster1->SetInput(smoothing->GetOutput());
writer1->SetInput(caster1->GetOutput());
writer1->SetFileName("CurvesImageFilterOutput1.png");
caster1->SetOutputMinimum(0);
caster1->SetOutputMaximum(255);
writer1->Update();
caster2->SetInput(gradientMagnitude->GetOutput());
writer2->SetInput(caster2->GetOutput());
writer2->SetFileName("CurvesImageFilterOutput2.png");
caster2->SetOutputMinimum(0);
caster2->SetOutputMaximum(255);
writer2->Update();
caster3->SetInput(sigmoid->GetOutput());
writer3->SetInput(caster3->GetOutput());
writer3->SetFileName("CurvesImageFilterOutput3.png");
caster3->SetOutputMinimum(0);
caster3->SetOutputMaximum(255);
writer3->Update();
caster4->SetInput(fastMarching->GetOutput());
writer4->SetInput(caster4->GetOutput());
writer4->SetFileName("CurvesImageFilterOutput4.png");
caster4->SetOutputMinimum(0);
caster4->SetOutputMaximum(255);
// 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(reader->GetOutput()->GetBufferedRegion().GetSize());
// 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 (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: " << geodesicActiveContour->GetNumberOfIterations()
<< std::endl;
std::cout << "Max. RMS error: " << geodesicActiveContour->GetMaximumRMSError()
<< std::endl;
std::cout << std::endl;
std::cout << "No. elpased iterations: "
<< geodesicActiveContour->GetElapsedIterations() << std::endl;
std::cout << "RMS change: " << geodesicActiveContour->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>;
InternalWriterType::Pointer mapWriter = InternalWriterType::New();
mapWriter->SetInput(fastMarching->GetOutput());
mapWriter->SetFileName("CurvesImageFilterOutput4.mha");
mapWriter->Update();
InternalWriterType::Pointer speedWriter = InternalWriterType::New();
speedWriter->SetInput(sigmoid->GetOutput());
speedWriter->SetFileName("CurvesImageFilterOutput3.mha");
speedWriter->Update();
InternalWriterType::Pointer gradientWriter = InternalWriterType::New();
gradientWriter->SetInput(gradientMagnitude->GetOutput());
gradientWriter->SetFileName("CurvesImageFilterOutput2.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:CurvesImageFilterOutput2} presents the
// parameters used for some structures.
//
// \begin{table}
// \begin{center}
// \begin{tabular}{|l|c|c|c|c|c|c|c|c|}
// \hline
// Structure & Seed Index & Distance & $\sigma$ &
// $\alpha$ & $\beta$ & Propag. & Output Image \\ \hline
// Left Ventricle & $(81,114)$ & 5.0 & 1.0 & -0.5 & 3.0 & 2.0 & First \\ \hline
// Right Ventricle & $(99,114)$ & 5.0 & 1.0 & -0.5 & 3.0 & 2.0 & Second \\ \hline
// White matter & $(56, 92)$ & 5.0 & 1.0 & -0.3 & 2.0 & 10.0 & Third \\ \hline
// Gray matter & $(40, 90)$ & 5.0 & 0.5 & -0.3 & 2.0 & 10.0 & Fourth \\ \hline
// \end{tabular}
// \end{center}
// \itkcaption[Curves segmentation example parameters]{Parameters used
// for segmenting some brain structures shown in
// Figure~\ref{fig:CurvesImageFilterOutput2} using the filter
// CurvesLevelSetImageFilter.
// \label{tab:CurvesImageFilterOutput2}}
// \end{table}
//
// Figure~\ref{fig:CurvesImageFilterOutput} presents the
// intermediate outputs of the pipeline illustrated in
// Figure~\ref{fig:CurvessCollaborationDiagram}. 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
// CurvesLevelSetImageFilter.
//
// \begin{figure} \center
// \includegraphics[height=0.40\textheight]{BrainProtonDensitySlice}
// \includegraphics[height=0.40\textheight]{CurvesImageFilterOutput1}
// \includegraphics[height=0.40\textheight]{CurvesImageFilterOutput2}
// \includegraphics[height=0.40\textheight]{CurvesImageFilterOutput3}
// \itkcaption[CurvesLevelSetImageFilter intermediate
// output]{Images generated by the segmentation process based on the
// CurvesLevelSetImageFilter. 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:CurvesImageFilterOutput} \end{figure}
//
// Segmentations of the main brain structures are presented in
// Figure~\ref{fig:CurvesImageFilterOutput2}. The results
// are quite similar to those obtained with the
// ShapeDetectionLevelSetImageFilter in
// Section~\ref{sec:ShapeDetectionLevelSetFilter}.
//
// Note that a relatively larger propagation scaling value was required to
// segment the white matter. This is due to two factors: the lower
// contrast at the border of the white matter and the complex shape of the
// structure. Unfortunately the optimal value of these scaling parameters
// can only be determined by experimentation. In a real application we
// could imagine an interactive mechanism by which a user supervises the
// contour evolution and adjusts these parameters accordingly.
//
// \begin{figure} \center
// \includegraphics[width=0.24\textwidth]{CurvesImageFilterOutput5}
// \includegraphics[width=0.24\textwidth]{CurvesImageFilterOutput6}
// \includegraphics[width=0.24\textwidth]{CurvesImageFilterOutput7}
// \includegraphics[width=0.24\textwidth]{CurvesImageFilterOutput8}
// \itkcaption[CurvesImageFilter segmentations]{Images generated by the
// segmentation process based on the CurvesImageFilter. 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:CurvesImageFilterOutput2}
// \end{figure}
//
// Software Guide : EndLatex
return EXIT_SUCCESS;
}
itkSigmoidImageFilter.h
itk::BinaryThresholdImageFilter
Binarize an input image by thresholding.
Definition: itkBinaryThresholdImageFilter.h:137
itk::CurvesLevelSetImageFilter
Segments structures in images based on user supplied edge potential map.
Definition: itkCurvesLevelSetImageFilter.h:102
itkImageFileReader.h
itk::SigmoidImageFilter
Computes the sigmoid function pixel-wise.
Definition: itkSigmoidImageFilter.h:148
itkCurvesLevelSetImageFilter.h
itkFastMarchingImageFilter.h
itk::ImageFileReader
Data source that reads image data from a single file.
Definition: itkImageFileReader.h:75
itk::GTest::TypedefsAndConstructors::Dimension2::IndexType
ImageBaseType::IndexType IndexType
Definition: itkGTestTypedefsAndConstructors.h:50
itk::CurvatureAnisotropicDiffusionImageFilter
Definition: itkCurvatureAnisotropicDiffusionImageFilter.h:58
itkCurvatureAnisotropicDiffusionImageFilter.h
itk::ImageFileWriter
Writes image data to a single file.
Definition: itkImageFileWriter.h:87
itk::FastMarchingImageFilter
Solve an Eikonal equation using Fast Marching.
Definition: itkFastMarchingImageFilter.h:107
itk::GradientMagnitudeRecursiveGaussianImageFilter
Computes the Magnitude of the Gradient of an image by convolution with the first derivative of a Gaus...
Definition: itkGradientMagnitudeRecursiveGaussianImageFilter.h:49
itkRescaleIntensityImageFilter.h
itkImageFileWriter.h
itk::RescaleIntensityImageFilter
Applies a linear transformation to the intensity levels of the input Image.
Definition: itkRescaleIntensityImageFilter.h:153
itkBinaryThresholdImageFilter.h
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:86
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44
itkGradientMagnitudeRecursiveGaussianImageFilter.h