From b6babebfd77364b12d81fb18d767f484ad1a0969 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Fri, 7 Jul 2017 15:48:01 -0500 Subject: [PATCH] ... --- lib/fpa/Base/Mori.h | 10 ++ lib/fpa/Base/Mori.hxx | 7 +- lib/fpa/CMakeLists.txt | 11 ++- lib/fpa/Generic/PeakDetector.cxx | 160 +++++++++++++++++++++++++++++++ lib/fpa/Generic/PeakDetector.h | 68 +++++++++++++ tests/image/MoriSegmentation.cxx | 1 + 6 files changed, 253 insertions(+), 4 deletions(-) create mode 100644 lib/fpa/Generic/PeakDetector.cxx create mode 100644 lib/fpa/Generic/PeakDetector.h diff --git a/lib/fpa/Base/Mori.h b/lib/fpa/Base/Mori.h index c2ea39c..d1774d7 100644 --- a/lib/fpa/Base/Mori.h +++ b/lib/fpa/Base/Mori.h @@ -14,6 +14,8 @@ #include #include #include + +#include #include 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 17415e7..7a1df9a 100644 --- 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 3dba4df..e57fdbd 100644 --- 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 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 +#include + +// ------------------------------------------------------------------------- +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 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 +#include + +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 f2dd92a..7a84a5e 100644 --- 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; -- 2.45.0