// -------------------------------------------------------------------------
// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
// -------------------------------------------------------------------------

#ifndef __cpExtensions__Algorithms__MultiScaleGaussianImageFilter__hxx__
#define __cpExtensions__Algorithms__MultiScaleGaussianImageFilter__hxx__

#include <vnl/vnl_vector.h>

#include <itkImageRegionIterator.h>
#include <itkImageRegionConstIterator.h>
#include <itkNumericTraits.h>
#include <itkProgressAccumulator.h>

#include <itkBinaryFunctorImageFilter.h>
#include <itkUnaryFunctorImageFilter.h>
#include <itkGradientRecursiveGaussianImageFilter.h>

// -------------------------------------------------------------------------
template< class I, class O >
cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::_Greater::
_Greater( )
{
}

// -------------------------------------------------------------------------
template< class I, class O >
cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::_Greater::
~_Greater( )
{
}

// -------------------------------------------------------------------------
template< class I, class O >
bool cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::_Greater::
operator!=( const _Greater& b ) const
{
  return( false );
}

// -------------------------------------------------------------------------
template< class I, class O >
bool cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::_Greater::
operator==( const _Greater& b ) const
{
  return !( *this != b );
}

// -------------------------------------------------------------------------
template< class I, class O >
typename cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::_Greater::
_T cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::_Greater::
operator()( const _T& a ) const
{
  return( a );
}

// -------------------------------------------------------------------------
template< class I, class O >
typename cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::_Greater::
_T cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::_Greater::
operator()( const _T& a, const _T& b ) const
{
  typedef itk::NumericTraits< _T >     _TTraits;
  typedef typename _TTraits::ValueType _TValue;
  typedef vnl_vector< _TValue >        _TVector;

  _TVector va( _TTraits::GetLength( ) );
  _TVector vb( _TTraits::GetLength( ) );

  _TTraits::AssignToArray( a, va );
  _TTraits::AssignToArray( b, vb );
  return( ( vb.magnitude( ) < va.magnitude( ) )? a: b );
}

// -------------------------------------------------------------------------
template< class I, class O >
void
cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::
AddScale( const double& s )
{
  if( this->m_Scales.insert( s ).second )
    this->Modified( );
}

// -------------------------------------------------------------------------
template< class I, class O >
unsigned long
cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::
GetNumberOfScales( ) const
{
  return( this->m_Scales.size( ) );
}

// -------------------------------------------------------------------------
template< class I, class O >
cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::
MultiScaleGaussianImageFilter( )
  : Superclass( )
{
}

// -------------------------------------------------------------------------
template< class I, class O >
cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::
~MultiScaleGaussianImageFilter( )
{
}

// -------------------------------------------------------------------------
template< class I, class O >
void
cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::
GenerateData( )
{
  typedef itk::GradientRecursiveGaussianImageFilter< I, O > _TGF;
  this->_GenerateData< _TGF >( );
}

// -------------------------------------------------------------------------
template< class I, class O >
template< class F >
void
cpExtensions::Algorithms::
MultiScaleGaussianImageFilter< I, O >::
_GenerateData( )
{
  // Some types
  typedef itk::BinaryFunctorImageFilter< O, O, O, _Greater > _TMaxFilter;
  typedef itk::UnaryFunctorImageFilter< O, O, _Greater >     _TCopyFilter;
  typedef itk::ImageRegionConstIterator< O >                 _TConstIt;
  typedef itk::ImageRegionIterator< O >                      _TIt;
  typedef itk::ProgressAccumulator                           _TProgress;

  // Some values
  typename I::ConstPointer input = this->GetInput( );
  typename O::Pointer output = this->GetOutput( );
  unsigned int nThreads = this->GetNumberOfThreads( );
  float fw = float( 1 ) / float( this->m_Scales.size( ) << 1 );

  // Progress accumulator
  _TProgress::Pointer pg = _TProgress::New( );
  pg->SetMiniPipelineFilter( this );

  // Copy image information
  output->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
  output->SetRequestedRegion( input->GetRequestedRegion( ) );
  output->SetBufferedRegion( input->GetBufferedRegion( ) );
  output->SetSpacing( input->GetSpacing( ) );
  output->SetOrigin( input->GetOrigin( ) );
  output->SetDirection( input->GetDirection( ) );
  output->Allocate( );

  // Intermediary buffer
  typename O::Pointer buffer = O::New( );
  buffer->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
  buffer->SetRequestedRegion( input->GetRequestedRegion( ) );
  buffer->SetBufferedRegion( input->GetBufferedRegion( ) );
  buffer->SetSpacing( input->GetSpacing( ) );
  buffer->SetOrigin( input->GetOrigin( ) );
  buffer->SetDirection( input->GetDirection( ) );
  buffer->Allocate( );

  // Perform all scales
  _TIt gIt( output, output->GetRequestedRegion( ) );
  TScalesContainer::const_iterator sIt = this->m_Scales.begin( );
  for( ; sIt != this->m_Scales.end( ); sIt++ )
  {
    // Single scale filter
    typename F::Pointer filter = F::New( );
    filter->SetInput( input );
    filter->SetNormalizeAcrossScale( true );
    filter->SetNumberOfThreads( nThreads );
    filter->SetSigma( *sIt );
    pg->RegisterInternalFilter( filter, fw );
    filter->GraftOutput( buffer );
    filter->Update( );
    buffer->Graft( filter->GetOutput( ) );

    // Get maximum response
    if( sIt == this->m_Scales.begin( ) )
    {
      // Copy first result
      typename _TCopyFilter::Pointer copy = _TCopyFilter::New( );
      copy->SetInput( buffer );
      copy->InPlaceOff( );
      copy->SetNumberOfThreads( nThreads );
      pg->RegisterInternalFilter( copy, fw );
      copy->GraftOutput( output );
      copy->Update( );
      output->Graft( copy->GetOutput( ) );
    }
    else
    {
      // Maximize results
      typename _TMaxFilter::Pointer max = _TMaxFilter::New( );
      max->SetInput1( output );
      max->SetInput2( buffer );
      max->InPlaceOn( );
      max->SetNumberOfThreads( nThreads );
      pg->RegisterInternalFilter( max, fw );
      max->GraftOutput( output );
      max->Update( );
      output->Graft( max->GetOutput( ) );

    } // fi

  } // rof
}

#endif // __cpExtensions__Algorithms__MultiScaleGaussianImageFilter__hxx__

// eof - $RCSfile$