X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=appli%2FCTBronchi%2FVesselness.cxx;h=3b0661824de9cb3841800bbd7f067c98310b8367;hb=9817c556a0b8b5e3b332d45f07faa84d91afb2d0;hp=ca7c99c358c1df14d63fd1db2794e34a1f452933;hpb=ae0e1b8916a0fb2188080b9134c1c2781c6c200f;p=FrontAlgorithms.git diff --git a/appli/CTBronchi/Vesselness.cxx b/appli/CTBronchi/Vesselness.cxx index ca7c99c..3b06618 100644 --- a/appli/CTBronchi/Vesselness.cxx +++ b/appli/CTBronchi/Vesselness.cxx @@ -3,153 +3,97 @@ // @email florez-l@javeriana.edu.co // ========================================================================= -#include -#include +#include +#include #include -#include -#include #include #include #include #include +#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 ) {