lib/cpExtensions/DataStructures/FourierSeriesContour.cxx

   1 /*
   2   #include <cmath>
   3   #include <limits>
   4   #include <map>
   5 */
   6 #include <vnl/vnl_math.h>
   7 #include <cpExtensions/DataStructures/FourierSeriesContour.h>
   8
   9 // -------------------------------------------------------------------------
  10 template< class _TScalar >
  11 cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
  12 FourierSeriesContour( unsigned int q )
  13 {
  14   this->resize( ( q << 1 ) + 1, TComplex( TScalar( 0 ), TScalar( 0 ) ) );
  15 }
  16
  17 // -------------------------------------------------------------------------
  18 template< class _TScalar >
  19 cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
  20 ~FourierSeriesContour( )
  21 {
  22 }
  23
  24 // -------------------------------------------------------------------------
  25 template< class _TScalar >
  26 _TScalar cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
  27 GetArea( ) const
  28 {
  29   TScalar a = TScalar( 0 );
  30   int q = this->GetNumberOfHarmonics( );
  31   typename Self::const_iterator i = this->begin( );
  32   for( int l = -q; i != this->end( ); ++i, ++l )
  33     a += TScalar( l ) * std::norm( *i );
  34   return( TScalar( vnl_math::pi ) * a );
  35 }
  36
  37 // -------------------------------------------------------------------------
  38 template< class _TScalar >
  39 _TScalar cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
  40 GetCurvature( const _TScalar& w ) const
  41 {
  42   TComplex d1 = this->_Z( w, 1 );
  43   TComplex d2 = this->_Z( w, 2 );
  44   TScalar d = std::abs( d1 );
  45   if( d > TScalar( 0 ) )
  46     return( std::imag( std::conj( d1 ) * d2 ) / ( d * d * d ) );
  47   else
  48     return( TScalar( 0 ) );
  49 }
  50
  51 // -------------------------------------------------------------------------
  52 template< class _TScalar >
  53 void cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
  54 SetOrdering( bool cw )
  55 {
  56   // Check if the area sign is different from the desired orientation
  57   bool ap = ( TScalar( 0 ) < this->GetArea( ) );
  58   if( ( !ap && cw ) || ( ap && !cw ) ) // XOR
  59     this->InvertOrdering( );
  60 }
  61
  62 // -------------------------------------------------------------------------
  63 template< class _TScalar >
  64 void cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
  65 InvertOrdering( )
  66 {
  67   int q = this->GetNumberOfHarmonics( );
  68   for( int l = 1; l <= q; l++ )
  69   {
  70     TComplex z = ( *this )[ l ];
  71     ( *this )[  l ] = ( *this )[ -l ];
  72     ( *this )[ -l ] = z;
  73
  74   } // rof
  75 }
  76
  77 // -------------------------------------------------------------------------
  78 template< class _TScalar >
  79 void cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
  80 SetOrderingToClockWise( )
  81 {
  82   this->SetOrdering( false );
  83 }
  84
  85 // -------------------------------------------------------------------------
  86 template< class _TScalar >
  87 void cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
  88 SetOrderingToCounterClockWise( )
  89 {
  90   this->SetOrdering( true );
  91 }
  92
  93 // -------------------------------------------------------------------------
  94 template< class _TScalar >
  95 int cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
  96 GetNumberOfHarmonics( ) const
  97 {
  98   return( ( this->size( ) - 1 ) >> 1 );
  99 }
 100
 101 // -------------------------------------------------------------------------
 102 template< class _TScalar >
 103 void cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 104 SetNumberOfHarmonics( const int& q )
 105 {
 106   int diff_q = this->GetNumberOfHarmonics( ) - q;
 107
 108   while( diff_q != 0 )
 109   {
 110     if( diff_q > 0 )
 111     {
 112       this->pop_front( );
 113       this->pop_back( );
 114       diff_q--;
 115     }
 116     else
 117     {
 118       this->push_front( TComplex( TScalar( 0 ), TScalar( 0 ) ) );
 119       this->push_back( TComplex( TScalar( 0 ), TScalar( 0 ) ) );
 120       diff_q++;
 121
 122     } // fi
 123
 124   } // elihw
 125 }
 126
 127 // -------------------------------------------------------------------------
 128 template< class _TScalar >
 129 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 130 TComplex& cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 131 operator[]( const int& l )
 132 {
 133   static TComplex _zero;
 134   int idx = ( ( this->size( ) - 1 ) >> 1 ) + l;
 135   if( idx < 0 || idx >= this->size( ) )
 136   {
 137     _zero = TComplex( TScalar( 0 ) );
 138     return( _zero );
 139   }
 140   else
 141     return( this->Superclass::operator[]( idx ) );
 142 }
 143
 144 // -------------------------------------------------------------------------
 145 template< class _TScalar >
 146 const typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 147 TComplex& cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 148 operator[]( const int& l ) const
 149 {
 150   static const TComplex _zero( TScalar( 0 ) );
 151   int idx = ( ( this->size( ) - 1 ) >> 1 ) + l;
 152   if( idx < 0 || idx >= this->size( ) )
 153     return( _zero );
 154   else
 155     return( this->Superclass::operator[]( idx ) );
 156 }
 157
 158 // -------------------------------------------------------------------------
 159 template< class _TScalar >
 160 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 161 TPoint cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 162 operator( )( const _TScalar& w ) const
 163 {
 164   TComplex z = this->_Z( w, 0 );
 165   TPoint p;
 166   p[ 0 ] = z.real( );
 167   p[ 1 ] = z.imag( );
 168   return( p );
 169 }
 170
 171 // -------------------------------------------------------------------------
 172 template< class _TScalar >
 173 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 174 TVector cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 175 operator( )( const _TScalar& w, const unsigned int& n ) const
 176 {
 177   TComplex z = this->_Z( w, n );
 178   TVector v;
 179   v[ 0 ] = z.real( );
 180   v[ 1 ] = z.imag( );
 181   return( v );
 182 }
 183
 184 // -------------------------------------------------------------------------
 185 template< class _TScalar >
 186 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 187 Self cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 188 operator+( const Self& o ) const
 189 {
 190   int q1 = o.GetNumberOfHarmonics( );
 191   int q2 = this->GetNumberOfHarmonics( );
 192   int q = ( q2 < q1 )? q1: q2;
 193
 194   Self res;
 195   res.SetNumberOfHarmonics( q );
 196   for( int l = -q; l <= q; ++l )
 197     res[ l ] = ( *this )[ l ] + o[ l ];
 198
 199   return( res );
 200 }
 201
 202 // -------------------------------------------------------------------------
 203 template< class _TScalar >
 204 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 205 Self& cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 206 operator+=( const Self& o )
 207 {
 208   int q1 = o.GetNumberOfHarmonics( );
 209   int q2 = this->GetNumberOfHarmonics( );
 210   int q = ( q2 < q1 )? q1: q2;
 211
 212   this->SetNumberOfHarmonics( q );
 213   for( int l = -q; l <= q; ++l )
 214     ( *this )[ l ] += o[ l ];
 215
 216   return( *this );
 217 }
 218
 219 // -------------------------------------------------------------------------
 220 template< class _TScalar >
 221 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 222 Self cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 223 operator-( const Self& o ) const
 224 {
 225   int q1 = o.GetNumberOfHarmonics( );
 226   int q2 = this->GetNumberOfHarmonics( );
 227   int q = ( q2 < q1 )? q1: q2;
 228
 229   Self res;
 230   res.SetNumberOfHarmonics( q );
 231   for( int l = -q; l <= q; ++l )
 232     res[ l ] = ( *this )[ l ] - o[ l ];
 233
 234   return( res );
 235 }
 236
 237 // -------------------------------------------------------------------------
 238 template< class _TScalar >
 239 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 240 Self& cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 241 operator-=( const Self& o )
 242 {
 243   int q1 = o.GetNumberOfHarmonics( );
 244   int q2 = this->GetNumberOfHarmonics( );
 245   int q = ( q2 < q1 )? q1: q2;
 246
 247   this->SetNumberOfHarmonics( q );
 248   for( int l = -q; l <= q; ++l )
 249     ( *this )[ l ] -= o[ l ];
 250
 251   return( *this );
 252 }
 253
 254 // -------------------------------------------------------------------------
 255 template< class _TScalar >
 256 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 257 Self cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 258 operator+( const TComplex& translation ) const
 259 {
 260   Self res = *this;
 261   res[ 0 ] += translation;
 262   return( res );
 263 }
 264
 265 // -------------------------------------------------------------------------
 266 template< class _TScalar >
 267 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 268 Self& cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 269 operator+=( const TComplex& translation )
 270 {
 271   ( *this )[ 0 ] += translation;
 272   return( *this );
 273 }
 274
 275 // -------------------------------------------------------------------------
 276 template< class _TScalar >
 277 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 278 Self cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 279 operator-( const TComplex& translation ) const
 280 {
 281   Self res = *this;
 282   res[ 0 ] -= translation;
 283   return( res );
 284 }
 285
 286 // -------------------------------------------------------------------------
 287 template< class _TScalar >
 288 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 289 Self& cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 290 operator-=( const TComplex& translation )
 291 {
 292   ( *this )[ 0 ] -= translation;
 293   return( *this );
 294 }
 295
 296 // -------------------------------------------------------------------------
 297 template< class _TScalar >
 298 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 299 Self cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 300 operator+( const TVector& translation ) const
 301 {
 302   return( ( *this ) + TComplex( translation[ 0 ], translation[ 1 ] ) );
 303 }
 304
 305 // -------------------------------------------------------------------------
 306 template< class _TScalar >
 307 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 308 Self& cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 309 operator+=( const TVector& translation )
 310 {
 311   ( *this ) += TComplex( translation[ 0 ], translation[ 1 ] );
 312   return( *this );
 313 }
 314
 315 // -------------------------------------------------------------------------
 316 template< class _TScalar >
 317 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 318 Self cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 319 operator-( const TVector& translation ) const
 320 {
 321   return( ( *this ) - TComplex( translation[ 0 ], translation[ 1 ] ) );
 322 }
 323
 324 // -------------------------------------------------------------------------
 325 template< class _TScalar >
 326 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 327 Self& cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 328 operator-=( const TVector& translation )
 329 {
 330   ( *this ) -= TComplex( translation[ 0 ], translation[ 1 ] );
 331   return( *this );
 332 }
 333
 334 // -------------------------------------------------------------------------
 335 template< class _TScalar >
 336 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 337 TPoint cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 338 GetCenter( ) const
 339 {
 340   TComplex z = ( *this )[ 0 ];
 341   TPoint p;
 342   p[ 0 ] = z.real( );
 343   p[ 1 ] = z.imag( );
 344   return( p );
 345 }
 346
 347 // -------------------------------------------------------------------------
 348 template< class _TScalar >
 349 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 350 Self cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 351 operator*( const TComplex& scale_or_rotation ) const
 352 {
 353   Self res;
 354   res.clear( );
 355   for( typename Self::const_iterator i = this->begin( ); i != this->end( ); ++i )
 356     res.push_back( *i * scale_or_rotation );
 357   return( res );
 358 }
 359
 360 // -------------------------------------------------------------------------
 361 template< class _TScalar >
 362 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 363 Self& cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 364 operator*=( const TComplex& scale_or_rotation )
 365 {
 366   for( typename Self::iterator i = this->begin( ); i != this->end( ); ++i )
 367     *i *= scale_or_rotation;
 368   return( *this );
 369 }
 370
 371 // -------------------------------------------------------------------------
 372 template< class _TScalar >
 373 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 374 Self cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 375 operator/( const _TScalar& scale ) const
 376 {
 377   Self res;
 378   res.clear( );
 379   for( typename Self::const_iterator i = this->begin( ); i != this->end( ); ++i )
 380     res.push_back( *i / scale );
 381   return( res );
 382 }
 383
 384 // -------------------------------------------------------------------------
 385 template< class _TScalar >
 386 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 387 Self& cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 388 operator/=( const _TScalar& scale )
 389 {
 390   for( typename Self::iterator i = this->begin( ); i != this->end( ); ++i )
 391     *i /= scale;
 392   return( *this );
 393 }
 394
 395 // -------------------------------------------------------------------------
 396 template< class _TScalar >
 397 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 398 Self cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 399 operator&( const _TScalar& phase ) const
 400 {
 401   Self res;
 402   res.clear( );
 403   int q = this->GetNumberOfHarmonics( );
 404   typename Self::const_iterator i = this->begin( );
 405   for( int l = -q; i != this->end( ); ++i, ++l )
 406     res.push_back( *i * std::polar( TScalar( 1 ), TScalar( l ) * phase ) );
 407   return( res );
 408 }
 409
 410 // -------------------------------------------------------------------------
 411 template< class _TScalar >
 412 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 413 Self& cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 414 operator&=( const _TScalar& phase )
 415 {
 416   int q = this->GetNumberOfHarmonics( );
 417   typename Self::iterator i = this->begin( );
 418   for( int l = -q; i != this->end( ); ++i, ++l )
 419     *i *= std::polar( TScalar( 1 ), TScalar( l ) * phase );
 420   return( *this );
 421 }
 422
 423 // -------------------------------------------------------------------------
 424 template< class _TScalar >
 425 _TScalar cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 426 GetPhase( const _TScalar& w ) const
 427 {
 428   // Some constant values
 429   static const TScalar _0 = TScalar( 0 );
 430   static const TScalar _1 = TScalar( 1 );
 431   static const TScalar _2 = TScalar( 2 );
 432   static const TScalar _2pi = _2 * TScalar( vnl_math::pi );
 433   static const TScalar _epsilon = TScalar( 1e-14 );
 434   static const TScalar _tolerance = TScalar( 1e-10 );
 435   static const unsigned int _samples = 100;
 436   static const unsigned int _maxIterations = 1000;
 437   static const TScalar _angleOff = _2pi / TScalar( _samples );
 438
 439   // Some other values
 440   int q = int( this->GetNumberOfHarmonics( ) );
 441   if( q == 0 )
 442     return( _0 );
 443
 444   // Get a centered copy
 445   Self contour = *this;
 446   contour[ 0 ] = TComplex( _0, _0 );
 447
 448   // Compute maximum coefficient norm
 449   TScalar max_norm = _0;
 450   for( int l = 1; l <= q; ++l )
 451   {
 452     TScalar np = std::abs( contour[  l ] );
 453     TScalar nn = std::abs( contour[ -l ] );
 454     TScalar n = ( np < nn )? nn: np;
 455     if( max_norm < n )
 456       max_norm = n;
 457
 458   } // rof
 459
 460   // Rotate to desired phase
 461   contour *= std::polar( _1, -w );
 462
 463   // Try to normalize contour: if not, malformed contour!
 464   if( max_norm > _0 )
 465   {
 466     contour /= max_norm;
 467
 468     // 1. Approximate initial guess by a simple dichotomy
 469     std::vector< std::pair< TScalar, TScalar > > function;
 470     TScalar minA = std::numeric_limits< TScalar >::max( );
 471     unsigned int minIdx = -1;
 472
 473     // 1.1. Roughly sample phases
 474     for( unsigned int s = 0; s < _samples; ++s )
 475     {
 476       TScalar w = TScalar( s ) * _angleOff;
 477       TScalar a = std::arg( contour._Z( w, 0 ) );
 478       function.push_back( std::pair< TScalar, TScalar >( w, a ) );
 479       if( a < minA )
 480       {
 481         minA = a;
 482         minIdx = s;
 483
 484       } // fi
 485
 486     } // rof
 487
 488     // 1.2. Get zero cuts by zero crossing analysis and keep the farthest
 489     //      point in the real axis (ie. the last data in "cuts")
 490     TScalar prevA = _1;
 491     std::map< TScalar, unsigned int > cuts;
 492     for( unsigned int s = 0; s < _samples; ++s )
 493     {
 494       unsigned int i = ( s + minIdx ) % _samples;
 495       if( function[ i ].second * prevA < _0 )
 496         cuts[ std::real( contour._Z( function[ i ].first, 0 ) ) ] = i;
 497       prevA = function[ i ].second;
 498
 499     } // rof
 500
 501     // 1.3. Get initial guess
 502     TScalar w0 = _0;
 503     if( cuts.size( ) > 0 )
 504     {
 505       unsigned int i = cuts.rbegin( )->second;
 506       if( i == 0 )
 507         i = function.size( );
 508       w0 = function[ i - 1 ].first;
 509
 510     } // fi
 511
 512     // 2. Refine by Newthon-Raphson
 513     bool stop = false;
 514     unsigned int i = 0;
 515     while( i < _maxIterations && !stop )
 516     {
 517       TComplex c = contour._Z( w0, 0 );
 518       TScalar d = std::imag( c * std::conj( contour._Z( w0, 1 ) ) );
 519
 520       // Emergency stop!!!
 521       if( !( std::fabs( d ) < _epsilon ) )
 522       {
 523         TScalar w1 = w0 + ( std::arg( c ) * ( std::norm( c ) / d ) );
 524         if( ( std::fabs( w1 - w0 ) / std::fabs( w1 ) ) < _tolerance )
 525           stop = true;
 526         else
 527           w0 = w1;
 528         i++;
 529       }
 530       else
 531         stop = true;
 532
 533     } // elihw
 534     return( w0 );
 535   }
 536   else
 537     return( _0 );
 538 }
 539
 540 // -------------------------------------------------------------------------
 541 template< class _TScalar >
 542 void cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 543 Sample( std::vector< TPoint >& p, const unsigned long& s ) const
 544 {
 545   static const TScalar _2pi = TScalar( 2 ) * TScalar( vnl_math::pi );
 546   TScalar off = _2pi / TScalar( s - 1 );
 547   for( unsigned int w = 0; w < s; ++w )
 548     p.push_back( ( *this )( TScalar( w ) * off ) );
 549 }
 550
 551 // -------------------------------------------------------------------------
 552 template< class _TScalar >
 553 void cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 554 GetEllipse( int l, _TScalar& a, _TScalar& b, _TScalar& t, _TScalar& p ) const
 555 {
 556   a = b = t = p = TScalar( 0 );
 557   if( l == 0 || l > this->GetNumberOfHarmonics( ) )
 558     return;
 559
 560   TComplex zp = ( *this )[  l ];
 561   TComplex zn = ( *this )[ -l ];
 562   TScalar np = std::abs( zp );
 563   TScalar nn = std::abs( zn );
 564
 565   // Radii
 566   a = np + nn;
 567   b = np - nn;
 568
 569   // Rotation and phase
 570   if( np > TScalar( 0 ) && nn > TScalar( 0 ) )
 571   {
 572     zp /= np;
 573     zn /= nn;
 574     t = std::real(
 575       std::log( zp * zn ) / TComplex( TScalar( 0 ), TScalar( 2 ) )
 576       );
 577     p = std::real(
 578       std::log( zp / zn ) / TComplex( TScalar( 0 ), TScalar( 2 * l ) )
 579       );
 580   }
 581   else
 582   {
 583     t = TScalar( 0 );
 584     p = ( np > TScalar( 0 ) )? std::arg( zp ): std::arg( zn );
 585
 586   } // fi
 587 }
 588
 589 // -------------------------------------------------------------------------
 590 template< class _TScalar >
 591 void cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 592 SetEllipse( int l, _TScalar& a, _TScalar& b, _TScalar& t, _TScalar& p )
 593 {
 594   TScalar nzp = ( a + b ) / TScalar( 2 );
 595   TScalar nzn = ( a - b ) / TScalar( 2 );
 596   TComplex et = std::polar( TScalar( 1 ), t );
 597   TComplex ep = std::polar( TScalar( 1 ), TScalar( l ) * p );
 598   ( *this )[  l ] = et * ep * nzp;
 599   ( *this )[ -l ] = et * std::conj( ep ) * nzn;
 600 }
 601
 602 // -------------------------------------------------------------------------
 603 template< class _TScalar >
 604 typename cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 605 TComplex cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 606 _Z( const _TScalar& w, const unsigned int& n ) const
 607 {
 608   static const TScalar _0 = TScalar( 0 );
 609   static const TScalar _1 = TScalar( 1 );
 610   TScalar _n = TScalar( n );
 611
 612   TComplex z( _0, _0 );
 613   int q = this->GetNumberOfHarmonics( );
 614   for( int l = -q; l <= q; ++l )
 615   {
 616     TComplex v = ( *this )[ l ] * std::polar( _1, w * TScalar( l ) );
 617     if( n > 0 )
 618       v *= std::pow( TComplex( _0, TScalar( l ) ), _n );
 619     z += v;
 620
 621   } // rof
 622   return( z );
 623 }
 624
 625 // -------------------------------------------------------------------------
 626 template< class _TScalar >
 627 void cpExtensions::DataStructures::FourierSeriesContour< _TScalar >::
 628 _DFT( const std::vector< TComplex >& p, unsigned int q )
 629 {
 630   static const TScalar _2pi = TScalar( 2 ) * TScalar( vnl_math::pi );
 631
 632   this->SetNumberOfHarmonics( q );
 633   *this *= TScalar( 0 );
 634   if( p.size( ) == 0 )
 635     return;
 636
 637   std::vector< TComplex > dft;
 638   for( long m = 0; m < p.size( ); ++m )
 639   {
 640     TComplex z( TScalar( 0 ), TScalar( 0 ) );
 641     for( long k = 0; k < p.size( ); ++k )
 642       z +=
 643         p[ k ] *
 644         std::polar(
 645           TScalar( 1 ),
 646           -( _2pi * TScalar( m * k ) ) / TScalar( p.size( ) )
 647           );
 648     z /= TScalar( p.size( ) );
 649     dft.push_back( z );
 650
 651   } // rof
 652
 653   ( *this )[ 0 ] = dft[ 0 ];
 654   for( int l = 1; l <= q; ++l )
 655   {
 656     ( *this )[  l ] = dft[ l ];
 657     ( *this )[ -l ] = dft[ p.size( ) - l ];
 658
 659   } // rof
 660 }
 661
 662 // -------------------------------------------------------------------------
 663 template class cpExtensions::DataStructures::FourierSeriesContour< float >;
 664 template class cpExtensions::DataStructures::FourierSeriesContour< double >;
 665
 666 // eof - $RCSfile$