#ifndef __cpExtensions__DataStructures__FourierSeriesContour__h__
#define __cpExtensions__DataStructures__FourierSeriesContour__h__

#include <cpExtensions/Config.h>
#include <complex>
#include <deque>
#include <itkMatrix.h>
#include <itkPoint.h>

namespace cpExtensions
{
  namespace DataStructures
  {
    /**
     * Planar series defined by its harmonic ellipses.
     *
     * The tuple [z(l),z(-l)], l != 0 defines an ellipse. The z(0)
     * harmonic defines the series's center.
     *
     */
    template< class _TScalar >
    class cpExtensions_EXPORT FourierSeries
      : public std::deque< std::complex< _TScalar > >
    {
    public:
      typedef FourierSeries               Self;
      typedef _TScalar                    TScalar;
      typedef itk::Matrix< S, 2, 2 >      TMatrix;
      typedef itk::Point< S, 2 >          TPoint;
      typedef typename TPoint::VectorType TVector;
      typedef std::complex< S >           TComplex;
      typedef std::deque< TComplex >      Superclass;

    public:
      FourierSeries( unsigned int q = 1 );

      template< class _TScalar2 >
      FourierSeries( const FourierSeries< _TScalar2 >& o );

      template< class _TIterator >
      FourierSeries(
        const _TIterator& b, const _TIterator& e, unsigned int q
        );

      virtual ~FourierSeries( );

      template< class _TScalar2 >
      Self& operator=( const FourierSeries< _TSCalar2 >& o );

      bool operator==( const Self& o ) const { return( false ); }
      bool operator!=( const Self& o ) const { return( true ); }

      /**
       * $A=\pi\sum_{l=-q}^{q}l|z_{l}|^{2}$
       */
      TScalar GetArea( ) const;

      /**
       * $\kappa(\omega)=\frac{|c'(\omega)\times c''(\omega)|}{|c'(\omega)|^{3}}$
       */
      TScalar GetCurvature( const TScalar& w ) const;

      /**
       * Sets the order of the ellipse. This order is defined as the
       * counterclock/clock wiseness of the series. Calculations are
       * done by computing the series's area (\see Area) and comparing
       * its sign.
       *
       * \param cw counter-clock
       */
      void SetOrdering( bool cw );
      void InvertOrdering( );
      void SetOrderingToClockWise( );
      void SetOrderingToCounterClockWise( );

      /**
       * Series manipulators.
       */
      int GetNumberOfHarmonics( ) const;
      void SetNumberOfHarmonics( const int& q );
      TComplex& operator[]( const int& l );
      const TComplex& operator[]( const int& l ) const;

      /**
       * $c(\omega)=\sum_{l=-q}^{q}z_{l}e^{jl\omega}$
       *
       * where $0\le\omega<2\pi$ is the angular curvilinear parameter, q is
       * the number of harmonics defining the series, $z_{l}\in\mathbb{C}$
       * is the l-th harmonic and $j$ is the imaginary unity.
       */
      TPoint operator( )( const TScalar& w ) const;

      /**
       * $\frac{d^{(n)}c(\omega)}{d\omega^{(n)}}=\sum_{l=-q}^{q}z_{l}e^{jl\omega}(jl)^{n}$
       *
       * where $n$ is the derivative number, $0\le\omega<2\pi$ is the
       * angular curvilinear parameter, q is the number of harmonics
       * defining the series, $z_{l}\in\mathbb{C}$ is the l-th harmonic
       * and $j$ is the imaginary unity.
       */
      TVector operator( )( const TScalar& w, const unsigned int& n ) const;

      /**
       * Coefficient-wise addition operators
       */
      Self operator+( const Self& o ) const;
      Self& operator+=( const Self& o );
      Self operator-( const Self& o ) const;
      Self& operator-=( const Self& o );

      /**
       * Translation
       */
      Self operator+( const TComplex& translation ) const;
      Self& operator+=( const TComplex& translation );
      Self operator-( const TComplex& translation ) const;
      Self& operator-=( const TComplex& translation );
      Self operator+( const TVector& translation ) const;
      Self& operator+=( const TVector& translation );
      Self operator-( const TVector& translation ) const;
      Self& operator-=( const TVector& translation );
      TPoint GetCenter( ) const;

      /**
       * Scale and/or rotate
       */
      Self operator*( const TComplex& scale_or_rotation ) const;
      Self& operator*=( const TComplex& scale_or_rotation );
      Self operator/( const TScalar& scale ) const;
      Self& operator/=( const TScalar& scale );

      /**
       * Phase
       */
      Self operator&( const TScalar& phase ) const;
      Self& operator&=( const TScalar& phase );
      TScalar GetPhase( const TScalar& w = TScalar( 0 ) ) const;

      /**
       * Sampling method
       */
      void Sample( std::vector< TPoint >& p, const unsigned long& s ) const;

      /**
       * Ellipse methods
       */
      void GetEllipse(
        int l, TScalar& a, TScalar& b, TScalar& t, TScalar& p
        ) const;
      void SetEllipse(
        int l, TScalar& a, TScalar& b, TScalar& t, TScalar& p
        );

      template< class _TScalar2 >
      void SetEllipses( const std::vector< _TScalar2 >& ellipses );

    protected:
      // Main method to compute all of series's values.
      TComplex _Z( const TScalar& w, const unsigned int& n ) const;

      // Fourier transform
      void _DFT( const std::vector< TComplex >& p, unsigned int q );

    public:
      friend std::ostream& operator<<( std::ostream& out, const Self& fs )
        {
          int q = fs.GetNumberOfHarmonics( );
          int l = -q;
          out << "{" << l << ":" << fs[ l ] << "}";
          for( ++l; l <= q; l++ )
            out << ", {" << l << ":" << fs[ l ] << "}";
          return( out );
        }
    };

  } // ecapseman

} // ecapseman

