#ifndef __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__
#define __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__

#include <limits>
#include <itkBinaryThresholdImageFilter.h>
#include <itkConstNeighborhoodIterator.h>

// -------------------------------------------------------------------------
template< class I >
void fpa::Image::RegionGrowWithMultipleThresholds< I >::
AddThreshold( const TPixel& v )
{
  this->m_Thresholds.insert( v );
  this->Modified( );
}

// -------------------------------------------------------------------------
template< class I >
void fpa::Image::RegionGrowWithMultipleThresholds< I >::
AddThresholds( const TPixel& t0, const TPixel& t1, const unsigned int& s )
{
  for( TPixel t = t0; t <= t1; t += TPixel( s ) )
    this->AddThreshold( t );
}

// -------------------------------------------------------------------------
template< class I >
fpa::Image::RegionGrowWithMultipleThresholds< I >::
RegionGrowWithMultipleThresholds( )
  : Superclass( ),
    m_InsideValue( TPixel( 1 ) ),
    m_OutsideValue( TPixel( 0 ) ),
    m_DifferenceThreshold( double( 3 ) ),
    m_TotalCount( 0 ),
    m_LastDiff( double( 0 ) ),
    m_StopForced( false ),
    m_StopThreshold( TPixel( 0 ) )
{
}

// -------------------------------------------------------------------------
template< class I >
fpa::Image::RegionGrowWithMultipleThresholds< I >::
~RegionGrowWithMultipleThresholds( )
{
}

// -------------------------------------------------------------------------
template< class I >
void fpa::Image::RegionGrowWithMultipleThresholds< I >::
_BeforeMainLoop( )
{
  const I* img = this->GetInput( );

  // Clear all states
  this->ClearMembershipFunctions( );
  this->m_Histogram.clear( );
  this->m_TotalCount = img->GetLargestPossibleRegion( ).GetNumberOfPixels( );
  this->m_LastDiff = double( -1 );
  this->m_StopForced = false;

  // Initialize thresholding functors
  typename TThresholds::const_iterator tIt = this->m_Thresholds.begin( );
  TPixel min_thr = std::numeric_limits< TPixel >::min( );
  for( ; tIt != this->m_Thresholds.end( ); ++tIt )
  {
    typename TFunction::Pointer function = TFunction::New( );
    function->SetInputImage( img );
    function->SetLowerThreshold( min_thr );
    function->SetUpperThreshold( *tIt );
    this->AddMembershipFunction( function );

  } // rof

  // Correct seeds
  typename I::SizeType radius;
  radius.Fill( 3 );
  itk::ConstNeighborhoodIterator< I > it(
    radius, img, img->GetRequestedRegion( )
    );
  for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s )
  {
    _TNode seed = this->m_Seeds[ s ];
    it.SetLocation( seed.Vertex );

    typename I::SizeType size = it.GetSize( );
    unsigned int nN = 1;
    for( unsigned int d = 0; d < I::ImageDimension; ++d )
      nN *= size[ d ];
    TPixel min_value = img->GetPixel( seed.Vertex );
    for( unsigned int i = 0; i < nN; ++i )
    {
      if( it.GetPixel( i ) < min_value )
      {
        seed.Vertex = it.GetIndex( i );
        seed.Parent = seed.Vertex;
        min_value = it.GetPixel( i );

      } // fi

    } // rof
    this->m_Seeds[ s ] = seed;

  } // rof

  // Continue all initializations
  this->Superclass::_BeforeMainLoop( );
}

// -------------------------------------------------------------------------
template< class I >
void fpa::Image::RegionGrowWithMultipleThresholds< I >::
_AfterMainLoop( )
{
  typedef itk::BinaryThresholdImageFilter< I, I > _TBinFilter;

  // Binarize, inplace, the grown region
  if( this->m_Histogram.size( ) > 1 )
  {
    typename _TBinFilter::Pointer bin = _TBinFilter::New( );
    bin->SetInput( this->GetOutput( ) );
    bin->SetInsideValue( this->m_InsideValue );
    bin->SetOutsideValue( this->m_OutsideValue );
    bin->InPlaceOn( );
    if( this->m_StopForced )
      bin->SetUpperThreshold( this->m_StopThreshold );
    else
      bin->SetUpperThreshold( this->m_Histogram.rbegin( )->first );
    bin->GraftOutput( this->GetOutput( ) );
    bin->Update( );
    this->GraftOutput( bin->GetOutput( ) );

  } // fi

  this->Superclass::_AfterMainLoop( );
}

