// @email florez-l@javeriana.edu.co
// =========================================================================
+// https://stackoverflow.com/questions/22583391/peak-signal-detection-in-realtime-timeseries-data
+
+
#ifndef __fpa__Base__Mori__h__
#define __fpa__Base__Mori__h__
#include <deque>
#include <set>
-#include <map>
+#include <vector>
#include <itkConceptChecking.h>
#include <itkFunctionBase.h>
#include <fpa/Base/Functors/RegionGrow/BinaryThreshold.h>
typedef std::deque< TNode > TQueue;
typedef std::set< TInputValue > TThresholds;
typedef std::pair< TInputValue, unsigned long > TSignalData;
- typedef std::map< TFrontId, TSignalData > TSignal;
+ typedef std::vector< TSignalData > TSignal;
typedef fpa::Base::Functors::RegionGrow::BinaryThreshold< TInputValue > TPredicate;
public:
itkGetConstMacro( InsideValue, TOutputValue );
itkSetMacro( InsideValue, TOutputValue );
+ itkGetConstReferenceMacro( Signal, TSignal );
+ itkGetConstMacro( SignalLag, unsigned long );
+ itkGetConstMacro( SignalThreshold, double );
+ itkGetConstMacro( SignalInfluence, double );
+
+ itkSetMacro( SignalLag, unsigned long );
+ itkSetMacro( SignalThreshold, double );
+ itkSetMacro( SignalInfluence, double );
+
public:
virtual itk::ModifiedTimeType GetMTime( ) const override;
const TInputValue& end,
const TInputValue& delta
);
+ unsigned long GetNumberOfEvaluatedThresholds( ) const;
+ TInputValue GetOptimumThreshold( ) const;
protected:
Mori( );
virtual ~Mori( );
virtual void _BeforeGenerateData( ) override;
- virtual void _AfterGenerateData( ) override;
virtual void _FinishOneLoop( ) override;
virtual void _ComputeOutputValue( TNode& n ) override;
virtual void _UpdateOutputValue( TNode& n ) override;
TQueue m_Queues[ 2 ];
unsigned int m_CurrentQueue;
unsigned long m_Count;
+
TSignal m_Signal;
+ unsigned long m_SignalLag;
+ double m_SignalThreshold;
+ double m_SignalInfluence;
+ std::vector< double > m_FilteredSignal;
+ std::vector< double > m_SignalAverages;
+ std::vector< double > m_SignalDeviations;
+ std::vector< short > m_SignalPeaks;
+ double m_CurrentAverage;
+ double m_CurrentVariance;
TOutputValue m_InsideValue;
};
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 ) );
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 );
}
// -------------------------------------------------------------------------
{
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( );
itkTypeMacro( fpa::Image::Mori, fpa::Base::Mori );
public:
- TMarks* GetOutputLevels( )
- {
- return( this->GetMarks( ) );
- }
- const TMarks* GetOutputLevels( ) const
- {
- return( this->GetMarks( ) );
- }
+ TOutputImage* GetThresholdedOutput( );
+ const TOutputImage* GetThresholdedOutput( ) const;
protected:
- Mori( )
- : Superclass( )
- {
- }
- virtual ~Mori( )
- {
- }
+ Mori( );
+ virtual ~Mori( );
+
+ virtual void _AfterGenerateData( ) override;
private:
// Purposely not implemented.
Mori( const Self& other );
Self& operator=( const Self& other );
+
+ protected:
+ unsigned long m_ThresholdedOutputIdx;
};
} // ecapseman
} // ecapseman
+#ifndef ITK_MANUAL_INSTANTIATION
+# include <fpa/Image/Mori.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+
#endif // __fpa__Image__Mori__h__
// eof - $RCSfile$
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#ifndef __fpa__Image__Mori__hxx__
+#define __fpa__Image__Mori__hxx__
+
+#include <itkBinaryThresholdImageFilter.h>
+
+// -------------------------------------------------------------------------
+template< class _TInputImage, class _TOutputImage >
+typename fpa::Image::Mori< _TInputImage, _TOutputImage >::
+TOutputImage* fpa::Image::Mori< _TInputImage, _TOutputImage >::
+GetThresholdedOutput( )
+{
+ return(
+ dynamic_cast< TOutputImage* >(
+ this->itk::ProcessObject::GetOutput( this->m_ThresholdedOutputIdx )
+ )
+ );
+}
+
+// -------------------------------------------------------------------------
+template< class _TInputImage, class _TOutputImage >
+const typename fpa::Image::Mori< _TInputImage, _TOutputImage >::
+TOutputImage* fpa::Image::Mori< _TInputImage, _TOutputImage >::
+GetThresholdedOutput( ) const
+{
+ return(
+ dynamic_cast< const TOutputImage* >(
+ this->itk::ProcessObject::GetOutput( this->m_ThresholdedOutputIdx )
+ )
+ );
+}
+
+// -------------------------------------------------------------------------
+template< class _TInputImage, class _TOutputImage >
+fpa::Image::Mori< _TInputImage, _TOutputImage >::
+Mori( )
+ : Superclass( )
+{
+ this->m_ThresholdedOutputIdx = this->GetNumberOfRequiredOutputs( );
+ this->itk::ProcessObject::SetNumberOfRequiredOutputs(
+ this->m_ThresholdedOutputIdx + 1
+ );
+ this->SetNthOutput( this->m_ThresholdedOutputIdx, TOutputImage::New( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TInputImage, class _TOutputImage >
+fpa::Image::Mori< _TInputImage, _TOutputImage >::
+~Mori( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class _TInputImage, class _TOutputImage >
+void fpa::Image::Mori< _TInputImage, _TOutputImage >::
+_AfterGenerateData( )
+{
+ typedef itk::BinaryThresholdImageFilter< TMarks, TOutputImage > _TFilter;
+
+ this->Superclass::_AfterGenerateData( );
+
+ typename _TFilter::Pointer filter = _TFilter::New( );
+ filter->SetInput( this->GetMarks( ) );
+ filter->SetInsideValue( this->GetInsideValue( ) );
+ filter->SetOutsideValue( this->GetOutsideValue( ) );
+ filter->SetLowerThreshold( 1 );
+ filter->SetUpperThreshold( this->GetNumberOfEvaluatedThresholds( ) - 1 );
+ filter->Update( );
+ this->GetThresholdedOutput( )->Graft( filter->GetOutput( ) );
+}
+
+#endif // __fpa__Image__Mori__hxx__
+
+// eof - $RCSfile$
#include "BaseFunctions.h"
+#include <chrono>
#include <itkImage.h>
#include <fpa/Image/Mori.h>
int main( int argc, char* argv[] )
{
// Get arguments
- if( argc < 8 + Dim )
+ if( argc < 9 + Dim )
{
std::cerr
<< "Usage: " << argv[ 0 ]
- << " input_image output_image output_levels"
+ << " input_image output_image output_marks output_signal"
<< " init_threshold end_threshold delta [index/point] seed"
<< std::endl;
return( 1 );
} // fi
std::string input_image_filename = argv[ 1 ];
std::string output_image_filename = argv[ 2 ];
- std::string output_levels_filename = argv[ 3 ];
- TPixel init_threshold = std::atoi( argv[ 4 ] );
- TPixel end_threshold = std::atoi( argv[ 5 ] );
- TPixel delta = std::atoi( argv[ 6 ] );
- std::string seed_type = argv[ 7 ];
+ std::string output_marks_filename = argv[ 3 ];
+ std::string output_signal_filename = argv[ 4 ];
+ TPixel init_threshold = std::atoi( argv[ 5 ] );
+ TPixel end_threshold = std::atoi( argv[ 6 ] );
+ TPixel delta = std::atoi( argv[ 7 ] );
+ std::string seed_type = argv[ 8 ];
TInputImage::IndexType iseed;
TInputImage::PointType pseed;
for( unsigned int i = 0; i < Dim; ++i )
{
if( seed_type == "index" )
- iseed[ i ] = std::atoi( argv[ 8 + i ] );
+ iseed[ i ] = std::atoi( argv[ 9 + i ] );
else
- pseed[ i ] = std::atof( argv[ 8 + i ] );
+ pseed[ i ] = std::atof( argv[ 9 + i ] );
} // rof
filter->SetThresholds( init_threshold, end_threshold, delta );
filter->SetInsideValue( 255 );
filter->SetOutsideValue( 0 );
+ filter->SetSignalLag( 20 );
+ filter->SetSignalThreshold( 500 );
+ filter->SetSignalInfluence( 0.5 );
// Execute filter
+ std::chrono::time_point< std::chrono::high_resolution_clock > tstart, tend;
+ tstart = std::chrono::high_resolution_clock::now( );
filter->Update( );
+ tend = std::chrono::high_resolution_clock::now( );
+ std::chrono::duration< double > telapsed = tend - tstart;
// Save results
std::string err1 =
- fpa::tests::image::Write( filter->GetOutput( ), output_image_filename );
+ fpa::tests::image::Write( filter->GetThresholdedOutput( ), output_image_filename );
std::string err2 =
- fpa::tests::image::Write( filter->GetOutputLevels( ), output_levels_filename );
+ fpa::tests::image::Write( filter->GetMarks( ), output_marks_filename );
if( err1 != "" ) std::cerr << err1 << std::endl;
if( err2 != "" ) std::cerr << err2 << std::endl;
+ std::ofstream osignal( output_signal_filename.c_str( ) );
+ const TFilter::TSignal& signal = filter->GetSignal( );
+ for( unsigned long i = 0; i < signal.size( ); ++i )
+ osignal << signal[ i ].first << " " << signal[ i ].second << std::endl;
+ osignal.close( );
+
+ std::cout
+ << "------------------------------------------------------" << std::endl
+ << "Elapsed time: " << telapsed.count( ) << " s" << std::endl
+ << "Optimum threshold: "
+ << filter->GetOptimumThreshold( ) << std::endl
+ << "Number of evaluated thresholds: "
+ << filter->GetNumberOfEvaluatedThresholds( ) << std::endl
+ << "------------------------------------------------------"
+ << std::endl;
+
return( 0 );
}
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff