1 #ifndef __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__
2 #define __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__
5 #include <itkBinaryThresholdImageFilter.h>
6 #include <itkConstNeighborhoodIterator.h>
8 // -------------------------------------------------------------------------
10 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
11 AddThreshold( const TPixel& v )
13 this->m_Thresholds.insert( v );
17 // -------------------------------------------------------------------------
19 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
20 AddThresholds( const TPixel& t0, const TPixel& t1, const unsigned int& s )
22 for( TPixel t = t0; t <= t1; t += TPixel( s ) )
23 this->AddThreshold( t );
26 // -------------------------------------------------------------------------
28 fpa::Image::RegionGrowWithMultipleThresholds< I >::
29 RegionGrowWithMultipleThresholds( )
31 m_InsideValue( TPixel( 1 ) ),
32 m_OutsideValue( TPixel( 0 ) ),
33 m_DifferenceThreshold( double( 3 ) ),
35 m_LastDiff( double( 0 ) ),
36 m_StopForced( false ),
37 m_StopThreshold( TPixel( 0 ) )
41 // -------------------------------------------------------------------------
43 fpa::Image::RegionGrowWithMultipleThresholds< I >::
44 ~RegionGrowWithMultipleThresholds( )
48 // -------------------------------------------------------------------------
50 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
53 const I* img = this->GetInput( );
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;
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 )
67 typename TFunction::Pointer function = TFunction::New( );
68 function->SetInputImage( img );
69 function->SetLowerThreshold( min_thr );
70 function->SetUpperThreshold( *tIt );
71 this->AddMembershipFunction( function );
76 typename I::SizeType radius;
78 itk::ConstNeighborhoodIterator< I > it(
79 radius, img, img->GetRequestedRegion( )
81 for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
83 _TNode seed = this->m_Seeds[ s ];
84 it.SetLocation( seed.Vertex );
86 typename I::SizeType size = it.GetSize( );
88 for( unsigned int d = 0; d < I::ImageDimension; ++d )
90 TPixel min_value = img->GetPixel( seed.Vertex );
91 for( unsigned int i = 0; i < nN; ++i )
93 if( it.GetPixel( i ) < min_value )
95 seed.Vertex = it.GetIndex( i );
96 seed.Parent = seed.Vertex;
97 min_value = it.GetPixel( i );
102 this->m_Seeds[ s ] = seed;
106 // Continue all initializations
107 this->Superclass::_BeforeMainLoop( );
110 // -------------------------------------------------------------------------
112 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
115 typedef itk::BinaryThresholdImageFilter< I, I > _TBinFilter;
117 // Binarize, inplace, the grown region
118 if( this->m_Histogram.size( ) > 1 )
120 typename _TBinFilter::Pointer bin = _TBinFilter::New( );
121 bin->SetInput( this->GetOutput( ) );
122 bin->SetInsideValue( this->m_InsideValue );
123 bin->SetOutsideValue( this->m_OutsideValue );
125 if( this->m_StopForced )
126 bin->SetUpperThreshold( this->m_StopThreshold );
128 bin->SetUpperThreshold( this->m_Histogram.rbegin( )->first );
129 bin->GraftOutput( this->GetOutput( ) );
131 this->GraftOutput( bin->GetOutput( ) );
135 this->Superclass::_AfterMainLoop( );
138 // -------------------------------------------------------------------------
140 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
143 if( this->m_ActualFunction != this->m_Functions.end( ) )
146 dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
147 this->m_Histogram[ f->GetUpperThreshold( ) ] = this->m_Marks.size( );
149 /* TODO: code to print the signal for demo purposes
151 << f->GetUpperThreshold( ) << " "
152 << this->m_Marks.size( )
156 // Get the two last complete count, if any, and compute last
158 if( this->m_Histogram.size( ) > 1 )
160 typename THistogram::const_reverse_iterator hIt =
161 this->m_Histogram.rbegin( );
162 typename THistogram::const_reverse_iterator gIt = hIt;
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;
170 this->m_LastDiff = double( -1 );
173 this->Superclass::_AfterLoop( );
176 // -------------------------------------------------------------------------
178 bool fpa::Image::RegionGrowWithMultipleThresholds< I >::
179 _UpdateResult( _TNode& n )
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( ) )
187 dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
189 this->GetOutput( )->SetPixel( n.Vertex, f->GetUpperThreshold( ) );
195 // -------------------------------------------------------------------------
197 void fpa::Image::RegionGrowWithMultipleThresholds< I >::
198 _Mark( const _TNode& n )
200 this->Superclass::_Mark( n );
202 // Check if the histogram's support is enough
203 if( this->m_Histogram.size( ) < 2 )
205 typename THistogram::const_reverse_iterator hIt =
206 this->m_Histogram.rbegin( );
208 // Get the actual function
210 dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
214 // Previous pixel count
215 unsigned long prev_aCount = hIt->second;
216 double prev_p = double( prev_aCount ) / double( this->m_TotalCount );
218 // Actual pixel count
219 unsigned long aCount = this->m_Marks.size( );
220 double p = double( aCount ) / double( this->m_TotalCount );
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 )
226 // Does the difference is worthy to be analyzed?
227 if( aCount > hIt->second )
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 )
234 /* TODO: comment this for demo purposes
236 this->m_StopForced = true;
237 this->m_StopThreshold = hIt->first;
246 // -------------------------------------------------------------------------
248 bool fpa::Image::RegionGrowWithMultipleThresholds< I >::
249 _CheckStopCondition( )
251 return( this->m_StopForced );
254 #endif // __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__