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

#ifndef __fpa__Image__MoriRegionGrow__hxx__
#define __fpa__Image__MoriRegionGrow__hxx__

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

// -------------------------------------------------------------------------
template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
typename
fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
TAuxiliaryImage*
fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
GetAuxiliaryImage( )
{
  return(
    dynamic_cast< TAuxiliaryImage* >(
      this->itk::ProcessObject::GetOutput( 1 )
      )
    );
}

// -------------------------------------------------------------------------
template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
const typename
fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
TAuxiliaryImage*
fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
GetAuxiliaryImage( ) const
{
  return(
    dynamic_cast< const TAuxiliaryImage* >(
      this->itk::ProcessObject::GetOutput( 1 )
      )
    );
}

// -------------------------------------------------------------------------
template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
void
fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
SetThresholdRange(
  const TInputPixel& lower, const TInputPixel& upper,
  const TInputPixel& delta
  )
{
  this->SetLowerThreshold( lower );
  this->SetUpperThreshold( upper );
  this->SetDeltaThreshold( delta );
}

// -------------------------------------------------------------------------
template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
MoriRegionGrow( )
  : Superclass( )
{
  this->m_InsideValue = TOutputPixel( 1 );
  this->m_OutsideValue = TOutputPixel( 0 );
  this->m_UpperThreshold = std::numeric_limits< TInputPixel >::max( );
  if( std::numeric_limits< TInputPixel >::is_integer )
    this->m_LowerThreshold = std::numeric_limits< TInputPixel >::min( );
  else
    this->m_LowerThreshold = -this->m_UpperThreshold;
  this->m_DeltaThreshold = TInputPixel( 1 );

  this->SetNumberOfRequiredOutputs( 2 );
  this->SetNthOutput( 0, TOutputImage::New( ) );
  this->SetNthOutput( 1, TAuxiliaryImage::New( ) );

  this->m_ThresholdFilter = TThresholdFilter::New( );
}

// -------------------------------------------------------------------------
template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
~MoriRegionGrow( )
{
}

// -------------------------------------------------------------------------
template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
void
fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
GenerateInputRequestedRegion( )
{
  this->Superclass::GenerateInputRequestedRegion( );
  if( this->GetInput( ) )
  {
    TInputImage* input =
      const_cast< TInputImage* >( this->GetInput( ) );
    input->SetRequestedRegionToLargestPossibleRegion( );

  } // fi
}

// -------------------------------------------------------------------------
template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
void
fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
EnlargeOutputRequestedRegion( itk::DataObject* output )
{
  this->Superclass::EnlargeOutputRequestedRegion( output );
  output->SetRequestedRegionToLargestPossibleRegion( );
}

// -------------------------------------------------------------------------
template< class _TInputImage, class _TOutputImage, class _TAuxiliaryPixel >
void fpa::Image::MoriRegionGrow< _TInputImage, _TOutputImage, _TAuxiliaryPixel >::
GenerateData( )
{
  const TInputImage* input = this->GetInput( );
  TAuxiliaryImage* auxiliary = this->GetAuxiliaryImage( );
  TRegion region = input->GetRequestedRegion( );

  // 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, 0.7 );

  // Configure outputs
  auxiliary->SetBufferedRegion( region );
  auxiliary->Allocate( );
  auxiliary->FillBuffer( 0 );

  // Init marks
  typedef itk::Image< bool, TInputImage::ImageDimension > _TMarks;
  typename _TMarks::Pointer marks = _TMarks::New( );
  marks->CopyInformation( input );
  marks->SetRequestedRegion( region );
  marks->SetBufferedRegion( input->GetBufferedRegion( ) );
  marks->Allocate( );
  marks->FillBuffer( false );

  // Prepare queues
  std::queue< TIndex > queues[ 2 ];
  queues[ 0 ].push( this->m_Seed );
  unsigned int cur_queue = 0;
  unsigned int aux_queue = 1;

  // Incremental strategy
  TInputPixel 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
      TIndex node = queues[ cur_queue ].front( );
      queues[ cur_queue ].pop( );
      if( marks->GetPixel( node ) )
        continue;
      marks->SetPixel( node, true );

      // Apply inclusion predicate
      TInputPixel value = input->GetPixel( node );
      bool in = ( ( this->m_LowerThreshold < value ) && ( value < upper ) );
      if( !in )
      {
        if( value < this->m_UpperThreshold )
          queues[ aux_queue ].push( node );
        marks->SetPixel( node, false );
        continue;

      } // fi

      // Ok, pixel lays inside region
      auxiliary->SetPixel( node, this->m_Curve.size( ) + 1 );
      count++;

      // Add neighborhood
      for( unsigned int d = 0; d < TInputImage::ImageDimension; ++d )
      {
        TIndex neigh = node;
        for( int i = -1; i <= 1; i += 2 )
        {
          neigh[ d ] = node[ d ] + i;
          if( region.IsInside( neigh ) )
            queues[ cur_queue ].push( neigh );

        } // rof

      } // 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( );
  int jmax = 0;
  for( int 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 = this->m_Curve[ jmax ].XValue;

  // Final threshold
  this->m_ThresholdFilter->SetInput( auxiliary );
  this->m_ThresholdFilter->SetInsideValue( this->m_InsideValue );
  this->m_ThresholdFilter->SetOutsideValue( this->m_OutsideValue );
  this->m_ThresholdFilter->SetLowerThreshold( 1 );
  this->m_ThresholdFilter->SetUpperThreshold( jmax );
  this->m_ThresholdFilter->Update( );
  this->GetOutput( )->Graft( this->m_ThresholdFilter->GetOutput( ) );
}

#endif // __fpa__Image__MoriRegionGrow__hxx__

// eof - $RCSfile$