KWHelloWorldExample-alex.cxx
From KitwarePublic
Jump to navigationJump to search
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#ifdef __BORLANDC__
#define ITK_LEAN_AND_MEAN
#endif
#include "itkImage.h"
#include "itkGeodesicActiveContourLevelSetImageFilter.h"
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
int main( int argc, char *argv[] )
{
if( argc < 12 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " inputImageFile outputImageFile";
std::cerr << " Sigma SigmoidAlpha SigmoidBeta";
std::cerr << " PropagationScaling CurvatureScaling";
std::cerr << " AdvectionScaling maxRMSError maxIterations";
std::cerr << " initialThreshold" << std::endl;
return 1;
}
// We 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.
typedef float InternalPixelType;
const unsigned char Dimension = 3;
typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
typedef unsigned char OutputPixelType;
typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
// This is the initial Threshold to determine the initial level set.
typedef itk::BinaryThresholdImageFilter<
InternalImageType,
InternalImageType > ThresholdingFilterType2;
ThresholdingFilterType2::Pointer inThresholder = ThresholdingFilterType2::New();
inThresholder->SetLowerThreshold( 0 );
inThresholder->SetUpperThreshold( atoi(argv[11]) );
inThresholder->SetOutsideValue( 0 );
inThresholder->SetInsideValue( 255 );
// Define the reader and writer.
typedef itk::ImageFileReader< InternalImageType > ReaderType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;
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.
typedef itk::RescaleIntensityImageFilter<
InternalImageType,
OutputImageType > CastFilterType;
// The types of the
// GradientMagnitudeRecursiveGaussianImageFilter and
// SigmoidImageFilter are instantiated using the internal image
// type.
//
typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<
InternalImageType,
InternalImageType > GradientFilterType;
GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();
// 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 = atof( argv[3] );
gradientMagnitude->SetSigma( sigma );
typedef itk::SigmoidImageFilter<
InternalImageType,
InternalImageType > SigmoidFilterType;
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 );
// 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 = atof( argv[4] );
const double beta = atof( argv[5] );
sigmoid->SetAlpha( alpha );
sigmoid->SetBeta( beta );
typedef itk::GeodesicActiveContourLevelSetImageFilter< InternalImageType,
InternalImageType > GeodesicActiveContourFilterType;
GeodesicActiveContourFilterType::Pointer geodesicActiveContour =
GeodesicActiveContourFilterType::New();
// For the GeodesicActiveContourLevelSetImageFilter, 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.
const double propagationScaling = atof( argv[6] );
const double curvatureScaling = atof( argv[7] );
const double advectionScaling = atof( argv[8] );
geodesicActiveContour->SetPropagationScaling( propagationScaling );
geodesicActiveContour->SetCurvatureScaling( curvatureScaling );
geodesicActiveContour->SetAdvectionScaling( advectionScaling );
// Once activated 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.
// !!!!!
const double maxRMSError = atof( argv[9] );
const unsigned int maxIter = atoi( argv[10] );
geodesicActiveContour->SetMaximumRMSError( maxRMSError );
geodesicActiveContour->SetNumberOfIterations( maxIter );
// The following lines instantiate the thresholding filter that will
// process the final level set at the output of the
// GeodesicActiveContourLevelSetImageFilter.
typedef itk::BinaryThresholdImageFilter<
InternalImageType,
OutputImageType > ThresholdingFilterType;
ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
thresholder->SetLowerThreshold( -1000.0 );
thresholder->SetUpperThreshold( 0.0 );
thresholder->SetOutsideValue( 0 );
thresholder->SetInsideValue( 255 );
// Configure the pipeline.
gradientMagnitude->SetInput( reader->GetOutput() );
sigmoid->SetInput( gradientMagnitude->GetOutput() );
inThresholder->SetInput( reader->GetOutput() );
geodesicActiveContour->SetInput( inThresholder->GetOutput() );
geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );
thresholder->SetInput( geodesicActiveContour->GetOutput() );
writer->SetInput( thresholder->GetOutput() );
// 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 caster2 = CastFilterType::New();
CastFilterType::Pointer caster3 = CastFilterType::New();
CastFilterType::Pointer caster4 = CastFilterType::New();
WriterType::Pointer writer2 = WriterType::New();
WriterType::Pointer writer3 = WriterType::New();
WriterType::Pointer writer4 = WriterType::New();
caster2->SetInput( gradientMagnitude->GetOutput() );
writer2->SetInput( caster2->GetOutput() );
writer2->SetFileName("GeodesicActiveContourImageFilterOutput-grad.tif");
caster2->SetOutputMinimum( 0 );
caster2->SetOutputMaximum( 255 );
writer2->Update();
caster2->ReleaseDataFlagOn();
writer2->ReleaseDataFlagOn();
caster3->SetInput( sigmoid->GetOutput() );
writer3->SetInput( caster3->GetOutput() );
writer3->SetFileName("GeodesicActiveContourImageFilterOutput-sigmoid.tif");
caster3->SetOutputMinimum( 0 );
caster3->SetOutputMaximum( 255 );
writer3->Update();
caster3->ReleaseDataFlagOn();
writer3->ReleaseDataFlagOn();
caster4->SetInput( inThresholder->GetOutput() );
writer4->SetInput( caster4->GetOutput() );
writer4->SetFileName("GeodesicActiveContourImageFilterOutput-inThresh.tif");
caster4->SetOutputMinimum( 0 );
caster4->SetOutputMaximum( 255 );
writer4->Update();
caster4->ReleaseDataFlagOn();
writer4->ReleaseDataFlagOn();
// 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;
}
// 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. elapsed iterations: " << geodesicActiveContour->GetElapsedIterations() << std::endl;
std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() << std::endl;
return 0;
}
