1 // =========================================================================
2 // @author Leonardo Florez Valencia
3 // @email florez-l@javeriana.edu.co
4 // =========================================================================
7 #include <tclap/CmdLine.h>
9 #include <itkMinimumMaximumImageCalculator.h>
10 #include <itkInvertIntensityImageFilter.h>
11 #include <itkHessianRecursiveGaussianImageFilter.h>
12 #include <itkHessian3DToVesselnessMeasureImageFilter.h>
13 #include "Functions.h"
15 // -------------------------------------------------------------------------
16 const unsigned int Dim = 3;
18 typedef itk::NumericTraits< TPixel >::RealType TScalar;
19 typedef itk::Image< TPixel, Dim > TImage;
20 typedef itk::Image< TScalar, Dim > TScalarImage;
22 // -------------------------------------------------------------------------
23 int main( int argc, char* argv[] )
25 typedef TCLAP::ValueArg< std::string > _TStringArg;
26 typedef TCLAP::ValueArg< TScalar > _TRealArg;
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)" );
36 TCLAP::CmdLine cmd( "Vesselness computation", ' ', "1.0.0" );
42 cmd.parse( argc, argv );
44 catch( TCLAP::ArgException& err )
47 << "===============================" << std::endl
48 << "Error caught: " << std::endl
49 << err.error( ) << " " << err.argId( )
50 << "===============================" << std::endl
59 TImage::Pointer input_image;
60 CTBronchi::ReadImage( input_image, in.getValue( ) );
63 typedef itk::MinimumMaximumImageCalculator< TImage > _TMinMax;
64 _TMinMax::Pointer minMax = _TMinMax::New( );
65 minMax->SetImage( input_image );
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;
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;
86 itk::Hessian3DToVesselnessMeasureImageFilter< TScalar >
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;
96 CTBronchi::WriteImage( vesselness->GetOutput( ), out.getValue( ) );
98 catch( std::exception& err )
101 << "===============================" << std::endl
102 << "Error caught: " << std::endl
104 << "===============================" << std::endl