// @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 )
{