// =========================================================================
// @author Leonardo Florez Valencia
// @email florez-l@javeriana.edu.co
// =========================================================================

#ifndef __fpa__Base__Mori__hxx__
#define __fpa__Base__Mori__hxx__

// -------------------------------------------------------------------------
template< class _TAlgorithm >
itk::ModifiedTimeType fpa::Base::Mori< _TAlgorithm >::
GetMTime( ) const
{
  itk::ModifiedTimeType t = this->Superclass::GetMTime( );
  itk::ModifiedTimeType q = this->m_Predicate->GetMTime( );
  return( ( q < t )? q: t );
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
typename fpa::Base::Mori< _TAlgorithm >::
TOutputValue fpa::Base::Mori< _TAlgorithm >::
GetOutsideValue( ) const
{
  return( this->GetInitValue( ) );
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
void fpa::Base::Mori< _TAlgorithm >::
SetOutsideValue( const TOutputValue& v )
{
  this->SetInitValue( v );
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
void fpa::Base::Mori< _TAlgorithm >::
ClearThresholds( )
{
  this->m_Thresholds.clear( );
  this->Modified( );
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
void fpa::Base::Mori< _TAlgorithm >::
AddThreshold( const TInputValue& thr )
{
  if( this->m_Thresholds.insert( thr ).second )
    this->Modified( );
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
void fpa::Base::Mori< _TAlgorithm >::
SetThresholds(
  const TInputValue& init,
  const TInputValue& end,
  const TInputValue& delta
  )
{
  for( TInputValue thr = init; thr <= end; thr += delta )
    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 = 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( );
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
fpa::Base::Mori< _TAlgorithm >::
~Mori( )
{
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
void fpa::Base::Mori< _TAlgorithm >::
_BeforeGenerateData( )
{
  this->Superclass::_BeforeGenerateData( );

  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 );
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
void fpa::Base::Mori< _TAlgorithm >::
_FinishOneLoop( )
{
  if( this->m_Queues[ this->m_CurrentQueue ].size( ) == 0 )
  {
    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->m_MinimumThreshold < *( this->m_CurrentThreshold ) )
        this->_QueueClear( );
    }
    else
      this->_QueueClear( );

  } // fi
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
void fpa::Base::Mori< _TAlgorithm >::
_ComputeOutputValue( TNode& n )
{
  // Do nothing!!!
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
void fpa::Base::Mori< _TAlgorithm >::
_UpdateOutputValue( TNode& n )
{
  TInputValue value = this->_GetInputValue( n.Vertex );
  bool inside = this->m_Predicate->Evaluate( value );
  if( !inside )
  {
    n.Value = this->m_InitValue;
    n.FrontId++;
    this->m_Queues[ ( this->m_CurrentQueue + 1 ) % 2 ].push_back( n );
    n.FrontId = 0;
  }
  else
  {
    n.Value = this->m_InsideValue;
    this->m_Count++;

  } // fi
  this->Superclass::_UpdateOutputValue( n );
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
void fpa::Base::Mori< _TAlgorithm >::
_QueueClear( )
{
  this->m_Queues[ 0 ].clear( );
  this->m_Queues[ 1 ].clear( );
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
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( );
  return( n );
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
void fpa::Base::Mori< _TAlgorithm >::
_QueuePush( const TNode& node )
{
  this->m_Queues[ this->m_CurrentQueue ].push_back( node );
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
unsigned long fpa::Base::Mori< _TAlgorithm >::
_QueueSize( ) const
{
  return( this->m_Queues[ this->m_CurrentQueue ].size( ) );
}

// -------------------------------------------------------------------------
template< class _TAlgorithm >
void fpa::Base::Mori< _TAlgorithm >::
_PrepareSeeds( TNodes& nodes )
{
  typename TNodes::iterator nIt = nodes.begin( );
  for( ; nIt != nodes.end( ); ++nIt )
    nIt->Value = this->m_InitValue;
}

#endif // __fpa__Base__Mori__hxx__

// eof - $RCSfile$