...

[cpPlugins.git] / lib / ivq / ITK / FourierSeries.cxx

diff --git a/lib/ivq/ITK/FourierSeries.cxx b/lib/ivq/ITK/FourierSeries.cxx
index 04c9c2b47829595e2ac722e80efe2eb4b9d61445..7fcadf0f98b1d2ef4a33f0c97c761b62b0e2e714 100644 (file)
--- a/lib/ivq/ITK/FourierSeries.cxx
+++ b/lib/ivq/ITK/FourierSeries.cxx
@@ -487,7 +487,7 @@ GetPhase( const TScalar& w ) const
   int q = int( this->GetNumberOfHarmonics( ) );
   if( q == 0 )
     return( _0 );
-  
+
   // Get a centered copy
   Self contour = *this;
   contour[ 0 ] = TComplex( _0, _0 );
@@ -670,57 +670,61 @@ template< class _TScalar, unsigned int _VDim >
 void ivq::ITK::FourierSeries< _TScalar, _VDim >::
 _DFT( const std::vector< TComplex >& p, unsigned int q )
 {
+  // Basic values and configuration
+  static const TScalar _0 = TScalar( 0 );
+  static const TScalar _1 = TScalar( 1 );
   static const TScalar _2pi = TScalar( 2 ) * TScalar( vnl_math::pi );
-
   this->SetNumberOfHarmonics( q );
   *this *= TScalar( 0 );
   if( p.size( ) == 0 )
     return;
+  TScalar N = TScalar( p.size( ) );
 
-  std::vector< TComplex > dft;
+  // Compute center
+  TComplex c( _0, _0 );
   for( long m = 0; m < p.size( ); ++m )
-  {
-    TComplex z( TScalar( 0 ), TScalar( 0 ) );
-    for( long k = 0; k < p.size( ); ++k )
-      z +=
-        p[ k ] *
-        std::polar( TScalar( 1 ), -( _2pi * TScalar( m * k ) ) / TScalar( p.size( ) ) );
-    z /= TScalar( p.size( ) );
-    dft.push_back( z );
+    c += p[ m ];
+  c /= N;
 
-  } // rof
-  
-  ( *this )[ 0 ] = dft[ 0 ];
-  for( int l = 1; l <= q; ++l )
+  // Minimize phase
+  long minIdx = 0;
+  TScalar minImag = std::fabs( std::imag( p[ 0 ] - c ) );
+  for( long m = 1; m < p.size( ); ++m )
   {
-    ( *this )[  l ] = dft[ l ];
-    ( *this )[ -l ] = dft[ p.size( ) - l ];
+    TScalar y = std::fabs( std::imag( p[ m ] - c ) );
+    if( y < minImag )
+    {
+      minIdx = m;
+      minImag = y;
+
+    } // fi
 
   } // rof
 
-  // Reset contour
-  /* TODO:
-  this->SetNumberOfHarmonics( q );
+  // Real DFT computation
+  std::vector< TComplex > dft;
+  dft.push_back( c );
+  TScalar _2piN = -_2pi / N;
+  for( long m = 1; m < p.size( ); ++m )
+  {
+    TComplex z( _0, _0 );
+    for( long i = 0; i < p.size( ); ++i )
+      z +=
+        p[ ( i + minIdx ) % p.size( ) ] *
+        std::polar( _1, _2piN * TScalar( m * i ) );
+    z /= N;
+    dft.push_back( z );
 
-  // Compute number of calculations K: cut harmonics to q
-  long N = long( p.size( ) );
-  if( N == 0 )
-    return;
-  long K = ( long( q ) << 1 ) + 1;
-  if( K < N - 1 )
-    K = N - 1;
+  } // rof
 
-  // Compute harmonics
-  for( long l = -K; l <= K; ++l )
+  // Assign clamped coefficients
+  ( *this )[ 0 ] = dft[ 0 ];
+  for( int x = 1; x <= q; ++x )
   {
-    TScalar k = ( ( TScalar( l ) + TScalar( ( l < 0 )? N: 0 ) ) * _2pi ) / TScalar( N );
-    ( *this )[ l ] = TComplex( TScalar( 0 ), TScalar( 0 ) );
-    for( long n = 0; n < N; ++n )
-      ( *this )[ l ] += p[ n ] * std::polar( TScalar( 1 ), k * TScalar( -n ) );
-    ( *this )[ l ] /= TScalar( N );
+    ( *this )[  x ] = dft[ x ];
+    ( *this )[ -x ] = dft[ p.size( ) - x ];
 
   } // rof
-  */
 }
 
 // -------------------------------------------------------------------------