]> Creatis software - FrontAlgorithms.git/commitdiff
...
authorLeonardo Flórez-Valencia <florez-l@javeriana.edu.co>
Mon, 10 Jul 2017 21:07:42 +0000 (16:07 -0500)
committerLeonardo Flórez-Valencia <florez-l@javeriana.edu.co>
Mon, 10 Jul 2017 21:07:42 +0000 (16:07 -0500)
CMakeLists.txt
appli/CMakeLists.txt [new file with mode: 0644]
appli/CTBronchi/CMakeLists.txt [new file with mode: 0644]
appli/CTBronchi/CTBronchi.cxx [new file with mode: 0644]
appli/CTBronchi/CTBronchi.ui [new file with mode: 0644]
appli/CTBronchi/MoriSegmentation.cxx [new file with mode: 0644]
lib/fpa/Base/Mori.h
lib/fpa/Base/Mori.hxx
tests/image/CMakeLists.txt
tests/image/MoriSegmentation.cxx

index d1ab775bbc258a20ba91c572d5110e6d14568e14..b098f7a5cebd7f702b6c0392b2edf4a34320e390 100644 (file)
@@ -91,11 +91,7 @@ mark_as_advanced(
 ## == Find needed packages and dependencies ==
 ## ===========================================
 
-find_package(ivq CONFIG QUIET)
-if(NOT ivq_FOUND)
-  find_package(ITK CONFIG REQUIRED)
-  include(${ITK_USE_FILE})
-endif(NOT ivq_FOUND)
+find_package(ivq CONFIG REQUIRED)
 
 ## =========================
 ## == Installation values ==
@@ -113,7 +109,7 @@ set(namespace "${PROJECT_NAME}::")
 ## == Build different parts ==
 ## ===========================
 
-subdirs(lib tests)
+subdirs(appli lib tests)
 
 ## ===============================
 ## == Global installation rules ==
diff --git a/appli/CMakeLists.txt b/appli/CMakeLists.txt
new file mode 100644 (file)
index 0000000..a82c437
--- /dev/null
@@ -0,0 +1,4 @@
+
+subdirs(CTBronchi)
+
+## eof - $RCSfile$
diff --git a/appli/CTBronchi/CMakeLists.txt b/appli/CTBronchi/CMakeLists.txt
new file mode 100644 (file)
index 0000000..446ff4c
--- /dev/null
@@ -0,0 +1,18 @@
+option(fpa_BUILD_CTBronchi "Build bronchi analysis from CT images applications?" OFF)
+if(fpa_BUILD_CTBronchi)
+  include_directories(
+    ${PROJECT_SOURCE_DIR}/lib
+    ${PROJECT_BINARY_DIR}/lib
+    )
+  set(_pfx fpa_CTBronchi_)
+  set(
+    _examples
+    MoriSegmentation
+    )
+  foreach(_e ${_examples})
+    add_executable(${_pfx}${_e} ${_e}.cxx)
+    target_link_libraries(${_pfx}${_e} ivq::ivq fpa)
+  endforeach(_e)
+endif(fpa_BUILD_CTBronchi)
+
+## eof - $RCSfile$
diff --git a/appli/CTBronchi/CTBronchi.cxx b/appli/CTBronchi/CTBronchi.cxx
new file mode 100644 (file)
index 0000000..6427532
--- /dev/null
@@ -0,0 +1,8 @@
+
+int main( int argc, char* argv[] )
+{
+  return( 0 );
+}
+
+
+// eof - $RCSfile$
diff --git a/appli/CTBronchi/CTBronchi.ui b/appli/CTBronchi/CTBronchi.ui
new file mode 100644 (file)
index 0000000..4fa5f30
--- /dev/null
@@ -0,0 +1,123 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<ui version="4.0">
+ <class>MainWindow</class>
+ <widget class="QMainWindow" name="MainWindow">
+  <property name="geometry">
+   <rect>
+    <x>0</x>
+    <y>0</y>
+    <width>800</width>
+    <height>600</height>
+   </rect>
+  </property>
+  <property name="windowTitle">
+   <string>MainWindow</string>
+  </property>
+  <widget class="QWidget" name="centralwidget">
+   <widget class="QWidget" name="Viewer" native="true">
+    <property name="geometry">
+     <rect>
+      <x>60</x>
+      <y>150</y>
+      <width>461</width>
+      <height>331</height>
+     </rect>
+    </property>
+   </widget>
+  </widget>
+  <widget class="QMenuBar" name="menubar">
+   <property name="geometry">
+    <rect>
+     <x>0</x>
+     <y>0</y>
+     <width>800</width>
+     <height>20</height>
+    </rect>
+   </property>
+   <widget class="QMenu" name="menuFile">
+    <property name="title">
+     <string>File</string>
+    </property>
+    <addaction name="aOpenRawImage"/>
+    <addaction name="separator"/>
+    <addaction name="aOpenMoriImage"/>
+    <addaction name="aSaveMoriImage"/>
+    <addaction name="separator"/>
+    <addaction name="aOpenRWImage"/>
+    <addaction name="aSaveRWImage"/>
+    <addaction name="separator"/>
+    <addaction name="aOpenSeedsFile"/>
+    <addaction name="aSaveSeedsFile"/>
+    <addaction name="separator"/>
+    <addaction name="aExit"/>
+   </widget>
+   <addaction name="menuFile"/>
+  </widget>
+  <widget class="QStatusBar" name="statusbar"/>
+  <widget class="QToolBar" name="toolBar">
+   <property name="windowTitle">
+    <string>toolBar</string>
+   </property>
+   <attribute name="toolBarArea">
+    <enum>TopToolBarArea</enum>
+   </attribute>
+   <attribute name="toolBarBreak">
+    <bool>false</bool>
+   </attribute>
+   <addaction name="aMori"/>
+   <addaction name="aRW"/>
+  </widget>
+  <action name="aOpenRawImage">
+   <property name="text">
+    <string>Open raw image</string>
+   </property>
+  </action>
+  <action name="aOpenMoriImage">
+   <property name="text">
+    <string>Open &quot;mori&quot; image</string>
+   </property>
+  </action>
+  <action name="aOpenRWImage">
+   <property name="text">
+    <string>Open &quot;random walker&quot; image</string>
+   </property>
+  </action>
+  <action name="aExit">
+   <property name="text">
+    <string>Exit</string>
+   </property>
+  </action>
+  <action name="aOpenSeedsFile">
+   <property name="text">
+    <string>Open seeds file</string>
+   </property>
+  </action>
+  <action name="aSaveSeedsFile">
+   <property name="text">
+    <string>Save seeds file</string>
+   </property>
+  </action>
+  <action name="aSaveMoriImage">
+   <property name="text">
+    <string>Save &quot;mori&quot; image</string>
+   </property>
+  </action>
+  <action name="aSaveRWImage">
+   <property name="text">
+    <string>Save &quot;random walker&quot; image</string>
+   </property>
+  </action>
+  <action name="aMori">
+   <property name="text">
+    <string>Mori</string>
+   </property>
+  </action>
+  <action name="aRW">
+   <property name="text">
+    <string>RandomWalker</string>
+   </property>
+  </action>
+ </widget>
+ <resources/>
+ <connections/>
+</ui>
diff --git a/appli/CTBronchi/MoriSegmentation.cxx b/appli/CTBronchi/MoriSegmentation.cxx
new file mode 100644 (file)
index 0000000..0b26649
--- /dev/null
@@ -0,0 +1,149 @@
+#include <chrono>
+#include <itkImage.h>
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+#include <fpa/Image/Mori.h>
+
+// -------------------------------------------------------------------------
+const unsigned int Dim = 3;
+typedef short          TPixel;
+typedef unsigned short TLabel;
+
+typedef itk::Image< TPixel, Dim > TInputImage;
+typedef itk::Image< TLabel, Dim > TLabelImage;
+typedef fpa::Image::Mori< TInputImage, TLabelImage > TFilter;
+
+// -------------------------------------------------------------------------
+int main( int argc, char* argv[] )
+{
+  // Get arguments
+  if( argc < 17 )
+  {
+    std::cerr
+      << "Usage: " << argv[ 0 ] << std::endl
+      << "   input_image output_image output_signal" << std::endl
+      << "   init_threshold(-1024) end_threshold(0) delta(1)" << std::endl
+      << "   minimum_threshold(-850)" << std::endl
+      << "   inside_value(255) outside_value(0)" << std::endl
+      << "   signal_kernel_size(20) signal_threshold(500) signal_influence(0.5)"
+      << std::endl
+      << "   [index/point] seed_x seed_y seed_z"
+      << std::endl;
+    return( 1 );
+
+  } // fi
+  std::string input_image_filename = argv[ 1 ];
+  std::string output_image_filename = argv[ 2 ];
+  std::string output_signal_filename = argv[ 3 ];
+  TPixel init_threshold = std::atoi( argv[ 4 ] );
+  TPixel end_threshold = std::atoi( argv[ 5 ] );
+  TPixel delta = std::atoi( argv[ 6 ] );
+  TPixel minimum_threshold = std::atoi( argv[ 7 ] );
+  TLabel inside_value = std::atoi( argv[ 8 ] );
+  TLabel outside_value = std::atoi( argv[ 9 ] );
+  unsigned long signal_kernel_size = std::atoi( argv[ 10 ] );
+  double signal_threshold = std::atof( argv[ 11 ] );
+  double signal_influence = std::atof( argv[ 12 ] );
+  std::string seed_type = argv[ 13 ];
+  TInputImage::IndexType iseed;
+  TInputImage::PointType pseed;
+  for( unsigned int i = 0; i < Dim; ++i )
+  {
+    if( seed_type == "index" )
+      iseed[ i ] = std::atoi( argv[ 14 + i ] );
+    else
+      pseed[ i ] = std::atof( argv[ 14 + i ] );
+
+  } // rof
+
+  // Read image
+  itk::ImageFileReader< TInputImage >::Pointer input_image_reader =
+    itk::ImageFileReader< TInputImage >::New( );
+  input_image_reader->SetFileName( input_image_filename );
+
+  // Prepare filter
+  TFilter::Pointer filter = TFilter::New( );
+  filter->SetInput( input_image_reader->GetOutput( ) );
+  if( seed_type == "index" )
+    filter->SetSeed( iseed );
+  else
+    filter->SetSeed( pseed );
+  filter->SetThresholds( init_threshold, end_threshold, delta );
+  filter->SetMinimumThreshold( minimum_threshold );
+  filter->SetInsideValue( inside_value );
+  filter->SetOutsideValue( outside_value );
+  filter->SetSignalKernelSize( signal_kernel_size );
+  filter->SetSignalThreshold( signal_threshold );
+  filter->SetSignalInfluence( signal_influence );
+
+  // Show some information
+  std::cout << "----------------------------------------------" << std::endl;
+  std::cout << "Image: " << input_image_filename << std::endl;
+
+  // Execute pipeline
+  std::chrono::time_point< std::chrono::high_resolution_clock > ts, te;
+  std::chrono::duration< double > tr;
+  try
+  {
+    ts = std::chrono::high_resolution_clock::now( );
+    input_image_reader->Update( );
+    te = std::chrono::high_resolution_clock::now( );
+    tr = te - ts;
+    std::cout << "Read time: " << tr.count( ) << " s" << std::endl;
+
+    ts = std::chrono::high_resolution_clock::now( );
+    filter->Update( );
+    te = std::chrono::high_resolution_clock::now( );
+    tr = te - ts;
+    std::cout
+      << "Mori time: " << tr.count( ) << " s" << std::endl
+      << "Optimum threshold: " << filter->GetOptimumThreshold( ) << std::endl
+      << "Number of thresholds: "
+      << filter->GetNumberOfEvaluatedThresholds( ) << std::endl;
+  }
+  catch( std::exception& err )
+  {
+    std::cerr << "Error caught: " << err.what( ) << std::endl;
+    return( 1 );
+
+  } // yrt
+
+  // Save output image
+  itk::ImageFileWriter< TLabelImage >::Pointer output_image_writer =
+    itk::ImageFileWriter< TLabelImage >::New( );
+  output_image_writer->SetInput( filter->GetThresholdedOutput( ) );
+  output_image_writer->SetFileName( output_image_filename );
+  try
+  {
+    ts = std::chrono::high_resolution_clock::now( );
+    output_image_writer->Update( );
+    te = std::chrono::high_resolution_clock::now( );
+    tr = te - ts;
+    std::cout << "Write time: " << tr.count( ) << " s" << std::endl;
+  }
+  catch( std::exception& err )
+  {
+    std::cerr << "Error caught: " << err.what( ) << std::endl;
+    return( 1 );
+
+  } // yrt
+
+  // Save output signal
+  std::ofstream osignal( output_signal_filename.c_str( ) );
+  for(
+    unsigned long i = 0; i < filter->GetNumberOfEvaluatedThresholds( ); ++i
+    )
+  {
+    double x, y;
+    TFilter::TPeak p;
+    filter->GetSignalValues( i, x, y, p );
+    osignal << x << " " << y << std::endl;
+
+  } // rof
+  osignal.close( );
+  std::cout << "----------------------------------------------" << std::endl;
+
+  return( 0 );
+}
+
+// eof - $RCSfile$
index d1774d708dcc3da8aca409fd4018666d691b6fbf..27f6915a671a8893805afc193973e9eb491278fc 100644 (file)
@@ -3,19 +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 <deque>
 #include <set>
-#include <vector>
 #include <itkConceptChecking.h>
 #include <itkFunctionBase.h>
 
-#include <fpa/Generic/PeakDetector.h>
+#include <ivq/ITK/PeakDetector.h>
 #include <fpa/Base/Functors/RegionGrow/BinaryThreshold.h>
 
 namespace fpa
@@ -43,10 +39,10 @@ namespace fpa
 
       typedef std::deque< TNode >     TQueue;
       typedef std::set< TInputValue > TThresholds;
-      typedef std::pair< TInputValue, unsigned long > TSignalData;
-      typedef std::vector< TSignalData >              TSignal;
 
-      typedef fpa::Generic::PeakDetector TPeakDetector;
+      typedef ivq::ITK::PeakDetector TPeakDetector;
+      typedef TPeakDetector::TPeak   TPeak;
+
       typedef fpa::Base::Functors::RegionGrow::BinaryThreshold< TInputValue > TPredicate;
 
     public:
@@ -56,19 +52,10 @@ namespace fpa
         );
 
       itkGetConstMacro( InsideValue, TOutputValue );
