...

diff --git a/lib/fpa/Base/Mori.h b/lib/fpa/Base/Mori.h
index c2ea39cff2c729c2f786a05c8f5cf5cee76a0c7e..d1774d708dcc3da8aca409fd4018666d691b6fbf 100644 (file)
--- a/lib/fpa/Base/Mori.h
+++ b/lib/fpa/Base/Mori.h
@@ -14,6 +14,8 @@
 #include <vector>
 #include <itkConceptChecking.h>
 #include <itkFunctionBase.h>
+
+#include <fpa/Generic/PeakDetector.h>
 #include <fpa/Base/Functors/RegionGrow/BinaryThreshold.h>
 
 namespace fpa
@@ -43,6 +45,8 @@ namespace fpa
       typedef std::set< TInputValue > TThresholds;
       typedef std::pair< TInputValue, unsigned long > TSignalData;
       typedef std::vector< TSignalData >              TSignal;
+
+      typedef fpa::Generic::PeakDetector TPeakDetector;
       typedef fpa::Base::Functors::RegionGrow::BinaryThreshold< TInputValue > TPredicate;
 
     public:
@@ -52,7 +56,10 @@ namespace fpa
         );
 
       itkGetConstMacro( InsideValue, TOutputValue );
+      itkGetConstMacro( MinimumThreshold, TInputValue );
+
       itkSetMacro( InsideValue, TOutputValue );
+      itkSetMacro( MinimumThreshold, TInputValue );
 
       itkGetConstReferenceMacro( Signal, TSignal );
       itkGetConstMacro( SignalLag, unsigned long );
@@ -106,6 +113,8 @@ namespace fpa
       unsigned int m_CurrentQueue;
       unsigned long m_Count;
 
+      TPeakDetector m_PeakDetector;
+
       TSignal m_Signal;
       unsigned long m_SignalLag;
       double m_SignalThreshold;
@@ -117,6 +126,7 @@ namespace fpa
       double m_CurrentAverage;
       double m_CurrentVariance;
 
+      TInputValue  m_MinimumThreshold;
       TOutputValue m_InsideValue;
     };
 

diff --git a/lib/fpa/Base/Mori.hxx b/lib/fpa/Base/Mori.hxx
index 17415e772bae2eba47f4fb777f2b4a6be1e1d7d4..7a1df9a854107f0f30e880edc926a5dc9345294d 100644 (file)
--- a/lib/fpa/Base/Mori.hxx
+++ b/lib/fpa/Base/Mori.hxx
@@ -97,6 +97,11 @@ Mori( )
   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( );
 }
 
 // -------------------------------------------------------------------------
@@ -218,7 +223,7 @@ _FinishOneLoop( )
       this->m_Count = 0;
 
       // Peak detected? -> stop!
-      if( this->m_SignalPeaks.back( ) == 1 )
+      if( this->m_SignalPeaks.back( ) == 1 && this->m_MinimumThreshold < *( this->m_CurrentThreshold ) )
         this->_QueueClear( );
     }
     else

diff --git a/lib/fpa/CMakeLists.txt b/lib/fpa/CMakeLists.txt
index 3dba4df8cba0cd95fbecce0e3b0713fac70813ac..e57fdbd16eb0ad742b9de293a0baeefcf44326c0 100644 (file)
--- a/lib/fpa/CMakeLists.txt
+++ b/lib/fpa/CMakeLists.txt
@@ -12,16 +12,21 @@ file(GLOB_RECURSE _image_src "${CMAKE_CURRENT_SOURCE_DIR}/Image/*.cxx")
 file(GLOB_RECURSE _image_hdr "${CMAKE_CURRENT_SOURCE_DIR}/Image/*.h")
 file(GLOB_RECURSE _image_hrc "${CMAKE_CURRENT_SOURCE_DIR}/Image/*.hxx")
 