// -------------------------------------------------------------------------
template< class I >
void fpa::Image::RegionGrowWithMultipleThresholds< I >::
_AfterLoop( )
{
  if( this->m_ActualFunction != this->m_Functions.end( ) )
  {
    TFunction* f =
      dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
    this->m_Histogram[ f->GetUpperThreshold( ) ] = this->m_Marks.size( );

    /* TODO: code to print the signal for demo purposes
       std::cout
       << f->GetUpperThreshold( ) << " "
       << this->m_Marks.size( )
       << std::endl;
    */

    // Get the two last complete count, if any, and compute last
    // finite difference
    if( this->m_Histogram.size( ) > 1 )
    {
      typename THistogram::const_reverse_iterator hIt =
        this->m_Histogram.rbegin( );
      typename THistogram::const_reverse_iterator gIt = hIt;
      gIt++;

      this->m_LastDiff  = double( hIt->second ) - double( gIt->second );
      this->m_LastDiff /= double( hIt->first )  - double( gIt->first );
      this->m_LastDiff *= m_DifferenceThreshold;
    }
    else
      this->m_LastDiff = double( -1 );
    
  } // fi
  this->Superclass::_AfterLoop( );
}

// -------------------------------------------------------------------------
template< class I >
bool fpa::Image::RegionGrowWithMultipleThresholds< I >::
_UpdateResult( _TNode& n )
{
  // Here, the output is changed to the threshold used to compute the
  // growing condition of the actual vertex n
  bool ret = this->Superclass::_UpdateResult( n );
  if( this->m_ActualFunction != this->m_Functions.end( ) )
  {
    TFunction* f =
      dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
    if( f != NULL )
      this->GetOutput( )->SetPixel( n.Vertex, f->GetUpperThreshold( ) );

  } // fi
  return( ret );
}

// -------------------------------------------------------------------------
template< class I >
void fpa::Image::RegionGrowWithMultipleThresholds< I >::
_Mark( const _TNode& n )
{
  this->Superclass::_Mark( n );

  // Check if the histogram's support is enough
  if( this->m_Histogram.size( ) < 2 )
    return;
  typename THistogram::const_reverse_iterator hIt =
    this->m_Histogram.rbegin( );

  // Get the actual function
  TFunction* f =
    dynamic_cast< TFunction* >( this->m_ActualFunction->GetPointer( ) );
  if( f == NULL )
    return;

  // Previous pixel count
  unsigned long prev_aCount = hIt->second;
  double prev_p = double( prev_aCount ) / double( this->m_TotalCount );

  // Actual pixel count
  unsigned long aCount = this->m_Marks.size( );
  double p = double( aCount ) / double( this->m_TotalCount );

  // Loop over, at least, 0.1% of the total number of pixels before
  // performing stop analysis 
  if( double( 1e-3 ) < p && double( 1e-3 ) < prev_p )
  {
    // Does the difference is worthy to be analyzed?
    if( aCount > hIt->second )
    {
      // Compute finite difference and compare it to the previous one
      double diff = double( aCount ) - double( hIt->second );
      diff /= double( f->GetUpperThreshold( ) ) - double( hIt->first );
      if( diff > this->m_LastDiff )
      {
        this->m_StopForced = true;
        this->m_StopThreshold = hIt->first;

      } // fi

    } // fi
    
  } // fi
}

// -------------------------------------------------------------------------
template< class I >
bool fpa::Image::RegionGrowWithMultipleThresholds< I >::
_CheckStopCondition( )
{
  return( this->m_StopForced );
}

#endif // __FPA__IMAGE__REGIONGROWWITHMULTIPLETHRESHOLDS__HXX__

// eof - $RCSfile$