From: Leonardo Flórez-Valencia Date: Thu, 6 Jul 2017 21:31:58 +0000 (-0500) Subject: ... X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=FrontAlgorithms.git;a=commitdiff_plain;h=53d56cb3d8fe139843d5b2308f821cc05e7593e1 ... --- diff --git a/lib/fpa/Base/Mori.h b/lib/fpa/Base/Mori.h index ef8d7cc..c2ea39c 100644 --- a/lib/fpa/Base/Mori.h +++ b/lib/fpa/Base/Mori.h @@ -3,12 +3,15 @@ // @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 #include -#include +#include #include #include #include @@ -39,7 +42,7 @@ namespace fpa 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: @@ -51,6 +54,15 @@ namespace fpa 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; @@ -64,13 +76,14 @@ namespace fpa 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; @@ -92,7 +105,17 @@ namespace fpa 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; }; 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( ); diff --git a/lib/fpa/Image/Mori.h b/lib/fpa/Image/Mori.h index a81016c..550864d 100644 --- a/lib/fpa/Image/Mori.h +++ b/lib/fpa/Image/Mori.h @@ -50,34 +50,32 @@ namespace fpa 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 +#endif // ITK_MANUAL_INSTANTIATION + #endif // __fpa__Image__Mori__h__ // eof - $RCSfile$ diff --git a/lib/fpa/Image/Mori.hxx b/lib/fpa/Image/Mori.hxx new file mode 100644 index 0000000..d33e5ec --- /dev/null +++ b/lib/fpa/Image/Mori.hxx @@ -0,0 +1,78 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#ifndef __fpa__Image__Mori__hxx__ +#define __fpa__Image__Mori__hxx__ + +#include + +// ------------------------------------------------------------------------- +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$ diff --git a/tests/image/MoriSegmentation.cxx b/tests/image/MoriSegmentation.cxx index d1b26f5..f2dd92a 100644 --- a/tests/image/MoriSegmentation.cxx +++ b/tests/image/MoriSegmentation.cxx @@ -1,4 +1,5 @@ #include "BaseFunctions.h" +#include #include #include @@ -15,11 +16,11 @@ typedef fpa::Image::Mori< TInputImage, TLabelImage > TFilter; 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 ); @@ -27,20 +28,21 @@ int main( int argc, char* argv[] ) } // 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 @@ -60,18 +62,41 @@ int main( int argc, char* argv[] ) 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 ); }