1 #ifndef __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__
2 #define __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__
5 #include <itkConstNeighborhoodIterator.h>
7 // -------------------------------------------------------------------------
9 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
10 AddThreshold( const TPixel& v )
12 this->m_Thresholds.insert( v );
16 // -------------------------------------------------------------------------
18 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
19 AddThresholds( const TPixel& t0, const TPixel& t1, const unsigned int& s )
21 for( TPixel t = t0; t <= t1; t += TPixel( s ) )
22 this->AddThreshold( t );
25 // -------------------------------------------------------------------------
27 fpa::Image::RegionGrowWithMultipleThresholds< I >::
28 RegionGrowWithMultipleThresholds( )
30 m_DifferenceThreshold( double( 3 ) ),
32 m_LastDiff( double( 0 ) ),
33 m_StopForced( false ),
34 m_StopThreshold( TPixel( 0 ) )
38 // -------------------------------------------------------------------------
40 fpa::Image::RegionGrowWithMultipleThresholds< I >::
41 ~RegionGrowWithMultipleThresholds( )
45 // -------------------------------------------------------------------------
47 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
50 const I* img = this->GetInput( );
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;
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 )
64 typename TFunction::Pointer function = TFunction::New( );
65 function->SetInputImage( img );
66 function->SetLowerThreshold( min_thr );
67 function->SetUpperThreshold( *tIt );
68 this->AddMembershipFunction( function );
73 typename I::SizeType radius;
75 itk::ConstNeighborhoodIterator< I > it(
76 radius, img, img->GetRequestedRegion( )
78 for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
80 _TNode seed = this->m_Seeds[ s ];
81 it.SetLocation( seed.Vertex );
83 typename I::SizeType size = it.GetSize( );
85 for( unsigned int d = 0; d < I::ImageDimension; ++d )
87 TPixel min_value = img->GetPixel( seed.Vertex );
88 for( unsigned int i = 0; i < nN; ++i )
90 if( it.GetPixel( i ) < min_value )
92 seed.Vertex = it.GetIndex( i );
93 min_value = it.GetPixel( i );
98 this->m_Seeds[ s ] = seed;
102 // Continue all initializations
103 this->Superclass::_BeforeMainLoop( );
106 // -------------------------------------------------------------------------
108 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
111 this->Superclass::_AfterMainLoop( );
114 // -------------------------------------------------------------------------
116 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
119 if( this->m_ActualFunction != this->m_Functions.end( ) )
122 dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
123 this->m_Histogram[ f->GetUpperThreshold( ) ] = this->m_Marks.size( );
125 /* TODO: code to print the signal for demo purposes
127 << f->GetUpperThreshold( ) << " "
128 << this->m_Marks.size( )
132 // Get the two last complete count, if any, and compute last
134 if( this->m_Histogram.size( ) > 1 )
136 typename THistogram::const_reverse_iterator hIt =
137 this->m_Histogram.rbegin( );
138 typename THistogram::const_reverse_iterator gIt = hIt;
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;
146 this->m_LastDiff = double( -1 );
149 this->Superclass::_AfterLoop( );
152 // -------------------------------------------------------------------------
154 bool fpa::Image::RegionGrowWithMultipleThresholds< I >::
155 _UpdateResult( _TNode& n )
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( ) )
163 dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
165 this->GetOutput( )->SetPixel( n.Vertex, f->GetUpperThreshold( ) );
171 // -------------------------------------------------------------------------
173 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
174 _Mark( const _TNode& n )
176 this->Superclass::_Mark( n );
178 // Check if the histogram's support is enough
179 if( this->m_Histogram.size( ) < 2 )
181 typename THistogram::const_reverse_iterator hIt =
182 this->m_Histogram.rbegin( );
184 // Get the actual function
186 dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
190 // Previous pixel count
191 unsigned long prev_aCount = hIt->second;
192 double prev_p = double( prev_aCount ) / double( this->m_TotalCount );
194 // Actual pixel count
195 unsigned long aCount = this->m_Marks.size( );
196 double p = double( aCount ) / double( this->m_TotalCount );
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 )
202 // Does the difference is worthy to be analyzed?
203 if( aCount > hIt->second )
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 )
210 this->m_StopForced = true;
211 this->m_StopThreshold = hIt->first;
220 // -------------------------------------------------------------------------
222 bool fpa::Image::RegionGrowWithMultipleThresholds< I >::
223 _CheckStopCondition( )
225 return( this->m_StopForced );
228 #endif // __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__