]> Creatis software - FrontAlgorithms.git/blobdiff - appli/CTBronchi/Vesselness.cxx
...
[FrontAlgorithms.git] / appli / CTBronchi / Vesselness.cxx
index ca7c99c358c1df14d63fd1db2794e34a1f452933..3b0661824de9cb3841800bbd7f067c98310b8367 100644 (file)
 // @email florez-l@javeriana.edu.co
 // =========================================================================
 
-#include <chrono>
-#include <map>
+#include <string>
+#include <tclap/CmdLine.h>
 #include <itkImage.h>
-#include <itkImageFileReader.h>
-#include <itkImageFileWriter.h>
 #include <itkMinimumMaximumImageCalculator.h>
 #include <itkInvertIntensityImageFilter.h>
 #include <itkHessianRecursiveGaussianImageFilter.h>
 #include <itkHessian3DToVesselnessMeasureImageFilter.h>
+#include "Functions.h"
 
 // -------------------------------------------------------------------------
 const unsigned int Dim = 3;
-typedef short TPixel;
+typedef short                                  TPixel;
 typedef itk::NumericTraits< TPixel >::RealType TScalar;
-typedef itk::Image< TPixel, Dim >  TImage;
-typedef itk::Image< TScalar, Dim > TScalarImage;
+typedef itk::Image< TPixel, Dim >              TImage;
+typedef itk::Image< TScalar, Dim >             TScalarImage;
 
 // -------------------------------------------------------------------------
-double MeasureTime( itk::ProcessObject* f )
-{
-  std::chrono::time_point< std::chrono::high_resolution_clock > s, e;
-  std::chrono::duration< double > t;
-  s = std::chrono::high_resolution_clock::now( );
-  f->Update( );
-  e = std::chrono::high_resolution_clock::now( );
-  t = e - s;
-  return( t.count( ) );
-}
-
-// -------------------------------------------------------------------------
-template< class _TImagePtr >
-void ReadImage( _TImagePtr& image, const std::string& fname )
-{
-  typedef typename _TImagePtr::ObjectType _TImage;
-  typedef itk::ImageFileReader< _TImage > _TReader;
-  typename _TReader::Pointer reader = _TReader::New( );
-  reader->SetFileName( fname );
-  double t = MeasureTime( reader );
-  std::cout << "Read " << fname << " in " << t << " s" << std::endl;
-  image = reader->GetOutput( );
-  image->DisconnectPipeline( );
-}
-
-// -------------------------------------------------------------------------
-template< class _TImage >
-void WriteImage( const _TImage* image, const std::string& fname )
-{
-  typedef itk::ImageFileWriter< _TImage > _TWriter;
-  typename _TWriter::Pointer writer = _TWriter::New( );
-  writer->SetFileName( fname );
-  writer->SetInput( image );
-  double t = MeasureTime( writer );
-  std::cout << "Wrote " << fname << " in " << t << " s" << std::endl;
-}
-
-// -------------------------------------------------------------------------
-bool ParseArgs(
-  std::map< std::string, std::string >& args, int argc, char* argv[]
-  )
+int main( int argc, char* argv[] )
 {
-  std::set< std::string > mandatory;
-  mandatory.insert( "in" );
-  mandatory.insert( "out" );
-
-  args[ "sigma" ] = "0.5";
-  args[ "alpha1" ] = "0.5";
-  args[ "alpha2" ] = "2";
-
-  int i = 1;
-  while( i < argc )
+  typedef TCLAP::ValueArg< std::string > _TStringArg;
+  typedef TCLAP::ValueArg< TScalar > _TRealArg;
+
+  // Parse input line
+  _TStringArg in( "i", "input", "Input image", true, "", "file" );
+  _TStringArg out( "o", "output", "Output image", true, "", "file" );
+  _TRealArg s( "s", "sigma", "Sigma", false, 0.5, "value (0.5)" );
+  _TRealArg a1( "a", "alpha1", "Alpha1", false, 0.5, "value (0.5)" );
+  _TRealArg a2( "b", "alpha2", "Alpha2", false, 2, "value (2)" );
+  try
   {
-    std::string cmd = argv[ i ] + 1;
-    args[ cmd ] = argv[ i + 1 ];
-    i += 2;
-
-  } // elihw
-
-  bool complete = true;
-  for( std::string t: mandatory )
-    complete &= ( args.find( t ) != args.end( ) );
-
-  if( !complete )
+    TCLAP::CmdLine cmd( "Vesselness computation", ' ', "1.0.0" );
+    cmd.add( a2 );
+    cmd.add( a1 );
+    cmd.add( s );
+    cmd.add( out );
+    cmd.add( in );
+    cmd.parse( argc, argv );
+  }
+  catch( TCLAP::ArgException& err )
   {
     std::cerr
-      << "Usage: " << argv[ 0 ] << std::endl
-      << "\t-in filename" << std::endl
-      << "\t-out filename" << std::endl
-      << "\t[-sigma val]" << std::endl
-      << "\t[-alpha1 val]" << std::endl
-      << "\t[-alpha2 val]" << std::endl;
-    return( false );
+      << "===============================" << std::endl
+      << "Error caught: " << std::endl
+      << err.error( ) << " " << err.argId( )
+      << "===============================" << std::endl
+      << std::endl;
+    return( 1 );
 
-  } // fi
-  return( true );
-}
+  } // yrt
 