// -------------------------------------------------------------------------
template< class _TScalar >
template< class _TScalar2 >
cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
FourierSeriesContour( const FourierSeriesContour< _TScalar2 >& o )
{
  int q = o.GetNumberOfHarmonics( );
  this->SetNumberOfHarmonics( q );
  for( int l = -q; l <= q; ++l )
    ( *this )[ l ] =
      TComplex(
        TScalar( std::real( o[ l ] ) ),
        TScalar( std::imag( o[ l ] ) )
        );
}

// -------------------------------------------------------------------------
template< class _TScalar >
template< class _TIterator >
cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
FourierSeriesContour(
  const _TIterator& b, const _TIterator& e, unsigned int q
  )
{
  std::vector< TComplex > aux;
  for( _TIterator it = b; it != e; ++it )
    aux.push_back(
      TComplex(
        TScalar( ( *it )[ 0 ] ),
        TScalar( ( *it )[ 1 ] )
        )
      );
  this->_DFT( aux, q );
}

// -------------------------------------------------------------------------
template< class _TScalar >
template< class _TScalar2 >
typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
Self& cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
operator=( const FourierSeriesContour< _TScalar2 >& o )
{
  int q = o.GetNumberOfHarmonics( );
  this->SetNumberOfHarmonics( q );
  for( int l = -q; l <= q; ++l )
    ( *this )[ l ] =
      TComplex(
        TScalar( std::real( o[ l ] ) ),
        TScalar( std::imag( o[ l ] ) )
        );
  return( *this );
}

// -------------------------------------------------------------------------
template< class _TScalar >
template< class _TScalar2 >
void cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
SetEllipses( const std::vector< _TScalar2 >& ellipses )
{
  int Q = ( ellipses.size( ) - 2 ) >> 2;
  Q = ( Q < 1 )? 1: Q;
  this->SetNumberOfHarmonics( Q );
  auto cIt = ellipses.begin( );
  TScalar cx = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
  if( cIt != ellipses.end( ) )
    ++cIt;
  TScalar cy = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
  if( cIt != ellipses.end( ) )
    ++cIt;
  ( *this )[ 0 ] = TComplex( cx, cy );
  for( int l = 1; l <= Q; ++l )
  {
    TScalar a = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
    if( cIt != ellipses.end( ) )
      ++cIt;
    TScalar b = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
    if( cIt != ellipses.end( ) )
      ++cIt;
    TScalar t = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
    if( cIt != ellipses.end( ) )
      ++cIt;
    TScalar p = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
    if( cIt != ellipses.end( ) )
      ++cIt;
    this->SetEllipse( l, a, b, t, p );

  } // rof
}

#endif // __cpExtensions__DataStructures__FourierSeriesContour__h__

// eof - $RCSfile$