-      itkGetConstMacro( MinimumThreshold, TInputValue );
-
       itkSetMacro( InsideValue, TOutputValue );
-      itkSetMacro( MinimumThreshold, TInputValue );
-
-      itkGetConstReferenceMacro( Signal, TSignal );
-      itkGetConstMacro( SignalLag, unsigned long );
-      itkGetConstMacro( SignalThreshold, double );
-      itkGetConstMacro( SignalInfluence, double );
 
-      itkSetMacro( SignalLag, unsigned long );
-      itkSetMacro( SignalThreshold, double );
-      itkSetMacro( SignalInfluence, double );
+      itkGetConstMacro( MinimumThreshold, TInputValue );
+      itkSetMacro( MinimumThreshold, TInputValue );
 
     public:
       virtual itk::ModifiedTimeType GetMTime( ) const override;
@@ -76,6 +63,14 @@ namespace fpa
       TOutputValue GetOutsideValue( ) const;
       void SetOutsideValue( const TOutputValue& v );
 
+      unsigned long GetSignalKernelSize( ) const;
+      double GetSignalThreshold( ) const;
+      double GetSignalInfluence( ) const;
+
+      void SetSignalKernelSize( unsigned long k );
+      void SetSignalThreshold( double t );
+      void SetSignalInfluence( double i );
+
       void ClearThresholds( );
       void AddThreshold( const TInputValue& thr );
       void SetThresholds(
@@ -83,8 +78,12 @@ namespace fpa
         const TInputValue& end,
         const TInputValue& delta
         );
