1 // =========================================================================
2 // @author Leonardo Florez Valencia
3 // @email florez-l@javeriana.edu.co
4 // =========================================================================
9 #include <itkImageFileReader.h>
10 #include <itkImageFileWriter.h>
11 #include <itkMinimumMaximumImageCalculator.h>
12 #include <itkInvertIntensityImageFilter.h>
13 #include <itkHessianRecursiveGaussianImageFilter.h>
14 #include <itkHessian3DToVesselnessMeasureImageFilter.h>
16 // -------------------------------------------------------------------------
17 const unsigned int Dim = 3;
19 typedef itk::NumericTraits< TPixel >::RealType TScalar;
20 typedef itk::Image< TPixel, Dim > TImage;
21 typedef itk::Image< TScalar, Dim > TScalarImage;
23 // -------------------------------------------------------------------------
24 double MeasureTime( itk::ProcessObject* f )
26 std::chrono::time_point< std::chrono::high_resolution_clock > s, e;
27 std::chrono::duration< double > t;
28 s = std::chrono::high_resolution_clock::now( );
30 e = std::chrono::high_resolution_clock::now( );
35 // -------------------------------------------------------------------------
36 template< class _TImagePtr >
37 void ReadImage( _TImagePtr& image, const std::string& fname )
39 typedef typename _TImagePtr::ObjectType _TImage;
40 typedef itk::ImageFileReader< _TImage > _TReader;
41 typename _TReader::Pointer reader = _TReader::New( );
42 reader->SetFileName( fname );
43 double t = MeasureTime( reader );
44 std::cout << "Read " << fname << " in " << t << " s" << std::endl;
45 image = reader->GetOutput( );
46 image->DisconnectPipeline( );
49 // -------------------------------------------------------------------------
50 template< class _TImage >
51 void WriteImage( const _TImage* image, const std::string& fname )
53 typedef itk::ImageFileWriter< _TImage > _TWriter;
54 typename _TWriter::Pointer writer = _TWriter::New( );
55 writer->SetFileName( fname );
56 writer->SetInput( image );
57 double t = MeasureTime( writer );
58 std::cout << "Wrote " << fname << " in " << t << " s" << std::endl;
61 // -------------------------------------------------------------------------
63 std::map< std::string, std::string >& args, int argc, char* argv[]
66 std::set< std::string > mandatory;
67 mandatory.insert( "in" );
68 mandatory.insert( "out" );
70 args[ "sigma" ] = "0.5";
71 args[ "alpha1" ] = "0.5";
72 args[ "alpha2" ] = "2";
77 std::string cmd = argv[ i ] + 1;
78 args[ cmd ] = argv[ i + 1 ];
84 for( std::string t: mandatory )
85 complete &= ( args.find( t ) != args.end( ) );
90 << "Usage: " << argv[ 0 ] << std::endl
91 << "\t-in filename" << std::endl
92 << "\t-out filename" << std::endl
93 << "\t[-sigma val]" << std::endl
94 << "\t[-alpha1 val]" << std::endl
95 << "\t[-alpha2 val]" << std::endl;
102 // -------------------------------------------------------------------------
103 int main( int argc, char* argv[] )
105 std::map< std::string, std::string > args;
108 if( ParseArgs( args, argc, argv ) )
111 TImage::Pointer input_image;
112 ReadImage( input_image, args[ "in" ] );
115 typedef itk::MinimumMaximumImageCalculator< TImage > TMinMax;
116 typedef itk::InvertIntensityImageFilter< TImage > TInverter;
117 typedef itk::HessianRecursiveGaussianImageFilter< TImage > THessian;
118 typedef itk::Hessian3DToVesselnessMeasureImageFilter< TScalar > TVesselness;
120 TMinMax::Pointer minMax = TMinMax::New( );
121 minMax->SetImage( input_image );
124 TInverter::Pointer inverter = TInverter::New( );
125 inverter->SetInput( input_image );
126 inverter->SetMaximum( minMax->GetMaximum( ) );
127 double t = MeasureTime( inverter );
128 std::cout << "Inversion executed in " << t << " s" << std::endl;
130 THessian::Pointer hessian = THessian::New( );
131 hessian->SetInput( inverter->GetOutput( ) );
133 std::atof( args[ "sigma" ].c_str( ) )
135 t = MeasureTime( hessian );
136 std::cout << "Hessian executed in " << t << " s" << std::endl;
138 TVesselness::Pointer vesselness = TVesselness::New( );
139 vesselness->SetInput( hessian->GetOutput( ) );
140 vesselness->SetAlpha1(
141 std::atof( args[ "alpha1" ].c_str( ) )
143 vesselness->SetAlpha2(
144 std::atof( args[ "alpha2" ].c_str( ) )
146 t = MeasureTime( vesselness );
147 std::cout << "Vesselness executed in " << t << " s" << std::endl;
149 WriteImage( vesselness->GetOutput( ), args[ "out" ] );
154 catch( std::exception& err )
157 << "===============================" << std::endl
158 << "Error caught: " << std::endl
160 << "===============================" << std::endl