-// -------------------------------------------------------------------------
-int main( int argc, char* argv[] )
-{
-  std::map< std::string, std::string > args;
   try
   {
-    if( ParseArgs( args, argc, argv ) )
-    {
-      // Read input image
-      TImage::Pointer input_image;
-      ReadImage( input_image, args[ "in" ] );
-
-      // Vesselness
-      typedef itk::MinimumMaximumImageCalculator< TImage >            TMinMax;
-      typedef itk::InvertIntensityImageFilter< TImage >               TInverter;
-      typedef itk::HessianRecursiveGaussianImageFilter< TImage >      THessian;
-      typedef itk::Hessian3DToVesselnessMeasureImageFilter< TScalar > TVesselness;
-
-      TMinMax::Pointer minMax = TMinMax::New( );
-      minMax->SetImage( input_image );
-      minMax->Compute( );
-
-      TInverter::Pointer inverter = TInverter::New( );
-      inverter->SetInput( input_image );
-      inverter->SetMaximum( minMax->GetMaximum( ) );
-      double t = MeasureTime( inverter );
-      std::cout << "Inversion executed in " << t << " s" << std::endl;
-
-      THessian::Pointer hessian = THessian::New( );
-      hessian->SetInput( inverter->GetOutput( ) );
-      hessian->SetSigma(
-        std::atof( args[ "sigma" ].c_str( ) )
-        );
-      t = MeasureTime( hessian );
-      std::cout << "Hessian executed in " << t << " s" << std::endl;
-
-      TVesselness::Pointer vesselness = TVesselness::New( );
-      vesselness->SetInput( hessian->GetOutput( ) );
-      vesselness->SetAlpha1(
-        std::atof( args[ "alpha1" ].c_str( ) )
-        );
-      vesselness->SetAlpha2(
-        std::atof( args[ "alpha2" ].c_str( ) )
-        );
-      t = MeasureTime( vesselness );
-      std::cout << "Vesselness executed in " << t << " s" << std::endl;
-
-      WriteImage( vesselness->GetOutput( ), args[ "out" ] );
-    }
-    else
-      return( 1 );
+    // Read image
+    TImage::Pointer input_image;
+    CTBronchi::ReadImage( input_image, in.getValue( ) );
+
+    // Min-max
+    typedef itk::MinimumMaximumImageCalculator< TImage > _TMinMax;
+    _TMinMax::Pointer minMax = _TMinMax::New( );
+    minMax->SetImage( input_image );
+    minMax->Compute( );
+
+    // Invert intensities
+    typedef itk::InvertIntensityImageFilter< TImage > _TInverter;
+    _TInverter::Pointer inverter = _TInverter::New( );
+    inverter->SetInput( input_image );
+    inverter->SetMaximum( minMax->GetMaximum( ) );
+    double t = CTBronchi::MeasureTime( inverter );
+    std::cout << "Inversion executed in " << t << " s" << std::endl;
+
+    // Compute hessian image
+    typedef itk::HessianRecursiveGaussianImageFilter< TImage > _THessian;
+    _THessian::Pointer hessian = _THessian::New( );
+    hessian->SetInput( inverter->GetOutput( ) );
+    hessian->SetSigma( s.getValue( ) );
+    t = CTBronchi::MeasureTime( hessian );
+    std::cout << "Hessian executed in " << t << " s" << std::endl;
+
+    // Vesselness
+    typedef
+      itk::Hessian3DToVesselnessMeasureImageFilter< TScalar >
+      _TVesselness;
+    _TVesselness::Pointer vesselness = _TVesselness::New( );
+    vesselness->SetInput( hessian->GetOutput( ) );
+    vesselness->SetAlpha1( a1.getValue( ) );
+    vesselness->SetAlpha2( a2.getValue( ) );
+    t = CTBronchi::MeasureTime( vesselness );
+    std::cout << "Vesselness executed in " << t << " s" << std::endl;
+
+    // Write result
+    CTBronchi::WriteImage( vesselness->GetOutput( ), out.getValue( ) );
   }
   catch( std::exception& err )
   {