+      const TThresholds& GetThresholds( ) const;
       unsigned long GetNumberOfEvaluatedThresholds( ) const;
       TInputValue GetOptimumThreshold( ) const;
+      void GetSignalValues(
+        unsigned long i, double& x, double& y, TPeak& p
+        ) const;
 
     protected:
       Mori( );
@@ -107,27 +106,17 @@ namespace fpa
 
     protected:
       typename TPredicate::Pointer m_Predicate;
-      TThresholds m_Thresholds;
-      typename TThresholds::const_iterator m_CurrentThreshold;
-      TQueue m_Queues[ 2 ];
-      unsigned int m_CurrentQueue;
-      unsigned long m_Count;
 
-      TPeakDetector m_PeakDetector;
+      TOutputValue m_InsideValue;
+      TInputValue  m_MinimumThreshold;
+      TThresholds  m_Thresholds;
+      typename TThresholds::const_iterator m_CurrThr;
 
-      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;
+      TQueue       m_Queues[ 2 ];
+      unsigned int m_CurrQueue;
+      double       m_CurrCount;
 
-      TInputValue  m_MinimumThreshold;
-      TOutputValue m_InsideValue;
+      TPeakDetector m_PeakDetector;
     };
 
   } // ecapseman
index 7a1df9a854107f0f30e880edc926a5dc9345294d..3d1a94dfa63e996cbf36aa580a8ae15bdd964660 100644 (file)
@@ -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( ) );
 }
 
 // -------------------------------------------------------------------------
