lib/cpExtensions/DataStructures/FourierSeriesContour.h

   1 #ifndef __cpExtensions__DataStructures__FourierSeriesContour__h__
   2 #define __cpExtensions__DataStructures__FourierSeriesContour__h__
   3
   4 #include <cpExtensions/Config.h>
   5 #include <complex>
   6 #include <deque>
   7 #include <itkMatrix.h>
   8 #include <itkPoint.h>
   9
  10 namespace cpExtensions
  11 {
  12   namespace DataStructures
  13   {
  14     /**
  15      * Planar series defined by its harmonic ellipses.
  16      *
  17      * The tuple [z(l),z(-l)], l != 0 defines an ellipse. The z(0)
  18      * harmonic defines the series's center.
  19      *
  20      */
  21     template< class _TScalar >
  22     class cpExtensions_EXPORT FourierSeriesContour
  23       : public std::deque< std::complex< _TScalar > >
  24     {
  25     public:
  26       typedef FourierSeriesContour         Self;
  27       typedef _TScalar                     TScalar;
  28       typedef itk::Matrix< TScalar, 2, 2 > TMatrix;
  29       typedef itk::Point< TScalar, 2 >     TPoint;
  30       typedef typename TPoint::VectorType  TVector;
  31       typedef std::complex< TScalar >      TComplex;
  32       typedef std::deque< TComplex >       Superclass;
  33
  34     public:
  35       FourierSeriesContour( unsigned int q = 1 );
  36
  37       template< class _TScalar2 >
  38       FourierSeriesContour( const FourierSeriesContour< _TScalar2 >& o );
  39
  40       template< class _TIterator >
  41       FourierSeriesContour(
  42         const _TIterator& b, const _TIterator& e, unsigned int q
  43         );
  44
  45       virtual ~FourierSeriesContour( );
  46
  47       template< class _TScalar2 >
  48       Self& operator=( const FourierSeriesContour< _TScalar2 >& o );
  49
  50       bool operator==( const Self& o ) const { return( false ); }
  51       bool operator!=( const Self& o ) const { return( true ); }
  52
  53       /**
  54        * $A=\pi\sum_{l=-q}^{q}l|z_{l}|^{2}$
  55        */
  56       TScalar GetArea( ) const;
  57
  58       /**
  59        * $\kappa(\omega)=\frac{|c'(\omega)\times c''(\omega)|}{|c'(\omega)|^{3}}$
  60        */
  61       TScalar GetCurvature( const TScalar& w ) const;
  62
  63       /**
  64        * Sets the order of the ellipse. This order is defined as the
  65        * counterclock/clock wiseness of the series. Calculations are
  66        * done by computing the series's area (\see Area) and comparing
  67        * its sign.
  68        *
  69        * \param cw counter-clock
  70        */
  71       void SetOrdering( bool cw );
  72       void InvertOrdering( );
  73       void SetOrderingToClockWise( );
  74       void SetOrderingToCounterClockWise( );
  75
  76       /**
  77        * Series manipulators.
  78        */
  79       int GetNumberOfHarmonics( ) const;
  80       void SetNumberOfHarmonics( const int& q );
  81       TComplex& operator[]( const int& l );
  82       const TComplex& operator[]( const int& l ) const;
  83
  84       /**
  85        * $c(\omega)=\sum_{l=-q}^{q}z_{l}e^{jl\omega}$
  86        *
  87        * where $0\le\omega<2\pi$ is the angular curvilinear parameter, q is
  88        * the number of harmonics defining the series, $z_{l}\in\mathbb{C}$
  89        * is the l-th harmonic and $j$ is the imaginary unity.
  90        */
  91       TPoint operator( )( const TScalar& w ) const;
  92
  93       /**
  94        * $\frac{d^{(n)}c(\omega)}{d\omega^{(n)}}=\sum_{l=-q}^{q}z_{l}e^{jl\omega}(jl)^{n}$
  95        *
  96        * where $n$ is the derivative number, $0\le\omega<2\pi$ is the
  97        * angular curvilinear parameter, q is the number of harmonics
  98        * defining the series, $z_{l}\in\mathbb{C}$ is the l-th harmonic
  99        * and $j$ is the imaginary unity.
 100        */
 101       TVector operator( )( const TScalar& w, const unsigned int& n ) const;
 102
 103       /**
 104        * Coefficient-wise addition operators
 105        */
 106       Self operator+( const Self& o ) const;
 107       Self& operator+=( const Self& o );
 108       Self operator-( const Self& o ) const;
 109       Self& operator-=( const Self& o );
 110
 111       /**
 112        * Translation
 113        */
 114       Self operator+( const TComplex& translation ) const;
 115       Self& operator+=( const TComplex& translation );
 116       Self operator-( const TComplex& translation ) const;
 117       Self& operator-=( const TComplex& translation );
 118       Self operator+( const TVector& translation ) const;
 119       Self& operator+=( const TVector& translation );
 120       Self operator-( const TVector& translation ) const;
 121       Self& operator-=( const TVector& translation );
 122       TPoint GetCenter( ) const;
 123
 124       /**
 125        * Scale and/or rotate
 126        */
 127       Self operator*( const TComplex& scale_or_rotation ) const;
 128       Self& operator*=( const TComplex& scale_or_rotation );
 129       Self operator/( const TScalar& scale ) const;
 130       Self& operator/=( const TScalar& scale );
 131
 132       /**
 133        * Phase
 134        */
 135       Self operator&( const TScalar& phase ) const;
 136       Self& operator&=( const TScalar& phase );
 137       TScalar GetPhase( const TScalar& w = TScalar( 0 ) ) const;
 138
 139       /**
 140        * Sampling method
 141        */
 142       void Sample( std::vector< TPoint >& p, const unsigned long& s ) const;
 143
 144       /**
 145        * Ellipse methods
 146        */
 147       void GetEllipse(
 148         int l, TScalar& a, TScalar& b, TScalar& t, TScalar& p
 149         ) const;
 150       void SetEllipse(
 151         int l, TScalar& a, TScalar& b, TScalar& t, TScalar& p
 152         );
 153
 154       template< class _TScalar2 >
 155       void SetEllipses( const std::vector< _TScalar2 >& ellipses );
 156
 157     protected:
 158       // Main method to compute all of series's values.
 159       TComplex _Z( const TScalar& w, const unsigned int& n ) const;
 160
 161       // Fourier transform
 162       void _DFT( const std::vector< TComplex >& p, unsigned int q );
 163
 164     public:
 165       friend std::ostream& operator<<( std::ostream& out, const Self& fs )
 166         {
 167           int q = fs.GetNumberOfHarmonics( );
 168           int l = -q;
 169           out << "{" << l << ":" << fs[ l ] << "}";
 170           for( ++l; l <= q; l++ )
 171             out << ", {" << l << ":" << fs[ l ] << "}";
 172           return( out );
 173         }
 174     };
 175
 176   } // ecapseman
 177
 178 } // ecapseman
 179
 180 // -------------------------------------------------------------------------
 181 template< class _TScalar >
 182 template< class _TScalar2 >
 183 cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 184 FourierSeriesContour( const FourierSeriesContour< _TScalar2 >& o )
 185 {
 186   int q = o.GetNumberOfHarmonics( );
 187   this->SetNumberOfHarmonics( q );
 188   for( int l = -q; l <= q; ++l )
 189     ( *this )[ l ] =
 190       TComplex(
 191         TScalar( std::real( o[ l ] ) ),
 192         TScalar( std::imag( o[ l ] ) )
 193         );
 194 }
 195
 196 // -------------------------------------------------------------------------
 197 template< class _TScalar >
 198 template< class _TIterator >
 199 cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 200 FourierSeriesContour(
 201   const _TIterator& b, const _TIterator& e, unsigned int q
 202   )
 203 {
 204   std::vector< TComplex > aux;
 205   for( _TIterator it = b; it != e; ++it )
 206     aux.push_back(
 207       TComplex(
 208         TScalar( ( *it )[ 0 ] ),
 209         TScalar( ( *it )[ 1 ] )
 210         )
 211       );
 212   this->_DFT( aux, q );
 213 }
 214
 215 // -------------------------------------------------------------------------
 216 template< class _TScalar >
 217 template< class _TScalar2 >
 218 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 219 Self& cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 220 operator=( const FourierSeriesContour< _TScalar2 >& o )
 221 {
 222   int q = o.GetNumberOfHarmonics( );
 223   this->SetNumberOfHarmonics( q );
 224   for( int l = -q; l <= q; ++l )
 225     ( *this )[ l ] =
 226       TComplex(
 227         TScalar( std::real( o[ l ] ) ),
 228         TScalar( std::imag( o[ l ] ) )
 229         );
 230   return( *this );
 231 }
 232
 233 // -------------------------------------------------------------------------
 234 template< class _TScalar >
 235 template< class _TScalar2 >
 236 void cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 237 SetEllipses( const std::vector< _TScalar2 >& ellipses )
 238 {
 239   int Q = ( ellipses.size( ) - 2 ) >> 2;
 240   Q = ( Q < 1 )? 1: Q;
 241   this->SetNumberOfHarmonics( Q );
 242   auto cIt = ellipses.begin( );
 243   TScalar cx = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
 244   if( cIt != ellipses.end( ) )
 245     ++cIt;
 246   TScalar cy = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
 247   if( cIt != ellipses.end( ) )
 248     ++cIt;
 249   ( *this )[ 0 ] = TComplex( cx, cy );
 250   for( int l = 1; l <= Q; ++l )
 251   {
 252     TScalar a = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
 253     if( cIt != ellipses.end( ) )
 254       ++cIt;
 255     TScalar b = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
 256     if( cIt != ellipses.end( ) )
 257       ++cIt;
 258     TScalar t = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
 259     if( cIt != ellipses.end( ) )
 260       ++cIt;
 261     TScalar p = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
 262     if( cIt != ellipses.end( ) )
 263       ++cIt;
 264     this->SetEllipse( l, a, b, t, p );
 265
 266   } // rof
 267 }
 268
 269 #endif // __cpExtensions__DataStructures__FourierSeriesContour__h__
 270
 271 // eof - $RCSfile$