From: Leonardo Florez-Valencia Date: Wed, 25 Oct 2017 15:19:21 +0000 (-0500) Subject: ... X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=1f1f72ac253fc553fe715ffd08f69c7526b92ce4;p=FrontAlgorithms.git ... --- diff --git a/appli/CTArteries/CTArteries.h b/appli/CTArteries/CTArteries.h index 0fd1b45..4052ab9 100644 --- a/appli/CTArteries/CTArteries.h +++ b/appli/CTArteries/CTArteries.h @@ -11,12 +11,10 @@ #include +#include #include - #include -#include "FourierSeries.h" - class vtkPolyData; class vtkSplineWidget; class QDialog; @@ -59,7 +57,7 @@ public: typedef itk::ImageToVTKImageFilter< TScalarImage > TVTKScalarImage; typedef ivq::ITK::ImagePath< Dim > TAxis; typedef ivq::ITK::Simple3DCurve< TScalar > TCurve; - typedef FourierSeries< TScalar, 2 > TFourier; + typedef ivq::ITK::FourierSeries< TScalar, 2 > TFourier; public: explicit CTArteries( diff --git a/appli/CTArteries/FourierSeries.cxx b/appli/CTArteries/FourierSeries.cxx deleted file mode 100644 index c7bf079..0000000 --- a/appli/CTArteries/FourierSeries.cxx +++ /dev/null @@ -1,728 +0,0 @@ -#include -#include -#include - -#include -#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$ diff --git a/appli/CTArteries/FourierSeries.h b/appli/CTArteries/FourierSeries.h deleted file mode 100644 index 9194bfa..0000000 --- a/appli/CTArteries/FourierSeries.h +++ /dev/null @@ -1,187 +0,0 @@ -/* ======================================================================= - * @author: Leonardo Florez-Valencia - * @email: florez-l@javeriana.edu.co - * ======================================================================= - */ -#ifndef __FourierSeries__h__ -#define __FourierSeries__h__ - -#include -#include -#include -#include - -#include -#include -#include - -/** - * Planar series defined by its harmonic ellipses. - * - * The tuple [z(l),z(-l)], l != 0 defines an ellipse. The z(0) - * harmonic defines the series's center. - * - */ -template< class S, unsigned int D > -class FourierSeries - : public std::deque< std::complex< S > > -{ -public: - typedef FourierSeries Self; - typedef S TScalar; - itkStaticConstMacro( Dimension, unsigned int, D ); - - typedef itk::Matrix< S, D, D > TMatrix; - typedef itk::Point< S, D > TPoint; - typedef typename TPoint::VectorType TVector; - typedef std::complex< S > TComplex; - typedef std::deque< TComplex > Superclass; - -public: - FourierSeries( unsigned int q = 1 ); - template< class S2, unsigned int D2 > - FourierSeries( const FourierSeries< S2, D2 >& o ); - template< class _TIt > - FourierSeries( const _TIt& b, const _TIt& e, unsigned int q ); - virtual ~FourierSeries( ); - - template< class S2, unsigned int D2 > - Self& operator=( const FourierSeries< S2, D2 >& o ); - - bool operator==( const Self& o ) const { return( false ); } - bool operator!=( const Self& o ) const { return( true ); } - - const unsigned int& GetRealAxis( ) const; - const unsigned int& GetImagAxis( ) const; - void SetRealAxis( const unsigned int& a ); - void SetImagAxis( const unsigned int& a ); - - /** - * $A=\pi\sum_{l=-q}^{q}l|z_{l}|^{2}$ - */ - S GetArea( ) const; - - /** - * $\kappa(\omega)=\frac{|c'(\omega)\times c''(\omega)|}{|c'(\omega)|^{3}}$ - */ - S GetCurvature( const S& w ) const; - - /** - * Sets the order of the ellipse. This order is defined as the - * counterclock/clock wiseness of the series. Calculations are - * done by computing the series's area (\see Area) and comparing - * its sign. - * - * \param cw counter-clock - */ - void SetOrdering( bool cw ); - void InvertOrdering( ); - void SetOrderingToClockWise( ); - void SetOrderingToCounterClockWise( ); - - /** - * Series manipulators. - */ - int GetNumberOfHarmonics( ) const; - void SetNumberOfHarmonics( const int& q ); - TComplex& operator[]( const int& l ); - const TComplex& operator[]( const int& l ) const; - - /** - * $c(\omega)=\sum_{l=-q}^{q}z_{l}e^{jl\omega}$ - * - * where $0\le\omega<2\pi$ is the angular curvilinear parameter, q is - * the number of harmonics defining the series, $z_{l}\in\mathbb{C}$ - * is the l-th harmonic and $j$ is the imaginary unity. - */ - TPoint operator( )( const S& w ) const; - - /** - * $\frac{d^{(n)}c(\omega)}{d\omega^{(n)}}=\sum_{l=-q}^{q}z_{l}e^{jl\omega}(jl)^{n}$ - * - * where $n$ is the derivative number, $0\le\omega<2\pi$ is the - * angular curvilinear parameter, q is the number of harmonics - * defining the series, $z_{l}\in\mathbb{C}$ is the l-th harmonic - * and $j$ is the imaginary unity. - */ - TVector operator( )( const S& w, const unsigned int& n ) const; - - /** - * Coefficient-wise addition operators - */ - Self operator+( const Self& o ) const; - Self& operator+=( const Self& o ); - Self operator-( const Self& o ) const; - Self& operator-=( const Self& o ); - - /** - * Translation - */ - Self operator+( const TComplex& translation ) const; - Self& operator+=( const TComplex& translation ); - Self operator-( const TComplex& translation ) const; - Self& operator-=( const TComplex& translation ); - Self operator+( const TVector& translation ) const; - Self& operator+=( const TVector& translation ); - Self operator-( const TVector& translation ) const; - Self& operator-=( const TVector& translation ); - TPoint GetCenter( ) const; - - /** - * Scale and/or rotate - */ - Self operator*( const TComplex& scale_or_rotation ) const; - Self& operator*=( const TComplex& scale_or_rotation ); - Self operator/( const S& scale ) const; - Self& operator/=( const S& scale ); - - /** - * Phase - */ - Self operator&( const S& phase ) const; - Self& operator&=( const S& phase ); - S GetPhase( const S& w = S( 0 ) ) const; - - /** - * Sampling method - */ - void Sample( std::vector< TPoint >& p, const unsigned long& s ) const; - - /** - * Ellipse methods - */ - void GetEllipse( int l, S& a, S& b, S& t, S& p ) const; - void SetEllipse( int l, S& a, S& b, S& t, S& p ); - - template< class _TOtherScalar > - void SetEllipses( const std::vector< _TOtherScalar >& ellipses ); - -protected: - // Main method to compute all of series's values. - TComplex _Z( const S& w, const unsigned int& n ) const; - - // Fourier transform - void _DFT( const std::vector< TComplex >& p, unsigned int q ); - -protected: - unsigned int m_RealAxis; - unsigned int m_ImagAxis; - -public: - friend std::ostream& operator<<( std::ostream& out, const Self& fs ) - { - int q = fs.GetNumberOfHarmonics( ); - int l = -q; - out << "{" << l << ":" << fs[ l ] << "}"; - for( ++l; l <= q; l++ ) - out << ", {" << l << ":" << fs[ l ] << "}"; - return( out ); - } -}; - -#ifndef ITK_MANUAL_INSTANTIATION -#include "FourierSeries.hxx" -#endif // ITK_MANUAL_INSTANTIATION -#endif // __FourierSeries__h__ - -// eof - $RCSfile$ diff --git a/appli/CTArteries/FourierSeries.hxx b/appli/CTArteries/FourierSeries.hxx deleted file mode 100644 index 6398da8..0000000 --- a/appli/CTArteries/FourierSeries.hxx +++ /dev/null @@ -1,102 +0,0 @@ -/* ======================================================================= - * @author: Leonardo Florez-Valencia - * @email: florez-l@javeriana.edu.co - * ======================================================================= - */ -#ifndef __FourierSeries__hxx__ -#define __FourierSeries__hxx__ - -// ------------------------------------------------------------------------- -template< class S, unsigned int D > -template< class S2, unsigned int D2 > -FourierSeries< S, D >:: -FourierSeries( const FourierSeries< S2, D2 >& o ) - : m_RealAxis( D - 2 ), - m_ImagAxis( D - 1 ) -{ - int q = o.GetNumberOfHarmonics( ); - this->SetNumberOfHarmonics( q ); - for( int l = -q; l <= q; ++l ) - ( *this )[ l ] = - TComplex( S( std::real( o[ l ] ) ), S( std::imag( o[ l ] ) ) ); -} - -// ------------------------------------------------------------------------- -template< class S, unsigned int D > -template< class _TIt > -FourierSeries< S, D >:: -FourierSeries( const _TIt& b, const _TIt& e, unsigned int q ) - : m_RealAxis( D - 2 ), - m_ImagAxis( D - 1 ) -{ - std::vector< TComplex > aux; - for( _TIt it = b; it != e; ++it ) - aux.push_back( TComplex( S( ( *it )[ 0 ] ), S( ( *it )[ 1 ] ) ) ); - this->_DFT( aux, q ); -} - -// ------------------------------------------------------------------------- -template< class S, unsigned int D > -template< class S2, unsigned int D2 > -typename FourierSeries< S, D >:: -Self& FourierSeries< S, D >:: -operator=( const FourierSeries< S2, D2 >& o ) -{ - if( D == D2 ) - { - this->m_RealAxis = o.GetRealAxis( ); - this->m_ImagAxis = o.GetImagAxis( ); - } - else - { - this->m_RealAxis = D - 2; - this->m_ImagAxis = D - 1; - - } // fi - int q = o.GetNumberOfHarmonics( ); - this->SetNumberOfHarmonics( q ); - for( int l = -q; l <= q; ++l ) - ( *this )[ l ] = - TComplex( S( std::real( o[ l ] ) ), S( std::imag( o[ l ] ) ) ); - return( *this ); -} - -// ------------------------------------------------------------------------- -template< class S, unsigned int D > -template< class _TOtherScalar > -void FourierSeries< S, D >:: -SetEllipses( const std::vector< _TOtherScalar >& ellipses ) -{ - int Q = ( ellipses.size( ) - 2 ) >> 2; - Q = ( Q < 1 )? 1: Q; - this->SetNumberOfHarmonics( Q ); - auto cIt = ellipses.begin( ); - S cx = ( cIt != ellipses.end( ) )? *cIt: S( 0 ); - if( cIt != ellipses.end( ) ) - ++cIt; - S cy = ( cIt != ellipses.end( ) )? *cIt: S( 0 ); - if( cIt != ellipses.end( ) ) - ++cIt; - ( *this )[ 0 ] = TComplex( cx, cy ); - for( int l = 1; l <= Q; ++l ) - { - S a = ( cIt != ellipses.end( ) )? *cIt: S( 0 ); - if( cIt != ellipses.end( ) ) - ++cIt; - S b = ( cIt != ellipses.end( ) )? *cIt: S( 0 ); - if( cIt != ellipses.end( ) ) - ++cIt; - S t = ( cIt != ellipses.end( ) )? *cIt: S( 0 ); - if( cIt != ellipses.end( ) ) - ++cIt; - S p = ( cIt != ellipses.end( ) )? *cIt: S( 0 ); - if( cIt != ellipses.end( ) ) - ++cIt; - this->SetEllipse( l, a, b, t, p ); - - } // rof -} - -#endif // __FourierSeries__hxx__ - -// eof - $RCSfile$