From: Leonardo Flórez-Valencia Date: Fri, 22 Sep 2017 09:36:11 +0000 (+0200) Subject: ... X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=30e4529b4758ee7b397c4e839ac0e3334adb4010;p=FrontAlgorithms.git ... --- diff --git a/appli/CTBronchi/CTBronchi.cxx b/appli/CTBronchi/CTBronchi.cxx deleted file mode 100644 index 6427532..0000000 --- a/appli/CTBronchi/CTBronchi.cxx +++ /dev/null @@ -1,8 +0,0 @@ - -int main( int argc, char* argv[] ) -{ - return( 0 ); -} - - -// eof - $RCSfile$ diff --git a/appli/CTBronchi/CTBronchi.ui b/appli/CTBronchi/CTBronchi.ui deleted file mode 100644 index 4fa5f30..0000000 --- a/appli/CTBronchi/CTBronchi.ui +++ /dev/null @@ -1,123 +0,0 @@ - - - MainWindow - - - - 0 - 0 - 800 - 600 - - - - MainWindow - - - - - - 60 - 150 - 461 - 331 - - - - - - - - 0 - 0 - 800 - 20 - - - - - File - - - - - - - - - - - - - - - - - - - - toolBar - - - TopToolBarArea - - - false - - - - - - - Open raw image - - - - - Open "mori" image - - - - - Open "random walker" image - - - - - Exit - - - - - Open seeds file - - - - - Save seeds file - - - - - Save "mori" image - - - - - Save "random walker" image - - - - - Mori - - - - - RandomWalker - - - - - - diff --git a/appli/CTBronchi/Mori.h b/appli/CTBronchi/Mori.h deleted file mode 100644 index 4693d19..0000000 --- a/appli/CTBronchi/Mori.h +++ /dev/null @@ -1,89 +0,0 @@ -// ========================================================================= -// @author Leonardo Florez Valencia -// @email florez-l@javeriana.edu.co -// ========================================================================= -#ifndef __CTBronchi__Mori__h__ -#define __CTBronchi__Mori__h__ - -#include -#include -#include - -namespace CTBronchi -{ - // ----------------------------------------------------------------------- - template< class _TInputPtr, class _TOutputPtr, class _TSeed > - typename _TInputPtr::ObjectType::PixelType Mori( - _TOutputPtr& output, const _TInputPtr& input, _TSeed& seed, - std::map< std::string, std::string >& args - ) - { - typedef typename _TInputPtr::ObjectType _TInput; - typedef typename _TOutputPtr::ObjectType _TOutput; - typedef fpa::Filters::Image::Mori< _TInput, _TOutput > _TMori; - typedef typename _TInput::PixelType _TInputPixel; - typedef typename _TOutput::PixelType _TOutputPixel; - - _TOutputPixel i_val = _TOutputPixel( 1 ); - _TOutputPixel o_val = _TOutputPixel( 0 ); - - typename _TMori::Pointer mori = _TMori::New( ); - mori->SetInput( input ); - mori->SetSeed( seed ); - mori->SetInsideValue( i_val ); - mori->SetOutsideValue( o_val ); - mori->SetMinimumThreshold( - _TInputPixel( std::atof( args[ "mori_minimum_threshold" ].c_str( ) ) ) - ); - mori->SetSignalKernelSize( - std::atoi( args[ "mori_signal_kernel_size" ].c_str( ) ) - ); - mori->SetSignalThreshold( - std::atof( args[ "mori_signal_threshold" ].c_str( ) ) - ); - mori->SetSignalInfluence( - std::atof( args[ "mori_signal_influence" ].c_str( ) ) - ); - mori->SetThresholds( - _TInputPixel( std::atof( args[ "mori_lower_threshold" ].c_str( ) ) ), - _TInputPixel( std::atof( args[ "mori_upper_threshold" ].c_str( ) ) ), - _TInputPixel( std::atof( args[ "mori_delta_threshold" ].c_str( ) ) ) - ); - double t = CTBronchi::MeasureTime( mori ); - std::cout << "Mori executed in " << t << " s" << std::endl; - output = mori->GetOutput( ); - - std::map< std::string, std::string >::const_iterator i = - args.find( "out_mori" ); - if( i != args.end( ) ) - CTBronchi::WriteImage( output, i->second ); - - 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 - output->DisconnectPipeline( ); - return( mori->GetOptimumThreshold( ) ); - } - -} // ecapseman - -#endif // __CTBronchi__Functions__h__ - -// eof - $RCSfile$ diff --git a/appli/CTBronchi/MoriLabelling.cxx b/appli/CTBronchi/MoriLabelling.cxx index 0ad694d..96575b4 100644 --- a/appli/CTBronchi/MoriLabelling.cxx +++ b/appli/CTBronchi/MoriLabelling.cxx @@ -5,10 +5,17 @@ #include #include + #include #include #include -#include +#include + +#include +#include +#include +#include +#include // ------------------------------------------------------------------------- const unsigned int Dim = 3; @@ -17,6 +24,132 @@ 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; + } + + 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; + + 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; +}; + // ------------------------------------------------------------------------- double MeasureTime( itk::ProcessObject* f ) { @@ -45,7 +178,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 +188,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[] @@ -87,6 +240,7 @@ bool ParseArgs( << "\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 ); @@ -110,57 +264,46 @@ int main( int argc, char* argv[] ) TLabels::Pointer input_labels; ReadImage( input_labels, args[ "labels" ] ); - // Mori segmentation + // 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 - 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 ); diff --git a/appli/CTBronchi/MoriLabelling.h b/appli/CTBronchi/MoriLabelling.h deleted file mode 100644 index d07968b..0000000 --- a/appli/CTBronchi/MoriLabelling.h +++ /dev/null @@ -1,91 +0,0 @@ -// ========================================================================= -// @author Leonardo Florez Valencia -// @email florez-l@javeriana.edu.co -// ========================================================================= -#ifndef __CTBronchi__MoriLabelling__h__ -#define __CTBronchi__MoriLabelling__h__ - -#include -#include -#include -#include -#include -#include - -namespace CTBronchi -{ - /** - */ - template< class _TInputImage, class _TLabelImage > - class MoriLabellingTraits - : public fpa::Filters::Image::DefaultTraits< _TInputImage, _TLabelImage, typename _TLabelImage::PixelType > - { - public: - typedef _TInputImage TInputImage; - typedef _TLabelImage TLabelImage; - typedef typename TLabelImage::PixelType TLabel; - typedef fpa::Filters::Image::DefaultTraits< TInputImage, TLabelImage, 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; - }; - - /** - */ - template< class _TInputImage, class _TLabelImage > - class MoriLabelling - : public fpa::Filters::Image::RegionGrow< _TInputImage, _TLabelImage, typename _TLabelImage::PixelType, CTBronchi::MoriLabellingTraits< _TInputImage, _TLabelImage > > - { - public: - typedef _TInputImage TInputImage; - typedef _TLabelImage TLabelImage; - typedef CTBronchi::MoriLabellingTraits< TInputImage, TLabelImage > TTraits; - fpaTraitsMacro( typename TTraits ); - - typedef fpa::Filters::Image::RegionGrow< TInputImage, TLabelImage, 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( CTBronchi::MoriLabelling, fpa::Filters::Image::RegionGrow ); - - itkGetConstMacro( LastThreshold, TInputValue ); - itkSetMacro( LastThreshold, TInputValue ); - - fpaFilterInputMacro( InputLabels, TLabelImage ); - - public: - TInputValue GetUpperThreshold( ) const; - void SetUpperThreshold( TInputValue t ); - - protected: - MoriLabelling( ); - virtual ~MoriLabelling( ); - - virtual const itk::DataObject* _GetReferenceInput( ) const override; - virtual void _PostComputeOutputValue( TNode& n ) override; - - private: - // Purposely not implemented. - MoriLabelling( const Self& other ); - Self& operator=( const Self& other ); - - protected: - typename TFunctor::Pointer m_Functor; - TInputValue m_LastThreshold; - }; - -} // ecapseman - -#ifndef ITK_MANUAL_INSTANTIATION -# include -#endif // ITK_MANUAL_INSTANTIATION -#endif // __CTBronchi__MoriLabelling__h__ - -// eof - $RCSfile$ diff --git a/appli/CTBronchi/MoriLabelling.hxx b/appli/CTBronchi/MoriLabelling.hxx deleted file mode 100644 index f786948..0000000 --- a/appli/CTBronchi/MoriLabelling.hxx +++ /dev/null @@ -1,111 +0,0 @@ -// ========================================================================= -// @author Leonardo Florez Valencia -// @email florez-l@javeriana.edu.co -// ========================================================================= -#ifndef __CTBronchi__MoriLabelling__hxx__ -#define __CTBronchi__MoriLabelling__hxx__ - -// ------------------------------------------------------------------------- -template< class _TInputImage, class _TLabelImage > -typename CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: -TInputValue CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: -GetUpperThreshold( ) const -{ - return( this->m_Functor->GetUpperThreshold( ) ); -} - -// ------------------------------------------------------------------------- -template< class _TInputImage, class _TLabelImage > -void CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: -SetUpperThreshold( TInputValue t ) -{ - this->m_Functor->SetUpperThreshold( t ); -} - -// ------------------------------------------------------------------------- -template< class _TInputImage, class _TLabelImage > -CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: -MoriLabelling( ) - : Superclass( ), - m_LastThreshold( TInputValue( 0 ) ) -{ - fpaFilterInputConfigureMacro( InputLabels, TLabelImage ); - this->m_Functor = TFunctor::New( ); - this->SetPredicate( this->m_Functor ); -} - -// ------------------------------------------------------------------------- -template< class _TInputImage, class _TLabelImage > -CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: -~MoriLabelling( ) -{ -} - -// ------------------------------------------------------------------------- -template< class _TInputImage, class _TLabelImage > -const itk::DataObject* -CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: -_GetReferenceInput( ) const -{ - return( this->GetInputLabels( ) ); -} - -// ------------------------------------------------------------------------- -template< class _TInputImage, class _TLabelImage > -void CTBronchi::MoriLabelling< _TInputImage, _TLabelImage >:: -_PostComputeOutputValue( TNode& n ) -{ - this->Superclass::_PostComputeOutputValue( n ); - if( n.Value == this->GetInsideValue( ) ) - { - const _TInputImage* input = this->GetInput( ); - const TLabelImage* 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; - - } // fi -} - -// ------------------------------------------------------------------------- -namespace CTBronchi -{ - template< class _TInputPtr, class _TOutputPtr, class _TInputValue > - void Label( - _TOutputPtr& output, - const _TInputPtr& input, const _TOutputPtr& input_labels, - const _TInputValue& last_thr, - std::map< std::string, std::string >& args - ) - { - typedef typename _TInputPtr::ObjectType _TInput; - typedef typename _TOutputPtr::ObjectType _TOutput; - typedef CTBronchi::MoriLabelling< _TInput, _TOutput > _TLabelling; - - typename _TLabelling::Pointer labelling = _TLabelling::New( ); - labelling->SetInput( input ); - labelling->SetInputLabels( input_labels ); - labelling->SetOutsideValue( 2 ); - labelling->SetInsideValue( 1 ); - labelling->SetUpperThreshold( - _TInputValue( std::atof( args[ "labelling_upper_threshold" ].c_str( ) ) ) - ); - labelling->SetLastThreshold( last_thr ); - double t = MeasureTime( labelling ); - std::cout << "Labelling executed in " << t << " s" << std::endl; - output = labelling->GetOutput( ); - std::map< std::string, std::string >::const_iterator i = - args.find( "out_labels" ); - if( i != args.end( ) ) - WriteImage( output, i->second ); - output->DisconnectPipeline( ); - } - -} // ecapseman -#endif // __CTBronchi__MoriLabelling__hxx__ - -// eof - $RCSfile$ diff --git a/appli/CTBronchi/MoriSegmentation.cxx b/appli/CTBronchi/MoriSegmentation.cxx index dcbae23..f5ca8c5 100644 --- a/appli/CTBronchi/MoriSegmentation.cxx +++ b/appli/CTBronchi/MoriSegmentation.cxx @@ -45,7 +45,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( ); diff --git a/appli/CTBronchi/RandomWalker.h b/appli/CTBronchi/RandomWalker.h deleted file mode 100644 index c284a54..0000000 --- a/appli/CTBronchi/RandomWalker.h +++ /dev/null @@ -1,47 +0,0 @@ -// ========================================================================= -// @author Leonardo Florez Valencia -// @email florez-l@javeriana.edu.co -// ========================================================================= -#ifndef __CTBronchi__RandomWalker__h__ -#define __CTBronchi__RandomWalker__h__ - -#include -#include -#include -#include - -namespace CTBronchi -{ - // ----------------------------------------------------------------------- - template< class _TInputPtr, class _TLabelPtr > - void RandomWalker( - _TLabelPtr& output, - const _TInputPtr& input, const _TLabelPtr& labels, - std::map< std::string, std::string >& args - ) - { - typedef typename _TInputPtr::ObjectType _TInput; - typedef typename _TLabelPtr::ObjectType _TLabel; - typedef float _TScalar; - typedef fpa::Filters::Image::RandomWalker< _TInput, _TLabel, _TScalar > _TRandomWalker; - typedef fpa::Functors::Dijkstra::Image::Gaussian< _TInput, _TScalar > _TWeight; - - typename _TWeight::Pointer weight = _TWeight::New( ); - weight->SetAlpha( std::atof( args[ "random_walker_alpha" ].c_str( ) ) ); - weight->SetBeta( std::atof( args[ "random_walker_beta" ].c_str( ) ) ); - - typename _TRandomWalker::Pointer rw = _TRandomWalker::New( ); - rw->SetInputImage( input ); - rw->SetInputLabels( labels ); - rw->SetWeightFunction( weight ); - double t = CTBronchi::MeasureTime( rw ); - std::cout << "RandomWalker executed in " << t << " s" << std::endl; - output = rw->GetOutputLabels( ); - output->DisconnectPipeline( ); - } - -} // ecapseman - -#endif // __CTBronchi__Functions__h__ - -// eof - $RCSfile$ diff --git a/lib/fpa/Common/OriginalRandomWalker.hxx b/lib/fpa/Common/OriginalRandomWalker.hxx index 46dd935..2d3dc73 100644 --- a/lib/fpa/Common/OriginalRandomWalker.hxx +++ b/lib/fpa/Common/OriginalRandomWalker.hxx @@ -58,57 +58,80 @@ GenerateData( ) // Allocate outputs this->AllocateOutputs( ); - // Build boundary triplets and count labels - _TTriplets St, Bt; - std::map< TLabel, unsigned long > labels; - this->_Boundary( St, labels ); - struct _TTripletsOrd - { - bool operator()( const _TTriplet& a, const _TTriplet& b ) - { - return( a.row( ) < b.row( ) ); - } - }; - std::sort( St.begin( ), St.end( ), _TTripletsOrd( ) ); - for( unsigned long i = 0; i < St.size( ); ++i ) - Bt.push_back( _TTriplet( i, labels[ St[ i ].col( ) ], St[ i ].value( ) ) ); - - // Laplacian triplets - _TTriplets At, Rt; - this->_Laplacian( At, Rt, St ); - - // Matrices - TRegion region = input->GetRequestedRegion( ); - unsigned long nSeeds = St.size( ); - unsigned long nLabels = labels.size( ); - unsigned long N = region.GetNumberOfPixels( ); - - std::vector< TLabel > invLabels( nLabels ); - for( typename std::map< TLabel, unsigned long >::value_type s: labels ) - invLabels[ s.second ] = s.first; - - _TMatrix B( nSeeds, nLabels ); - B.setFromTriplets( Bt.begin( ), Bt.end( ) ); - B.makeCompressed( ); - - _TMatrix R( N - nSeeds, nSeeds ); - R.setFromTriplets( Rt.begin( ), Rt.end( ) ); - R.makeCompressed( ); - - _TMatrix A( N - nSeeds, N - nSeeds ); - A.setFromTriplets( At.begin( ), At.end( ) ); - A.makeCompressed( ); + // Persisting objects + _TMatrix A( 1, 1 ), C( 1, 1 ); + _TTriplets St; + std::vector< TLabel > invLabels; + + { // begin + // Build boundary triplets and count labels + _TTriplets Bt; + std::map< TLabel, unsigned long > labels; + itkDebugMacro( << "Building boundary matrix..." ); + this->_Boundary( St, labels ); + struct _TTripletsOrd + { + bool operator()( const _TTriplet& a, const _TTriplet& b ) + { + return( a.row( ) < b.row( ) ); + } + }; + itkDebugMacro( << "Sorting boundary pixels..." ); + std::sort( St.begin( ), St.end( ), _TTripletsOrd( ) ); + itkDebugMacro( << "Assigning boundary pixels..." ); + for( unsigned long i = 0; i < St.size( ); ++i ) + Bt.push_back( + _TTriplet( i, labels[ St[ i ].col( ) ], St[ i ].value( ) ) + ); + + // Laplacian triplets + itkDebugMacro( << "Building laplacian matrix..." ); + _TTriplets At, Rt; + this->_Laplacian( At, Rt, St ); + + // Matrices + TRegion region = input->GetRequestedRegion( ); + unsigned long nSeeds = St.size( ); + unsigned long nLabels = labels.size( ); + unsigned long N = region.GetNumberOfPixels( ); + + itkDebugMacro( << "Creating inverse labels..." ); + invLabels.resize( nLabels ); + for( typename std::map< TLabel, unsigned long >::value_type s: labels ) + invLabels[ s.second ] = s.first; + + itkDebugMacro( << "Creating B matrix..." ); + _TMatrix B( nSeeds, nLabels ); + B.setFromTriplets( Bt.begin( ), Bt.end( ) ); + B.makeCompressed( ); + + itkDebugMacro( << "Creating R matrix..." ); + _TMatrix R( N - nSeeds, nSeeds ); + R.setFromTriplets( Rt.begin( ), Rt.end( ) ); + R.makeCompressed( ); + + itkDebugMacro( << "Creating C matrix..." ); + C = R * B; + + itkDebugMacro( << "Creating A matrix..." ); + A.resize( N - nSeeds, N - nSeeds ); + A.setFromTriplets( At.begin( ), At.end( ) ); + A.makeCompressed( ); + } // end // Solve dirichlet problem _TSolver solver; + itkDebugMacro( << "Factorizing problem..." ); solver.compute( A ); if( solver.info( ) != Eigen::Success ) itkExceptionMacro( << "Error decomposing matrix." ); - _TMatrix x = solver.solve( R * B ); + itkDebugMacro( << "Solving problem..." ); + _TMatrix x = solver.solve( C ); if( solver.info( ) != Eigen::Success ) itkExceptionMacro( << "Error solving system." ); // Fill outputs + itkDebugMacro( << "Filling output..." ); this->_Output( x, St, invLabels ); } diff --git a/tests/image/RandomWalker/CompareRandomWalkers.cxx b/tests/image/RandomWalker/CompareRandomWalkers.cxx index 67afb48..b3177a9 100644 --- a/tests/image/RandomWalker/CompareRandomWalkers.cxx +++ b/tests/image/RandomWalker/CompareRandomWalkers.cxx @@ -95,6 +95,7 @@ void Original( filter->SetInput( input ); filter->SetInputLabels( labels ); filter->SetEdgeFunction( edge ); + filter->DebugOn( ); double t = MeasureTime( filter ); output_labels = filter->GetOutput( ); output_probabilities = filter->GetOutputProbabilities( ); @@ -183,8 +184,8 @@ int main( int argc, char* argv[] ) beta, epsilon ); - WriteImage( original_labels, "orig.png" ); - WriteImage( fp_labels, "fp.png" ); + WriteImage( original_labels, "orig.mhd" ); + WriteImage( fp_labels, "fp.mhd" ); return( 0 ); }