lib/fpa/Base/Mori.hxx

   1 // =========================================================================
   2 // @author Leonardo Florez Valencia
   3 // @email florez-l@javeriana.edu.co
   4 // =========================================================================
   5
   6 #ifndef __fpa__Base__Mori__hxx__
   7 #define __fpa__Base__Mori__hxx__
   8
   9 // -------------------------------------------------------------------------
  10 template< class _TAlgorithm >
  11 itk::ModifiedTimeType fpa::Base::Mori< _TAlgorithm >::
  12 GetMTime( ) const
  13 {
  14   itk::ModifiedTimeType t = this->Superclass::GetMTime( );
  15   itk::ModifiedTimeType q = this->m_Predicate->GetMTime( );
  16   return( ( q < t )? q: t );
  17 }
  18
  19 // -------------------------------------------------------------------------
  20 template< class _TAlgorithm >
  21 typename fpa::Base::Mori< _TAlgorithm >::
  22 TOutputValue fpa::Base::Mori< _TAlgorithm >::
  23 GetOutsideValue( ) const
  24 {
  25   return( this->GetInitValue( ) );
  26 }
  27
  28 // -------------------------------------------------------------------------
  29 template< class _TAlgorithm >
  30 void fpa::Base::Mori< _TAlgorithm >::
  31 SetOutsideValue( const TOutputValue& v )
  32 {
  33   this->SetInitValue( v );
  34 }
  35
  36 // -------------------------------------------------------------------------
  37 template< class _TAlgorithm >
  38 void fpa::Base::Mori< _TAlgorithm >::
  39 ClearThresholds( )
  40 {
  41   this->m_Thresholds.clear( );
  42   this->Modified( );
  43 }
  44
  45 // -------------------------------------------------------------------------
  46 template< class _TAlgorithm >
  47 void fpa::Base::Mori< _TAlgorithm >::
  48 AddThreshold( const TInputValue& thr )
  49 {
  50   if( this->m_Thresholds.insert( thr ).second )
  51     this->Modified( );
  52 }
  53
  54 // -------------------------------------------------------------------------
  55 template< class _TAlgorithm >
  56 void fpa::Base::Mori< _TAlgorithm >::
  57 SetThresholds(
  58   const TInputValue& init,
  59   const TInputValue& end,
  60   const TInputValue& delta
  61   )
  62 {
  63   for( TInputValue thr = init; thr <= end; thr += delta )
  64     this->AddThreshold( thr );
  65 }
  66
  67 // -------------------------------------------------------------------------
  68 template< class _TAlgorithm >
  69 unsigned long fpa::Base::Mori< _TAlgorithm >::
  70 GetNumberOfEvaluatedThresholds( ) const
  71 {
  72   return( this->m_Signal.size( ) );
  73 }
  74
  75 // -------------------------------------------------------------------------
  76 template< class _TAlgorithm >
  77 typename fpa::Base::Mori< _TAlgorithm >::
  78 TInputValue fpa::Base::Mori< _TAlgorithm >::
  79 GetOptimumThreshold( ) const
  80 {
  81   TInputValue thr = TInputValue( 0 );
  82   if( this->m_Signal.size( ) > 1 )
  83     thr = this->m_Signal[ this->m_Signal.size( ) - 2 ].first;
  84   return( thr );
  85 }
  86
  87 // -------------------------------------------------------------------------
  88 template< class _TAlgorithm >
  89 fpa::Base::Mori< _TAlgorithm >::
  90 Mori( )
  91   : Superclass( ),
  92     m_SignalLag( 20 ),
  93     m_SignalThreshold( 500 ),
  94     m_SignalInfluence( 0.5 ),
  95     m_InsideValue( TOutputValue( 1 ) )
  96 {
  97   this->SetInitValue( TOutputValue( 0 ) );
  98   this->m_Predicate = TPredicate::New( );
  99   this->m_Predicate->StrictOff( );
 100 }
 101
 102 // -------------------------------------------------------------------------
 103 template< class _TAlgorithm >
 104 fpa::Base::Mori< _TAlgorithm >::
 105 ~Mori( )
 106 {
 107 }
 108
 109 // -------------------------------------------------------------------------
 110 template< class _TAlgorithm >
 111 void fpa::Base::Mori< _TAlgorithm >::
 112 _BeforeGenerateData( )
 113 {
 114   this->Superclass::_BeforeGenerateData( );
 115
 116   this->_QueueClear( );
 117   this->m_CurrentQueue = 0;
 118   this->m_CurrentThreshold = this->m_Thresholds.begin( );
 119   this->m_Predicate->SetLower( *( this->m_CurrentThreshold ) );
 120   this->m_CurrentThreshold++;
 121   this->m_Predicate->SetUpper( *( this->m_CurrentThreshold ) );
 122   this->m_Count = 0;
 123   this->m_Signal.clear( );
 124   this->m_FilteredSignal.clear( );
 125   this->m_SignalAverages.clear( );
 126   this->m_SignalDeviations.clear( );
 127   this->m_SignalPeaks.clear( );
 128   this->m_CurrentAverage = double( 0 );
 129   this->m_CurrentVariance = double( 0 );
 130 }
 131
 132 // -------------------------------------------------------------------------
 133 template< class _TAlgorithm >
 134 void fpa::Base::Mori< _TAlgorithm >::
 135 _FinishOneLoop( )
 136 {
 137   if( this->m_Queues[ this->m_CurrentQueue ].size( ) == 0 )
 138   {
 139     this->m_Signal.push_back(
 140       TSignalData( *this->m_CurrentThreshold, this->m_Count )
 141       );
 142     this->m_CurrentThreshold++;
 143     this->m_CurrentQueue = ( this->m_CurrentQueue + 1 ) % 2;
 144     if( this->m_CurrentThreshold != this->m_Thresholds.end( ) )
 145     {
 146       if( this->m_FilteredSignal.size( ) < this->m_SignalLag )
 147       {
 148         double v = double( this->m_Count );
 149         double n = double( this->m_FilteredSignal.size( ) +  1 );
 150         this->m_FilteredSignal.push_back( v );
 151         this->m_SignalAverages.push_back( double( 0 ) );
 152         this->m_SignalDeviations.push_back( double( 0 ) );
 153         this->m_SignalPeaks.push_back( 0 );
 154         if( n > double( 1 ) )
 155           this->m_CurrentVariance =
 156             ( ( ( n - 2.0 ) / ( n - 1.0 ) ) * this->m_CurrentVariance ) +
 157             ( (
 158                 ( v - this->m_CurrentAverage ) *
 159                 ( v - this->m_CurrentAverage )
 160               ) / n );
 161         this->m_CurrentAverage += ( v - this->m_CurrentAverage ) / n;
 162         if( this->m_FilteredSignal.size( ) == this->m_SignalLag )
 163         {
 164           this->m_SignalAverages.push_back( this->m_CurrentAverage );
 165           this->m_SignalDeviations.push_back(
 166             std::sqrt( this->m_CurrentVariance )
 167             );
 168
 169         } // fi
 170       }
 171       else
 172       {
 173         unsigned long i = this->m_Signal.size( ) - 1;
 174         double v = double( this->m_Count );
 175         if(
 176           ( std::fabs( v - this->m_SignalAverages[ i - 1 ] ) ) >
 177           ( this->m_SignalThreshold * this->m_SignalDeviations[ i - 1 ] )
 178           )
 179         {
 180           if( v > this->m_SignalAverages[ i - 1 ] )
 181             this->m_SignalPeaks.push_back( 1 );
 182           else
 183             this->m_SignalPeaks.push_back( -1 );
 184           this->m_FilteredSignal.push_back(
 185             ( this->m_SignalInfluence * v ) +
 186             (
 187               ( 1.0 - this->m_SignalInfluence ) *
 188               this->m_FilteredSignal[ i - 1 ]
 189               )
 190             );
 191         }
 192         else
 193         {
 194           this->m_SignalPeaks.push_back( 0 );
 195           this->m_FilteredSignal.push_back( v );
 196
 197         } // fi
 198
 199         double avg = double( 0 );
 200         double var = double( 0 );
 201         unsigned long k = 0;
 202         for( unsigned long j = i - this->m_SignalLag; j <= i; ++j, ++k )
 203         {
 204           double v = this->m_FilteredSignal[ j ];
 205           double n = double( k + 1 );
 206           if( k > 1 )
 207             var =
 208               ( ( ( n - 2.0 ) / ( n - 1.0 ) ) * var ) +
 209               ( ( ( v - avg ) * ( v - avg ) ) / n );
 210           avg += ( v - avg ) / n;
 211
 212         } // rof
 213         this->m_SignalAverages.push_back( avg );
 214         this->m_SignalDeviations.push_back( std::sqrt( var ) );
 215
 216       } // fi
 217       this->m_Predicate->SetUpper( *( this->m_CurrentThreshold ) );
 218       this->m_Count = 0;
 219
 220       // Peak detected? -> stop!
 221       if( this->m_SignalPeaks.back( ) == 1 )
 222         this->_QueueClear( );
 223     }
 224     else
 225       this->_QueueClear( );
 226
 227   } // fi
 228 }
 229
 230 // -------------------------------------------------------------------------
 231 template< class _TAlgorithm >
 232 void fpa::Base::Mori< _TAlgorithm >::
 233 _ComputeOutputValue( TNode& n )
 234 {
 235   // Do nothing!!!
 236 }
 237
 238 // -------------------------------------------------------------------------
 239 template< class _TAlgorithm >
 240 void fpa::Base::Mori< _TAlgorithm >::
 241 _UpdateOutputValue( TNode& n )
 242 {
 243   TInputValue value = this->_GetInputValue( n.Vertex );
 244   bool inside = this->m_Predicate->Evaluate( value );
 245   if( !inside )
 246   {
 247     n.Value = this->m_InitValue;
 248     n.FrontId++;
 249     this->m_Queues[ ( this->m_CurrentQueue + 1 ) % 2 ].push_back( n );
 250     n.FrontId = 0;
 251   }
 252   else
 253   {
 254     n.Value = this->m_InsideValue;
 255     this->m_Count++;
 256
 257   } // fi
 258   this->Superclass::_UpdateOutputValue( n );
 259 }
 260
 261 // -------------------------------------------------------------------------
 262 template< class _TAlgorithm >
 263 void fpa::Base::Mori< _TAlgorithm >::
 264 _QueueClear( )
 265 {
 266   this->m_Queues[ 0 ].clear( );
 267   this->m_Queues[ 1 ].clear( );
 268 }
 269
 270 // -------------------------------------------------------------------------
 271 template< class _TAlgorithm >
 272 typename fpa::Base::Mori< _TAlgorithm >::
 273 TNode fpa::Base::Mori< _TAlgorithm >::
 274 _QueuePop( )
 275 {
 276   TNode n = this->m_Queues[ this->m_CurrentQueue ].front( );
 277   this->m_Queues[ this->m_CurrentQueue ].pop_front( );
 278   return( n );
 279 }
 280
 281 // -------------------------------------------------------------------------
 282 template< class _TAlgorithm >
 283 void fpa::Base::Mori< _TAlgorithm >::
 284 _QueuePush( const TNode& node )
 285 {
 286   this->m_Queues[ this->m_CurrentQueue ].push_back( node );
 287 }
 288
 289 // -------------------------------------------------------------------------
 290 template< class _TAlgorithm >
 291 unsigned long fpa::Base::Mori< _TAlgorithm >::
 292 _QueueSize( ) const
 293 {
 294   return( this->m_Queues[ this->m_CurrentQueue ].size( ) );
 295 }
 296
 297 // -------------------------------------------------------------------------
 298 template< class _TAlgorithm >
 299 void fpa::Base::Mori< _TAlgorithm >::
 300 _PrepareSeeds( TNodes& nodes )
 301 {
 302   typename TNodes::iterator nIt = nodes.begin( );
 303   for( ; nIt != nodes.end( ); ++nIt )
 304     nIt->Value = this->m_InitValue;
 305 }
 306
 307 #endif // __fpa__Base__Mori__hxx__
 308
 309 // eof - $RCSfile$