1 // =========================================================================
2 // @author Leonardo Florez Valencia (florez-l@javeriana.edu.co)
3 // =========================================================================
4 #ifndef __CTBronchi__Process__h__
5 #define __CTBronchi__Process__h__
13 // Some types and values
14 static const unsigned int Dim = 3;
16 typedef unsigned char TLabel;
17 typedef float TScalar;
19 double MeasureTime( itk::ProcessObject* f );
21 template< class _TImagePtr >
22 void ReadImage( _TImagePtr& image, const std::string& fname );
24 template< class _TImage >
25 void WriteImage( const _TImage* image, const std::string& fname );
32 #include <tclap/CmdLine.h>
34 #include <itkMinimumMaximumImageCalculator.h>
35 #include <itkInvertIntensityImageFilter.h>
36 #include <itkHessianRecursiveGaussianImageFilter.h>
37 #include <itkHessian3DToVesselnessMeasureImageFilter.h>
38 #include "Functions.h"
40 // -------------------------------------------------------------------------
41 const unsigned int Dim = 3;
43 typedef itk::NumericTraits< TPixel >::RealType TScalar;
44 typedef itk::Image< TPixel, Dim > TImage;
45 typedef itk::Image< TScalar, Dim > TScalarImage;
47 // -------------------------------------------------------------------------
48 int main( int argc, char* argv[] )
50 typedef TCLAP::ValueArg< std::string > _TStringArg;
51 typedef TCLAP::ValueArg< TScalar > _TRealArg;
54 _TStringArg in( "i", "input", "Input image", true, "", "file" );
55 _TStringArg out( "o", "output", "Output image", true, "", "file" );
56 _TRealArg s( "s", "sigma", "Sigma", false, 0.5, "value (0.5)" );
57 _TRealArg a1( "a", "alpha1", "Alpha1", false, 0.5, "value (0.5)" );
58 _TRealArg a2( "b", "alpha2", "Alpha2", false, 2, "value (2)" );
61 TCLAP::CmdLine cmd( "Vesselness computation", ' ', "1.0.0" );
67 cmd.parse( argc, argv );
69 catch( TCLAP::ArgException& err )
72 << "===============================" << std::endl
73 << "Error caught: " << std::endl
74 << err.error( ) << " " << err.argId( )
75 << "===============================" << std::endl
84 TImage::Pointer input_image;
85 CTBronchi::ReadImage( input_image, in.getValue( ) );
88 typedef itk::MinimumMaximumImageCalculator< TImage > _TMinMax;
89 _TMinMax::Pointer minMax = _TMinMax::New( );
90 minMax->SetImage( input_image );
94 typedef itk::InvertIntensityImageFilter< TImage > _TInverter;
95 _TInverter::Pointer inverter = _TInverter::New( );
96 inverter->SetInput( input_image );
97 inverter->SetMaximum( minMax->GetMaximum( ) );
98 double t = CTBronchi::MeasureTime( inverter );
99 std::cout << "Inversion executed in " << t << " s" << std::endl;
101 // Compute hessian image
102 typedef itk::HessianRecursiveGaussianImageFilter< TImage > _THessian;
103 _THessian::Pointer hessian = _THessian::New( );
104 hessian->SetInput( inverter->GetOutput( ) );
105 hessian->SetSigma( s.getValue( ) );
106 t = CTBronchi::MeasureTime( hessian );
107 std::cout << "Hessian executed in " << t << " s" << std::endl;
111 itk::Hessian3DToVesselnessMeasureImageFilter< TScalar >
113 _TVesselness::Pointer vesselness = _TVesselness::New( );
114 vesselness->SetInput( hessian->GetOutput( ) );
115 vesselness->SetAlpha1( a1.getValue( ) );
116 vesselness->SetAlpha2( a2.getValue( ) );
117 t = CTBronchi::MeasureTime( vesselness );
118 std::cout << "Vesselness executed in " << t << " s" << std::endl;
121 CTBronchi::WriteImage( vesselness->GetOutput( ), out.getValue( ) );
123 catch( std::exception& err )
126 << "===============================" << std::endl
127 << "Error caught: " << std::endl
129 << "===============================" << std::endl
137 #endif // __CTBronchi__Process__h__