+file(GLOB_RECURSE _generic_src "${CMAKE_CURRENT_SOURCE_DIR}/Generic/*.cxx")
+file(GLOB_RECURSE _generic_hdr "${CMAKE_CURRENT_SOURCE_DIR}/Generic/*.h")
+file(GLOB_RECURSE _generic_hrc "${CMAKE_CURRENT_SOURCE_DIR}/Generic/*.hxx")
+
+
 set(_src
-  ${_base_src} ${_image_src}
+  ${_base_src} ${_image_src} ${_generic_src}
   "${CMAKE_CURRENT_BINARY_DIR}/Version.cxx"
   )
 set(
   _hdr
-  ${_base_hdr} ${_image_hdr}
+  ${_base_hdr} ${_image_hdr} ${_generic_hdr}
   "${CMAKE_CURRENT_BINARY_DIR}/Config.h"
   )
-set(_hrc ${_base_hrc} ${_image_hrc})
+set(_hrc ${_base_hrc} ${_image_hrc} ${_generic_hrc})
 
 ## =====================
 ## == Compile library ==

diff --git a/lib/fpa/Generic/PeakDetector.cxx b/lib/fpa/Generic/PeakDetector.cxx
new file mode 100644 (file)
index 0000000..5e1043e
--- /dev/null
+++ b/lib/fpa/Generic/PeakDetector.cxx
@@ -0,0 +1,160 @@
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#include <fpa/Generic/PeakDetector.h>
+#include <cmath>
+
+// -------------------------------------------------------------------------
+fpa::Generic::PeakDetector::
+PeakDetector( )
+  : m_K( 3 ),
+    m_T( 3.5 ),
+    m_I( 0.5 )
+{
+}
+
+// -------------------------------------------------------------------------
+fpa::Generic::PeakDetector::
+~PeakDetector( )
+{
+}
+
+// -------------------------------------------------------------------------
+unsigned long fpa::Generic::PeakDetector::
+GetKernelSize( ) const
+{
+  return( this->m_K );
+}
+
+// -------------------------------------------------------------------------
+double fpa::Generic::PeakDetector::
+GetThreshold( ) const
+{
+  return( this->m_T );
+}
+
+// -------------------------------------------------------------------------
+double fpa::Generic::PeakDetector::
+GetInfluence( ) const
+{
+  return( this->m_I );
+}
+
+// -------------------------------------------------------------------------
+void fpa::Generic::PeakDetector::
+SetKernelSize( unsigned long k )
+{
+  this->m_K = k;
+  this->Clear( );
+}
+
+// -------------------------------------------------------------------------
+void fpa::Generic::PeakDetector::
+SetThreshold( double t )
+{
+  this->m_T = t;
+  this->Clear( );
+}
+
+// -------------------------------------------------------------------------
+void fpa::Generic::PeakDetector::
+SetInfluence( double i )
+{
+  this->m_I = i;
+  this->Clear( );
+}
+
+// -------------------------------------------------------------------------
+void fpa::Generic::PeakDetector::
+Clear( )
+{
+  this->m_X.clear( );
+  this->m_Y.clear( );
+  this->m_YF.clear( );
+  this->m_Avg.clear( );
+  this->m_STD.clear( );
+  this->m_Peaks.clear( );
+  this->m_Mean = 0;
+  this->m_Vari = 0;
+}
+
+// -------------------------------------------------------------------------
+unsigned long fpa::Generic::PeakDetector::
+GetNumberOfSamples( ) const
+{
+  return( this->m_X.size( ) );
+}
+
+// -------------------------------------------------------------------------
+fpa::Generic::PeakDetector::
+TPeak fpa::Generic::PeakDetector::
+AddValue( double x, double y )
+{
+  this->m_X.push_back( x );
+  this->m_Y.push_back( y );
+
+  if( this->m_YF.size( ) < this->m_K )
+  {
+    double n = double( this->m_YF.size( ) + 1 );
+    this->m_YF.push_back( y );
+    this->m_Avg.push_back( double( 0 ) );
+    this->m_STD.push_back( double( 0 ) );
+    this->m_Peaks.push_back( Self::NoPeak );
+    if( n > double( 1 ) )
+      this->m_Vari =
+        ( ( ( n - 2.0 ) / ( n - 1.0 ) ) * this->m_Vari ) +
+        ( ( ( y - this->m_Mean ) * ( y - this->m_Mean ) ) / n );
+    this->m_Mean += ( y - this->m_Mean ) / n;
+    if( this->m_YF.size( ) == this->m_K )
+    {
+      this->m_Avg.push_back( this->m_Mean );
+      this->m_STD.push_back( std::sqrt( this->m_Vari ) );
+
+    } // fi
+  }
+  else
+  {
+    unsigned long i = this->m_X.size( ) - 1;
+    if(
+      ( std::fabs( y - this->m_Avg[ i - 1 ] ) ) >
+      ( this->m_T * this->m_STD[ i - 1 ] )
+      )
+    {
+      this->m_Peaks.push_back(
+        ( y > this->m_Avg[ i - 1 ] )? Self::PosPeak: Self::NegPeak
+        );
+      this->m_YF.push_back(
+        ( this->m_I * y ) + ( ( 1.0 - this->m_I ) * this->m_YF[ i - 1 ] )
+        );
+    }
+    else
+    {
+      this->m_Peaks.push_back( Self::NoPeak );
+      this->m_YF.push_back( y );
+
+    } // fi
+
+    double avg = double( 0 );
+    double var = double( 0 );
+    unsigned long k = 0;
+    for( unsigned long j = i - this->m_K; j <= i; ++j, ++k )
+    {
+      double v = this->m_YF[ 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_Avg.push_back( avg );
+    this->m_STD.push_back( std::sqrt( var ) );
+
+  } // fi
+  return( this->m_Peaks.back( ) );
+}
+
+// eof - $RCSfile$

diff --git a/lib/fpa/Generic/PeakDetector.h b/lib/fpa/Generic/PeakDetector.h
new file mode 100644 (file)
index 0000000..cc8308b
--- /dev/null
+++ b/lib/fpa/Generic/PeakDetector.h
@@ -0,0 +1,68 @@
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Generic__PeakDetector__h__
+#define __fpa__Generic__PeakDetector__h__
+
+#include <fpa/fpa_export.h>
+#include <vector>
+
+namespace fpa
+{
+  namespace Generic
+  {
+    /**
+     * https://stackoverflow.com/questions/22583391/peak-signal-detection-in-realtime-timeseries-data
+     */
+    class FPA_EXPORT PeakDetector
+    {
+    public:
+      typedef PeakDetector Self;
+
+      enum TPeak
+      {
+        NoPeak = 0,
+        PosPeak,
+        NegPeak
+      };
+
+    public:
+      PeakDetector( );
+      virtual ~PeakDetector( );
+
+      unsigned long GetKernelSize( ) const;
+      double GetThreshold( ) const;
+      double GetInfluence( ) const;
+
+      void SetKernelSize( unsigned long k );
+      void SetThreshold( double t );
+      void SetInfluence( double i );
+
+      void Clear( );
+      unsigned long GetNumberOfSamples( ) const;
+      TPeak AddValue( double x, double y );
+
+    protected:
+      unsigned long m_K;
+      double m_T;
+      double m_I;
+
+      std::vector< double > m_X;
+      std::vector< double > m_Y;
+      std::vector< double > m_YF;
+      std::vector< double > m_Avg;
+      std::vector< double > m_STD;
+      std::vector< TPeak >  m_Peaks;
+      double m_Mean;
+      double m_Vari;
+    };
+
+  } // ecapseman
+
+} // ecapseman
+
+#endif // __fpa__Generic__PeakDetector__h__
+
+// eof - $RCSfile$

diff --git a/tests/image/MoriSegmentation.cxx b/tests/image/MoriSegmentation.cxx
index f2dd92a0e8147e22436acf13e8fbcc78d93d68d3..7a84a5e467cab0dec20ec7701a2e9aee968c69b4 100644 (file)
--- a/tests/image/MoriSegmentation.cxx
+++ b/tests/image/MoriSegmentation.cxx
@@ -65,6 +65,7 @@ int main( int argc, char* argv[] )
   filter->SetSignalLag( 20 );
   filter->SetSignalThreshold( 500 );
   filter->SetSignalInfluence( 0.5 );
+  filter->SetMinimumThreshold( -850 );
 
   // Execute filter
   std::chrono::time_point< std::chrono::high_resolution_clock > tstart, tend;