X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=lib%2Ffpa%2FImage%2FRegionGrowWithMultipleThresholds.hxx;h=2e8c995c20d26ece1a8992b09c2c703f0d604a9c;hb=c7563ffe89c76b52736c3bd6bb2d5e6fac009cdd;hp=0969b0cd59e082ccdf8d33c30778e75d0c171697;hpb=22051707794a9153acbfc489ec9a933471f13406;p=FrontAlgorithms.git diff --git a/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx b/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx index 0969b0c..2e8c995 100644 --- a/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx +++ b/lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx @@ -2,43 +2,16 @@ #define __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__ #include -#include +#include +#include // ------------------------------------------------------------------------- template< class I > void fpa::Image::RegionGrowWithMultipleThresholds< I >:: AddThreshold( const TPixel& v ) { - typedef - fpa::Image::Functors::RegionGrowThresholdFunction< I > - TFunction; - typename TFunction::Pointer function = TFunction::New( ); - - if( this->GetNumberOfMembershipFunctions( ) > 0 ) - { - } - else - function->SetLowerThreshold( std::numeric_limits< TPixel >::min( ) ); - std::cout - << typeid( TPixel ).name( ) << " <<<----->>> " - << function->GetLowerThreshold( ) - << std::endl; - std::exit( 1 ); - function->SetUpperThreshold( v ); - this->AddMembershipFunction( function ); - - /* TODO - this->m_Histogram[ v ] = 0; - - TFunction* function = - dynamic_cast< TFunction* >( this->GetMembershipFunction( ) ); - if( function != NULL ) - { - function->SetLowerThreshold( this->m_Histogram.begin( )->first ); - - } // fi - this->Modified( ); - */ + this->m_Thresholds.insert( v ); + this->Modified( ); } // ------------------------------------------------------------------------- @@ -46,9 +19,8 @@ template< class I > void fpa::Image::RegionGrowWithMultipleThresholds< I >:: AddThresholds( const TPixel& t0, const TPixel& t1, const unsigned int& s ) { - double o = ( double( t1 ) - double( t0 ) ) / ( double( s ) - double( 1 ) ); - for( unsigned int i = 0; i < s; i++ ) - this->AddThreshold( TPixel( ( double( i ) * o ) + double( t0 ) ) ); + for( TPixel t = t0; t <= t1; t += TPixel( s ) ) + this->AddThreshold( t ); } // ------------------------------------------------------------------------- @@ -56,15 +28,14 @@ template< class I > fpa::Image::RegionGrowWithMultipleThresholds< I >:: RegionGrowWithMultipleThresholds( ) : Superclass( ), - m_DerivativeThreshold( double( 3 ) ) + m_InsideValue( TPixel( 1 ) ), + m_OutsideValue( TPixel( 0 ) ), + m_DifferenceThreshold( double( 3 ) ), + m_TotalCount( 0 ), + m_LastDiff( double( 0 ) ), + m_StopForced( false ), + m_StopThreshold( TPixel( 0 ) ) { - /* TODO - typedef - fpa::Image::Functors::RegionGrowThresholdFunction< I > - TFunction; - typename TFunction::Pointer function = TFunction::New( ); - this->SetMembershipFunction( function ); - */ } // ------------------------------------------------------------------------- @@ -76,53 +47,92 @@ fpa::Image::RegionGrowWithMultipleThresholds< I >:: // ------------------------------------------------------------------------- template< class I > -bool fpa::Image::RegionGrowWithMultipleThresholds< I >:: -_UpdateResult( _TNode& n ) +void fpa::Image::RegionGrowWithMultipleThresholds< I >:: +_BeforeMainLoop( ) { - bool ret = this->Superclass::_UpdateResult( n ); + const I* img = this->GetInput( ); - /* TODO - if( ret ) - { - TPixel v = TPixel( this->_Cost( n.Vertex, n.Vertex ) ); - - typename THistogram::reverse_iterator hIt = this->m_Histogram.rbegin( ); - while( hIt != this->m_Histogram.rend( ) ) - { - if( v <= hIt->first ) - { - hIt->second += 1; - hIt++; - } - else - hIt = this->m_Histogram.rend( ); - - } // elihw - this->GetOutput( )->SetPixel( n.Vertex, v ); - - } // fi - */ - return( ret ); + // Clear all states + this->ClearMembershipFunctions( ); + this->m_Histogram.clear( ); + this->m_TotalCount = img->GetLargestPossibleRegion( ).GetNumberOfPixels( ); + this->m_LastDiff = double( -1 ); + this->m_StopForced = false; + + // Initialize thresholding functors + typename TThresholds::const_iterator tIt = this->m_Thresholds.begin( ); + TPixel min_thr = std::numeric_limits< TPixel >::min( ); + for( ; tIt != this->m_Thresholds.end( ); ++tIt ) + { + typename TFunction::Pointer function = TFunction::New( ); + function->SetInputImage( img ); + function->SetLowerThreshold( min_thr ); + function->SetUpperThreshold( *tIt ); + this->AddMembershipFunction( function ); + + } // rof + + // Correct seeds + typename I::SizeType radius; + radius.Fill( 3 ); + itk::ConstNeighborhoodIterator< I > it( + radius, img, img->GetRequestedRegion( ) + ); + for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s ) + { + _TNode seed = this->m_Seeds[ s ]; + it.SetLocation( seed.Vertex ); + + typename I::SizeType size = it.GetSize( ); + unsigned int nN = 1; + for( unsigned int d = 0; d < I::ImageDimension; ++d ) + nN *= size[ d ]; + TPixel min_value = img->GetPixel( seed.Vertex ); + for( unsigned int i = 0; i < nN; ++i ) + { + if( it.GetPixel( i ) < min_value ) + { + seed.Vertex = it.GetIndex( i ); + seed.Parent = seed.Vertex; + min_value = it.GetPixel( i ); + + } // fi + + } // rof + this->m_Seeds[ s ] = seed; + + } // rof + + // Continue all initializations + this->Superclass::_BeforeMainLoop( ); } // ------------------------------------------------------------------------- template< class I > void fpa::Image::RegionGrowWithMultipleThresholds< I >:: -_BeforeLoop( ) +_AfterMainLoop( ) { - std::cout << "**1" << std::endl; - const I* img = this->GetInput( ); - std::cout << "**2" << std::endl; - typename TFunctions::iterator fIt = this->m_Functions.begin( ); - for( ; fIt != this->m_Functions.end( ); ++fIt ) - { - TMembershipFunction* f = - dynamic_cast< TMembershipFunction* >( fIt->GetPointer( ) ); - if( f != NULL ) - f->SetInputImage( this->GetInput( ) ); + typedef itk::BinaryThresholdImageFilter< I, I > _TBinFilter; - } // rof - this->Superclass::_BeforeLoop( ); + // Binarize, inplace, the grown region + if( this->m_Histogram.size( ) > 1 ) + { + typename _TBinFilter::Pointer bin = _TBinFilter::New( ); + bin->SetInput( this->GetOutput( ) ); + bin->SetInsideValue( this->m_InsideValue ); + bin->SetOutsideValue( this->m_OutsideValue ); + bin->InPlaceOn( ); + if( this->m_StopForced ) + bin->SetUpperThreshold( this->m_StopThreshold ); + else + bin->SetUpperThreshold( this->m_Histogram.rbegin( )->first ); + bin->GraftOutput( this->GetOutput( ) ); + bin->Update( ); + this->GraftOutput( bin->GetOutput( ) ); + + } // fi + + this->Superclass::_AfterMainLoop( ); } // ------------------------------------------------------------------------- @@ -130,43 +140,113 @@ template< class I > void fpa::Image::RegionGrowWithMultipleThresholds< I >:: _AfterLoop( ) { - /* - typename THistogram::iterator prevIt = this->m_Histogram.begin( ); - typename THistogram::iterator currIt = prevIt; currIt++; - typename THistogram::iterator nextIt = currIt; nextIt++; - double prev_d1 = double( 0 ); - for( ; nextIt != this->m_Histogram.end( ); ++prevIt, ++currIt, ++nextIt ) + if( this->m_ActualFunction != this->m_Functions.end( ) ) + { + TFunction* f = + dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) ); + this->m_Histogram[ f->GetUpperThreshold( ) ] = this->m_Marks.size( ); + + /* TODO: code to print the signal for demo purposes + std::cout + << f->GetUpperThreshold( ) << " " + << this->m_Marks.size( ) + << std::endl; + */ + + // Get the two last complete count, if any, and compute last + // finite difference + if( this->m_Histogram.size( ) > 1 ) { - double d1 = double( nextIt->second ) - double( prevIt->second ); - d1 /= double( nextIt->first ) - double( prevIt->first ); - - std::cout - << currIt->first << " " - << currIt->second << " " - << d1 << " " - << ( d1 - prev_d1 ) << std::endl; + typename THistogram::const_reverse_iterator hIt = + this->m_Histogram.rbegin( ); + typename THistogram::const_reverse_iterator gIt = hIt; + gIt++; + + this->m_LastDiff = double( hIt->second ) - double( gIt->second ); + this->m_LastDiff /= double( hIt->first ) - double( gIt->first ); + this->m_LastDiff *= m_DifferenceThreshold; + } + else + this->m_LastDiff = double( -1 ); + + } // fi + this->Superclass::_AfterLoop( ); +} - prev_d1 = d1; +// ------------------------------------------------------------------------- +template< class I > +bool fpa::Image::RegionGrowWithMultipleThresholds< I >:: +_UpdateResult( _TNode& n ) +{ + // Here, the output is changed to the threshold used to compute the + // growing condition of the actual vertex n + bool ret = this->Superclass::_UpdateResult( n ); + if( this->m_ActualFunction != this->m_Functions.end( ) ) + { + TFunction* f = + dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) ); + if( f != NULL ) + this->GetOutput( )->SetPixel( n.Vertex, f->GetUpperThreshold( ) ); - } // rof - */ + } // fi + return( ret ); +} - /* - double prev = double( 0 ); - while( hIt != this->m_Histogram.end( ) ) +// ------------------------------------------------------------------------- +template< class I > +void fpa::Image::RegionGrowWithMultipleThresholds< I >:: +_Mark( const _TNode& n ) +{ + this->Superclass::_Mark( n ); + + // Check if the histogram's support is enough + if( this->m_Histogram.size( ) < 2 ) + return; + typename THistogram::const_reverse_iterator hIt = + this->m_Histogram.rbegin( ); + + // Get the actual function + TFunction* f = + dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) ); + if( f == NULL ) + return; + + // Previous pixel count + unsigned long prev_aCount = hIt->second; + double prev_p = double( prev_aCount ) / double( this->m_TotalCount ); + + // Actual pixel count + unsigned long aCount = this->m_Marks.size( ); + double p = double( aCount ) / double( this->m_TotalCount ); + + // Loop over, at least, 0.1% of the total number of pixels before + // performing stop analysis + if( double( 1e-3 ) < p && double( 1e-3 ) < prev_p ) + { + // Does the difference is worthy to be analyzed? + if( aCount > hIt->second ) { - double curr = double( hIt->second ); - double deri = curr - prev; - - if( hIt != this->m_Histogram.begin( ) ) - std::cout << hIt->first << " " << curr << " " << deri << std::endl; - - prev = curr; - hIt++; + // Compute finite difference and compare it to the previous one + double diff = double( aCount ) - double( hIt->second ); + diff /= double( f->GetUpperThreshold( ) ) - double( hIt->first ); + if( diff > this->m_LastDiff ) + { + this->m_StopForced = true; + this->m_StopThreshold = hIt->first; + + } // fi + + } // fi + + } // fi +} - } // rof - */ - this->Superclass::_AfterLoop( ); +// ------------------------------------------------------------------------- +template< class I > +bool fpa::Image::RegionGrowWithMultipleThresholds< I >:: +_CheckStopCondition( ) +{ + return( this->m_StopForced ); } #endif // __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__