appli/CTBronchi/Vesselness.cxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5
   6 #include <chrono>
   7 #include <map>
   8 #include <itkImage.h>
   9 #include <itkImageFileReader.h>
  10 #include <itkImageFileWriter.h>
  11 #include <itkMinimumMaximumImageCalculator.h>
  12 #include <itkInvertIntensityImageFilter.h>
  13 #include <itkHessianRecursiveGaussianImageFilter.h>
  14 #include <itkHessian3DToVesselnessMeasureImageFilter.h>
  15
  16 // -------------------------------------------------------------------------
  17 const unsigned int Dim = 3;
  18 typedef short TPixel;
  19 typedef itk::NumericTraits< TPixel >::RealType TScalar;
  20 typedef itk::Image< TPixel, Dim >  TImage;
  21 typedef itk::Image< TScalar, Dim > TScalarImage;
  22
  23 // -------------------------------------------------------------------------
  24 double MeasureTime( itk::ProcessObject* f )
  25 {
  26   std::chrono::time_point< std::chrono::high_resolution_clock > s, e;
  27   std::chrono::duration< double > t;
  28   s = std::chrono::high_resolution_clock::now( );
  29   f->Update( );
  30   e = std::chrono::high_resolution_clock::now( );
  31   t = e - s;
  32   return( t.count( ) );
  33 }
  34
  35 // -------------------------------------------------------------------------
  36 template< class _TImagePtr >
  37 void ReadImage( _TImagePtr& image, const std::string& fname )
  38 {
  39   typedef typename _TImagePtr::ObjectType _TImage;
  40   typedef itk::ImageFileReader< _TImage > _TReader;
  41   typename _TReader::Pointer reader = _TReader::New( );
  42   reader->SetFileName( fname );
  43   double t = MeasureTime( reader );
  44   std::cout << "Read " << fname << " in " << t << " s" << std::endl;
  45   image = reader->GetOutput( );
  46   image->DisconnectPipeline( );
  47 }
  48
  49 // -------------------------------------------------------------------------
  50 template< class _TImage >
  51 void WriteImage( const _TImage* image, const std::string& fname )
  52 {
  53   typedef itk::ImageFileWriter< _TImage > _TWriter;
  54   typename _TWriter::Pointer writer = _TWriter::New( );
  55   writer->SetFileName( fname );
  56   writer->SetInput( image );
  57   double t = MeasureTime( writer );
  58   std::cout << "Wrote " << fname << " in " << t << " s" << std::endl;
  59 }
  60
  61 // -------------------------------------------------------------------------
  62 bool ParseArgs(
  63   std::map< std::string, std::string >& args, int argc, char* argv[]
  64   )
  65 {
  66   std::set< std::string > mandatory;
  67   mandatory.insert( "in" );
  68   mandatory.insert( "out" );
  69
  70   args[ "sigma" ] = "0.5";
  71   args[ "alpha1" ] = "0.5";
  72   args[ "alpha2" ] = "2";
  73
  74   int i = 1;
  75   while( i < argc )
  76   {
  77     std::string cmd = argv[ i ] + 1;
  78     args[ cmd ] = argv[ i + 1 ];
  79     i += 2;
  80
  81   } // elihw
  82
  83   bool complete = true;
  84   for( std::string t: mandatory )
  85     complete &= ( args.find( t ) != args.end( ) );
  86
  87   if( !complete )
  88   {
  89     std::cerr
  90       << "Usage: " << argv[ 0 ] << std::endl
  91       << "\t-in filename" << std::endl
  92       << "\t-out filename" << std::endl
  93       << "\t[-sigma val]" << std::endl
  94       << "\t[-alpha1 val]" << std::endl
  95       << "\t[-alpha2 val]" << std::endl;
  96     return( false );
  97
  98   } // fi
  99   return( true );
 100 }
 101
 102 // -------------------------------------------------------------------------
 103 int main( int argc, char* argv[] )
 104 {
 105   std::map< std::string, std::string > args;
 106   try
 107   {
 108     if( ParseArgs( args, argc, argv ) )
 109     {
 110       // Read input image
 111       TImage::Pointer input_image;
 112       ReadImage( input_image, args[ "in" ] );
 113
 114       // Vesselness
 115       typedef itk::MinimumMaximumImageCalculator< TImage >            TMinMax;
 116       typedef itk::InvertIntensityImageFilter< TImage >               TInverter;
 117       typedef itk::HessianRecursiveGaussianImageFilter< TImage >      THessian;
 118       typedef itk::Hessian3DToVesselnessMeasureImageFilter< TScalar > TVesselness;
 119
 120       TMinMax::Pointer minMax = TMinMax::New( );
 121       minMax->SetImage( input_image );
 122       minMax->Compute( );
 123
 124       TInverter::Pointer inverter = TInverter::New( );
 125       inverter->SetInput( input_image );
 126       inverter->SetMaximum( minMax->GetMaximum( ) );
 127       double t = MeasureTime( inverter );
 128       std::cout << "Inversion executed in " << t << " s" << std::endl;
 129
 130       THessian::Pointer hessian = THessian::New( );
 131       hessian->SetInput( inverter->GetOutput( ) );
 132       hessian->SetSigma(
 133         std::atof( args[ "sigma" ].c_str( ) )
 134         );
 135       t = MeasureTime( hessian );
 136       std::cout << "Hessian executed in " << t << " s" << std::endl;
 137
 138       TVesselness::Pointer vesselness = TVesselness::New( );
 139       vesselness->SetInput( hessian->GetOutput( ) );
 140       vesselness->SetAlpha1(
 141         std::atof( args[ "alpha1" ].c_str( ) )
 142         );
 143       vesselness->SetAlpha2(
 144         std::atof( args[ "alpha2" ].c_str( ) )
 145         );
 146       t = MeasureTime( vesselness );
 147       std::cout << "Vesselness executed in " << t << " s" << std::endl;
 148
 149       WriteImage( vesselness->GetOutput( ), args[ "out" ] );
 150     }
 151     else
 152       return( 1 );
 153   }
 154   catch( std::exception& err )
 155   {
 156     std::cerr
 157       << "===============================" << std::endl
 158       << "Error caught: " << std::endl
 159       << err.what( )
 160       << "===============================" << std::endl
 161       << std::endl;
 162     return( 1 );
 163
 164   } // yrt
 165   return( 0 );
 166 }
 167
 168 // eof - $RCSfile$