lib/ivq/ITK/FourierSeries.h

   1 // =========================================================================
   2 // @author: Leonardo Florez-Valencia
   3 // @email: florez-l@javeriana.edu.co
   4 // =========================================================================
   5 #ifndef __ivq__ITK__FourierSeries__h__
   6 #define __ivq__ITK__FourierSeries__h__
   7
   8 #include <complex>
   9 #include <deque>
  10 #include <vector>
  11 #include <utility>
  12
  13 #include <itkMatrix.h>
  14 #include <itkPoint.h>
  15 #include <itkVector.h>
  16
  17 namespace ivq
  18 {
  19   namespace ITK
  20   {
  21     /**
  22      * Planar series defined by its harmonic ellipses.
  23      *
  24      * The tuple [z(l),z(-l)], l != 0 defines an ellipse. The z(0)
  25      * harmonic defines the series's center.
  26      *
  27      */
  28     template< class _TScalar, unsigned int _VDim >
  29     class FourierSeries
  30       : public std::deque< std::complex< _TScalar > >
  31     {
  32     public:
  33       typedef FourierSeries Self;
  34       typedef _TScalar      TScalar;
  35       itkStaticConstMacro( Dimension, unsigned int, _VDim );
  36
  37       typedef itk::Matrix< _TScalar, _VDim, _VDim > TMatrix;
  38       typedef itk::Point< _TScalar, _VDim >         TPoint;
  39       typedef typename TPoint::VectorType           TVector;
  40       typedef std::complex< _TScalar >              TComplex;
  41       typedef std::deque< TComplex >                Superclass;
  42
  43     public:
  44       FourierSeries( unsigned int q = 1 );
  45       template< class _TScalar2, unsigned int _VDim2 >
  46       FourierSeries( const FourierSeries< _TScalar2, _VDim2 >& o );
  47       template< class _TIt >
  48       FourierSeries( const _TIt& b, const _TIt& e, unsigned int q );
  49       virtual ~FourierSeries( );
  50
  51       template< class _TScalar2, unsigned int _VDim2 >
  52       Self& operator=( const FourierSeries< _TScalar2, _VDim2 >& o );
  53
  54       bool operator==( const Self& o ) const { return( false ); }
  55       bool operator!=( const Self& o ) const { return( true ); }
  56
  57       const unsigned int& GetRealAxis( ) const;
  58       const unsigned int& GetImagAxis( ) const;
  59       void SetRealAxis( const unsigned int& a );
  60       void SetImagAxis( const unsigned int& a );
  61
  62       /**
  63        * $A=\pi\sum_{l=-q}^{q}l|z_{l}|^{2}$
  64        */
  65       TScalar GetArea( ) const;
  66
  67       /**
  68        * $\kappa(\omega)=\frac{|c'(\omega)\times c''(\omega)|}{|c'(\omega)|^{3}}$
  69        */
  70       TScalar GetCurvature( const TScalar& w ) const;
  71
  72       /**
  73        * Sets the order of the ellipse. This order is defined as the
  74        * counterclock/clock wiseness of the series. Calculations are
  75        * done by computing the series's area (\see Area) and comparing
  76        * its sign.
  77        *
  78        * \param cw counter-clock
  79        */
  80       void SetOrdering( bool cw );
  81       void InvertOrdering( );
  82       void SetOrderingToClockWise( );
  83       void SetOrderingToCounterClockWise( );
  84
  85       /**
  86        * Series manipulators.
  87        */
  88       int GetNumberOfHarmonics( ) const;
  89       void SetNumberOfHarmonics( const int& q );
  90       TComplex& operator[]( const int& l );
  91       const TComplex& operator[]( const int& l ) const;
  92
  93       /**
  94        * $c(\omega)=\sum_{l=-q}^{q}z_{l}e^{jl\omega}$
  95        *
  96        * where $0\le\omega<2\pi$ is the angular curvilinear parameter, q is
  97        * the number of harmonics defining the series, $z_{l}\in\mathbb{C}$
  98        * is the l-th harmonic and $j$ is the imaginary unity.
  99        */
 100       TPoint operator( )( const TScalar& w ) const;
 101
 102       /**
 103        * $\frac{d^{(n)}c(\omega)}{d\omega^{(n)}}=\sum_{l=-q}^{q}z_{l}e^{jl\omega}(jl)^{n}$
 104        *
 105        * where $n$ is the derivative number, $0\le\omega<2\pi$ is the
 106        * angular curvilinear parameter, q is the number of harmonics
 107        * defining the series, $z_{l}\in\mathbb{C}$ is the l-th harmonic
 108        * and $j$ is the imaginary unity.
 109        */
 110       TVector operator( )( const TScalar& w, const unsigned int& n ) const;
 111
 112       /**
 113        * Coefficient-wise addition operators
 114        */
 115       Self operator+( const Self& o ) const;
 116       Self& operator+=( const Self& o );
 117       Self operator-( const Self& o ) const;
 118       Self& operator-=( const Self& o );
 119
 120       /**
 121        * Translation
 122        */
 123       Self operator+( const TComplex& translation ) const;
 124       Self& operator+=( const TComplex& translation );
 125       Self operator-( const TComplex& translation ) const;
 126       Self& operator-=( const TComplex& translation );
 127       Self operator+( const TVector& translation ) const;
 128       Self& operator+=( const TVector& translation );
 129       Self operator-( const TVector& translation ) const;
 130       Self& operator-=( const TVector& translation );
 131       TPoint GetCenter( ) const;
 132
 133       /**
 134        * Scale and/or rotate
 135        */
 136       Self operator*( const TComplex& scale_or_rotation ) const;
 137       Self& operator*=( const TComplex& scale_or_rotation );
 138       Self operator/( const TScalar& scale ) const;
 139       Self& operator/=( const TScalar& scale );
 140
 141       /**
 142        * Phase
 143        */
 144       Self operator&( const TScalar& phase ) const;
 145       Self& operator&=( const TScalar& phase );
 146       TScalar GetPhase( const TScalar& w = TScalar( 0 ) ) const;
 147
 148       /**
 149        * Sampling method
 150        */
 151       void Sample( std::vector< TPoint >& p, const unsigned long& s ) const;
 152
 153       /**
 154        * Ellipse methods
 155        */
 156       void GetEllipse( int l, TScalar& a, TScalar& b, TScalar& t, TScalar& p ) const;
 157       void SetEllipse( int l, TScalar& a, TScalar& b, TScalar& t, TScalar& p );
 158
 159       template< class _TOtherScalar >
 160       void SetEllipses( const std::vector< _TOtherScalar >& ellipses );
 161
 162     protected:
 163       // Main method to compute all of series's values.
 164       TComplex _Z( const TScalar& w, const unsigned int& n ) const;
 165
 166       // Fourier transform
 167       void _DFT( const std::vector< TComplex >& p, unsigned int q );
 168
 169     protected:
 170       unsigned int m_RealAxis;
 171       unsigned int m_ImagAxis;
 172
 173     public:
 174       friend std::ostream& operator<<( std::ostream& out, const Self& fs )
 175         {
 176           int q = fs.GetNumberOfHarmonics( );
 177           int l = -q;
 178           out << "{" << l << ":" << fs[ l ] << "}";
 179           for( ++l; l <= q; l++ )
 180             out << ", {" << l << ":" << fs[ l ] << "}";
 181           return( out );
 182         }
 183     };
 184
 185   } // ecapseman
 186
 187 } // ecapseman
 188
 189 #include <ivq/ITK/FourierSeries.hxx>
 190 #endif // __ivq__ITK__FourierSeries__h__
 191
 192 // eof - $RCSfile$