X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=lib%2Ffpa%2FBase%2FMori.hxx;h=17415e772bae2eba47f4fb777f2b4a6be1e1d7d4;hb=53d56cb3d8fe139843d5b2308f821cc05e7593e1;hp=548d755983b58756b471339c92f842444e10814a;hpb=89393f2267e42e921571c0184320d6c6382f34ab;p=FrontAlgorithms.git

diff --git a/lib/fpa/Base/Mori.hxx b/lib/fpa/Base/Mori.hxx
index 548d755..17415e7 100644
--- a/lib/fpa/Base/Mori.hxx
+++ b/lib/fpa/Base/Mori.hxx
@@ -64,11 +64,34 @@ SetThresholds(
     this->AddThreshold( thr );
 }
 
+// -------------------------------------------------------------------------
+template< class _TAlgorithm >
+unsigned long fpa::Base::Mori< _TAlgorithm >::
+GetNumberOfEvaluatedThresholds( ) const
+{
+  return( this->m_Signal.size( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TAlgorithm >
+typename fpa::Base::Mori< _TAlgorithm >::
+TInputValue fpa::Base::Mori< _TAlgorithm >::
+GetOptimumThreshold( ) const
+{
+  TInputValue thr = TInputValue( 0 );
+  if( this->m_Signal.size( ) > 1 )
+    thr = this->m_Signal[ this->m_Signal.size( ) - 2 ].first;
+  return( thr );
+}
+
 // -------------------------------------------------------------------------
 template< class _TAlgorithm >
 fpa::Base::Mori< _TAlgorithm >::
 Mori( )
   : Superclass( ),
+    m_SignalLag( 20 ),
+    m_SignalThreshold( 500 ),
+    m_SignalInfluence( 0.5 ),
     m_InsideValue( TOutputValue( 1 ) )
 {
   this->SetInitValue( TOutputValue( 0 ) );
@@ -98,22 +121,12 @@ _BeforeGenerateData( )
   this->m_Predicate->SetUpper( *( this->m_CurrentThreshold ) );
   this->m_Count = 0;
   this->m_Signal.clear( );
-}
-
-// -------------------------------------------------------------------------
-template< class _TAlgorithm >
-void fpa::Base::Mori< _TAlgorithm >::
-_AfterGenerateData( )
-{
-  // https://stackoverflow.com/questions/22583391/peak-signal-detection-in-realtime-timeseries-data
-  this->Superclass::_AfterGenerateData( );
-
-  typename TSignal::const_iterator sIt = this->m_Signal.begin( );
-  for( ; sIt != this->m_Signal.end( ); ++sIt )
-  {
-    std::cout << int( sIt->first ) << " " << int( sIt->second.first ) << " " << sIt->second.second << std::endl;
-
-  } // rof
+  this->m_FilteredSignal.clear( );
+  this->m_SignalAverages.clear( );
+  this->m_SignalDeviations.clear( );
+  this->m_SignalPeaks.clear( );
+  this->m_CurrentAverage = double( 0 );
+  this->m_CurrentVariance = double( 0 );
 }
 
 // -------------------------------------------------------------------------
@@ -123,14 +136,90 @@ _FinishOneLoop( )
 {
   if( this->m_Queues[ this->m_CurrentQueue ].size( ) == 0 )
   {
-    this->m_Signal[ this->m_Signal.size( ) + 1 ] =
-      TSignalData( *this->m_CurrentThreshold, this->m_Count );
+    this->m_Signal.push_back(
+      TSignalData( *this->m_CurrentThreshold, this->m_Count )
+      );
     this->m_CurrentThreshold++;
     this->m_CurrentQueue = ( this->m_CurrentQueue + 1 ) % 2;
     if( this->m_CurrentThreshold != this->m_Thresholds.end( ) )
     {
+      if( this->m_FilteredSignal.size( ) < this->m_SignalLag )
+      {
+        double v = double( this->m_Count );
+        double n = double( this->m_FilteredSignal.size( ) +  1 );
+        this->m_FilteredSignal.push_back( v );
+        this->m_SignalAverages.push_back( double( 0 ) );
+        this->m_SignalDeviations.push_back( double( 0 ) );
+        this->m_SignalPeaks.push_back( 0 );
+        if( n > double( 1 ) )
+          this->m_CurrentVariance =
+            ( ( ( n - 2.0 ) / ( n - 1.0 ) ) * this->m_CurrentVariance ) +
+            ( (
+                ( v - this->m_CurrentAverage ) *
+                ( v - this->m_CurrentAverage )
+              ) / n );
+        this->m_CurrentAverage += ( v - this->m_CurrentAverage ) / n;
+        if( this->m_FilteredSignal.size( ) == this->m_SignalLag )
+        {
+          this->m_SignalAverages.push_back( this->m_CurrentAverage );
+          this->m_SignalDeviations.push_back(
+            std::sqrt( this->m_CurrentVariance )
+            );
+
+        } // fi
+      }
+      else
+      {
+        unsigned long i = this->m_Signal.size( ) - 1;
+        double v = double( this->m_Count );
+        if(
+          ( std::fabs( v - this->m_SignalAverages[ i - 1 ] ) ) >
+          ( this->m_SignalThreshold * this->m_SignalDeviations[ i - 1 ] )
+          )
+        {
+          if( v > this->m_SignalAverages[ i - 1 ] )
+            this->m_SignalPeaks.push_back( 1 );
+          else
+            this->m_SignalPeaks.push_back( -1 );
+          this->m_FilteredSignal.push_back(
+            ( this->m_SignalInfluence * v ) +
+            (
+              ( 1.0 - this->m_SignalInfluence ) *
+              this->m_FilteredSignal[ i - 1 ]
+              )
+            );
+        }
+        else
+        {
+          this->m_SignalPeaks.push_back( 0 );
+          this->m_FilteredSignal.push_back( v );
+
+        } // fi
+
+        double avg = double( 0 );
+        double var = double( 0 );
+        unsigned long k = 0;
+        for( unsigned long j = i - this->m_SignalLag; j <= i; ++j, ++k )
+        {
+          double v = this->m_FilteredSignal[ j ];
+          double n = double( k + 1 );
+          if( k > 1 )
+            var =
+              ( ( ( n - 2.0 ) / ( n - 1.0 ) ) * var ) +
+              ( ( ( v - avg ) * ( v - avg ) ) / n );
+          avg += ( v - avg ) / n;
+
+        } // rof
+        this->m_SignalAverages.push_back( avg );
+        this->m_SignalDeviations.push_back( std::sqrt( var ) );
+
+      } // fi
       this->m_Predicate->SetUpper( *( this->m_CurrentThreshold ) );
       this->m_Count = 0;
+
+      // Peak detected? -> stop!
+      if( this->m_SignalPeaks.back( ) == 1 )
+        this->_QueueClear( );
     }
     else
       this->_QueueClear( );