KWHelloWorldExample-alex.cxx

#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;
}

---

Complete Setup: [Welcome | Site Map]