X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=lib%2Ffpa%2FBase%2FMori.hxx;h=3d1a94dfa63e996cbf36aa580a8ae15bdd964660;hb=38eec70e80acde7adfc59e3eae666c848a437abd;hp=7a1df9a854107f0f30e880edc926a5dc9345294d;hpb=5c78aecb0f5a207ff020e24f99d1e9bd8c388ad1;p=FrontAlgorithms.git

diff --git a/lib/fpa/Base/Mori.hxx b/lib/fpa/Base/Mori.hxx
index 7a1df9a..3d1a94d 100644
--- a/lib/fpa/Base/Mori.hxx
+++ b/lib/fpa/Base/Mori.hxx
@@ -33,6 +33,57 @@ SetOutsideValue( const TOutputValue& v )
   this->SetInitValue( v );
 }
 
+// -------------------------------------------------------------------------
+template< class _TAlgorithm >
+unsigned long fpa::Base::Mori< _TAlgorithm >::
+GetSignalKernelSize( ) const
+{
+  return( this->m_PeakDetector.GetKernelSize( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TAlgorithm >
+double fpa::Base::Mori< _TAlgorithm >::
+GetSignalThreshold( ) const
+{
+  return( this->m_PeakDetector.GetThreshold( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TAlgorithm >
+double fpa::Base::Mori< _TAlgorithm >::
+GetSignalInfluence( ) const
+{
+  return( this->m_PeakDetector.GetInfluence( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TAlgorithm >
+void fpa::Base::Mori< _TAlgorithm >::
+SetSignalKernelSize( unsigned long k )
+{
+  this->m_PeakDetector.SetKernelSize( k );
+  this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TAlgorithm >
+void fpa::Base::Mori< _TAlgorithm >::
+SetSignalThreshold( double t )
+{
+  this->m_PeakDetector.SetThreshold( t );
+  this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class _TAlgorithm >
+void fpa::Base::Mori< _TAlgorithm >::
+SetSignalInfluence( double i )
+{
+  this->m_PeakDetector.SetInfluence( i );
+  this->Modified( );
+}
+
 // -------------------------------------------------------------------------
 template< class _TAlgorithm >
 void fpa::Base::Mori< _TAlgorithm >::
@@ -64,12 +115,21 @@ SetThresholds(
     this->AddThreshold( thr );
 }
 
+// -------------------------------------------------------------------------
+template< class _TAlgorithm >
+const typename fpa::Base::Mori< _TAlgorithm >::
+TThresholds& fpa::Base::Mori< _TAlgorithm >::
+GetThresholds( ) const
+{
+  return( this->m_Thresholds );
+}
+
 // -------------------------------------------------------------------------
 template< class _TAlgorithm >
 unsigned long fpa::Base::Mori< _TAlgorithm >::
 GetNumberOfEvaluatedThresholds( ) const
 {
-  return( this->m_Signal.size( ) );
+  return( this->m_PeakDetector.GetNumberOfSamples( ) );
 }
 
 // -------------------------------------------------------------------------
@@ -79,29 +139,43 @@ 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;
+  unsigned long n = this->m_PeakDetector.GetNumberOfSamples( );
+  if( n > 1 )
+    thr = TInputValue( this->m_PeakDetector.GetXValues( )[ n - 2 ] );
   return( thr );
 }
 
+// -------------------------------------------------------------------------
+template< class _TAlgorithm >
+void fpa::Base::Mori< _TAlgorithm >::
+GetSignalValues( unsigned long i, double& x, double& y, TPeak& p ) const
+{
+  if( i < this->m_PeakDetector.GetNumberOfSamples( ) )
+  {
+    x = this->m_PeakDetector.GetXValues( )[ i ];
+    y = this->m_PeakDetector.GetYValues( )[ i ];
+    p = this->m_PeakDetector.GetPeaks( )[ i ];
+
+  } // fi
+}
+
 // -------------------------------------------------------------------------
 template< class _TAlgorithm >
 fpa::Base::Mori< _TAlgorithm >::
 Mori( )
   : Superclass( ),
-    m_SignalLag( 20 ),
-    m_SignalThreshold( 500 ),
-    m_SignalInfluence( 0.5 ),
-    m_InsideValue( TOutputValue( 1 ) )
+     m_InsideValue( TOutputValue( 1 ) )
 {
   this->SetInitValue( TOutputValue( 0 ) );
   this->m_Predicate = TPredicate::New( );
   this->m_Predicate->StrictOff( );
-
   if( std::numeric_limits< TInputValue >::is_integer )
     this->m_MinimumThreshold = std::numeric_limits< TInputValue >::min( );
   else
     this->m_MinimumThreshold = -std::numeric_limits< TInputValue >::max( );
+  this->m_PeakDetector.SetKernelSize( 20 );
+  this->m_PeakDetector.SetThreshold( 500 );
+  this->m_PeakDetector.SetInfluence( 0.5 );
 }
 
 // -------------------------------------------------------------------------
@@ -118,20 +192,19 @@ _BeforeGenerateData( )
 {
   this->Superclass::_BeforeGenerateData( );
 
+  // Prepare queues
   this->_QueueClear( );
-  this->m_CurrentQueue = 0;
-  this->m_CurrentThreshold = this->m_Thresholds.begin( );
-  this->m_Predicate->SetLower( *( this->m_CurrentThreshold ) );
-  this->m_CurrentThreshold++;
-  this->m_Predicate->SetUpper( *( this->m_CurrentThreshold ) );
-  this->m_Count = 0;
-  this->m_Signal.clear( );
-  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 );
+  this->m_CurrQueue = 0;
+
+  // Prepare iteration over all thresholds
+  this->m_CurrThr = this->m_Thresholds.begin( );
+  this->m_Predicate->SetLower( *( this->m_CurrThr ) );
+  this->m_CurrThr++;
+  this->m_Predicate->SetUpper( *( this->m_CurrThr ) );
+
+  // Prepare counting signal
+  this->m_CurrCount = double( 0 );
+  this->m_PeakDetector.Clear( );
 }
 
 // -------------------------------------------------------------------------
@@ -139,91 +212,25 @@ template< class _TAlgorithm >
 void fpa::Base::Mori< _TAlgorithm >::
 _FinishOneLoop( )
 {
-  if( this->m_Queues[ this->m_CurrentQueue ].size( ) == 0 )
+  if( this->m_Queues[ this->m_CurrQueue ].size( ) == 0 )
   {
-    this->m_Signal.push_back(
-      TSignalData( *this->m_CurrentThreshold, this->m_Count )
+    // Update peak detector
+    TPeak p = this->m_PeakDetector.AddValue(
+      *this->m_CurrThr, this->m_CurrCount
       );
-    this->m_CurrentThreshold++;
-    this->m_CurrentQueue = ( this->m_CurrentQueue + 1 ) % 2;
-    if( this->m_CurrentThreshold != this->m_Thresholds.end( ) )
+    std::cout << *( this->m_CurrThr ) << " " << this->m_CurrCount << std::endl;
+    this->m_CurrThr++;
+    if( this->m_CurrThr != 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;
+      // Update predicate and counting value
+      this->m_Predicate->SetUpper( *( this->m_CurrThr ) );
+      this->m_CurrCount = double( 0 );
 
       // Peak detected? -> stop!
-      if( this->m_SignalPeaks.back( ) == 1 && this->m_MinimumThreshold < *( this->m_CurrentThreshold ) )
+      if(
+        p == TPeakDetector::PosPeak &&
+        this->m_MinimumThreshold < *( this->m_CurrThr )
+        )
         this->_QueueClear( );
     }
     else
@@ -251,13 +258,13 @@ _UpdateOutputValue( TNode& n )
   {
     n.Value = this->m_InitValue;
     n.FrontId++;
-    this->m_Queues[ ( this->m_CurrentQueue + 1 ) % 2 ].push_back( n );
+    this->m_Queues[ ( this->m_CurrQueue + 1 ) % 2 ].push_back( n );
     n.FrontId = 0;
   }
   else
   {
     n.Value = this->m_InsideValue;
-    this->m_Count++;
+    this->m_CurrCount += double( 1 );
 
   } // fi
   this->Superclass::_UpdateOutputValue( n );
@@ -278,8 +285,8 @@ typename fpa::Base::Mori< _TAlgorithm >::
 TNode fpa::Base::Mori< _TAlgorithm >::
 _QueuePop( )
 {
-  TNode n = this->m_Queues[ this->m_CurrentQueue ].front( );
-  this->m_Queues[ this->m_CurrentQueue ].pop_front( );
+  TNode n = this->m_Queues[ this->m_CurrQueue ].front( );
+  this->m_Queues[ this->m_CurrQueue ].pop_front( );
   return( n );
 }
 
@@ -288,7 +295,7 @@ template< class _TAlgorithm >
 void fpa::Base::Mori< _TAlgorithm >::
 _QueuePush( const TNode& node )
 {
-  this->m_Queues[ this->m_CurrentQueue ].push_back( node );
+  this->m_Queues[ this->m_CurrQueue ].push_back( node );
 }
 
 // -------------------------------------------------------------------------
@@ -296,7 +303,7 @@ template< class _TAlgorithm >
 unsigned long fpa::Base::Mori< _TAlgorithm >::
 _QueueSize( ) const
 {
-  return( this->m_Queues[ this->m_CurrentQueue ].size( ) );
+  return( this->m_Queues[ this->m_CurrQueue ].size( ) );
 }
 
 // -------------------------------------------------------------------------