appli/CTBronchi/Vesselness.cxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5
   6 #include <string>
   7 #include <tclap/CmdLine.h>
   8 #include <itkImage.h>
   9 #include <itkMinimumMaximumImageCalculator.h>
  10 #include <itkInvertIntensityImageFilter.h>
  11 #include <itkHessianRecursiveGaussianImageFilter.h>
  12 #include <itkHessian3DToVesselnessMeasureImageFilter.h>
  13 #include "Functions.h"
  14
  15 // -------------------------------------------------------------------------
  16 const unsigned int Dim = 3;
  17 typedef short                                  TPixel;
  18 typedef itk::NumericTraits< TPixel >::RealType TScalar;
  19 typedef itk::Image< TPixel, Dim >              TImage;
  20 typedef itk::Image< TScalar, Dim >             TScalarImage;
  21
  22 // -------------------------------------------------------------------------
  23 int main( int argc, char* argv[] )
  24 {
  25   typedef TCLAP::ValueArg< std::string > _TStringArg;
  26   typedef TCLAP::ValueArg< TScalar > _TRealArg;
  27
  28   // Parse input line
  29   _TStringArg in( "i", "input", "Input image", true, "", "file" );
  30   _TStringArg out( "o", "output", "Output image", true, "", "file" );
  31   _TRealArg s( "s", "sigma", "Sigma", false, 0.5, "value (0.5)" );
  32   _TRealArg a1( "a", "alpha1", "Alpha1", false, 0.5, "value (0.5)" );
  33   _TRealArg a2( "b", "alpha2", "Alpha2", false, 2, "value (2)" );
  34   try
  35   {
  36     TCLAP::CmdLine cmd( "Vesselness computation", ' ', "1.0.0" );
  37     cmd.add( a2 );
  38     cmd.add( a1 );
  39     cmd.add( s );
  40     cmd.add( out );
  41     cmd.add( in );
  42     cmd.parse( argc, argv );
  43   }
  44   catch( TCLAP::ArgException& err )
  45   {
  46     std::cerr
  47       << "===============================" << std::endl
  48       << "Error caught: " << std::endl
  49       << err.error( ) << " " << err.argId( )
  50       << "===============================" << std::endl
  51       << std::endl;
  52     return( 1 );
  53
  54   } // yrt
  55
  56   try
  57   {
  58     // Read image
  59     TImage::Pointer input_image;
  60     CTBronchi::ReadImage( input_image, in.getValue( ) );
  61
  62     // Min-max
  63     typedef itk::MinimumMaximumImageCalculator< TImage > _TMinMax;
  64     _TMinMax::Pointer minMax = _TMinMax::New( );
  65     minMax->SetImage( input_image );
  66     minMax->Compute( );
  67
  68     // Invert intensities
  69     typedef itk::InvertIntensityImageFilter< TImage > _TInverter;
  70     _TInverter::Pointer inverter = _TInverter::New( );
  71     inverter->SetInput( input_image );
  72     inverter->SetMaximum( minMax->GetMaximum( ) );
  73     double t = CTBronchi::MeasureTime( inverter );
  74     std::cout << "Inversion executed in " << t << " s" << std::endl;
  75
  76     // Compute hessian image
  77     typedef itk::HessianRecursiveGaussianImageFilter< TImage > _THessian;
  78     _THessian::Pointer hessian = _THessian::New( );
  79     hessian->SetInput( inverter->GetOutput( ) );
  80     hessian->SetSigma( s.getValue( ) );
  81     t = CTBronchi::MeasureTime( hessian );
  82     std::cout << "Hessian executed in " << t << " s" << std::endl;
  83
  84     // Vesselness
  85     typedef
  86       itk::Hessian3DToVesselnessMeasureImageFilter< TScalar >
  87       _TVesselness;
  88     _TVesselness::Pointer vesselness = _TVesselness::New( );
  89     vesselness->SetInput( hessian->GetOutput( ) );
  90     vesselness->SetAlpha1( a1.getValue( ) );
  91     vesselness->SetAlpha2( a2.getValue( ) );
  92     t = CTBronchi::MeasureTime( vesselness );
  93     std::cout << "Vesselness executed in " << t << " s" << std::endl;
  94
  95     // Write result
  96     CTBronchi::WriteImage( vesselness->GetOutput( ), out.getValue( ) );
  97   }
  98   catch( std::exception& err )
  99   {
 100     std::cerr
 101       << "===============================" << std::endl
 102       << "Error caught: " << std::endl
 103       << err.what( )
 104       << "===============================" << std::endl
 105       << std::endl;
 106     return( 1 );
 107
 108   } // yrt
 109   return( 0 );
 110 }
 111
 112 // eof - $RCSfile$