From: Leonardo Flórez-Valencia Date: Wed, 11 Oct 2017 20:59:00 +0000 (-0500) Subject: ... X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=926cb50be2ef54fce0d9493119dc89626e9b7526;p=FrontAlgorithms.git ... --- diff --git a/appli/CTBronchi/CMakeLists.txt b/appli/CTBronchi/CMakeLists.txt index d81b5e8..a5c2b61 100644 --- a/appli/CTBronchi/CMakeLists.txt +++ b/appli/CTBronchi/CMakeLists.txt @@ -17,7 +17,7 @@ if(fpa_BUILD_CTBronchi) SOURCE ${_e}.cxx INSTALL RECURRENT - LINKS fpa + LINKS fpa cpPlugins::tclap ) endforeach(_e) endif(fpa_BUILD_CTBronchi) diff --git a/appli/CTBronchi/Functions.h b/appli/CTBronchi/Functions.h index 7f20083..2e9732c 100644 --- a/appli/CTBronchi/Functions.h +++ b/appli/CTBronchi/Functions.h @@ -6,8 +6,6 @@ #define __CTBronchi__Functions__h__ #include -#include -#include #include #include @@ -34,97 +32,21 @@ namespace CTBronchi typename _TReader::Pointer reader = _TReader::New( ); reader->SetFileName( fname ); double t = MeasureTime( reader ); - std::cout << "Read " << fname << " in " << t << " s" << std::endl; + std::cout << "Read \"" << fname << "\" in " << t << " s" << std::endl; image = reader->GetOutput( ); image->DisconnectPipeline( ); } // ----------------------------------------------------------------------- - template< class _TImagePtr > - void WriteImage( const _TImagePtr& image, const std::string& fname ) + template< class _TImage > + void WriteImage( const _TImage* image, const std::string& fname ) { - typedef typename _TImagePtr::ObjectType _TImage; 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[] - ) - { - std::set< std::string > mandatory; - mandatory.insert( "in" ); - mandatory.insert( "out" ); - mandatory.insert( "seed" ); - - args[ "mori_minimum_threshold" ] = "-850"; - args[ "mori_signal_kernel_size" ] = "20"; - args[ "mori_signal_threshold" ] = "100"; - args[ "mori_signal_influence" ] = "0.5"; - args[ "mori_lower_threshold" ] = "-1024"; - args[ "mori_upper_threshold" ] = "0"; - args[ "mori_delta_threshold" ] = "1"; - args[ "labelling_upper_threshold" ] = "-600"; - args[ "random_walker_alpha" ] = "1"; - args[ "random_walker_beta" ] = "100"; - - int i = 1; - while( i < argc ) - { - std::string cmd = argv[ i ] + 1; - if( cmd == "seed" ) - { - std::stringstream seed; - seed - << argv[ i + 1 ] << ";" - << argv[ i + 2 ] << ";" - << argv[ i + 3 ]; - args[ cmd ] = seed.str( ); - i += 4; - } - else - { - args[ cmd ] = argv[ i + 1 ]; - i += 2; - - } // fi - - } // elihw - - bool complete = true; - for( std::string t: mandatory ) - complete &= ( args.find( t ) != args.end( ) ); - - if( !complete ) - { - std::cerr - << "Usage: " << argv[ 0 ] << std::endl - << "\t-in filename" << std::endl - << "\t-out filename" << std::endl - << "\t-seed x y z" << std::endl - << "\t[-out_mori filename]" << std::endl - << "\t[-out_signal filename]" << std::endl - << "\t[-out_labels filename]" << std::endl - << "\t[-out_diff filename]" << std::endl - << "\t[-mori_minimum_threshold value]" << std::endl - << "\t[-mori_signal_kernel_size value]" << std::endl - << "\t[-mori_signal_threshold value]" << std::endl - << "\t[-mori_signal_influence value]" << std::endl - << "\t[-mori_lower_threshold value]" << std::endl - << "\t[-mori_upper_threshold value]" << std::endl - << "\t[-mori_delta_threshold value]" << std::endl - << "\t[-labelling_upper_threshold value]" << std::endl - << "\t[-random_walker_alpha value]" << std::endl - << "\t[-random_walker_beta value]" << std::endl; - return( false ); - - } // fi - return( true ); + std::cout << "Wrote \"" << fname << "\" in " << t << " s" << std::endl; } } // ecapseman diff --git a/appli/CTBronchi/MoriLabelling.cxx b/appli/CTBronchi/MoriLabelling.cxx index 848bf6d..96f4dc4 100644 --- a/appli/CTBronchi/MoriLabelling.cxx +++ b/appli/CTBronchi/MoriLabelling.cxx @@ -3,20 +3,12 @@ // @email florez-l@javeriana.edu.co // ========================================================================= -#include -#include - +#include +#include +#include #include -#include -#include -#include -#include - -#include -#include -#include -#include -#include +#include "MoriLabelling.h" +#include "Functions.h" // ------------------------------------------------------------------------- const unsigned int Dim = 3; @@ -27,323 +19,80 @@ typedef itk::Image< TPixel, Dim > TImage; typedef itk::Image< TLabel, Dim > TLabels; typedef itk::Image< TScalar, Dim > TScalarImage; -/** - */ -class MoriLabellingTraits - : public fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel > -{ -public: - typedef fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel > Superclass; - typedef typename Superclass::TInternalTraits TInternalTraits; - typedef typename Superclass::TMarksImage TMarksImage; - typedef typename Superclass::TFilterInterface TFilterInterface; - - typedef fpa::Filters::BaseMarksInterface< TInternalTraits > TMarksInterface; - typedef fpa::Filters::Image::SeedsFromLabelsInterface< TInternalTraits > TSeedsInterface; -}; - -/** - */ -class MoriLabelling - : public fpa::Filters::Image::RegionGrow< TImage, TLabels, TLabel, MoriLabellingTraits > -{ -public: - typedef MoriLabellingTraits TTraits; - fpaTraitsMacro( TTraits ); - - typedef fpa::Filters::Image::RegionGrow< TImage, TLabels, TMark, TTraits > Superclass; - typedef MoriLabelling Self; - typedef itk::SmartPointer< Self > Pointer; - typedef itk::SmartPointer< const Self > ConstPointer; - typedef fpa::Functors::RegionGrow::BinaryThreshold< TInputValue > TFunctor; - -public: - itkNewMacro( Self ); - itkTypeMacro( MoriLabelling, fpa::Filters::Image::RegionGrow ); - - itkGetConstMacro( LastThreshold, TInputValue ); - itkSetMacro( LastThreshold, TInputValue ); - - itkGetConstMacro( MinVertex, TVertex ); - itkGetConstMacro( MaxVertex, TVertex ); - - ivqITKInputMacro( InputLabels, TLabels ); - ivqITKInputMacro( InputVesselness, TScalarImage ); - -public: - TInputValue GetUpperThreshold( ) const - { - return( this->m_Functor->GetUpperThreshold( ) ); - } - void SetUpperThreshold( TInputValue t ) - { - this->m_Functor->SetUpperThreshold( t ); - } - -protected: - MoriLabelling( ) - : Superclass( ), - m_LastThreshold( TInputValue( 0 ) ), - m_VesselnessThr( TScalar( 0.05 ) ) - { - ivqITKInputConfigureMacro( InputLabels, TLabels ); - ivqITKInputConfigureMacro( InputVesselness, TScalarImage ); - this->m_Functor = TFunctor::New( ); - this->SetPredicate( this->m_Functor ); - } - - virtual ~MoriLabelling( ) - { - } - - virtual const itk::DataObject* _GetReferenceInput( ) const override - { - return( this->GetInputLabels( ) ); - } - virtual void _BeforeGenerateData( ) override - { - this->Superclass::_BeforeGenerateData( ); - this->m_FirstVertex = true; - - typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax; - _TMinMax::Pointer minMax = _TMinMax::New( ); - minMax->SetImage( this->GetInputVesselness( ) ); - minMax->Compute( ); - this->m_MaxVesselness = ( 1.0 - this->m_VesselnessThr ) * minMax->GetMaximum( ); - } - - virtual void _PostComputeOutputValue( TNode& n ) override - { - this->Superclass::_PostComputeOutputValue( n ); - if( n.Value == this->GetInsideValue( ) ) - { - const TImage* input = this->GetInput( ); - const TLabels* labels = this->GetInputLabels( ); - const TScalarImage* vesselness = this->GetInputVesselness( ); - double x = input->GetPixel( n.Vertex ); - /* TODO - double a = std::fabs( x - double( this->m_Functor->GetUpperThreshold( ) ) ); - double b = std::fabs( x - double( this->m_LastThreshold ) ); - */ - if( labels->GetPixel( n.Vertex ) == 0 ) - { - if( vesselness->GetPixel( n.Vertex ) > this->m_MaxVesselness ) - n.Value = this->GetInsideValue( ); - else - n.Value = 0; - - } // fi - - if( !( this->m_FirstVertex ) ) - { - for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) - { - if( n.Vertex[ d ] < this->m_MinVertex[ d ] ) - this->m_MinVertex[ d ] = n.Vertex[ d ]; - if( this->m_MaxVertex[ d ] < n.Vertex[ d ] ) - this->m_MaxVertex[ d ] = n.Vertex[ d ]; - - } // rof - } - else - { - this->m_MinVertex = n.Vertex; - this->m_MaxVertex = n.Vertex; - this->m_FirstVertex = false; - - } // fi - - } // fi - } - -private: - // Purposely not implemented. - MoriLabelling( const Self& other ); - Self& operator=( const Self& other ); - -protected: - TFunctor::Pointer m_Functor; - TInputValue m_LastThreshold; - bool m_FirstVertex; - TVertex m_MinVertex; - TVertex m_MaxVertex; - - TScalar m_MaxVesselness; - TScalar m_VesselnessThr; -}; - -// ------------------------------------------------------------------------- -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; -} - -// ------------------------------------------------------------------------- -template< class _TImagePtr, class _TRegion > -void ROI( - _TImagePtr& output, - const typename _TImagePtr::ObjectType* input, - const _TRegion& roi - ) -{ - typedef typename _TImagePtr::ObjectType _TImage; - typedef itk::RegionOfInterestImageFilter< _TImage, _TImage > _TFilter; - - typename _TFilter::Pointer filter = _TFilter::New( ); - filter->SetInput( input ); - filter->SetRegionOfInterest( roi ); - double t = MeasureTime( filter ); - std::cout << "ROI computed in " << t << " s" << std::endl; - output = filter->GetOutput( ); - output->DisconnectPipeline( ); -} - // ------------------------------------------------------------------------- -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" ); - mandatory.insert( "labels" ); - mandatory.insert( "vesselness" ); - - args[ "upper_threshold" ] = "-650"; - - int i = 1; - while( i < argc ) + typedef TCLAP::ValueArg< std::string > _TStringArg; + typedef TCLAP::ValueArg< double > _TRealArg; + typedef TCLAP::ValueArg< TPixel > _TPixelArg; + + // Parse input line + _TStringArg in( "i", "input", "Input image", true, "", "file" ); + _TStringArg labels( "l", "labels", "Input labels", true, "", "file" ); + _TStringArg vesselness( "v", "vesselness", "Input vesselness", true, "", "file" ); + _TStringArg out( "o", "output", "Output image", true, "", "file" ); + _TRealArg vThr( "a", "vesselness_thr", "Vesselness threshold", false, 0.05, "value (0.05)" ); + _TPixelArg uThr( "u", "upper_thr", "Upper threshold", false, -650, "value (-650)" ); + 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( "Labelling", ' ', "1.0.0" ); + cmd.add( uThr ); + cmd.add( vThr ); + cmd.add( out ); + cmd.add( vesselness ); + cmd.add( labels ); + 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-labels filename" << std::endl - << "\t-vesselness filename" << std::endl - << "\t[-roi] filename" << std::endl - << "\t[-upper_threshold value]" << std::endl; - return( false ); - - } // fi - return( true ); -} + << "===============================" << std::endl + << "Error caught: " << std::endl + << err.error( ) << " " << err.argId( ) << std::endl + << "===============================" << std::endl + << std::endl; + return( 1 ); + + } // 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" ] ); - - // Read labels image - TLabels::Pointer input_labels; - ReadImage( input_labels, args[ "labels" ] ); - - // Read vesselness image - TScalarImage::Pointer input_vesselness; - ReadImage( input_vesselness, args[ "vesselness" ] ); - - // Mori labelling - MoriLabelling::Pointer labelling = MoriLabelling::New( ); - labelling->SetInput( input_image ); - labelling->SetInputLabels( input_labels ); - labelling->SetInputVesselness( input_vesselness ); - labelling->SetOutsideValue( 2 ); - labelling->SetFillValue( 2 ); - labelling->SetInsideValue( 1 ); - labelling->SetUpperThreshold( - TPixel( std::atof( args[ "upper_threshold" ].c_str( ) ) ) - ); - /* TODO - labelling->SetLastThreshold( last_thr ); - */ - double t = MeasureTime( labelling ); - std::cout << "Labelling executed in " << t << " s" << std::endl; - - std::map< std::string, std::string >::const_iterator i = - args.find( "roi" ); - if( i != args.end( ) ) - { - TImage::IndexType minV = labelling->GetMinVertex( ); - TImage::IndexType maxV = labelling->GetMaxVertex( ); - TImage::SizeType roiS; - for( unsigned d = 0; d < TImage::ImageDimension; ++d ) - roiS[ d ] = maxV[ d ] - minV[ d ] + 1; - TImage::RegionType roi; - roi.SetIndex( minV ); - roi.SetSize( roiS ); - - // ROI input image - TImage::Pointer input_image_roi; - ROI( input_image_roi, input_image.GetPointer( ), roi ); - WriteImage( input_image_roi.GetPointer( ), i->second ); - - // ROI output image - TLabels::Pointer output_roi; - ROI( output_roi, labelling->GetOutput( ), roi ); - WriteImage( output_roi.GetPointer( ), args[ "out" ] ); - } - else - WriteImage( labelling->GetOutput( ), args[ "out" ] ); - } - else - return( 1 ); + // Read input image + TImage::Pointer input_image; + CTBronchi::ReadImage( input_image, in.getValue( ) ); + + // Read input labels + TLabels::Pointer input_labels; + CTBronchi::ReadImage( input_labels, labels.getValue( ) ); + + // Read input vesselness + TScalarImage::Pointer input_vesselness; + CTBronchi::ReadImage( input_vesselness, vesselness.getValue( ) ); + + // Create labels + typedef CTBronchi::MoriLabelling< TImage, TLabels, TScalarImage > _TLabelling; + _TLabelling::Pointer labelling = _TLabelling::New( ); + labelling->SetInput( input_image ); + labelling->SetInputLabels( input_labels ); + labelling->SetInputVesselness( input_vesselness ); + labelling->SetVesselnessThreshold( vThr.getValue( ) ); + labelling->SetUpperThreshold( uThr.getValue( ) ); + labelling->SetInsideValue( 1 ); + labelling->SetOutsideValue( 2 ); + labelling->SetFillValue( 2 ); + double t = CTBronchi::MeasureTime( labelling ); + std::cout << "Labelling executed in " << t << " s" << std::endl; + + // Write result + CTBronchi::WriteImage( labelling->GetOutput( ), out.getValue( ) ); } catch( std::exception& err ) { std::cerr << "===============================" << std::endl << "Error caught: " << std::endl - << err.what( ) + << err.what( ) << std::endl << "===============================" << std::endl << std::endl; return( 1 ); @@ -352,4 +101,344 @@ int main( int argc, char* argv[] ) return( 0 ); } +/* TODO + #include + #include + + #include + #include + #include + #include + #include + + #include + #include + #include + #include + #include + + // ------------------------------------------------------------------------- + const unsigned int Dim = 3; + typedef short TPixel; + typedef unsigned char TLabel; + typedef itk::NumericTraits< TPixel >::RealType TScalar; + typedef itk::Image< TPixel, Dim > TImage; + typedef itk::Image< TLabel, Dim > TLabels; + typedef itk::Image< TScalar, Dim > TScalarImage; + + class MoriLabellingTraits + : public fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel > + { + public: + typedef fpa::Filters::Image::DefaultTraits< TImage, TLabels, TLabel > Superclass; + typedef typename Superclass::TInternalTraits TInternalTraits; + typedef typename Superclass::TMarksImage TMarksImage; + typedef typename Superclass::TFilterInterface TFilterInterface; + + typedef fpa::Filters::BaseMarksInterface< TInternalTraits > TMarksInterface; + typedef fpa::Filters::Image::SeedsFromLabelsInterface< TInternalTraits > TSeedsInterface; + }; + + class MoriLabelling + : public fpa::Filters::Image::RegionGrow< TImage, TLabels, TLabel, MoriLabellingTraits > + { + public: + typedef MoriLabellingTraits TTraits; + fpaTraitsMacro( TTraits ); + + typedef fpa::Filters::Image::RegionGrow< TImage, TLabels, TMark, TTraits > Superclass; + typedef MoriLabelling Self; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + typedef fpa::Functors::RegionGrow::BinaryThreshold< TInputValue > TFunctor; + + public: + itkNewMacro( Self ); + itkTypeMacro( MoriLabelling, fpa::Filters::Image::RegionGrow ); + + itkGetConstMacro( LastThreshold, TInputValue ); + itkSetMacro( LastThreshold, TInputValue ); + + itkGetConstMacro( MinVertex, TVertex ); + itkGetConstMacro( MaxVertex, TVertex ); + + ivqITKInputMacro( InputLabels, TLabels ); + ivqITKInputMacro( InputVesselness, TScalarImage ); + + public: + TInputValue GetUpperThreshold( ) const + { + return( this->m_Functor->GetUpperThreshold( ) ); + } + void SetUpperThreshold( TInputValue t ) + { + this->m_Functor->SetUpperThreshold( t ); + } + + protected: + MoriLabelling( ) + : Superclass( ), + m_LastThreshold( TInputValue( 0 ) ), + m_VesselnessThr( TScalar( 0.05 ) ) + { + ivqITKInputConfigureMacro( InputLabels, TLabels ); + ivqITKInputConfigureMacro( InputVesselness, TScalarImage ); + this->m_Functor = TFunctor::New( ); + this->SetPredicate( this->m_Functor ); + } + + virtual ~MoriLabelling( ) + { + } + + virtual const itk::DataObject* _GetReferenceInput( ) const override + { + return( this->GetInputLabels( ) ); + } + virtual void _BeforeGenerateData( ) override + { + this->Superclass::_BeforeGenerateData( ); + this->m_FirstVertex = true; + + typedef itk::MinimumMaximumImageCalculator< TScalarImage > _TMinMax; + _TMinMax::Pointer minMax = _TMinMax::New( ); + minMax->SetImage( this->GetInputVesselness( ) ); + minMax->Compute( ); + this->m_MaxVesselness = ( 1.0 - this->m_VesselnessThr ) * minMax->GetMaximum( ); + } + + virtual void _PostComputeOutputValue( TNode& n ) override + { + this->Superclass::_PostComputeOutputValue( n ); + if( n.Value == this->GetInsideValue( ) ) + { + const TImage* input = this->GetInput( ); + const TLabels* labels = this->GetInputLabels( ); + const TScalarImage* vesselness = this->GetInputVesselness( ); + double x = input->GetPixel( n.Vertex ); + if( labels->GetPixel( n.Vertex ) == 0 ) + { + if( vesselness->GetPixel( n.Vertex ) > this->m_MaxVesselness ) + n.Value = this->GetInsideValue( ); + else + n.Value = 0; + + } // fi + + if( !( this->m_FirstVertex ) ) + { + for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) + { + if( n.Vertex[ d ] < this->m_MinVertex[ d ] ) + this->m_MinVertex[ d ] = n.Vertex[ d ]; + if( this->m_MaxVertex[ d ] < n.Vertex[ d ] ) + this->m_MaxVertex[ d ] = n.Vertex[ d ]; + + } // rof + } + else + { + this->m_MinVertex = n.Vertex; + this->m_MaxVertex = n.Vertex; + this->m_FirstVertex = false; + + } // fi + + } // fi + } + + private: + // Purposely not implemented. + MoriLabelling( const Self& other ); + Self& operator=( const Self& other ); + + protected: + TFunctor::Pointer m_Functor; + TInputValue m_LastThreshold; + bool m_FirstVertex; + TVertex m_MinVertex; + TVertex m_MaxVertex; + + TScalar m_MaxVesselness; + TScalar m_VesselnessThr; + }; + + // ------------------------------------------------------------------------- + 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; + } + + // ------------------------------------------------------------------------- + template< class _TImagePtr, class _TRegion > + void ROI( + _TImagePtr& output, + const typename _TImagePtr::ObjectType* input, + const _TRegion& roi + ) + { + typedef typename _TImagePtr::ObjectType _TImage; + typedef itk::RegionOfInterestImageFilter< _TImage, _TImage > _TFilter; + + typename _TFilter::Pointer filter = _TFilter::New( ); + filter->SetInput( input ); + filter->SetRegionOfInterest( roi ); + double t = MeasureTime( filter ); + std::cout << "ROI computed in " << t << " s" << std::endl; + output = filter->GetOutput( ); + output->DisconnectPipeline( ); + } + + // ------------------------------------------------------------------------- + bool ParseArgs( + std::map< std::string, std::string >& args, int argc, char* argv[] + ) + { + std::set< std::string > mandatory; + mandatory.insert( "in" ); + mandatory.insert( "out" ); + mandatory.insert( "labels" ); + mandatory.insert( "vesselness" ); + + args[ "upper_threshold" ] = "-650"; + + int i = 1; + while( i < argc ) + { + 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 ) + { + std::cerr + << "Usage: " << argv[ 0 ] << std::endl + << "\t-in filename" << std::endl + << "\t-out filename" << std::endl + << "\t-labels filename" << std::endl + << "\t-vesselness filename" << std::endl + << "\t[-roi] filename" << std::endl + << "\t[-upper_threshold value]" << std::endl; + return( false ); + + } // fi + return( true ); + } + + // ------------------------------------------------------------------------- + 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" ] ); + + // Read labels image + TLabels::Pointer input_labels; + ReadImage( input_labels, args[ "labels" ] ); + + // Read vesselness image + TScalarImage::Pointer input_vesselness; + ReadImage( input_vesselness, args[ "vesselness" ] ); + + // Mori labelling + MoriLabelling::Pointer labelling = MoriLabelling::New( ); + labelling->SetInput( input_image ); + labelling->SetInputLabels( input_labels ); + labelling->SetInputVesselness( input_vesselness ); + labelling->SetOutsideValue( 2 ); + labelling->SetFillValue( 2 ); + labelling->SetInsideValue( 1 ); + labelling->SetUpperThreshold( + TPixel( std::atof( args[ "upper_threshold" ].c_str( ) ) ) + ); + double t = MeasureTime( labelling ); + std::cout << "Labelling executed in " << t << " s" << std::endl; + + std::map< std::string, std::string >::const_iterator i = + args.find( "roi" ); + if( i != args.end( ) ) + { + TImage::IndexType minV = labelling->GetMinVertex( ); + TImage::IndexType maxV = labelling->GetMaxVertex( ); + TImage::SizeType roiS; + for( unsigned d = 0; d < TImage::ImageDimension; ++d ) + roiS[ d ] = maxV[ d ] - minV[ d ] + 1; + TImage::RegionType roi; + roi.SetIndex( minV ); + roi.SetSize( roiS ); + + // ROI input image + TImage::Pointer input_image_roi; + ROI( input_image_roi, input_image.GetPointer( ), roi ); + WriteImage( input_image_roi.GetPointer( ), i->second ); + + // ROI output image + TLabels::Pointer output_roi; + ROI( output_roi, labelling->GetOutput( ), roi ); + WriteImage( output_roi.GetPointer( ), args[ "out" ] ); + } + else + WriteImage( labelling->GetOutput( ), args[ "out" ] ); + } + else + return( 1 ); + } + catch( std::exception& err ) + { + std::cerr + << "===============================" << std::endl + << "Error caught: " << std::endl + << err.what( ) + << "===============================" << std::endl + << std::endl; + return( 1 ); + + } // yrt + return( 0 ); + } +*/ + // eof - $RCSfile$ diff --git a/appli/CTBronchi/MoriLabelling.h b/appli/CTBronchi/MoriLabelling.h new file mode 100644 index 0000000..89ce001 --- /dev/null +++ b/appli/CTBronchi/MoriLabelling.h @@ -0,0 +1,93 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= +#ifndef __CTBronchi__MoriLabelling__h__ +#define __CTBronchi__MoriLabelling__h__ + +#include +#include +#include +#include +#include + +namespace CTBronchi +{ + /** + */ + template< class _TImage, class _TLabels > + class MoriLabellingTraits + : public fpa::Filters::Image::DefaultTraits< _TImage, _TLabels, typename _TLabels::PixelType > + { + public: + typedef fpa::Filters::Image::DefaultTraits< _TImage, _TLabels, typename _TLabels::PixelType > Superclass; + typedef typename Superclass::TInternalTraits TInternalTraits; + typedef typename Superclass::TMarksImage TMarksImage; + typedef typename Superclass::TFilterInterface TFilterInterface; + + typedef fpa::Filters::BaseMarksInterface< TInternalTraits > TMarksInterface; + typedef fpa::Filters::Image::SeedsFromLabelsInterface< TInternalTraits > TSeedsInterface; + }; + + /** + */ + template< class _TImage, class _TLabels, class _TScalarImage > + class MoriLabelling + : public fpa::Filters::Image::RegionGrow< _TImage, _TLabels, typename _TLabels::PixelType, CTBronchi::MoriLabellingTraits< _TImage, _TLabels > > + { + public: + typedef CTBronchi::MoriLabellingTraits< _TImage, _TLabels > TTraits; + fpaTraitsMacro( typename TTraits ); + + typedef fpa::Filters::Image::RegionGrow< _TImage, _TLabels, TMark, TTraits > Superclass; + typedef MoriLabelling Self; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef fpa::Functors::RegionGrow::BinaryThreshold< TInputValue > TFunctor; + + public: + itkNewMacro( Self ); + itkTypeMacro( MoriLabelling, fpa::Filters::Image::RegionGrow ); + + itkGetConstMacro( VesselnessThreshold, double ); + itkSetMacro( VesselnessThreshold, double ); + + itkGetConstMacro( UpperThreshold, TInputValue ); + itkSetMacro( UpperThreshold, TInputValue ); + + ivqITKInputMacro( InputLabels, _TLabels ); + ivqITKInputMacro( InputVesselness, _TScalarImage ); + + protected: + MoriLabelling( ) + : Superclass( ), + m_VesselnessThreshold( 0.05 ), + m_UpperThreshold( -650 ) + { + ivqITKInputConfigureMacro( InputLabels, _TLabels ); + ivqITKInputConfigureMacro( InputVesselness, _TScalarImage ); + this->m_Functor = TFunctor::New( ); + this->SetPredicate( this->m_Functor ); + } + + virtual ~MoriLabelling( ) + { + } + + private: + // Purposely not implemented. + MoriLabelling( const Self& other ); + Self& operator=( const Self& other ); + + protected: + typename TFunctor::Pointer m_Functor; + double m_VesselnessThreshold; + TInputValue m_UpperThreshold; + }; + +} // ecapseman + +#endif // __CTBronchi__Functions__h__ + +// eof - $RCSfile$ diff --git a/appli/CTBronchi/MoriSegmentation.cxx b/appli/CTBronchi/MoriSegmentation.cxx index f5ca8c5..90a1bf3 100644 --- a/appli/CTBronchi/MoriSegmentation.cxx +++ b/appli/CTBronchi/MoriSegmentation.cxx @@ -3,196 +3,119 @@ // @email florez-l@javeriana.edu.co // ========================================================================= -#include -#include +#include +#include +#include #include -#include -#include #include +#include "Functions.h" // ------------------------------------------------------------------------- const unsigned int Dim = 3; -typedef short TPixel; -typedef unsigned char TLabel; +typedef short TPixel; +typedef unsigned char TLabel; typedef itk::Image< TPixel, Dim > TImage; typedef itk::Image< TLabel, Dim > TLabels; // ------------------------------------------------------------------------- -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" ); - mandatory.insert( "seed" ); - - args[ "minimum_threshold" ] = "-850"; - args[ "signal_kernel_size" ] = "20"; - args[ "signal_threshold" ] = "100"; - args[ "signal_influence" ] = "0.5"; - args[ "lower_threshold" ] = "-1024"; - args[ "upper_threshold" ] = "0"; - args[ "delta_threshold" ] = "1"; - - int i = 1; - while( i < argc ) + typedef TCLAP::ValueArg< std::string > _TStringArg; + typedef TCLAP::ValueArg< double > _TRealArg; + typedef TCLAP::ValueArg< unsigned long > _TUIntArg; + typedef TCLAP::ValueArg< TPixel > _TPixelArg; + + // Parse input line + _TStringArg in( "i", "input", "Input image", true, "", "file" ); + _TStringArg out( "o", "output", "Output image", true, "", "file" ); + _TStringArg signal( "s", "signal", "Output signal", false, "", "file" ); + _TStringArg aSeed( "p", "seed", "Input seed", true, "", "\"x y z\"" ); + _TPixelArg mThr( "t", "min_thr", "Minimum threshold", false, -850, "value (-850)" ); + _TUIntArg sKernel( "k", "kernel", "Kernel size", false, 20, "value (20)" ); + _TRealArg sThr( "r", "signal_thr", "Signal threshold", false, 100, "value (100)" ); + _TRealArg sInfluence( "f", "signal_influence", "Signal influence", false, 0.5, "value (0.5)" ); + _TPixelArg lower( "l", "lower", "Lower threshold", false, -1024, "value (-1024)" ); + _TPixelArg upper( "u", "upper", "Upper threshold", false, 0, "value (0)" ); + _TPixelArg delta( "d", "delta", "Delta threshold", false, 1, "value (1)" ); + try { - std::string cmd = argv[ i ] + 1; - if( cmd == "seed" ) - { - std::stringstream seed; - seed - << argv[ i + 1 ] << ";" - << argv[ i + 2 ] << ";" - << argv[ i + 3 ]; - args[ cmd ] = seed.str( ); - i += 4; - } - else - { - args[ cmd ] = argv[ i + 1 ]; - i += 2; - - } // fi - - } // elihw - - bool complete = true; - for( std::string t: mandatory ) - complete &= ( args.find( t ) != args.end( ) ); - - if( !complete ) + TCLAP::CmdLine cmd( "Mori segmentation", ' ', "1.0.0" ); + cmd.add( delta ); + cmd.add( upper ); + cmd.add( lower ); + cmd.add( sInfluence ); + cmd.add( sThr ); + cmd.add( sKernel ); + cmd.add( mThr ); + cmd.add( aSeed ); + cmd.add( signal ); + 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-seed x y z" << std::endl - << "\t[-out_signal filename]" << std::endl - << "\t[-minimum_threshold value]" << std::endl - << "\t[-signal_kernel_size value]" << std::endl - << "\t[-signal_threshold value]" << std::endl - << "\t[-signal_influence value]" << std::endl - << "\t[-lower_threshold value]" << std::endl - << "\t[-upper_threshold value]" << std::endl - << "\t[-delta_threshold value]" << std::endl; - return( false ); + << "===============================" << std::endl + << "Error caught: " << std::endl + << err.error( ) << " " << err.argId( ) + << "===============================" << std::endl + << std::endl; + return( 1 ); - } // fi - return( true ); -} + } // yrt + + // Get seed + std::istringstream sSeed( aSeed.getValue( ) ); + TImage::PointType seed; + sSeed >> seed[ 0 ] >> seed[ 1 ] >> seed[ 2 ]; -// ------------------------------------------------------------------------- -int main( int argc, char* argv[] ) -{ - std::map< std::string, std::string > args; try { - if( ParseArgs( args, argc, argv ) ) + // Read image + TImage::Pointer input_image; + CTBronchi::ReadImage( input_image, in.getValue( ) ); + + // Mori segmentation + typedef fpa::Filters::Image::Mori< TImage, TLabels > _TMori; + _TMori::Pointer mori = _TMori::New( ); + mori->SetInput( input_image ); + mori->SetSeed( seed ); + mori->SetInsideValue( 1 ); + mori->SetOutsideValue( 0 ); + mori->SetMinimumThreshold( mThr.getValue( ) ); + mori->SetSignalKernelSize( sKernel.getValue( ) ); + mori->SetSignalThreshold( sThr.getValue( ) ); + mori->SetSignalInfluence( sInfluence.getValue( ) ); + mori->SetThresholds( + lower.getValue( ), upper.getValue( ), delta.getValue( ) + ); + double t = CTBronchi::MeasureTime( mori ); + std::cout << "Mori executed in " << t << " s" << std::endl; + + // Write result + CTBronchi::WriteImage( mori->GetOutput( ), out.getValue( ) ); + + // Write signal + if( signal.getValue( ) != "" ) { - // Parse seed - TImage::PointType seed; - char* str = new char[ args[ "seed" ].size( ) + 1 ]; - strcpy( str, args[ "seed" ].c_str( ) ); - seed[ 0 ] = std::atof( strtok( str, ";" ) ); - seed[ 1 ] = std::atof( strtok( NULL, ";" ) ); - seed[ 2 ] = std::atof( strtok( NULL, ";" ) ); - - // Read input image - TImage::Pointer input_image; - ReadImage( input_image, args[ "in" ] ); - - // Mori segmentation - typedef fpa::Filters::Image::Mori< TImage, TLabels > TMori; - TMori::Pointer mori = TMori::New( ); - mori->SetInput( input_image ); - mori->SetSeed( seed ); - mori->SetInsideValue( 1 ); - mori->SetOutsideValue( 0 ); - mori->SetMinimumThreshold( - TPixel( std::atof( args[ "minimum_threshold" ].c_str( ) ) ) - ); - mori->SetSignalKernelSize( - std::atoi( args[ "signal_kernel_size" ].c_str( ) ) - ); - mori->SetSignalThreshold( - std::atof( args[ "signal_threshold" ].c_str( ) ) - ); - mori->SetSignalInfluence( - std::atof( args[ "signal_influence" ].c_str( ) ) - ); - mori->SetThresholds( - TPixel( std::atof( args[ "lower_threshold" ].c_str( ) ) ), - TPixel( std::atof( args[ "upper_threshold" ].c_str( ) ) ), - TPixel( std::atof( args[ "delta_threshold" ].c_str( ) ) ) - ); - double t = MeasureTime( mori ); - std::cout << "Mori executed in " << t << " s" << std::endl; - WriteImage( mori->GetOutput( ), args[ "out" ] ); - - std::map< std::string, std::string >::const_iterator i = - args.find( "out_signal" ); - if( i != args.end( ) ) + std::stringstream signalStr; + unsigned long nthr = mori->GetNumberOfEvaluatedThresholds( ); + signalStr << "# nThr = " << nthr << std::endl; + signalStr << "# Opt = " << mori->GetOptimumThreshold( ) << std::endl; + for( unsigned long j = 0; j < nthr; ++j ) { - std::stringstream signal; - unsigned long nthr = mori->GetNumberOfEvaluatedThresholds( ); - signal << "# nThr = " << nthr << std::endl; - signal << "# Opt = " << mori->GetOptimumThreshold( ) << std::endl; - for( unsigned long j = 0; j < nthr; ++j ) - { - typename TMori::TPeak p; - double x, y; - mori->GetSignalValues( j, x, y, p ); - signal << x << " " << y << std::endl; + typename _TMori::TPeak p; + double x, y; + mori->GetSignalValues( j, x, y, p ); + signalStr << x << " " << y << std::endl; - } // rof - std::ofstream signals_str( i->second.c_str( ) ); - signals_str << signal.str( ); - signals_str.close( ); + } // rof + std::ofstream oStrSignals( signal.getValue( ).c_str( ) ); + oStrSignals << signalStr.str( ); + oStrSignals.close( ); - } // fi - } - else - return( 1 ); + } // fi } catch( std::exception& err ) { 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 ) {