index b62186cfa093c990d305d9faec19816e3dda6d84..fa58c65a7f2188d4c0836fca668bcd410a506fec 100644 (file)
@@ -12,7 +12,7 @@ set(
   )
 foreach(_e ${_examples})
   add_executable(${_pfx}${_e} ${_e}.cxx)
-  target_link_libraries(${_pfx}${_e} fpa)
+  target_link_libraries(${_pfx}${_e} ivq::ivq fpa)
 endforeach(_e)
 
 ## eof - $RCSfile$
index 7a84a5e467cab0dec20ec7701a2e9aee968c69b4..5a5a472757121fc7ed1425198a9630565a3437d7 100644 (file)
@@ -62,7 +62,7 @@ int main( int argc, char* argv[] )
   filter->SetThresholds( init_threshold, end_threshold, delta );
   filter->SetInsideValue( 255 );
   filter->SetOutsideValue( 0 );
-  filter->SetSignalLag( 20 );
+  filter->SetSignalKernelSize( 20 );
   filter->SetSignalThreshold( 500 );
   filter->SetSignalInfluence( 0.5 );
   filter->SetMinimumThreshold( -850 );
@@ -83,18 +83,22 @@ int main( int argc, char* argv[] )
   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;
+  unsigned long nThr = filter->GetNumberOfEvaluatedThresholds( );
+  for( unsigned long i = 0; i < nThr; ++i )
+  {
+    double x, y;
+    TFilter::TPeak p;
+    filter->GetSignalValues( i, x, y, p );
+    osignal << x << " " << y << std::endl;
+
+  } // rof
   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
+    << "Optimum threshold: " << filter->GetOptimumThreshold( ) << std::endl
+    << "Number of evaluated thresholds: " << nThr << std::endl
     << "------------------------------------------------------"
     << std::endl;