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

#ifndef __fpa__Base__MoriRegionGrow__hxx__
#define __fpa__Base__MoriRegionGrow__hxx__

#include <queue>
#include <itkProgressReporter.h>

// -------------------------------------------------------------------------
template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
void
fpa::Base::MoriRegionGrow< _TFilter, _TMarksInterface, _TSeedsInterface >::
SetThresholdRange(
  const TInputValue& lower, const TInputValue& upper,
  const TInputValue& delta
  )
{
  this->SetLowerThreshold( lower );
  this->SetUpperThreshold( upper );
  this->SetDeltaThreshold( delta );
}

// -------------------------------------------------------------------------
template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
fpa::Base::MoriRegionGrow< _TFilter, _TMarksInterface, _TSeedsInterface >::
MoriRegionGrow( )
  : Superclass( ),
    _TMarksInterface( this ),
    _TSeedsInterface( this )
{
  this->m_UpperThreshold = std::numeric_limits< TInputValue >::max( );
  if( std::numeric_limits< TInputValue >::is_integer )
    this->m_LowerThreshold = std::numeric_limits< TInputValue >::min( );
  else
    this->m_LowerThreshold = -this->m_UpperThreshold;
  this->m_DeltaThreshold = TInputValue( 1 );
}

// -------------------------------------------------------------------------
template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
fpa::Base::MoriRegionGrow< _TFilter, _TMarksInterface, _TSeedsInterface >::
~MoriRegionGrow( )
{
}

// -------------------------------------------------------------------------
template< class _TFilter, class _TMarksInterface, class _TSeedsInterface >
void
fpa::Base::MoriRegionGrow< _TFilter, _TMarksInterface, _TSeedsInterface >::
GenerateData( )
{
  // Prepare progress counter
  double prog_size = double( this->m_UpperThreshold - this->m_LowerThreshold );
  prog_size /= double( this->m_DeltaThreshold );
  itk::ProgressReporter progress( this, 0, prog_size, 100, 0, 1 );

  // Init objects
  this->_ConfigureOutputs( std::numeric_limits< TOutputValue >::max( ) );
  this->_InitMarks( this->GetNumberOfSeeds( ) );

  // Init queues
  typedef std::pair< TVertex, unsigned long > _TNode;
  std::queue< _TNode > queues[ 2 ];
  unsigned long frontId = 1;
  typename TSeedsInterface::TSeeds::const_iterator sIt = this->BeginSeeds( );
  for( ; sIt != this->EndSeeds( ); ++sIt )
    queues[ 0 ].push( _TNode( *sIt, frontId++ ) );
  unsigned int cur_queue = 0;
  unsigned int aux_queue = 1;

  // Incremental strategy
  TInputValue upper = this->m_LowerThreshold + this->m_DeltaThreshold;
  this->m_Curve.clear( );
  unsigned long count = 0;
  while( upper <= this->m_UpperThreshold )
  {
    // Main growing strategy
    while( queues[ cur_queue ].size( ) > 0 )
    {
      // Get next candidate
      _TNode node = queues[ cur_queue ].front( );
      queues[ cur_queue ].pop( );
      if( this->_IsMarked( node.first ) )
        continue;
      this->_Mark( node.first, node.second );

      // Apply inclusion predicate
      TInputValue value = this->_GetInputValue( node.first );
      bool in = ( ( this->m_LowerThreshold < value ) && ( value < upper ) );
      if( !in )
      {
        if( value < this->m_UpperThreshold )
          queues[ aux_queue ].push( node );
        this->_Mark( node.first, 0 );
        continue;

      } // fi

      // Ok, value lays inside region
      this->_SetOutputValue( node.first, this->m_Curve.size( ) + 1 );
      count++;

      // Add neighborhood
      TVertices neighbors = this->_GetNeighbors( node.first );
      typename TVertices::const_iterator neighIt = neighbors.begin( );
      for( ; neighIt != neighbors.end( ); ++neighIt )
      {
        TVertex neigh = *neighIt;
        if( this->_IsMarked( neigh ) )
        {
          // Invoke stop at collisions
          unsigned long nColl = this->_Collisions( node.first, neigh );
          if(
            this->StopAtOneFront( ) &&
            this->GetNumberOfSeeds( ) > 1 &&
            nColl == 1
            )
          {
            while( queues[ 0 ].size( ) > 0 )
              queues[ 0 ].pop( );
            while( queues[ 1 ].size( ) > 0 )
              queues[ 1 ].pop( );

          } // fi
        }
        else
          queues[ cur_queue ].push( _TNode( neigh, node.second ) );

      } // rof

    } // elihw

    // Update curve
    if( this->m_Curve.size( ) > 0 )
    {
      if( this->m_Curve.back( ).YValue < count )
        this->m_Curve.push_back( TCurveData( upper, count ) );
      if( this->m_Curve.size( ) > 2 )
      {
        long j = this->m_Curve.size( ) - 2;
        double yp = double( this->m_Curve[ j + 1 ].YValue );
        double yn = double( this->m_Curve[ j - 1 ].YValue );
        double xp = double( this->m_Curve[ j + 1 ].XValue );
        double xn = double( this->m_Curve[ j - 1 ].XValue );
        this->m_Curve[ j ].Diff1 = ( yp - yn ) / ( xp - xn );

      } // fi
    }
    else
      this->m_Curve.push_back( TCurveData( upper, count ) );

    // Update queue
    cur_queue = aux_queue;
    aux_queue = ( cur_queue + 1 ) % 2;

    // Update threshold
    upper += this->m_DeltaThreshold;
    progress.CompletedPixel( );

  } // elihw

  // Compute optimum threshold
  double dmax = -std::numeric_limits< double >::max( );
  long jmax = 0;
  for( long j = 1; j < this->m_Curve.size( ) - 1; ++j )
  {
    double dp = this->m_Curve[ j + 1 ].Diff1;
    double dn = this->m_Curve[ j - 1 ].Diff1;
    double xp = double( this->m_Curve[ j + 1 ].XValue );
    double xn = double( this->m_Curve[ j - 1 ].XValue );
    double d2 = ( dp - dn ) / ( xp - xn );
    if( d2 > dmax )
    {
      dmax = d2;
      jmax = j;

    } // fi

  } // rof
  this->m_OptimumThreshold = TOutputValue( jmax );
  this->_FreeMarks( );
}

#endif // __fpa__Base__MoriRegionGrow__hxx__

// eof - $RCSfile$