appli/CTArteries/FourierSeries.h

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