lib/fpa/Image/RegionGrowWithMultipleThresholds.hxx

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