#ifndef __FPA__IMAGE__ALGORITHM__HXX__
#define __FPA__IMAGE__ALGORITHM__HXX__

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

// -------------------------------------------------------------------------
template< class I, class A, class CC >
fpa::Image::Algorithm< I, A, CC >::
Algorithm( )
  : Superclass( ),
    m_NeighborhoodOrder( 1 )
{
  this->m_CostConversion = TCostConversionFunction::New( );
}

// -------------------------------------------------------------------------
template< class I, class A, class CC >
fpa::Image::Algorithm< I, A, CC >::
~Algorithm( )
{
}

// -------------------------------------------------------------------------
template< class I, class A, class CC >
bool fpa::Image::Algorithm< I, A, CC >::
_UpdateResult( _TNode& n )
{
  bool ret = this->Superclass::_UpdateResult( n );
  this->GetOutput( )->SetPixel( n.Vertex, n.Result );
  return( ret );
}

// -------------------------------------------------------------------------
template< class I, class A, class CC >
unsigned long fpa::Image::Algorithm< I, A, CC >::
_NumberOfVertices( ) const
{
  return(
    this->GetInput( )->GetLargestPossibleRegion( ).GetNumberOfPixels( )
    );
}

// -------------------------------------------------------------------------
template< class I, class A, class CC >
typename fpa::Image::Algorithm< I, A, CC >::
TVertexValue fpa::Image::Algorithm< I, A, CC >::
_Value( const TVertex& v ) const
{
  return( this->GetInput( )->GetPixel( v ) );
}

// -------------------------------------------------------------------------
template< class I, class A, class CC >
typename fpa::Image::Algorithm< I, A, CC >::
TResult fpa::Image::Algorithm< I, A, CC >::
_Result( const TVertex& v ) const
{
  return( this->GetOutput( )->GetPixel( v ) );
}

// -------------------------------------------------------------------------
template< class I, class A, class CC >
double fpa::Image::Algorithm< I, A, CC >::
_Norm( const TVertex& a, const TVertex& b ) const
{
  typename I::PointType pa, pb;
  this->GetInput( )->TransformIndexToPhysicalPoint( a, pa );
  this->GetInput( )->TransformIndexToPhysicalPoint( b, pb );
  return( double( pa.EuclideanDistanceTo( pb ) ) );
}

// -------------------------------------------------------------------------
template< class I, class A, class CC >
bool fpa::Image::Algorithm< I, A, CC >::
_Edge( const TVertex& a, const TVertex& b ) const
{
  unsigned long dist = 0;
  for( unsigned int d = 0; d < I::ImageDimension; d++ )
    dist += std::abs( long( a[ d ] ) - long( b[ d ] ) );
  if( this->m_NeighborhoodOrder == 1 )
    return( dist <= 1 );
  else
    return( dist <= I::ImageDimension );
}

// -------------------------------------------------------------------------
template< class I, class A, class CC >
typename fpa::Image::Algorithm< I, A, CC >::
TCost fpa::Image::Algorithm< I, A, CC >::
_Cost( const TVertex& a, const TVertex& b ) const
{
  static const TCost INF_COST = std::numeric_limits< TCost >::max( );
  if( this->_Edge( a, b ) )
  {
    return(
      this->m_CostConversion->Evaluate( this->GetInput( )->GetPixel( b ) )
      );
  }
  else
    return( INF_COST );
}

// -------------------------------------------------------------------------
template< class I, class A, class CC >
void fpa::Image::Algorithm< I, A, CC >::
_Neighs( const _TNode& n, _TNodes& N ) const
{
  typename I::RegionType reg = this->GetInput( )->GetRequestedRegion( );

  N.clear( );
  if( this->m_NeighborhoodOrder == 1 )
  {
    for( unsigned int d = 0; d < I::ImageDimension; d++ )
    {
      for( int i = -1; i <= 1; i += 2 )
      {
        TVertex v = n.Vertex;
        v[ d ] += i;
        if( reg.IsInside( v ) )
          N.push_back( _TNode( v, n.FrontId ) );
        else
          N.push_back( n );

      } // rof

    } // rof
  }
  else
  {
    typedef itk::ConstNeighborhoodIterator< I > TNeighIt;
    typename I::SizeType nSize;
    nSize.Fill( 1 );

    TNeighIt nIt( nSize, this->GetInput( ), reg );
    nIt.SetLocation( n.Vertex );
    for( unsigned int i = 0; i < nIt.Size( ); i++ )
    {
      TVertex idxN = nIt.GetIndex( i );
      if( idxN == n.Vertex )
        continue;
      if( !reg.IsInside( idxN ) )
        continue;
      N.push_back( _TNode( idxN, n.FrontId ) );

    } // rof

  } // fi
}

// -------------------------------------------------------------------------
template< class I, class A, class CC >
void fpa::Image::Algorithm< I, A, CC >::
_NeighsInDim( const _TNode& n, const unsigned int& d, _TNodes& N )
{
  typename I::RegionType reg = this->GetInput( )->GetRequestedRegion( );

  N.clear( );
  for( int i = -1; i <= 1; i += 2 )
  {
    TVertex v = n.Vertex;
    v[ d ] += i;
    if( reg.IsInside( v ) )
      N.push_back( _TNode( v, n.FrontId ) );

  } // rof
}

// -------------------------------------------------------------------------
template< class I, class A, class CC >
void fpa::Image::Algorithm< I, A, CC >::
_InitializeResults( )
{
}

#endif // __FPA__IMAGE__ALGORITHM__HXX__

// eof - $RCSfile$