From 4daffc93df1565a88ca3eae48567475a6fba8319 Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Wed, 11 Oct 2017 20:44:56 -0500 Subject: [PATCH] ... --- appli/CTBronchi/CMakeLists.txt | 1 + appli/CTBronchi/FastRandomWalker.cxx | 103 ++++++++ appli/CTBronchi/MoriLabelling.cxx | 343 +-------------------------- appli/CTBronchi/MoriLabelling.h | 40 ++++ 4 files changed, 145 insertions(+), 342 deletions(-) create mode 100644 appli/CTBronchi/FastRandomWalker.cxx diff --git a/appli/CTBronchi/CMakeLists.txt b/appli/CTBronchi/CMakeLists.txt index a5c2b61..c88add4 100644 --- a/appli/CTBronchi/CMakeLists.txt +++ b/appli/CTBronchi/CMakeLists.txt @@ -10,6 +10,7 @@ if(fpa_BUILD_CTBronchi) Vesselness MoriSegmentation MoriLabelling + FastRandomWalker ) foreach(_e ${_examples}) BuildApplication( diff --git a/appli/CTBronchi/FastRandomWalker.cxx b/appli/CTBronchi/FastRandomWalker.cxx new file mode 100644 index 0000000..b2c901b --- /dev/null +++ b/appli/CTBronchi/FastRandomWalker.cxx @@ -0,0 +1,103 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#include +#include +#include +#include "MoriLabelling.h" +#include "Functions.h" + +// ------------------------------------------------------------------------- +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; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + 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, -400, "value (-400)" ); + try + { + 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 + << "===============================" << std::endl + << "Error caught: " << std::endl + << err.error( ) << " " << err.argId( ) << std::endl + << "===============================" << std::endl + << std::endl; + return( 1 ); + + } // yrt + + try + { + // 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( ) << std::endl + << "===============================" << std::endl + << std::endl; + return( 1 ); + + } // yrt + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/appli/CTBronchi/MoriLabelling.cxx b/appli/CTBronchi/MoriLabelling.cxx index 96f4dc4..b2c901b 100644 --- a/appli/CTBronchi/MoriLabelling.cxx +++ b/appli/CTBronchi/MoriLabelling.cxx @@ -3,7 +3,6 @@ // @email florez-l@javeriana.edu.co // ========================================================================= -#include #include #include #include @@ -32,7 +31,7 @@ int main( int argc, char* argv[] ) _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)" ); + _TPixelArg uThr( "u", "upper_thr", "Upper threshold", false, -400, "value (-400)" ); try { TCLAP::CmdLine cmd( "Labelling", ' ', "1.0.0" ); @@ -101,344 +100,4 @@ 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 index 89ce001..ff6553e 100644 --- a/appli/CTBronchi/MoriLabelling.h +++ b/appli/CTBronchi/MoriLabelling.h @@ -5,6 +5,7 @@ #ifndef __CTBronchi__MoriLabelling__h__ #define __CTBronchi__MoriLabelling__h__ +#include #include #include #include @@ -75,6 +76,44 @@ namespace CTBronchi { } + virtual const itk::DataObject* _GetReferenceInput( ) const override + { + return( this->GetInputLabels( ) ); + } + virtual void _BeforeGenerateData( ) override + { + this->Superclass::_BeforeGenerateData( ); + + this->m_Functor->SetUpperThreshold( this->m_UpperThreshold ); + + typedef itk::MinimumMaximumImageCalculator< _TScalarImage > _TMinMax; + typename _TMinMax::Pointer minMax = _TMinMax::New( ); + minMax->SetImage( this->GetInputVesselness( ) ); + minMax->Compute( ); + this->m_MinVesselness = + ( 1.0 - this->m_VesselnessThreshold ) * + double( minMax->GetMaximum( ) ); + } + + virtual void _PostComputeOutputValue( TNode& n ) override + { + this->Superclass::_PostComputeOutputValue( n ); + if( n.Value == this->GetInsideValue( ) ) + { + const _TLabels* labels = this->GetInputLabels( ); + const _TScalarImage* vesselness = this->GetInputVesselness( ); + if( labels->GetPixel( n.Vertex ) == 0 ) + { + if( this->m_MinVesselness < vesselness->GetPixel( n.Vertex ) ) + n.Value = this->GetInsideValue( ); + else + n.Value = 0; + + } // fi + + } // fi + } + private: // Purposely not implemented. MoriLabelling( const Self& other ); @@ -83,6 +122,7 @@ namespace CTBronchi protected: typename TFunctor::Pointer m_Functor; double m_VesselnessThreshold; + double m_MinVesselness; TInputValue m_UpperThreshold; }; -- 2.45.1