#ifndef __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
#define __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__

#include <cmath>
#include <vnl/vnl_math.h>
#include <itkLineConstIterator.h>

// -------------------------------------------------------------------------
template< class _TGradient >
cpExtensions::Algorithms::GulsunTekMedialness< _TGradient >::
GulsunTekMedialness( )
  : Superclass( ),
    m_MinRadius( double( 0 ) ),
    m_MaxRadius( double( 1 ) ),
    m_ProfileSampling( 4 ),
    m_RadialSampling( 10 )
{
}

// -------------------------------------------------------------------------
template< class _TGradient >
cpExtensions::Algorithms::GulsunTekMedialness< _TGradient >::
~GulsunTekMedialness( )
{
}

// -------------------------------------------------------------------------
template< class _TGradient >
typename cpExtensions::Algorithms::GulsunTekMedialness< _TGradient >::
TOutput cpExtensions::Algorithms::GulsunTekMedialness< _TGradient >::
_Evaluate( const TIndex& i ) const
{
  itk::Object::GlobalWarningDisplayOff( );

  // Various values
  const _TGradient* in = this->GetInputImage( );
  double pi2n =
    double( 2 ) * double( vnl_math::pi ) /
    double( this->m_ProfileSampling );
  double rOff = this->m_MaxRadius / double( this->m_RadialSampling - 1 );
  double optR = double( 0 );
  TPoint pnt;
  in->TransformIndexToPhysicalPoint( i, pnt );

  // Main loop
  for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ )
  {
    for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ )
    {
      TProfile maxProfile( this->m_RadialSampling, double( 0 ) );
      for( unsigned int p = 0; p < this->m_ProfileSampling; p++ )
      {
        double a = pi2n * double( p );

        // Direction of this profile
        TVector dir;
        dir.Fill( TScalar( 0 ) );
        dir[ cx ] = TScalar( std::cos( a ) );
        dir[ cy ] = TScalar( std::sin( a ) );

        double maxrise = double( 0 );
        double maxfall = double( -1 );
        TProfile profile;
        for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
        {
          double radius = double( r ) * rOff;
          TIndex idx;
          typename TPoint::VectorType aux;
          aux.SetVnlVector( dir.GetVnlVector( ) );
          if(
            in->TransformPhysicalPointToIndex( pnt + ( aux * radius ), idx )
            )
          {
            TVector g = in->GetPixel( idx );
            double b = double( g.GetNorm( ) );
            if( double( g * dir ) < double( 0 ) )
              b *= double( -1 );
            maxrise = ( b > maxrise )? b: maxrise;
            if( radius >= this->m_MinRadius )
              maxfall = ( b < maxfall )? b: maxfall;
            profile.push_back( -b - maxrise );
          }
          else
            profile.push_back( double( 0 ) );

        } // rof

        for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
        {
          double E = profile[ r ] / -maxfall;
          E = ( E < double( 0 ) )? double( 0 ): E;
          E = ( E > double( 1 ) )? double( 1 ): E;
          maxProfile[ r ] += E;

        } // rof

      } // rof

      for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
      {
        double E = maxProfile[ r ] / double( this->m_RadialSampling );
        optR = ( E > optR )? E: optR;

      } // rof

    } // rof

  } // rof
  return( TScalar( optR ) );
}

#endif // __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__

// eof - $RCSfile$