#ifndef __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__HXX__
#define __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__HXX__

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

// -------------------------------------------------------------------------
template< class _TGradient, class _TMask >
cpExtensions::Algorithms::FluxMedialness< _TGradient, _TMask >::
FluxMedialness( )
  : Superclass( ),
    m_MinRadius( double( 0 ) ),
    m_MaxRadius( double( 1 ) ),
    m_RadialSampling( 4 ),
    m_RadiusStep( double( 1 ) )
{
}

// -------------------------------------------------------------------------
template< class _TGradient, class _TMask >
cpExtensions::Algorithms::FluxMedialness< _TGradient, _TMask >::
~FluxMedialness( )
{
}

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

  double pi2n = double( 2 ) * double( vnl_math::pi );
  pi2n /= int( this->m_RadialSampling );
  const _TGradient* img = this->GetInputImage( );
  // Gradient in central pixel
  //TVector grad_idx = img->GetPixel( i );

  double Flux1 = 0;
  double Flux = 0;

  TRCandidates FluxFinal;
  TRCandidates radiusGenerated;
  double dR = double( 0 );
  double optR = double( 0 );
  TPoint center;
  img->TransformIndexToPhysicalPoint( i, center );
  double radius = double(0);

  for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ )
  {
    for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ )
    {
      dR = double( 0 );
      FluxFinal.clear();
      radiusGenerated.clear();
      radius = this->m_MinRadius;
      while( radius <= this->m_MaxRadius )
      {
        Flux = 0;
        for( unsigned int I_radial = 0; I_radial < this->m_RadialSampling ; I_radial++ )
        {
          Flux1 = 0;

          // Direction of first profile
          typename TPoint::VectorType dir1;
          dir1.Fill( double( 0 ) );
          dir1[ cx ] = std::cos( pi2n * double( I_radial ) );
          dir1[ cy ] = std::sin( pi2n * double( I_radial ) );
          dir1 *= (radius);
          TIndex rIdx;
          if (img->TransformPhysicalPointToIndex( center + dir1, rIdx ))
          {
            TVector grad_rIdx = img->GetPixel( rIdx );
            TVector u_i1;
            u_i1.SetVnlVector( ( center - ( center + dir1 ) ).GetVnlVector( ) );
            u_i1.Normalize( );
            // dot product
            Flux1 = grad_rIdx * u_i1;
          }

          Flux += Flux1;

        } // rof
        //std::cout<<" radius:"<<radius<<std::endl;
        //std::cout<<"Center:"<<center[0]<<" "<<center[1]<<std::endl;
        //std::cout<<"edge:"<<center[0]+radius*std::cos( pi2n * double( 0 ) )<<std::endl;
        Flux /= this->m_RadialSampling;
        FluxFinal.push_back(Flux);
        radiusGenerated.push_back(radius);
        radius += this->m_RadiusStep;

      }     //elihw

      dR= *( std::max_element( FluxFinal.begin(), FluxFinal.end() ) );
      optR= (dR>optR)? dR:optR;

    } // rof

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

#endif // __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__HXX__

// eof - $RCSfile$