1 #ifndef __cpExtensions__DataStructures__FourierSeriesContour__h__
2 #define __cpExtensions__DataStructures__FourierSeriesContour__h__
4 #include <cpExtensions/Config.h>
10 namespace cpExtensions
12 namespace DataStructures
15 * Planar series defined by its harmonic ellipses.
17 * The tuple [z(l),z(-l)], l != 0 defines an ellipse. The z(0)
18 * harmonic defines the series's center.
21 template< class _TScalar >
22 class cpExtensions_EXPORT FourierSeriesContour
23 : public std::deque< std::complex< _TScalar > >
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;
35 FourierSeriesContour( unsigned int q = 1 );
37 template< class _TScalar2 >
38 FourierSeriesContour( const FourierSeriesContour< _TScalar2 >& o );
40 template< class _TIterator >
42 const _TIterator& b, const _TIterator& e, unsigned int q
45 virtual ~FourierSeriesContour( );
47 template< class _TScalar2 >
48 Self& operator=( const FourierSeriesContour< _TScalar2 >& o );
50 bool operator==( const Self& o ) const { return( false ); }
51 bool operator!=( const Self& o ) const { return( true ); }
54 * $A=\pi\sum_{l=-q}^{q}l|z_{l}|^{2}$
56 TScalar GetArea( ) const;
59 * $\kappa(\omega)=\frac{|c'(\omega)\times c''(\omega)|}{|c'(\omega)|^{3}}$
61 TScalar GetCurvature( const TScalar& w ) const;
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
69 * \param cw counter-clock
71 void SetOrdering( bool cw );
72 void InvertOrdering( );
73 void SetOrderingToClockWise( );
74 void SetOrderingToCounterClockWise( );
77 * Series manipulators.
79 int GetNumberOfHarmonics( ) const;
80 void SetNumberOfHarmonics( const int& q );
81 TComplex& operator[]( const int& l );
82 const TComplex& operator[]( const int& l ) const;
85 * $c(\omega)=\sum_{l=-q}^{q}z_{l}e^{jl\omega}$
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.
91 TPoint operator( )( const TScalar& w ) const;
94 * $\frac{d^{(n)}c(\omega)}{d\omega^{(n)}}=\sum_{l=-q}^{q}z_{l}e^{jl\omega}(jl)^{n}$
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.
101 TVector operator( )( const TScalar& w, const unsigned int& n ) const;
104 * Coefficient-wise addition operators
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 );
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;
125 * Scale and/or rotate
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 );
135 Self operator&( const TScalar& phase ) const;
136 Self& operator&=( const TScalar& phase );
137 TScalar GetPhase( const TScalar& w = TScalar( 0 ) ) const;
142 void Sample( std::vector< TPoint >& p, const unsigned long& s ) const;
148 int l, TScalar& a, TScalar& b, TScalar& t, TScalar& p
151 int l, TScalar& a, TScalar& b, TScalar& t, TScalar& p
154 template< class _TScalar2 >
155 void SetEllipses( const std::vector< _TScalar2 >& ellipses );
158 // Main method to compute all of series's values.
159 TComplex _Z( const TScalar& w, const unsigned int& n ) const;
162 void _DFT( const std::vector< TComplex >& p, unsigned int q );
165 friend std::ostream& operator<<( std::ostream& out, const Self& fs )
167 int q = fs.GetNumberOfHarmonics( );
169 out << "{" << l << ":" << fs[ l ] << "}";
170 for( ++l; l <= q; l++ )
171 out << ", {" << l << ":" << fs[ l ] << "}";
180 // -------------------------------------------------------------------------
181 template< class _TScalar >
182 template< class _TScalar2 >
183 cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
184 FourierSeriesContour( const FourierSeriesContour< _TScalar2 >& o )
186 int q = o.GetNumberOfHarmonics( );
187 this->SetNumberOfHarmonics( q );
188 for( int l = -q; l <= q; ++l )
191 TScalar( std::real( o[ l ] ) ),
192 TScalar( std::imag( o[ l ] ) )
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
204 std::vector< TComplex > aux;
205 for( _TIterator it = b; it != e; ++it )
208 TScalar( ( *it )[ 0 ] ),
209 TScalar( ( *it )[ 1 ] )
212 this->_DFT( aux, q );
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 )
222 int q = o.GetNumberOfHarmonics( );
223 this->SetNumberOfHarmonics( q );
224 for( int l = -q; l <= q; ++l )
227 TScalar( std::real( o[ l ] ) ),
228 TScalar( std::imag( o[ l ] ) )
233 // -------------------------------------------------------------------------
234 template< class _TScalar >
235 template< class _TScalar2 >
236 void cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
237 SetEllipses( const std::vector< _TScalar2 >& ellipses )
239 int Q = ( ellipses.size( ) - 2 ) >> 2;
241 this->SetNumberOfHarmonics( Q );
242 auto cIt = ellipses.begin( );
243 TScalar cx = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
244 if( cIt != ellipses.end( ) )
246 TScalar cy = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
247 if( cIt != ellipses.end( ) )
249 ( *this )[ 0 ] = TComplex( cx, cy );
250 for( int l = 1; l <= Q; ++l )
252 TScalar a = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
253 if( cIt != ellipses.end( ) )
255 TScalar b = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
256 if( cIt != ellipses.end( ) )
258 TScalar t = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
259 if( cIt != ellipses.end( ) )
261 TScalar p = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 );
262 if( cIt != ellipses.end( ) )
264 this->SetEllipse( l, a, b, t, p );
269 #endif // __cpExtensions__DataStructures__FourierSeriesContour__h__