lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx

   1 #ifndef __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__
   2 #define __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__
   3
   4 #include <limits>
   5 #include <itkBinaryThresholdImageFilter.h>
   6 #include <itkConstNeighborhoodIterator.h>
   7
   8 // -------------------------------------------------------------------------
   9 template< class I >
  10 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
  11 AddThreshold( const TPixel& v )
  12 {
  13   this->m_Thresholds.insert( v );
  14   this->Modified( );
  15 }
  16
  17 // -------------------------------------------------------------------------
  18 template< class I >
  19 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
  20 AddThresholds( const TPixel& t0, const TPixel& t1, const unsigned int& s )
  21 {
  22   for( TPixel t = t0; t <= t1; t += TPixel( s ) )
  23     this->AddThreshold( t );
  24 }
  25
  26 // -------------------------------------------------------------------------
  27 template< class I >
  28 fpa::Image::RegionGrowWithMultipleThresholds< I >::
  29 RegionGrowWithMultipleThresholds( )
  30   : Superclass( ),
  31     m_InsideValue( TPixel( 1 ) ),
  32     m_OutsideValue( TPixel( 0 ) ),
  33     m_DifferenceThreshold( double( 3 ) ),
  34     m_TotalCount( 0 ),
  35     m_LastDiff( double( 0 ) ),
  36     m_StopForced( false ),
  37     m_StopThreshold( TPixel( 0 ) )
  38 {
  39 }
  40
  41 // -------------------------------------------------------------------------
  42 template< class I >
  43 fpa::Image::RegionGrowWithMultipleThresholds< I >::
  44 ~RegionGrowWithMultipleThresholds( )
  45 {
  46 }
  47
  48 // -------------------------------------------------------------------------
  49 template< class I >
  50 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
  51 _BeforeMainLoop( )
  52 {
  53   const I* img = this->GetInput( );
  54
  55   // Clear all states
  56   this->ClearMembershipFunctions( );
  57   this->m_Histogram.clear( );
  58   this->m_TotalCount = img->GetLargestPossibleRegion( ).GetNumberOfPixels( );
  59   this->m_LastDiff = double( -1 );
  60   this->m_StopForced = false;
  61
  62   // Initialize thresholding functors
  63   typename TThresholds::const_iterator tIt = this->m_Thresholds.begin( );
  64   TPixel min_thr = std::numeric_limits< TPixel >::min( );
  65   for( ; tIt != this->m_Thresholds.end( ); ++tIt )
  66   {
  67     typename TFunction::Pointer function = TFunction::New( );
  68     function->SetInputImage( img );
  69     function->SetLowerThreshold( min_thr );
  70     function->SetUpperThreshold( *tIt );
  71     this->AddMembershipFunction( function );
  72
  73   } // rof
  74
  75   // Correct seeds
  76   typename I::SizeType radius;
  77   radius.Fill( 3 );
  78   itk::ConstNeighborhoodIterator< I > it(
  79     radius, img, img->GetRequestedRegion( )
  80     );
  81   for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
  82   {
  83     _TNode seed = this->m_Seeds[ s ];
  84     it.SetLocation( seed.Vertex );
  85
  86     typename I::SizeType size = it.GetSize( );
  87     unsigned int nN = 1;
  88     for( unsigned int d = 0; d < I::ImageDimension; ++d )
  89       nN *= size[ d ];
  90     TPixel min_value = img->GetPixel( seed.Vertex );
  91     for( unsigned int i = 0; i < nN; ++i )
  92     {
  93       if( it.GetPixel( i ) < min_value )
  94       {
  95         seed.Vertex = it.GetIndex( i );
  96         seed.Parent = seed.Vertex;
  97         min_value = it.GetPixel( i );
  98
  99       } // fi
 100
 101     } // rof
 102     this->m_Seeds[ s ] = seed;
 103
 104   } // rof
 105
 106   // Continue all initializations
 107   this->Superclass::_BeforeMainLoop( );
 108 }
 109
 110 // -------------------------------------------------------------------------
 111 template< class I >
 112 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
 113 _AfterMainLoop( )
 114 {
 115   typedef itk::BinaryThresholdImageFilter< I, I > _TBinFilter;
 116
 117   // Binarize, inplace, the grown region
 118   if( this->m_Histogram.size( ) > 1 )
 119   {
 120     typename _TBinFilter::Pointer bin = _TBinFilter::New( );
 121     bin->SetInput( this->GetOutput( ) );
 122     bin->SetInsideValue( this->m_InsideValue );
 123     bin->SetOutsideValue( this->m_OutsideValue );
 124     bin->InPlaceOn( );
 125     if( this->m_StopForced )
 126       bin->SetUpperThreshold( this->m_StopThreshold );
 127     else
 128       bin->SetUpperThreshold( this->m_Histogram.rbegin( )->first );
 129     bin->GraftOutput( this->GetOutput( ) );
 130     bin->Update( );
 131     this->GraftOutput( bin->GetOutput( ) );
 132
 133   } // fi
 134
 135   this->Superclass::_AfterMainLoop( );
 136 }
 137
 138 // -------------------------------------------------------------------------
 139 template< class I >
 140 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
 141 _AfterLoop( )
 142 {
 143   if( this->m_ActualFunction != this->m_Functions.end( ) )
 144   {
 145     TFunction* f =
 146       dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
 147     this->m_Histogram[ f->GetUpperThreshold( ) ] = this->m_Marks.size( );
 148
 149     /* TODO: code to print the signal for demo purposes
 150        std::cout
 151        << f->GetUpperThreshold( ) << " "
 152        << this->m_Marks.size( )
 153        << std::endl;
 154     */
 155
 156     // Get the two last complete count, if any, and compute last
 157     // finite difference
 158     if( this->m_Histogram.size( ) > 1 )
 159     {
 160       typename THistogram::const_reverse_iterator hIt =
 161         this->m_Histogram.rbegin( );
 162       typename THistogram::const_reverse_iterator gIt = hIt;
 163       gIt++;
 164
 165       this->m_LastDiff  = double( hIt->second ) - double( gIt->second );
 166       this->m_LastDiff /= double( hIt->first )  - double( gIt->first );
 167       this->m_LastDiff *= m_DifferenceThreshold;
 168     }
 169     else
 170       this->m_LastDiff = double( -1 );
 171
 172   } // fi
 173   this->Superclass::_AfterLoop( );
 174 }
 175
 176 // -------------------------------------------------------------------------
 177 template< class I >
 178 bool fpa::Image::RegionGrowWithMultipleThresholds< I >::
 179 _UpdateResult( _TNode& n )
 180 {
 181   // Here, the output is changed to the threshold used to compute the
 182   // growing condition of the actual vertex n
 183   bool ret = this->Superclass::_UpdateResult( n );
 184   if( this->m_ActualFunction != this->m_Functions.end( ) )
 185   {
 186     TFunction* f =
 187       dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
 188     if( f != NULL )
 189       this->GetOutput( )->SetPixel( n.Vertex, f->GetUpperThreshold( ) );
 190
 191   } // fi
 192   return( ret );
 193 }
 194
 195 // -------------------------------------------------------------------------
 196 template< class I >
 197 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
 198 _Mark( const _TNode& n )
 199 {
 200   this->Superclass::_Mark( n );
 201
 202   // Check if the histogram's support is enough
 203   if( this->m_Histogram.size( ) < 2 )
 204     return;
 205   typename THistogram::const_reverse_iterator hIt =
 206     this->m_Histogram.rbegin( );
 207
 208   // Get the actual function
 209   TFunction* f =
 210     dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
 211   if( f == NULL )
 212     return;
 213
 214   // Previous pixel count
 215   unsigned long prev_aCount = hIt->second;
 216   double prev_p = double( prev_aCount ) / double( this->m_TotalCount );
 217
 218   // Actual pixel count
 219   unsigned long aCount = this->m_Marks.size( );
 220   double p = double( aCount ) / double( this->m_TotalCount );
 221
 222   // Loop over, at least, 0.1% of the total number of pixels before
 223   // performing stop analysis
 224   if( double( 1e-3 ) < p && double( 1e-3 ) < prev_p )
 225   {
 226     // Does the difference is worthy to be analyzed?
 227     if( aCount > hIt->second )
 228     {
 229       // Compute finite difference and compare it to the previous one
 230       double diff = double( aCount ) - double( hIt->second );
 231       diff /= double( f->GetUpperThreshold( ) ) - double( hIt->first );
 232       if( diff > this->m_LastDiff )
 233       {
 234         this->m_StopForced = true;
 235         this->m_StopThreshold = hIt->first;
 236
 237       } // fi
 238
 239     } // fi
 240
 241   } // fi
 242 }
 243
 244 // -------------------------------------------------------------------------
 245 template< class I >
 246 bool fpa::Image::RegionGrowWithMultipleThresholds< I >::
 247 _CheckStopCondition( )
 248 {
 249   return( this->m_StopForced );
 250 }
 251
 252 #endif // __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__
 253
 254 // eof - $RCSfile$