#ifndef __FPA__IMAGE__ALGORITHM__HXX__
#define __FPA__IMAGE__ALGORITHM__HXX__

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

// -------------------------------------------------------------------------
template< class I, class O, class A >
typename fpa::Image::Algorithm< I, O, A >::
TMinimumSpanningTree* fpa::Image::Algorithm< I, O, A >::
GetMinimumSpanningTree( )
{
  return(
    dynamic_cast< TMinimumSpanningTree* >(
      this->itk::ProcessObject::GetOutput( 1 )
      )
    );
}

// -------------------------------------------------------------------------
template< class I, class O, class A >
const typename fpa::Image::Algorithm< I, O, A >::
TMinimumSpanningTree* fpa::Image::Algorithm< I, O, A >::
GetMinimumSpanningTree( ) const
{
  return(
    dynamic_cast< const TMinimumSpanningTree* >(
      this->itk::ProcessObject::GetOutput( 1 )
      )
    );
}

// -------------------------------------------------------------------------
template< class I, class O, class A >
void fpa::Image::Algorithm< I, O, A >::
GraftMinimumSpanningTree( itk::DataObject* obj )
{
  TMinimumSpanningTree* mst = dynamic_cast< TMinimumSpanningTree* >( obj );
  if( mst != NULL )
    this->GraftNthOutput( 1, mst );
}

// -------------------------------------------------------------------------
template< class I, class O, class A >
fpa::Image::Algorithm< I, O, A >::
Algorithm( )
  : Superclass( ),
    m_NeighborhoodOrder( 1 )
{
  this->itk::ProcessObject::SetNumberOfRequiredOutputs( 2 );
  this->itk::ProcessObject::SetNthOutput( 0, O::New( ) );
  this->itk::ProcessObject::SetNthOutput( 1, TMinimumSpanningTree::New( ) );
}

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

// -------------------------------------------------------------------------
template< class I, class O, class A >
void fpa::Image::Algorithm< I, O, A >::
_BeforeGenerateData( )
{
  this->Superclass::_BeforeGenerateData( );
  this->AllocateOutputs( );
}

// -------------------------------------------------------------------------
template< class I, class O, class A >
void fpa::Image::Algorithm< I, O, A >::
_AfterGenerateData( )
{
  this->Superclass::_AfterGenerateData( );
  this->GetMinimumSpanningTree( )->SetCollisions( this->m_Collisions );
}

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

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

// -------------------------------------------------------------------------
template< class I, class O, class A >
double fpa::Image::Algorithm< I, O, A >::
_Distance( 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 O, class A >
bool fpa::Image::Algorithm< I, O, A >::
_HasEdge( 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 O, class A >
void fpa::Image::Algorithm< I, O, A >::
_Neighborhood( _TVertices& neighborhood, const TVertex& v ) const
{
  typename I::RegionType reg = this->GetInput( )->GetRequestedRegion( );

  neighborhood.clear( );
  if( this->m_NeighborhoodOrder == 1 )
  {
    for( unsigned int d = 0; d < I::ImageDimension; d++ )
    {
      for( int i = -1; i <= 1; i += 2 )
      {
        TVertex n = v;
        n[ d ] += i;
        if( reg.IsInside( n ) )
          neighborhood.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( v );
    for( unsigned int i = 0; i < nIt.Size( ); i++ )
    {
      TVertex n = nIt.GetIndex( i );
      if( n == v )
        continue;
      if( reg.IsInside( n ) )
        neighborhood.push_back( n );

    } // rof

  } // fi
}

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

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

// -------------------------------------------------------------------------
template< class I, class O, class A >
void fpa::Image::Algorithm< I, O, A >::
_SetResult( const TVertex& v, const TResult& r )
{
  this->GetOutput( )->SetPixel( v, r );
}

// -------------------------------------------------------------------------
template< class I, class O, class A >
const typename fpa::Image::Algorithm< I, O, A >::
_TNode& fpa::Image::Algorithm< I, O, A >::
_Node( const TVertex& v ) const
{
  return( this->GetMinimumSpanningTree( )->GetPixel( v ) );
}

// -------------------------------------------------------------------------
template< class I, class O, class A >
void fpa::Image::Algorithm< I, O, A >::
_InitMarks( )
{
  _TNode far_node;
  far_node.Label = Self::FarLabel;
  this->GetMinimumSpanningTree( )->FillBuffer( far_node );
}

// -------------------------------------------------------------------------
template< class I, class O, class A >
void fpa::Image::Algorithm< I, O, A >::
_Mark( const _TNode& node )
{
  this->GetMinimumSpanningTree( )->SetPixel( node.Vertex, node );
}

#endif // __FPA__IMAGE__ALGORITHM__HXX__

// eof - $RCSfile$