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   if( std::numeric_limits< TInputValue >::is_integer )
 102     this->m_MinimumThreshold = std::numeric_limits< TInputValue >::min( );
 103   else
 104     this->m_MinimumThreshold = -std::numeric_limits< TInputValue >::max( );
 105 }
 106
 107 // -------------------------------------------------------------------------
 108 template< class _TAlgorithm >
 109 fpa::Base::Mori< _TAlgorithm >::
 110 ~Mori( )
 111 {
 112 }
 113
 114 // -------------------------------------------------------------------------
 115 template< class _TAlgorithm >
 116 void fpa::Base::Mori< _TAlgorithm >::
 117 _BeforeGenerateData( )
 118 {
 119   this->Superclass::_BeforeGenerateData( );
 120
 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 ) );
 127   this->m_Count = 0;
 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 );
 135 }
 136
 137 // -------------------------------------------------------------------------
 138 template< class _TAlgorithm >
 139 void fpa::Base::Mori< _TAlgorithm >::
 140 _FinishOneLoop( )
 141 {
 142   if( this->m_Queues[ this->m_CurrentQueue ].size( ) == 0 )
 143   {
 144     this->m_Signal.push_back(
 145       TSignalData( *this->m_CurrentThreshold, this->m_Count )
 146       );
 147     this->m_CurrentThreshold++;
 148     this->m_CurrentQueue = ( this->m_CurrentQueue + 1 ) % 2;
 149     if( this->m_CurrentThreshold != this->m_Thresholds.end( ) )
 150     {
 151       if( this->m_FilteredSignal.size( ) < this->m_SignalLag )
 152       {
 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 ) +
 162             ( (
 163                 ( v - this->m_CurrentAverage ) *
 164                 ( v - this->m_CurrentAverage )
 165               ) / n );
 166         this->m_CurrentAverage += ( v - this->m_CurrentAverage ) / n;
 167         if( this->m_FilteredSignal.size( ) == this->m_SignalLag )
 168         {
 169           this->m_SignalAverages.push_back( this->m_CurrentAverage );
 170           this->m_SignalDeviations.push_back(
 171             std::sqrt( this->m_CurrentVariance )
 172             );
 173
 174         } // fi
 175       }
 176       else
 177       {
 178         unsigned long i = this->m_Signal.size( ) - 1;
 179         double v = double( this->m_Count );
 180         if(
 181           ( std::fabs( v - this->m_SignalAverages[ i - 1 ] ) ) >
 182           ( this->m_SignalThreshold * this->m_SignalDeviations[ i - 1 ] )
 183           )
 184         {
 185           if( v > this->m_SignalAverages[ i - 1 ] )
 186             this->m_SignalPeaks.push_back( 1 );
 187           else
 188             this->m_SignalPeaks.push_back( -1 );
 189           this->m_FilteredSignal.push_back(
 190             ( this->m_SignalInfluence * v ) +
 191             (
 192               ( 1.0 - this->m_SignalInfluence ) *
 193               this->m_FilteredSignal[ i - 1 ]
 194               )
 195             );
 196         }
 197         else
 198         {
 199           this->m_SignalPeaks.push_back( 0 );
 200           this->m_FilteredSignal.push_back( v );
 201
 202         } // fi
 203
 204         double avg = double( 0 );
 205         double var = double( 0 );
 206         unsigned long k = 0;
 207         for( unsigned long j = i - this->m_SignalLag; j <= i; ++j, ++k )
 208         {
 209           double v = this->m_FilteredSignal[ j ];
 210           double n = double( k + 1 );
 211           if( k > 1 )
 212             var =
 213               ( ( ( n - 2.0 ) / ( n - 1.0 ) ) * var ) +
 214               ( ( ( v - avg ) * ( v - avg ) ) / n );
 215           avg += ( v - avg ) / n;
 216
 217         } // rof
 218         this->m_SignalAverages.push_back( avg );
 219         this->m_SignalDeviations.push_back( std::sqrt( var ) );
 220
 221       } // fi
 222       this->m_Predicate->SetUpper( *( this->m_CurrentThreshold ) );
 223       this->m_Count = 0;
 224
 225       // Peak detected? -> stop!
 226       if( this->m_SignalPeaks.back( ) == 1 && this->m_MinimumThreshold < *( this->m_CurrentThreshold ) )
 227         this->_QueueClear( );
 228     }
 229     else
 230       this->_QueueClear( );
 231
 232   } // fi
 233 }
 234
 235 // -------------------------------------------------------------------------
 236 template< class _TAlgorithm >
 237 void fpa::Base::Mori< _TAlgorithm >::
 238 _ComputeOutputValue( TNode& n )
 239 {
 240   // Do nothing!!!
 241 }
 242
 243 // -------------------------------------------------------------------------
 244 template< class _TAlgorithm >
 245 void fpa::Base::Mori< _TAlgorithm >::
 246 _UpdateOutputValue( TNode& n )
 247 {
 248   TInputValue value = this->_GetInputValue( n.Vertex );
 249   bool inside = this->m_Predicate->Evaluate( value );
 250   if( !inside )
 251   {
 252     n.Value = this->m_InitValue;
 253     n.FrontId++;
 254     this->m_Queues[ ( this->m_CurrentQueue + 1 ) % 2 ].push_back( n );
 255     n.FrontId = 0;
 256   }
 257   else
 258   {
 259     n.Value = this->m_InsideValue;
 260     this->m_Count++;
 261
 262   } // fi
 263   this->Superclass::_UpdateOutputValue( n );
 264 }
 265
 266 // -------------------------------------------------------------------------
 267 template< class _TAlgorithm >
 268 void fpa::Base::Mori< _TAlgorithm >::
 269 _QueueClear( )
 270 {
 271   this->m_Queues[ 0 ].clear( );
 272   this->m_Queues[ 1 ].clear( );
 273 }
 274
 275 // -------------------------------------------------------------------------
 276 template< class _TAlgorithm >
 277 typename fpa::Base::Mori< _TAlgorithm >::
 278 TNode fpa::Base::Mori< _TAlgorithm >::
 279 _QueuePop( )
 280 {
 281   TNode n = this->m_Queues[ this->m_CurrentQueue ].front( );
 282   this->m_Queues[ this->m_CurrentQueue ].pop_front( );
 283   return( n );
 284 }
 285
 286 // -------------------------------------------------------------------------
 287 template< class _TAlgorithm >
 288 void fpa::Base::Mori< _TAlgorithm >::
 289 _QueuePush( const TNode& node )
 290 {
 291   this->m_Queues[ this->m_CurrentQueue ].push_back( node );
 292 }
 293
 294 // -------------------------------------------------------------------------
 295 template< class _TAlgorithm >
 296 unsigned long fpa::Base::Mori< _TAlgorithm >::
 297 _QueueSize( ) const
 298 {
 299   return( this->m_Queues[ this->m_CurrentQueue ].size( ) );
 300 }
 301
 302 // -------------------------------------------------------------------------
 303 template< class _TAlgorithm >
 304 void fpa::Base::Mori< _TAlgorithm >::
 305 _PrepareSeeds( TNodes& nodes )
 306 {
 307   typename TNodes::iterator nIt = nodes.begin( );
 308   for( ; nIt != nodes.end( ); ++nIt )
 309     nIt->Value = this->m_InitValue;
 310 }
 311
 312 #endif // __fpa__Base__Mori__hxx__
 313
 314 // eof - $RCSfile$