X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=appli%2FCTBronchi%2FMoriLabelling.cxx;h=78544b5ed98b094493c71e849f4c26e09721e9bb;hb=ae0e1b8916a0fb2188080b9134c1c2781c6c200f;hp=0ad694d9027a1249f56f76c6b760c447fa499841;hpb=d2f5533aef455b7d862dbfbe72c42c5497eaf98b;p=FrontAlgorithms.git diff --git a/appli/CTBronchi/MoriLabelling.cxx b/appli/CTBronchi/MoriLabelling.cxx index 0ad694d..78544b5 100644 --- a/appli/CTBronchi/MoriLabelling.cxx +++ b/appli/CTBronchi/MoriLabelling.cxx @@ -5,17 +5,172 @@ #include #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::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; + +/** + */ +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 ); + fpaFilterInputMacro( 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 ) ) + { + fpaFilterInputConfigureMacro( InputLabels, TLabels ); + fpaFilterInputConfigureMacro( 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 ) @@ -45,7 +200,7 @@ void ReadImage( _TImagePtr& image, const std::string& fname ) // ------------------------------------------------------------------------- template< class _TImage > -void WriteImage( const _TImage*image, const std::string& fname ) +void WriteImage( const _TImage* image, const std::string& fname ) { typedef itk::ImageFileWriter< _TImage > _TWriter; typename _TWriter::Pointer writer = _TWriter::New( ); @@ -55,6 +210,26 @@ void WriteImage( const _TImage*image, const std::string& fname ) 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[] @@ -64,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 ) @@ -87,6 +263,8 @@ 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 ); @@ -110,57 +288,52 @@ int main( int argc, char* argv[] ) TLabels::Pointer input_labels; ReadImage( input_labels, args[ "labels" ] ); - // Mori segmentation + // 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 - 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 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; - - } // rof - std::ofstream signals_str( i->second.c_str( ) ); - signals_str << signal.str( ); - signals_str.close( ); - - } // fi + 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 );