+++ /dev/null
-#include <cmath>
-#include <limits>
-#include <map>
-
-#include <vnl/vnl_math.h>
-#include "FourierSeries.h"
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-FourierSeries< S, D >::
-FourierSeries( unsigned int q )
- : m_RealAxis( D - 2 ),
- m_ImagAxis( D - 1 )
-{
- this->resize( ( q << 1 ) + 1, TComplex( S( 0 ), S( 0 ) ) );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-FourierSeries< S, D >::
-~FourierSeries( )
-{
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-const unsigned int& FourierSeries< S, D >::
-GetRealAxis( ) const
-{
- return( this->m_RealAxis );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-const unsigned int& FourierSeries< S, D >::
-GetImagAxis( ) const
-{
- return( this->m_ImagAxis );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-void FourierSeries< S, D >::
-SetRealAxis( const unsigned int& a )
-{
- this->m_RealAxis = a;
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-void FourierSeries< S, D >::
-SetImagAxis( const unsigned int& a )
-{
- this->m_ImagAxis = a;
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-S FourierSeries< S, D >::
-GetArea( ) const
-{
- S a = S( 0 );
- int q = this->GetNumberOfHarmonics( );
- typename Self::const_iterator i = this->begin( );
- for( int l = -q; i != this->end( ); ++i, ++l )
- a += S( l ) * std::norm( *i );
- return( S( vnl_math::pi ) * a );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-S FourierSeries< S, D >::
-GetCurvature( const S& w ) const
-{
- TComplex d1 = this->_Z( w, 1 );
- TComplex d2 = this->_Z( w, 2 );
- S d = std::abs( d1 );
- if( d > S( 0 ) )
- return( std::imag( std::conj( d1 ) * d2 ) / ( d * d * d ) );
- else
- return( S( 0 ) );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-void FourierSeries< S, D >::
-SetOrdering( bool cw )
-{
- // Check if the area sign is different from the desired orientation
- bool ap = ( S( 0 ) < this->GetArea( ) );
- if( ( !ap && cw ) || ( ap && !cw ) ) // XOR
- this->InvertOrdering( );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-void FourierSeries< S, D >::
-InvertOrdering( )
-{
- int q = this->GetNumberOfHarmonics( );
- for( int l = 1; l <= q; l++ )
- {
- TComplex z = ( *this )[ l ];
- ( *this )[ l ] = ( *this )[ -l ];
- ( *this )[ -l ] = z;
-
- } // rof
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-void FourierSeries< S, D >::
-SetOrderingToClockWise( )
-{
- this->SetOrdering( false );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-void FourierSeries< S, D >::
-SetOrderingToCounterClockWise( )
-{
- this->SetOrdering( true );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-int FourierSeries< S, D >::
-GetNumberOfHarmonics( ) const
-{
- return( ( this->size( ) - 1 ) >> 1 );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-void FourierSeries< S, D >::
-SetNumberOfHarmonics( const int& q )
-{
- int diff_q = this->GetNumberOfHarmonics( ) - q;
-
- while( diff_q != 0 )
- {
- if( diff_q > 0 )
- {
- this->pop_front( );
- this->pop_back( );
- diff_q--;
- }
- else
- {
- this->push_front( TComplex( S( 0 ), S( 0 ) ) );
- this->push_back( TComplex( S( 0 ), S( 0 ) ) );
- diff_q++;
-
- } // fi
-
- } // elihw
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-TComplex& FourierSeries< S, D >::
-operator[]( const int& l )
-{
- static TComplex _zero;
- int idx = ( ( this->size( ) - 1 ) >> 1 ) + l;
- if( idx < 0 || idx >= this->size( ) )
- {
- _zero = TComplex( S( 0 ) );
- return( _zero );
- }
- else
- return( this->Superclass::operator[]( idx ) );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-const typename FourierSeries< S, D >::
-TComplex& FourierSeries< S, D >::
-operator[]( const int& l ) const
-{
- static const TComplex _zero( S( 0 ) );
- int idx = ( ( this->size( ) - 1 ) >> 1 ) + l;
- if( idx < 0 || idx >= this->size( ) )
- return( _zero );
- else
- return( this->Superclass::operator[]( idx ) );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-TPoint FourierSeries< S, D >::
-operator( )( const S& w ) const
-{
- TComplex z = this->_Z( w, 0 );
- TPoint p;
- p.Fill( S( 0 ) );
- p[ this->m_RealAxis ] = z.real( );
- p[ this->m_ImagAxis ] = z.imag( );
- return( p );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-TVector FourierSeries< S, D >::
-operator( )( const S& w, const unsigned int& n ) const
-{
- TComplex z = this->_Z( w, n );
- TVector v( S( 0 ) );
- v[ this->m_RealAxis ] = z.real( );
- v[ this->m_ImagAxis ] = z.imag( );
- return( v );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self FourierSeries< S, D >::
-operator+( const Self& o ) const
-{
- int q1 = o.GetNumberOfHarmonics( );
- int q2 = this->GetNumberOfHarmonics( );
- int q = ( q2 < q1 )? q1: q2;
-
- Self res;
- res.SetNumberOfHarmonics( q );
- for( int l = -q; l <= q; ++l )
- res[ l ] = ( *this )[ l ] + o[ l ];
-
- return( res );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self& FourierSeries< S, D >::
-operator+=( const Self& o )
-{
- int q1 = o.GetNumberOfHarmonics( );
- int q2 = this->GetNumberOfHarmonics( );
- int q = ( q2 < q1 )? q1: q2;
-
- this->SetNumberOfHarmonics( q );
- for( int l = -q; l <= q; ++l )
- ( *this )[ l ] += o[ l ];
-
- return( *this );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self FourierSeries< S, D >::
-operator-( const Self& o ) const
-{
- int q1 = o.GetNumberOfHarmonics( );
- int q2 = this->GetNumberOfHarmonics( );
- int q = ( q2 < q1 )? q1: q2;
-
- Self res;
- res.SetNumberOfHarmonics( q );
- for( int l = -q; l <= q; ++l )
- res[ l ] = ( *this )[ l ] - o[ l ];
-
- return( res );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self& FourierSeries< S, D >::
-operator-=( const Self& o )
-{
- int q1 = o.GetNumberOfHarmonics( );
- int q2 = this->GetNumberOfHarmonics( );
- int q = ( q2 < q1 )? q1: q2;
-
- this->SetNumberOfHarmonics( q );
- for( int l = -q; l <= q; ++l )
- ( *this )[ l ] -= o[ l ];
-
- return( *this );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self FourierSeries< S, D >::
-operator+( const TComplex& translation ) const
-{
- Self res = *this;
- res[ 0 ] += translation;
- return( res );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self& FourierSeries< S, D >::
-operator+=( const TComplex& translation )
-{
- ( *this )[ 0 ] += translation;
- return( *this );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self FourierSeries< S, D >::
-operator-( const TComplex& translation ) const
-{
- Self res = *this;
- res[ 0 ] -= translation;
- return( res );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self& FourierSeries< S, D >::
-operator-=( const TComplex& translation )
-{
- ( *this )[ 0 ] -= translation;
- return( *this );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self FourierSeries< S, D >::
-operator+( const TVector& translation ) const
-{
- return(
- ( *this ) +
- TComplex( translation[ this->m_RealAxis ], translation[ this->m_ImagAxis ] )
- );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self& FourierSeries< S, D >::
-operator+=( const TVector& translation )
-{
- ( *this ) +=
- TComplex( translation[ this->m_RealAxis ], translation[ this->m_ImagAxis ] );
- return( *this );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self FourierSeries< S, D >::
-operator-( const TVector& translation ) const
-{
- return(
- ( *this ) -
- TComplex( translation[ this->m_RealAxis ], translation[ this->m_ImagAxis ] )
- );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self& FourierSeries< S, D >::
-operator-=( const TVector& translation )
-{
- ( *this ) -=
- TComplex( translation[ this->m_RealAxis ], translation[ this->m_ImagAxis ] );
- return( *this );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-TPoint FourierSeries< S, D >::
-GetCenter( ) const
-{
- TComplex z = ( *this )[ 0 ];
- TPoint p;
- p.Fill( S( 0 ) );
- p[ this->m_RealAxis ] = z.real( );
- p[ this->m_ImagAxis ] = z.imag( );
- return( p );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self FourierSeries< S, D >::
-operator*( const TComplex& scale_or_rotation ) const
-{
- Self res;
- res.clear( );
- for( typename Self::const_iterator i = this->begin( ); i != this->end( ); ++i )
- res.push_back( *i * scale_or_rotation );
- return( res );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self& FourierSeries< S, D >::
-operator*=( const TComplex& scale_or_rotation )
-{
- for( typename Self::iterator i = this->begin( ); i != this->end( ); ++i )
- *i *= scale_or_rotation;
- return( *this );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self FourierSeries< S, D >::
-operator/( const S& scale ) const
-{
- Self res;
- res.clear( );
- for( typename Self::const_iterator i = this->begin( ); i != this->end( ); ++i )
- res.push_back( *i / scale );
- return( res );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self& FourierSeries< S, D >::
-operator/=( const S& scale )
-{
- for( typename Self::iterator i = this->begin( ); i != this->end( ); ++i )
- *i /= scale;
- return( *this );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self FourierSeries< S, D >::
-operator&( const S& phase ) const
-{
- Self res;
- res.clear( );
- int q = this->GetNumberOfHarmonics( );
- typename Self::const_iterator i = this->begin( );
- for( int l = -q; i != this->end( ); ++i, ++l )
- res.push_back( *i * std::polar( S( 1 ), S( l ) * phase ) );
- return( res );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-Self& FourierSeries< S, D >::
-operator&=( const S& phase )
-{
- int q = this->GetNumberOfHarmonics( );
- typename Self::iterator i = this->begin( );
- for( int l = -q; i != this->end( ); ++i, ++l )
- *i *= std::polar( S( 1 ), S( l ) * phase );
- return( *this );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-S FourierSeries< S, D >::
-GetPhase( const S& w ) const
-{
- // Some constant values
- static const S _0 = S( 0 );
- static const S _1 = S( 1 );
- static const S _2 = S( 2 );
- static const S _2pi = _2 * S( vnl_math::pi );
- static const S _epsilon = S( 1e-14 );
- static const S _tolerance = S( 1e-10 );
- static const unsigned int _samples = 100;
- static const unsigned int _maxIterations = 1000;
- static const S _angleOff = _2pi / S( _samples );
-
- // Some other values
- int q = int( this->GetNumberOfHarmonics( ) );
- if( q == 0 )
- return( _0 );
-
- // Get a centered copy
- Self contour = *this;
- contour[ 0 ] = TComplex( _0, _0 );
-
- // Compute maximum coefficient norm
- S max_norm = _0;
- for( int l = 1; l <= q; ++l )
- {
- S np = std::abs( contour[ l ] );
- S nn = std::abs( contour[ -l ] );
- S n = ( np < nn )? nn: np;
- if( max_norm < n )
- max_norm = n;
-
- } // rof
-
- // Rotate to desired phase
- contour *= std::polar( _1, -w );
-
- // Try to normalize contour: if not, malformed contour!
- if( max_norm > _0 )
- {
- contour /= max_norm;
-
- // 1. Approximate initial guess by a simple dichotomy
- std::vector< std::pair< S, S > > function;
- S minA = std::numeric_limits< S >::max( );
- unsigned int minIdx = -1;
-
- // 1.1. Roughly sample phases
- for( unsigned int s = 0; s < _samples; ++s )
- {
- S w = S( s ) * _angleOff;
- S a = std::arg( contour._Z( w, 0 ) );
- function.push_back( std::pair< S, S >( w, a ) );
- if( a < minA )
- {
- minA = a;
- minIdx = s;
-
- } // fi
-
- } // rof
-
- // 1.2. Get zero cuts by zero crossing analysis and keep the farthest
- // point in the real axis (ie. the last data in "cuts")
- S prevA = _1;
- std::map< S, unsigned int > cuts;
- for( unsigned int s = 0; s < _samples; ++s )
- {
- unsigned int i = ( s + minIdx ) % _samples;
- if( function[ i ].second * prevA < _0 )
- cuts[ std::real( contour._Z( function[ i ].first, 0 ) ) ] = i;
- prevA = function[ i ].second;
-
- } // rof
-
- // 1.3. Get initial guess
- S w0 = _0;
- if( cuts.size( ) > 0 )
- {
- unsigned int i = cuts.rbegin( )->second;
- if( i == 0 )
- i = function.size( );
- w0 = function[ i - 1 ].first;
-
- } // fi
-
- // 2. Refine by Newthon-Raphson
- bool stop = false;
- unsigned int i = 0;
- while( i < _maxIterations && !stop )
- {
- TComplex c = contour._Z( w0, 0 );
- S d = std::imag( c * std::conj( contour._Z( w0, 1 ) ) );
-
- // Emergency stop!!!
- if( !( std::fabs( d ) < _epsilon ) )
- {
- S w1 = w0 + ( std::arg( c ) * ( std::norm( c ) / d ) );
- if( ( std::fabs( w1 - w0 ) / std::fabs( w1 ) ) < _tolerance )
- stop = true;
- else
- w0 = w1;
- i++;
- }
- else
- stop = true;
-
- } // elihw
- return( w0 );
- }
- else
- return( _0 );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-void FourierSeries< S, D >::
-Sample( std::vector< TPoint >& p, const unsigned long& s ) const
-{
- static const S _2pi = S( 2 ) * S( vnl_math::pi );
- S off = _2pi / S( s - 1 );
- for( unsigned int w = 0; w < s; ++w )
- p.push_back( ( *this )( S( w ) * off ) );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-void FourierSeries< S, D >::
-GetEllipse( int l, S& a, S& b, S& t, S& p ) const
-{
- a = b = t = p = S( 0 );
- if( l == 0 || l > this->GetNumberOfHarmonics( ) )
- return;
-
- TComplex zp = ( *this )[ l ];
- TComplex zn = ( *this )[ -l ];
- S np = std::abs( zp );
- S nn = std::abs( zn );
-
- // Radii
- a = np + nn;
- b = np - nn;
-
- // Rotation and phase
- if( np > S( 0 ) && nn > S( 0 ) )
- {
- zp /= np;
- zn /= nn;
- t = std::real( std::log( zp * zn ) / TComplex( S( 0 ), S( 2 ) ) );
- p = std::real( std::log( zp / zn ) / TComplex( S( 0 ), S( 2 * l ) ) );
- }
- else
- {
- t = S( 0 );
- p = ( np > S( 0 ) )? std::arg( zp ): std::arg( zn );
-
- } // fi
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-void FourierSeries< S, D >::
-SetEllipse( int l, S& a, S& b, S& t, S& p )
-{
- S nzp = ( a + b ) / S( 2 );
- S nzn = ( a - b ) / S( 2 );
- TComplex et = std::polar( S( 1 ), t );
- TComplex ep = std::polar( S( 1 ), S( l ) * p );
- ( *this )[ l ] = et * ep * nzp;
- ( *this )[ -l ] = et * std::conj( ep ) * nzn;
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-typename FourierSeries< S, D >::
-TComplex FourierSeries< S, D >::
-_Z( const S& w, const unsigned int& n ) const
-{
- static const S _0 = S( 0 );
- static const S _1 = S( 1 );
- S _n = S( n );
-
- TComplex z( _0, _0 );
- int q = this->GetNumberOfHarmonics( );
- for( int l = -q; l <= q; ++l )
- {
- TComplex v = ( *this )[ l ] * std::polar( _1, w * S( l ) );
- if( n > 0 )
- v *= std::pow( TComplex( _0, S( l ) ), _n );
- z += v;
-
- } // rof
- return( z );
-}
-
-// -------------------------------------------------------------------------
-template< class S, unsigned int D >
-void FourierSeries< S, D >::
-_DFT( const std::vector< TComplex >& p, unsigned int q )
-{
- static const S _2pi = S( 2 ) * S( vnl_math::pi );
-
- this->SetNumberOfHarmonics( q );
- *this *= S( 0 );
- if( p.size( ) == 0 )
- return;
-
- std::vector< TComplex > dft;
- for( long m = 0; m < p.size( ); ++m )
- {
- TComplex z( S( 0 ), S( 0 ) );
- for( long k = 0; k < p.size( ); ++k )
- z +=
- p[ k ] *
- std::polar( S( 1 ), -( _2pi * S( m * k ) ) / S( p.size( ) ) );
- z /= S( p.size( ) );
- dft.push_back( z );
-
- } // rof
-
- ( *this )[ 0 ] = dft[ 0 ];
- for( int l = 1; l <= q; ++l )
- {
- ( *this )[ l ] = dft[ l ];
- ( *this )[ -l ] = dft[ p.size( ) - l ];
-
- } // rof
-
- // Reset contour
- /* TODO:
- this->SetNumberOfHarmonics( q );
-
- // 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;
-
- // Compute harmonics
- for( long l = -K; l <= K; ++l )
- {
- S k = ( ( S( l ) + S( ( l < 0 )? N: 0 ) ) * _2pi ) / S( N );
- ( *this )[ l ] = TComplex( S( 0 ), S( 0 ) );
- for( long n = 0; n < N; ++n )
- ( *this )[ l ] += p[ n ] * std::polar( S( 1 ), k * S( -n ) );
- ( *this )[ l ] /= S( N );
-
- } // rof
- */
-}
-
-// -------------------------------------------------------------------------
-template class FourierSeries< float, 2 >;
-template class FourierSeries< float, 3 >;
-template class FourierSeries< double, 2 >;
-template class FourierSeries< double, 3 >;
-
-// eof - $RCSfile$
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / blobdiff