From 5be07a28290fed315829547221848a0ee90558d1 Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Thu, 5 Oct 2017 22:38:14 -0500 Subject: [PATCH] ... --- appli/CMakeLists.txt | 2 +- appli/CTBronchi/CMakeLists.txt | 39 +++-- appli/CTBronchi/MoriLabelling.cxx | 46 +++++- appli/CTBronchi/Vesselness.cxx | 168 +++++++++++++++++++++ examples/image/RandomWalker/Gaussian.cxx | 6 +- lib/fpa/Functors/Dijkstra/Image/Gaussian.h | 5 +- 6 files changed, 231 insertions(+), 35 deletions(-) create mode 100644 appli/CTBronchi/Vesselness.cxx diff --git a/appli/CMakeLists.txt b/appli/CMakeLists.txt index 898efa2..7e5406b 100644 --- a/appli/CMakeLists.txt +++ b/appli/CMakeLists.txt @@ -6,6 +6,6 @@ include_directories( ${PROJECT_SOURCE_DIR}/lib ${PROJECT_BINARY_DIR}/lib ) -subdirs(CTArteries) +subdirs(CTArteries CTBronchi) ## eof - $RCSfile$ diff --git a/appli/CTBronchi/CMakeLists.txt b/appli/CTBronchi/CMakeLists.txt index fccfb97..d81b5e8 100644 --- a/appli/CTBronchi/CMakeLists.txt +++ b/appli/CTBronchi/CMakeLists.txt @@ -1,26 +1,25 @@ +## ========================================================================= +## @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co) +## ========================================================================= + option(fpa_BUILD_CTBronchi "Build bronchi analysis from CT images applications?" OFF) if(fpa_BUILD_CTBronchi) - configure_file( - CTBronchi_process.sh - ${PROJECT_BINARY_DIR}/CTBronchi_process.sh - COPYONLY + set(_pfx fpa_CTBronchi_) + set( + _examples + Vesselness + MoriSegmentation + MoriLabelling ) - include_directories( - ${PROJECT_SOURCE_DIR}/lib - ${PROJECT_BINARY_DIR}/lib - ${PROJECT_SOURCE_DIR}/appli - ${PROJECT_BINARY_DIR}/appli - ) -set(_pfx fpa_CTBronchi_) -set( - _examples - MoriSegmentation - MoriLabelling - ) -foreach(_e ${_examples}) - add_executable(${_pfx}${_e} ${_e}.cxx) - target_link_libraries(${_pfx}${_e} fpa) -endforeach(_e) + foreach(_e ${_examples}) + BuildApplication( + ${_pfx}${_e} + SOURCE ${_e}.cxx + INSTALL + RECURRENT + LINKS fpa + ) + endforeach(_e) endif(fpa_BUILD_CTBronchi) ## eof - $RCSfile$ diff --git a/appli/CTBronchi/MoriLabelling.cxx b/appli/CTBronchi/MoriLabelling.cxx index 96575b4..78544b5 100644 --- a/appli/CTBronchi/MoriLabelling.cxx +++ b/appli/CTBronchi/MoriLabelling.cxx @@ -10,6 +10,7 @@ #include #include #include +#include #include #include @@ -19,10 +20,12 @@ // ------------------------------------------------------------------------- const unsigned int Dim = 3; -typedef short TPixel; -typedef unsigned char TLabel; -typedef itk::Image< TPixel, Dim > TImage; -typedef itk::Image< TLabel, Dim > TLabels; +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; /** */ @@ -65,6 +68,7 @@ public: itkGetConstMacro( MaxVertex, TVertex ); fpaFilterInputMacro( InputLabels, TLabels ); + fpaFilterInputMacro( InputVesselness, TScalarImage ); public: TInputValue GetUpperThreshold( ) const @@ -79,9 +83,11 @@ public: protected: MoriLabelling( ) : Superclass( ), - m_LastThreshold( TInputValue( 0 ) ) + m_LastThreshold( TInputValue( 0 ) ), + m_VesselnessThr( TScalar( 0.05 ) ) { fpaFilterInputConfigureMacro( InputLabels, TLabels ); + fpaFilterInputConfigureMacro( InputVesselness, TScalarImage ); this->m_Functor = TFunctor::New( ); this->SetPredicate( this->m_Functor ); } @@ -98,6 +104,12 @@ protected: { 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 @@ -107,13 +119,20 @@ protected: { 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 /* && b < a*/ ) - n.Value = 0; + 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 ) ) { @@ -148,6 +167,9 @@ protected: bool m_FirstVertex; TVertex m_MinVertex; TVertex m_MaxVertex; + + TScalar m_MaxVesselness; + TScalar m_VesselnessThr; }; // ------------------------------------------------------------------------- @@ -217,8 +239,9 @@ bool ParseArgs( mandatory.insert( "in" ); mandatory.insert( "out" ); mandatory.insert( "labels" ); + mandatory.insert( "vesselness" ); - args[ "upper_threshold" ] = "-600"; + args[ "upper_threshold" ] = "-650"; int i = 1; while( i < argc ) @@ -240,6 +263,7 @@ bool ParseArgs( << "\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 ); @@ -264,11 +288,17 @@ int main( int argc, char* argv[] ) 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( ) ) ) diff --git a/appli/CTBronchi/Vesselness.cxx b/appli/CTBronchi/Vesselness.cxx new file mode 100644 index 0000000..ca7c99c --- /dev/null +++ b/appli/CTBronchi/Vesselness.cxx @@ -0,0 +1,168 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 3; +typedef short TPixel; +typedef itk::NumericTraits< TPixel >::RealType TScalar; +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[] + ) +{ + 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 ) + { + 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[-sigma val]" << std::endl + << "\t[-alpha1 val]" << std::endl + << "\t[-alpha2 val]" << 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" ] ); + + // 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 ); + } + 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/examples/image/RandomWalker/Gaussian.cxx b/examples/image/RandomWalker/Gaussian.cxx index 5c2c798..9a012e9 100644 --- a/examples/image/RandomWalker/Gaussian.cxx +++ b/examples/image/RandomWalker/Gaussian.cxx @@ -10,10 +10,10 @@ #include // ------------------------------------------------------------------------- -const unsigned int Dim = 2; -typedef unsigned char TPixel; +const unsigned int Dim = 3; +typedef short TPixel; typedef unsigned char TLabel; -typedef float TCost; +typedef double TCost; typedef itk::Image< TPixel, Dim > TImage; typedef itk::Image< TLabel, Dim > TLabelImage; diff --git a/lib/fpa/Functors/Dijkstra/Image/Gaussian.h b/lib/fpa/Functors/Dijkstra/Image/Gaussian.h index b819f01..8fa04be 100644 --- a/lib/fpa/Functors/Dijkstra/Image/Gaussian.h +++ b/lib/fpa/Functors/Dijkstra/Image/Gaussian.h @@ -62,9 +62,8 @@ namespace fpa TValue d = TValue( image->GetPixel( v ) ); d -= TValue( image->GetPixel( p ) ); d /= this->m_Beta; - d *= d; - if( this->m_TreatAsWeight ) d = std::exp( d ) - TValue( 1 ); - else d = std::exp( -d ); + if( this->m_TreatAsWeight ) d = std::exp( d * d ) - TValue( 1 ); + else d = std::exp( -std::fabs( d ) ); if( d < this->m_Epsilon ) return( this->m_Epsilon ); else return( d ); -- 2.45.1