X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=appli%2FCTBronchi%2FMoriLabelling.cxx;h=96575b420d5935bfa6aa8c529216b7cf3a860639;hb=3bcd064e000912302f2f0d3556d518ff46c146aa;hp=0915cd5f1cc0448a0b2dda4816cd96ed1cc68fda;hpb=6f716a6d2ce76a7bc9abb37406ce3802aa1e4451;p=FrontAlgorithms.git diff --git a/appli/CTBronchi/MoriLabelling.cxx b/appli/CTBronchi/MoriLabelling.cxx index 0915cd5..96575b4 100644 --- a/appli/CTBronchi/MoriLabelling.cxx +++ b/appli/CTBronchi/MoriLabelling.cxx @@ -1,117 +1,324 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + #include +#include + #include #include #include -#include +#include + +#include +#include +#include +#include +#include // ------------------------------------------------------------------------- const unsigned int Dim = 3; -typedef short TPixel; -typedef unsigned short TLabel; +typedef short TPixel; +typedef unsigned char TLabel; +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::Image< TLabel, Dim > TLabels; + +/** + */ +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 ); + + fpaFilterInputMacro( InputLabels, TLabels ); + +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 ) ) + { + fpaFilterInputConfigureMacro( InputLabels, TLabels ); + 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::Image< TPixel, Dim > TInputImage; -typedef itk::Image< TLabel, Dim > TLabelImage; + 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( ); + 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; -typedef CTBronchi::MoriLabelling< TInputImage, TLabelImage > TFilter; + 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; +}; // ------------------------------------------------------------------------- -int main( int argc, char* argv[] ) +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 + ) { - // Get arguments - if( argc < 8 ) + 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" ); + + args[ "upper_threshold" ] = "-600"; + + 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 - << " input_image label_image output_image" << std::endl - << " upper_threshold(-400)" << std::endl - << " inside_value(255)" << std::endl - << " inside_label(1) outside_label(2)" - << std::endl; - return( 1 ); + << "\t-in filename" << std::endl + << "\t-out filename" << std::endl + << "\t-labels filename" << std::endl + << "\t[-roi] filename" << std::endl + << "\t[-upper_threshold value]" << std::endl; + return( false ); } // fi - std::string input_image_filename = argv[ 1 ]; - std::string label_image_filename = argv[ 2 ]; - std::string output_image_filename = argv[ 3 ]; - TLabel upper_threshold = std::atoi( argv[ 4 ] ); - TLabel inside_value = std::atoi( argv[ 5 ] ); - TLabel inside_label = std::atoi( argv[ 6 ] ); - TLabel outside_label = std::atoi( argv[ 7 ] ); - - // Read images - itk::ImageFileReader< TInputImage >::Pointer input_image_reader = - itk::ImageFileReader< TInputImage >::New( ); - input_image_reader->SetFileName( input_image_filename ); - - itk::ImageFileReader< TLabelImage >::Pointer label_image_reader = - itk::ImageFileReader< TLabelImage >::New( ); - label_image_reader->SetFileName( label_image_filename ); - - // Prepare filter - TFilter::Pointer filter = TFilter::New( ); - filter->SetInputLabelImage( label_image_reader->GetOutput( ) ); - filter->SetInputRawImage( input_image_reader->GetOutput( ) ); - filter->SetUpperThreshold( upper_threshold ); - filter->SetInsideValue( inside_value ); - filter->SetInsideLabel( inside_label ); - filter->SetOutsideLabel( outside_label ); - - // Show some information - std::cout << "----------------------------------------------" << std::endl; - std::cout << "Image: " << input_image_filename << std::endl; - - // Execute pipeline - std::chrono::time_point< std::chrono::high_resolution_clock > ts, te; - std::chrono::duration< double > tr; + return( true ); +} + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + std::map< std::string, std::string > args; try { - ts = std::chrono::high_resolution_clock::now( ); - input_image_reader->Update( ); - te = std::chrono::high_resolution_clock::now( ); - tr = te - ts; - std::cout << "Raw read time: " << tr.count( ) << " s" << std::endl; - - ts = std::chrono::high_resolution_clock::now( ); - label_image_reader->Update( ); - te = std::chrono::high_resolution_clock::now( ); - tr = te - ts; - std::cout << "Label read time: " << tr.count( ) << " s" << std::endl; - - ts = std::chrono::high_resolution_clock::now( ); - filter->Update( ); - te = std::chrono::high_resolution_clock::now( ); - tr = te - ts; - std::cout - << "Labelling time: " << tr.count( ) << " s" << std::endl; - } - catch( std::exception& err ) - { - std::cerr << "Error caught: " << err.what( ) << std::endl; - return( 1 ); + if( ParseArgs( args, argc, argv ) ) + { + // Read input image + TImage::Pointer input_image; + ReadImage( input_image, args[ "in" ] ); - } // yrt + // Read labels image + TLabels::Pointer input_labels; + ReadImage( input_labels, args[ "labels" ] ); - // Save output image - itk::ImageFileWriter< TLabelImage >::Pointer output_image_writer = - itk::ImageFileWriter< TLabelImage >::New( ); - output_image_writer->SetInput( filter->GetOutput( ) ); - output_image_writer->SetFileName( output_image_filename ); - try - { - ts = std::chrono::high_resolution_clock::now( ); - output_image_writer->Update( ); - te = std::chrono::high_resolution_clock::now( ); - tr = te - ts; - std::cout << "Write time: " << tr.count( ) << " s" << std::endl; + // Mori labelling + MoriLabelling::Pointer labelling = MoriLabelling::New( ); + labelling->SetInput( input_image ); + labelling->SetInputLabels( input_labels ); + labelling->SetOutsideValue( 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 ); } catch( std::exception& err ) { - std::cerr << "Error caught: " << err.what( ) << std::endl; + std::cerr + << "===============================" << std::endl + << "Error caught: " << std::endl + << err.what( ) + << "===============================" << std::endl + << std::endl; return( 1 ); } // yrt - std::cout << "----------------------------------------------" << std::endl; - return( 0 ); }