1 // =========================================================================
2 // @author Leonardo Florez Valencia
3 // @email florez-l@javeriana.edu.co
4 // =========================================================================
6 #ifndef __fpa__Base__Mori__hxx__
7 #define __fpa__Base__Mori__hxx__
9 // -------------------------------------------------------------------------
10 template< class _TAlgorithm >
11 itk::ModifiedTimeType fpa::Base::Mori< _TAlgorithm >::
14 itk::ModifiedTimeType t = this->Superclass::GetMTime( );
15 itk::ModifiedTimeType q = this->m_Predicate->GetMTime( );
16 return( ( q < t )? q: t );
19 // -------------------------------------------------------------------------
20 template< class _TAlgorithm >
21 typename fpa::Base::Mori< _TAlgorithm >::
22 TOutputValue fpa::Base::Mori< _TAlgorithm >::
23 GetOutsideValue( ) const
25 return( this->GetInitValue( ) );
28 // -------------------------------------------------------------------------
29 template< class _TAlgorithm >
30 void fpa::Base::Mori< _TAlgorithm >::
31 SetOutsideValue( const TOutputValue& v )
33 this->SetInitValue( v );
36 // -------------------------------------------------------------------------
37 template< class _TAlgorithm >
38 void fpa::Base::Mori< _TAlgorithm >::
41 this->m_Thresholds.clear( );
45 // -------------------------------------------------------------------------
46 template< class _TAlgorithm >
47 void fpa::Base::Mori< _TAlgorithm >::
48 AddThreshold( const TInputValue& thr )
50 if( this->m_Thresholds.insert( thr ).second )
54 // -------------------------------------------------------------------------
55 template< class _TAlgorithm >
56 void fpa::Base::Mori< _TAlgorithm >::
58 const TInputValue& init,
59 const TInputValue& end,
60 const TInputValue& delta
63 for( TInputValue thr = init; thr <= end; thr += delta )
64 this->AddThreshold( thr );
67 // -------------------------------------------------------------------------
68 template< class _TAlgorithm >
69 unsigned long fpa::Base::Mori< _TAlgorithm >::
70 GetNumberOfEvaluatedThresholds( ) const
72 return( this->m_Signal.size( ) );
75 // -------------------------------------------------------------------------
76 template< class _TAlgorithm >
77 typename fpa::Base::Mori< _TAlgorithm >::
78 TInputValue fpa::Base::Mori< _TAlgorithm >::
79 GetOptimumThreshold( ) const
81 TInputValue thr = TInputValue( 0 );
82 if( this->m_Signal.size( ) > 1 )
83 thr = this->m_Signal[ this->m_Signal.size( ) - 2 ].first;
87 // -------------------------------------------------------------------------
88 template< class _TAlgorithm >
89 fpa::Base::Mori< _TAlgorithm >::
93 m_SignalThreshold( 500 ),
94 m_SignalInfluence( 0.5 ),
95 m_InsideValue( TOutputValue( 1 ) )
97 this->SetInitValue( TOutputValue( 0 ) );
98 this->m_Predicate = TPredicate::New( );
99 this->m_Predicate->StrictOff( );
101 if( std::numeric_limits< TInputValue >::is_integer )
102 this->m_MinimumThreshold = std::numeric_limits< TInputValue >::min( );
104 this->m_MinimumThreshold = -std::numeric_limits< TInputValue >::max( );
107 // -------------------------------------------------------------------------
108 template< class _TAlgorithm >
109 fpa::Base::Mori< _TAlgorithm >::
114 // -------------------------------------------------------------------------
115 template< class _TAlgorithm >
116 void fpa::Base::Mori< _TAlgorithm >::
117 _BeforeGenerateData( )
119 this->Superclass::_BeforeGenerateData( );
121 this->_QueueClear( );
122 this->m_CurrentQueue = 0;
123 this->m_CurrentThreshold = this->m_Thresholds.begin( );
124 this->m_Predicate->SetLower( *( this->m_CurrentThreshold ) );
125 this->m_CurrentThreshold++;
126 this->m_Predicate->SetUpper( *( this->m_CurrentThreshold ) );
128 this->m_Signal.clear( );
129 this->m_FilteredSignal.clear( );
130 this->m_SignalAverages.clear( );
131 this->m_SignalDeviations.clear( );
132 this->m_SignalPeaks.clear( );
133 this->m_CurrentAverage = double( 0 );
134 this->m_CurrentVariance = double( 0 );
137 // -------------------------------------------------------------------------
138 template< class _TAlgorithm >
139 void fpa::Base::Mori< _TAlgorithm >::
142 if( this->m_Queues[ this->m_CurrentQueue ].size( ) == 0 )
144 this->m_Signal.push_back(
145 TSignalData( *this->m_CurrentThreshold, this->m_Count )
147 this->m_CurrentThreshold++;
148 this->m_CurrentQueue = ( this->m_CurrentQueue + 1 ) % 2;
149 if( this->m_CurrentThreshold != this->m_Thresholds.end( ) )
151 if( this->m_FilteredSignal.size( ) < this->m_SignalLag )
153 double v = double( this->m_Count );
154 double n = double( this->m_FilteredSignal.size( ) + 1 );
155 this->m_FilteredSignal.push_back( v );
156 this->m_SignalAverages.push_back( double( 0 ) );
157 this->m_SignalDeviations.push_back( double( 0 ) );
158 this->m_SignalPeaks.push_back( 0 );
159 if( n > double( 1 ) )
160 this->m_CurrentVariance =
161 ( ( ( n - 2.0 ) / ( n - 1.0 ) ) * this->m_CurrentVariance ) +
163 ( v - this->m_CurrentAverage ) *
164 ( v - this->m_CurrentAverage )
166 this->m_CurrentAverage += ( v - this->m_CurrentAverage ) / n;
167 if( this->m_FilteredSignal.size( ) == this->m_SignalLag )
169 this->m_SignalAverages.push_back( this->m_CurrentAverage );
170 this->m_SignalDeviations.push_back(
171 std::sqrt( this->m_CurrentVariance )
178 unsigned long i = this->m_Signal.size( ) - 1;
179 double v = double( this->m_Count );
181 ( std::fabs( v - this->m_SignalAverages[ i - 1 ] ) ) >
182 ( this->m_SignalThreshold * this->m_SignalDeviations[ i - 1 ] )
185 if( v > this->m_SignalAverages[ i - 1 ] )
186 this->m_SignalPeaks.push_back( 1 );
188 this->m_SignalPeaks.push_back( -1 );
189 this->m_FilteredSignal.push_back(
190 ( this->m_SignalInfluence * v ) +
192 ( 1.0 - this->m_SignalInfluence ) *
193 this->m_FilteredSignal[ i - 1 ]
199 this->m_SignalPeaks.push_back( 0 );
200 this->m_FilteredSignal.push_back( v );
204 double avg = double( 0 );
205 double var = double( 0 );
207 for( unsigned long j = i - this->m_SignalLag; j <= i; ++j, ++k )
209 double v = this->m_FilteredSignal[ j ];
210 double n = double( k + 1 );
213 ( ( ( n - 2.0 ) / ( n - 1.0 ) ) * var ) +
214 ( ( ( v - avg ) * ( v - avg ) ) / n );
215 avg += ( v - avg ) / n;
218 this->m_SignalAverages.push_back( avg );
219 this->m_SignalDeviations.push_back( std::sqrt( var ) );
222 this->m_Predicate->SetUpper( *( this->m_CurrentThreshold ) );
225 // Peak detected? -> stop!
226 if( this->m_SignalPeaks.back( ) == 1 && this->m_MinimumThreshold < *( this->m_CurrentThreshold ) )
227 this->_QueueClear( );
230 this->_QueueClear( );
235 // -------------------------------------------------------------------------
236 template< class _TAlgorithm >
237 void fpa::Base::Mori< _TAlgorithm >::
238 _ComputeOutputValue( TNode& n )
243 // -------------------------------------------------------------------------
244 template< class _TAlgorithm >
245 void fpa::Base::Mori< _TAlgorithm >::
246 _UpdateOutputValue( TNode& n )
248 TInputValue value = this->_GetInputValue( n.Vertex );
249 bool inside = this->m_Predicate->Evaluate( value );
252 n.Value = this->m_InitValue;
254 this->m_Queues[ ( this->m_CurrentQueue + 1 ) % 2 ].push_back( n );
259 n.Value = this->m_InsideValue;
263 this->Superclass::_UpdateOutputValue( n );
266 // -------------------------------------------------------------------------
267 template< class _TAlgorithm >
268 void fpa::Base::Mori< _TAlgorithm >::
271 this->m_Queues[ 0 ].clear( );
272 this->m_Queues[ 1 ].clear( );
275 // -------------------------------------------------------------------------
276 template< class _TAlgorithm >
277 typename fpa::Base::Mori< _TAlgorithm >::
278 TNode fpa::Base::Mori< _TAlgorithm >::
281 TNode n = this->m_Queues[ this->m_CurrentQueue ].front( );
282 this->m_Queues[ this->m_CurrentQueue ].pop_front( );
286 // -------------------------------------------------------------------------
287 template< class _TAlgorithm >
288 void fpa::Base::Mori< _TAlgorithm >::
289 _QueuePush( const TNode& node )
291 this->m_Queues[ this->m_CurrentQueue ].push_back( node );
294 // -------------------------------------------------------------------------
295 template< class _TAlgorithm >
296 unsigned long fpa::Base::Mori< _TAlgorithm >::
299 return( this->m_Queues[ this->m_CurrentQueue ].size( ) );
302 // -------------------------------------------------------------------------
303 template< class _TAlgorithm >
304 void fpa::Base::Mori< _TAlgorithm >::
305 _PrepareSeeds( TNodes& nodes )
307 typename TNodes::iterator nIt = nodes.begin( );
308 for( ; nIt != nodes.end( ); ++nIt )
309 nIt->Value = this->m_InitValue;
312 #endif // __fpa__Base__Mori__hxx__