ITK  5.2.0
Insight Toolkit
Examples/Segmentation/ThresholdSegmentationLevelSetImageFilter.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: {ThresholdSegmentationLevelSetImageFilterWhiteMatter.png}
// ARGUMENTS: 60 116 5 150 180
// Software Guide : EndCommandLineArgs
// Software Guide : BeginCommandLineArgs
// INPUTS: {BrainProtonDensitySlice.png}
// OUTPUTS: {ThresholdSegmentationLevelSetImageFilterVentricle.png}
// ARGUMENTS: 81 112 5 210 250
// Software Guide : EndCommandLineArgs
// Software Guide : BeginCommandLineArgs
// INPUTS: {BrainProtonDensitySlice.png}
// OUTPUTS: {ThresholdSegmentationLevelSetImageFilterGrayMatter.png}
// ARGUMENTS: 107 69 5 180 210
// Software Guide : EndCommandLineArgs
// Software Guide : BeginLatex
//
// \index{itk::Threshold\-Segmentation\-Level\-Set\-Image\-Filter}
//
// The \doxygen{ThresholdSegmentationLevelSetImageFilter} is an extension of
// the threshold connected-component segmentation to the level set framework.
// The goal is to define a range of intensity values that classify the tissue
// type of interest and then base the propagation term on the level set
// equation for that intensity range. Using the level set approach, the
// smoothness of the evolving surface can be constrained to prevent some of
// the ``leaking'' that is common in connected-component schemes.
//
// The propagation term $P$ from Equation~\ref{eqn:LevelSetEquation} is
// calculated from the \code{FeatureImage} input $g$ with
// \code{UpperThreshold} $U$ and \code{LowerThreshold} $L$ according to the
// following formula.
//
// \begin{equation}
// \label{eqn:ThresholdSegmentationLevelSetImageFilterPropagationTerm}
// P(\mathbf{x}) = \left\{ \begin{array}{ll} g(\mathbf{x}) - L &
// \mbox{if $g(\mathbf{x}) < (U-L)/2 + L$} \\ U - g(\mathbf{x}) &
// \mbox{otherwise} \end{array} \right. \end{equation}
//
// Figure~\ref{fig:ThresholdSegmentationSpeedTerm} illustrates the propagation
// term function. Intensity values in $g$ between $L$ and $H$ yield positive
// values in $P$, while outside intensities yield negative values in $P$.
//
// \begin{figure} \center
// \includegraphics[width=0.8\textwidth]{ThresholdSegmentationLevelSetImageFilterCollaborationDiagram1}
// \itkcaption[ThresholdSegmentationLevelSetImageFilter collaboration
// diagram]{Collaboration diagram for the
// ThresholdSegmentationLevelSetImageFilter applied to a segmentation task.}
// \label{fig:ThresholdSegmentationLevelSetImageFilterDiagram}
// \end{figure}
//
// \begin{figure} \center
// \includegraphics[width=6.5cm]{ThresholdSegmentationLevelSetImageFilterFigure1}
// \itkcaption[Propagation term for threshold-based level set segmentation]
// {Propagation term for threshold-based level set segmentation.
// From
// Equation~\ref{eqn:ThresholdSegmentationLevelSetImageFilterPropagationTerm}.
// \label{fig:ThresholdSegmentationSpeedTerm}}
// \end{figure}
//
// The threshold segmentation filter expects two inputs. The first is an
// initial level set in the form of an \doxygen{Image}. The second input is
// the feature image $g$. For many applications, this filter requires little
// or no preprocessing of its input. Smoothing the input image is not
// usually required to produce reasonable solutions, though it may still be
// warranted in some cases.
//
// Figure~\ref{fig:ThresholdSegmentationLevelSetImageFilterDiagram} shows how
// the image processing pipeline is constructed. The initial surface is
// generated using the fast marching filter. The output of the segmentation
// filter is passed to a \doxygen{BinaryThresholdImageFilter} to create a
// binary representation of the segmented object. Let's start by including
// the appropriate header file.
//
// Software Guide : EndLatex
#include "itkImage.h"
// Software Guide : BeginCodeSnippet
// Software Guide : EndCodeSnippet
int
main(int argc, char * argv[])
{
if (argc < 8)
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " inputImage outputImage";
std::cerr << " seedX seedY InitialDistance";
std::cerr << " LowerThreshold";
std::cerr << " UpperThreshold";
std::cerr << " [CurvatureScaling == 1.0]";
std::cerr << std::endl;
return EXIT_FAILURE;
}
// Software Guide : BeginLatex
//
// We define the image type using a particular pixel type and
// dimension. In this case we will use 2D \code{float} images.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using InternalPixelType = float;
constexpr unsigned int Dimension = 2;
using InternalImageType = itk::Image<InternalPixelType, Dimension>;
// Software Guide : EndCodeSnippet
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);
ReaderType::Pointer reader = ReaderType::New();
WriterType::Pointer writer = WriterType::New();
reader->SetFileName(argv[1]);
writer->SetFileName(argv[2]);
// We now declare the type of the \doxygen{FastMarchingImageFilter} that
// will be used to generate the initial level set in the form of a distance
// map.
//
using FastMarchingFilterType =
FastMarchingFilterType::Pointer fastMarching =
FastMarchingFilterType::New();
// Software Guide : BeginLatex
//
// The following lines instantiate a
// ThresholdSegmentationLevelSetImageFilter using the \code{New()} method.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
using ThresholdSegmentationLevelSetImageFilterType =
InternalImageType>;
ThresholdSegmentationLevelSetImageFilterType::Pointer
thresholdSegmentation =
ThresholdSegmentationLevelSetImageFilterType::New();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// For the ThresholdSegmentationLevelSetImageFilter, scaling
// parameters are used to balance the influence of the propagation
// (inflation) and the curvature (surface smoothing) terms from
// Equation~\ref{eqn:LevelSetEquation}. The advection term is not used in
// this filter. Set the terms with methods \code{SetPropagationScaling()}
// and \code{SetCurvatureScaling()}. Both terms are set to 1.0 in this
// example.
//
// \index{itk::Threshold\-Segmentation\-Level\-Set\-Image\-Filter!SetPropagationScaling()}
// \index{itk::Segmentation\-Level\-Set\-Image\-Filter!SetPropagationScaling()}
// \index{itk::Threshold\-Segmentation\-Level\-Set\-Image\-Filter!SetCurvatureScaling()}
// \index{itk::Segmentation\-Level\-Set\-Image\-Filter!SetCurvatureScaling()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
thresholdSegmentation->SetPropagationScaling(1.0);
if (argc > 8)
{
thresholdSegmentation->SetCurvatureScaling(std::stod(argv[8]));
}
else
{
thresholdSegmentation->SetCurvatureScaling(1.0);
}
// Software Guide : EndCodeSnippet
// The level set solver will stop if the convergence criteria has been
// reached or if the maximum number of iterations has elasped. The
// convergence criteria is defined in terms of the root mean squared (RMS)
// change in the level set function. When RMS change for an iteration is
// below a user-specified threshold, the solution is considered to have
// converged.
thresholdSegmentation->SetMaximumRMSError(0.02);
thresholdSegmentation->SetNumberOfIterations(1200);
// thresholdSegmentation->SetMaximumRMSError( std::stod(argv[8]) );
// thresholdSegmentation->SetNumberOfIterations( std::stoi(argv[9]) );
// Software Guide : BeginLatex
//
// The convergence criteria \code{MaximumRMSError} and
// \code{MaximumIterations} are set as in previous examples. We now set
// the upper and lower threshold values $U$ and $L$, and the isosurface
// value to use in the initial model.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
thresholdSegmentation->SetUpperThreshold(::std::stod(argv[7]));
thresholdSegmentation->SetLowerThreshold(::std::stod(argv[6]));
thresholdSegmentation->SetIsoSurfaceValue(0.0);
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The filters are now connected in a pipeline indicated in
// Figure~\ref{fig:ThresholdSegmentationLevelSetImageFilterDiagram}.
// Remember that before calling \code{Update()} on the file writer object,
// the fast marching filter must be initialized with the seed points and
// the output from the reader object. See previous examples and the
// source code for this section for details.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
thresholdSegmentation->SetInput(fastMarching->GetOutput());
thresholdSegmentation->SetFeatureImage(reader->GetOutput());
thresholder->SetInput(thresholdSegmentation->GetOutput());
writer->SetInput(thresholder->GetOutput());
// Software Guide : EndCodeSnippet
//
// 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 he 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);
// The FastMarchingImageFilter requires the user to specify the size of the
// image to be produced as output. This is done using the
// \code{SetOutputRegion()} method. 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. Other image parameters such as Origin,
// Spacing and Direction are set in a similar manner.
// Software Guide : BeginLatex
//
// Invoking 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
{
reader->Update();
const InternalImageType * inputImage = reader->GetOutput();
fastMarching->SetOutputRegion(inputImage->GetBufferedRegion());
fastMarching->SetOutputSpacing(inputImage->GetSpacing());
fastMarching->SetOutputOrigin(inputImage->GetOrigin());
fastMarching->SetOutputDirection(inputImage->GetDirection());
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: "
<< thresholdSegmentation->GetNumberOfIterations() << std::endl;
std::cout << "Max. RMS error: "
<< thresholdSegmentation->GetMaximumRMSError() << std::endl;
std::cout << std::endl;
std::cout << "No. elpased iterations: "
<< thresholdSegmentation->GetElapsedIterations() << std::endl;
std::cout << "RMS change: " << thresholdSegmentation->GetRMSChange()
<< std::endl;
// We write out some intermediate images for debugging. These images can
// help tune parameters.
//
using InternalWriterType = itk::ImageFileWriter<InternalImageType>;
InternalWriterType::Pointer mapWriter = InternalWriterType::New();
mapWriter->SetInput(fastMarching->GetOutput());
mapWriter->SetFileName("fastMarchingImage.mha");
mapWriter->Update();
InternalWriterType::Pointer speedWriter = InternalWriterType::New();
speedWriter->SetInput(thresholdSegmentation->GetSpeedImage());
speedWriter->SetFileName("speedTermImage.mha");
speedWriter->Update();
// Software Guide : BeginLatex
//
// Let's run this application with the same data and parameters as the
// example given for \doxygen{ConnectedThresholdImageFilter} in
// Section~\ref{sec:ConnectedThreshold}. We will use a value of 5 as the
// initial distance of the surface from the seed points. The algorithm is
// relatively insensitive to this initialization. Compare the results in
// Figure~\ref{fig:ThresholdSegmentationLevelSetImageFilter} with those in
// Figure~\ref{fig:ConnectedThresholdOutput}. Notice how the smoothness
// constraint on the surface prevents leakage of the segmentation into
// both ventricles, but also localizes the segmentation to a smaller
// portion of the gray matter.
//
// \begin{figure}
// \includegraphics[width=0.24\textwidth]{BrainProtonDensitySlice}
// \includegraphics[width=0.24\textwidth]{ThresholdSegmentationLevelSetImageFilterWhiteMatter}
// \includegraphics[width=0.24\textwidth]{ThresholdSegmentationLevelSetImageFilterVentricle}
// \includegraphics[width=0.24\textwidth]{ThresholdSegmentationLevelSetImageFilterGrayMatter}
// \itkcaption[ThresholdSegmentationLevelSet segmentations]{Images
// generated by the segmentation process based on the
// ThresholdSegmentationLevelSetImageFilter. 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. The parameters used in this segmentations are presented in
// Table~\ref{tab:ThresholdSegmentationLevelSetImageFilter}.}
// \label{fig:ThresholdSegmentationLevelSetImageFilter} \end{figure}
//
// \begin{table}
// \begin{center}
// \begin{tabular}{|l|c|c|c|c|c|}
// \hline
// Structure & Seed Index & Lower & Upper & Output Image \\ \hline
// White matter & $(60,116)$ & 150 & 180 & Second from left \\ \hline
// Ventricle & $(81,112)$ & 210 & 250 & Third from left \\ \hline
// Gray matter & $(107,69)$ & 180 & 210 & Fourth from left \\ \hline
// \end{tabular}
// \itkcaption[ThresholdSegmentationLevelSet segmentation parameters]
// {Segmentation results using the
// ThresholdSegmentationLevelSetImageFilter for various seed points. The
// resulting images are shown in
// Figure~\ref{fig:ThresholdSegmentationLevelSetImageFilter}
// \label{tab:ThresholdSegmentationLevelSetImageFilter}.}\end{center}
// \end{table}
//
// Software Guide : EndLatex
return EXIT_SUCCESS;
}
itk::DataObject::Update
virtual void Update()
itk::ThresholdSegmentationLevelSetImageFilter
Segments structures in images based on intensity values.
Definition: itkThresholdSegmentationLevelSetImageFilter.h:84
itk::BinaryThresholdImageFilter
Binarize an input image by thresholding.
Definition: itkBinaryThresholdImageFilter.h:139
itkImageFileReader.h
itkImage.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::ImageFileWriter
Writes image data to a single file.
Definition: itkImageFileWriter.h:87
itk::FastMarchingImageFilter
Solve an Eikonal equation using Fast Marching.
Definition: itkFastMarchingImageFilter.h:136
itkImageFileWriter.h
itkThresholdSegmentationLevelSetImageFilter.h
itkZeroCrossingImageFilter.h
itkBinaryThresholdImageFilter.h
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:86
itk::ImageSource::GetOutput
OutputImageType * GetOutput()
itk::GTest::TypedefsAndConstructors::Dimension2::Dimension
constexpr unsigned int Dimension
Definition: itkGTestTypedefsAndConstructors.h:44