/* =======================================================================
 * @author: Leonardo Florez-Valencia
 * @email: florez-l@javeriana.edu.co
 * =======================================================================
 */
#ifndef __FourierSeries__h__
#define __FourierSeries__h__

#include <complex>
#include <deque>
#include <vector>
#include <utility>

#include <itkMatrix.h>
#include <itkPoint.h>
#include <itkVector.h>

/**
 * 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 S, unsigned int D >
class FourierSeries
  : public std::deque< std::complex< S > >
{
public:
  typedef FourierSeries Self;
  typedef S TScalar;
  itkStaticConstMacro( Dimension, unsigned int, D );

  typedef itk::Matrix< S, D, D >      TMatrix;
  typedef itk::Point< S, D >          TPoint;
  typedef typename TPoint::VectorType TVector;
  typedef std::complex< S >           TComplex;
  typedef std::deque< TComplex >      Superclass;

public:
  FourierSeries( unsigned int q = 1 );
  template< class S2, unsigned int D2 >
  FourierSeries( const FourierSeries< S2, D2 >& o );
  template< class _TIt >
  FourierSeries( const _TIt& b, const _TIt& e, unsigned int q );
  virtual ~FourierSeries( );

  template< class S2, unsigned int D2 >
  Self& operator=( const FourierSeries< S2, D2 >& o );

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

  const unsigned int& GetRealAxis( ) const;
  const unsigned int& GetImagAxis( ) const;
  void SetRealAxis( const unsigned int& a );
  void SetImagAxis( const unsigned int& a );

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

  /**
   * $\kappa(\omega)=\frac{|c'(\omega)\times c''(\omega)|}{|c'(\omega)|^{3}}$
   */
  S GetCurvature( const S& 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 S& 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 S& 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 S& scale ) const;
  Self& operator/=( const S& scale );

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

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

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

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

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

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

protected:
  unsigned int m_RealAxis;
  unsigned int m_ImagAxis;

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 );
    }
};

#ifndef ITK_MANUAL_INSTANTIATION
#include "FourierSeries.hxx"
#endif // ITK_MANUAL_INSTANTIATION
#endif // __FourierSeries__h__

// eof - $RCSfile$