From 3944412129d5da673a93005bdfb0d22043cbaf28 Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Sun, 22 Oct 2017 22:06:52 -0500 Subject: [PATCH] ... --- lib/ivq/ITK/ExtractNativeSlicesImageFilter.h | 70 ++ .../ITK/ExtractNativeSlicesImageFilter.hxx | 70 ++ lib/ivq/ITK/FourierSeries.cxx | 732 ++++++++++++++++++ lib/ivq/ITK/FourierSeries.h | 192 +++++ lib/ivq/ITK/FourierSeries.hxx | 101 +++ 5 files changed, 1165 insertions(+) create mode 100644 lib/ivq/ITK/ExtractNativeSlicesImageFilter.h create mode 100644 lib/ivq/ITK/ExtractNativeSlicesImageFilter.hxx create mode 100644 lib/ivq/ITK/FourierSeries.cxx create mode 100644 lib/ivq/ITK/FourierSeries.h create mode 100644 lib/ivq/ITK/FourierSeries.hxx diff --git a/lib/ivq/ITK/ExtractNativeSlicesImageFilter.h b/lib/ivq/ITK/ExtractNativeSlicesImageFilter.h new file mode 100644 index 0000000..df242cc --- /dev/null +++ b/lib/ivq/ITK/ExtractNativeSlicesImageFilter.h @@ -0,0 +1,70 @@ +// ======================================================================= +// @author: Leonardo Florez-Valencia +// @email: florez-l@javeriana.edu.co +// ======================================================================= +#ifndef __ivq__ITK__ExtractNativeSlicesImageFilter__h__ +#define __ivq__ITK__ExtractNativeSlicesImageFilter__h__ + +#include +#include +#include +#include +#include + +namespace ivq +{ + namespace ITK + { + /** + */ + template< class _TImage > + class ExtractNativeSlicesImageFilter + : public itk::ProcessObject + { + public: + typedef ExtractNativeSlicesImageFilter Self; + typedef itk::ProcessObject Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TImage TImage; + typedef typename TImage::PixelType TPixel; + itkStaticConstMacro( + ImageDimension, unsigned int, TImage::ImageDimension + ); + itkStaticConstMacro( SliceDimension, unsigned int, ImageDimension - 1 ); + typedef itk::Image< TPixel, SliceDimension > TSlice; + + typedef std::vector< typename TSlice::Pointer > TSlices; + typedef itk::SimpleDataObjectDecorator< TSlices > TSlicesVector; + + public: + itkNewMacro( Self ); + itkTypeMacro( ExtractNativeSlicesImageFilter, itk::ImageToImageFilter ); + + itkGetConstMacro( Direction, unsigned int ); + itkSetMacro( Direction, unsigned int ); + + ivqITKInputMacro( Input, TImage ); + ivqITKOutputMacro( Output, TSlicesVector ); + + protected: + ExtractNativeSlicesImageFilter( ); + virtual ~ExtractNativeSlicesImageFilter( ); + + virtual void GenerateData( ) override; + + protected: + unsigned int m_Direction; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION +#endif // __ivq__ITK__ExtractNativeSlicesImageFilter__h__ + +// eof - $RCSfile$ diff --git a/lib/ivq/ITK/ExtractNativeSlicesImageFilter.hxx b/lib/ivq/ITK/ExtractNativeSlicesImageFilter.hxx new file mode 100644 index 0000000..f669558 --- /dev/null +++ b/lib/ivq/ITK/ExtractNativeSlicesImageFilter.hxx @@ -0,0 +1,70 @@ +// ======================================================================= +// @author: Leonardo Florez-Valencia +// @email: florez-l@javeriana.edu.co +// ======================================================================= +#ifndef __ivq__ITK__ExtractNativeSlicesImageFilter__hxx__ +#define __ivq__ITK__ExtractNativeSlicesImageFilter__hxx__ + +#include + +// ------------------------------------------------------------------------- +template< class _TImage > +ivq::ITK::ExtractNativeSlicesImageFilter< _TImage >:: +ExtractNativeSlicesImageFilter( ) + : Superclass( ), + m_Direction( 0 ) +{ + ivqITKInputConfigureMacro( Input, TImage ); + ivqITKOutputConfigureMacro( Output, TSlicesVector ); +} + +// ------------------------------------------------------------------------- +template< class _TImage > +ivq::ITK::ExtractNativeSlicesImageFilter< _TImage >:: +~ExtractNativeSlicesImageFilter( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TImage > +void ivq::ITK::ExtractNativeSlicesImageFilter< _TImage >:: +GenerateData( ) +{ + typedef itk::ExtractImageFilter< TImabe, TSlice > _TExtract; + + if( this->m_Direction >= Self::ImageDimension ) + itkExceptionMacro( << "Invalid direction: " << this->m_Direction ); + + // Some values + const TImage* input = this->GetInput( ); + TSlicesVector* output = this->GetOutput( ); + typename TImage::RegionType region = input->GetRequestedRegion( ); + typename TImage::IndexType regionIndex = region.GetIndex( ); + typename TImage::SizeType sliceSize = region.GetSize( ); + unsigned int numSlices = sliceSize[ this->m_Direction ]; + sliceSize[ this->m_Direction ] = 0; + + // Extract all slices + output->Get( ).clear( ); + for( unsigned int i = 0; i < numSlices; ++i ) + { + // Prepare extraction region + typename TImage::IndexType sliceIndex = regionIndex; + sliceIndex[ this->m_Direction ] += i; + typename TImage::RegionType sliceRegion( sliceIndex, sliceSize ); + + // Extract slice + typename _TExtract::Pointer extract = _TExtract::New( ); + extract->SetInput( input ); + extract->SetExtractionRegion( sliceRegion ); + extract->SetDirectionCollapseToIdentity( ); + extract->Update( ); + typename TSlice::Pointer slice = extract->GetOutput( ); + slice->DisconnectPipeline( ); + output->Get( ).push_back( slice ); + + } // rof +} + +#endif // __ivq__ITK__ExtractNativeSlicesImageFilter__hxx__ +// eof - $RCSfile$ diff --git a/lib/ivq/ITK/FourierSeries.cxx b/lib/ivq/ITK/FourierSeries.cxx new file mode 100644 index 0000000..04c9c2b --- /dev/null +++ b/lib/ivq/ITK/FourierSeries.cxx @@ -0,0 +1,732 @@ +// ========================================================================= +// @author: Leonardo Florez-Valencia +// @email: florez-l@javeriana.edu.co +// ========================================================================= +#include +#include +#include + +#include +#include + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +ivq::ITK::FourierSeries< _TScalar, _VDim >:: +FourierSeries( unsigned int q ) + : m_RealAxis( _VDim - 2 ), + m_ImagAxis( _VDim - 1 ) +{ + this->resize( ( q << 1 ) + 1, TComplex( TScalar( 0 ), TScalar( 0 ) ) ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +ivq::ITK::FourierSeries< _TScalar, _VDim >:: +~FourierSeries( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +const unsigned int& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +GetRealAxis( ) const +{ + return( this->m_RealAxis ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +const unsigned int& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +GetImagAxis( ) const +{ + return( this->m_ImagAxis ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +void ivq::ITK::FourierSeries< _TScalar, _VDim >:: +SetRealAxis( const unsigned int& a ) +{ + this->m_RealAxis = a; +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +void ivq::ITK::FourierSeries< _TScalar, _VDim >:: +SetImagAxis( const unsigned int& a ) +{ + this->m_ImagAxis = a; +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +_TScalar ivq::ITK::FourierSeries< _TScalar, _VDim >:: +GetArea( ) const +{ + TScalar a = TScalar( 0 ); + int q = this->GetNumberOfHarmonics( ); + typename Self::const_iterator i = this->begin( ); + for( int l = -q; i != this->end( ); ++i, ++l ) + a += TScalar( l ) * std::norm( *i ); + return( TScalar( vnl_math::pi ) * a ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +_TScalar ivq::ITK::FourierSeries< _TScalar, _VDim >:: +GetCurvature( const TScalar& w ) const +{ + TComplex d1 = this->_Z( w, 1 ); + TComplex d2 = this->_Z( w, 2 ); + TScalar d = std::abs( d1 ); + if( d > TScalar( 0 ) ) + return( std::imag( std::conj( d1 ) * d2 ) / ( d * d * d ) ); + else + return( TScalar( 0 ) ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +void ivq::ITK::FourierSeries< _TScalar, _VDim >:: +SetOrdering( bool cw ) +{ + // Check if the area sign is different from the desired orientation + bool ap = ( TScalar( 0 ) < this->GetArea( ) ); + if( ( !ap && cw ) || ( ap && !cw ) ) // XOR + this->InvertOrdering( ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +void ivq::ITK::FourierSeries< _TScalar, _VDim >:: +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 _TScalar, unsigned int _VDim > +void ivq::ITK::FourierSeries< _TScalar, _VDim >:: +SetOrderingToClockWise( ) +{ + this->SetOrdering( false ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +void ivq::ITK::FourierSeries< _TScalar, _VDim >:: +SetOrderingToCounterClockWise( ) +{ + this->SetOrdering( true ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +int ivq::ITK::FourierSeries< _TScalar, _VDim >:: +GetNumberOfHarmonics( ) const +{ + return( ( this->size( ) - 1 ) >> 1 ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +void ivq::ITK::FourierSeries< _TScalar, _VDim >:: +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( TScalar( 0 ), TScalar( 0 ) ) ); + this->push_back( TComplex( TScalar( 0 ), TScalar( 0 ) ) ); + diff_q++; + + } // fi + + } // elihw +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +TComplex& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator[]( const int& l ) +{ + static TComplex _zero; + int idx = ( ( this->size( ) - 1 ) >> 1 ) + l; + if( idx < 0 || idx >= this->size( ) ) + { + _zero = TComplex( TScalar( 0 ) ); + return( _zero ); + } + else + return( this->Superclass::operator[]( idx ) ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +const typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +TComplex& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator[]( const int& l ) const +{ + static const TComplex _zero( TScalar( 0 ) ); + int idx = ( ( this->size( ) - 1 ) >> 1 ) + l; + if( idx < 0 || idx >= this->size( ) ) + return( _zero ); + else + return( this->Superclass::operator[]( idx ) ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +TPoint ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator( )( const TScalar& w ) const +{ + TComplex z = this->_Z( w, 0 ); + TPoint p; + p.Fill( TScalar( 0 ) ); + p[ this->m_RealAxis ] = z.real( ); + p[ this->m_ImagAxis ] = z.imag( ); + return( p ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +TVector ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator( )( const TScalar& w, const unsigned int& n ) const +{ + TComplex z = this->_Z( w, n ); + TVector v( TScalar( 0 ) ); + v[ this->m_RealAxis ] = z.real( ); + v[ this->m_ImagAxis ] = z.imag( ); + return( v ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self ivq::ITK::FourierSeries< _TScalar, _VDim >:: +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 _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +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 _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self ivq::ITK::FourierSeries< _TScalar, _VDim >:: +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 _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +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 _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator+( const TComplex& translation ) const +{ + Self res = *this; + res[ 0 ] += translation; + return( res ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator+=( const TComplex& translation ) +{ + ( *this )[ 0 ] += translation; + return( *this ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator-( const TComplex& translation ) const +{ + Self res = *this; + res[ 0 ] -= translation; + return( res ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator-=( const TComplex& translation ) +{ + ( *this )[ 0 ] -= translation; + return( *this ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator+( const TVector& translation ) const +{ + return( + ( *this ) + + TComplex( translation[ this->m_RealAxis ], translation[ this->m_ImagAxis ] ) + ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator+=( const TVector& translation ) +{ + ( *this ) += + TComplex( translation[ this->m_RealAxis ], translation[ this->m_ImagAxis ] ); + return( *this ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator-( const TVector& translation ) const +{ + return( + ( *this ) - + TComplex( translation[ this->m_RealAxis ], translation[ this->m_ImagAxis ] ) + ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator-=( const TVector& translation ) +{ + ( *this ) -= + TComplex( translation[ this->m_RealAxis ], translation[ this->m_ImagAxis ] ); + return( *this ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +TPoint ivq::ITK::FourierSeries< _TScalar, _VDim >:: +GetCenter( ) const +{ + TComplex z = ( *this )[ 0 ]; + TPoint p; + p.Fill( TScalar( 0 ) ); + p[ this->m_RealAxis ] = z.real( ); + p[ this->m_ImagAxis ] = z.imag( ); + return( p ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self ivq::ITK::FourierSeries< _TScalar, _VDim >:: +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 _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +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 _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator/( const TScalar& 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 _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator/=( const TScalar& scale ) +{ + for( typename Self::iterator i = this->begin( ); i != this->end( ); ++i ) + *i /= scale; + return( *this ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator&( const TScalar& 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( TScalar( 1 ), TScalar( l ) * phase ) ); + return( res ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator&=( const TScalar& phase ) +{ + int q = this->GetNumberOfHarmonics( ); + typename Self::iterator i = this->begin( ); + for( int l = -q; i != this->end( ); ++i, ++l ) + *i *= std::polar( TScalar( 1 ), TScalar( l ) * phase ); + return( *this ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +_TScalar ivq::ITK::FourierSeries< _TScalar, _VDim >:: +GetPhase( const TScalar& w ) const +{ + // Some constant values + static const TScalar _0 = TScalar( 0 ); + static const TScalar _1 = TScalar( 1 ); + static const TScalar _2 = TScalar( 2 ); + static const TScalar _2pi = _2 * TScalar( vnl_math::pi ); + static const TScalar _epsilon = TScalar( 1e-14 ); + static const TScalar _tolerance = TScalar( 1e-10 ); + static const unsigned int _samples = 100; + static const unsigned int _maxIterations = 1000; + static const TScalar _angleOff = _2pi / TScalar( _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 + TScalar max_norm = _0; + for( int l = 1; l <= q; ++l ) + { + TScalar np = std::abs( contour[ l ] ); + TScalar nn = std::abs( contour[ -l ] ); + TScalar 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< TScalar, TScalar > > function; + TScalar minA = std::numeric_limits< TScalar >::max( ); + unsigned int minIdx = -1; + + // 1.1. Roughly sample phases + for( unsigned int s = 0; s < _samples; ++s ) + { + TScalar w = TScalar( s ) * _angleOff; + TScalar a = std::arg( contour._Z( w, 0 ) ); + function.push_back( std::pair< TScalar, TScalar >( 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") + TScalar prevA = _1; + std::map< TScalar, 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 + TScalar 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 ); + TScalar d = std::imag( c * std::conj( contour._Z( w0, 1 ) ) ); + + // Emergency stop!!! + if( !( std::fabs( d ) < _epsilon ) ) + { + TScalar 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 _TScalar, unsigned int _VDim > +void ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Sample( std::vector< TPoint >& p, const unsigned long& s ) const +{ + static const TScalar _2pi = TScalar( 2 ) * TScalar( vnl_math::pi ); + TScalar off = _2pi / TScalar( s - 1 ); + for( unsigned int w = 0; w < s; ++w ) + p.push_back( ( *this )( TScalar( w ) * off ) ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +void ivq::ITK::FourierSeries< _TScalar, _VDim >:: +GetEllipse( int l, TScalar& a, TScalar& b, TScalar& t, TScalar& p ) const +{ + a = b = t = p = TScalar( 0 ); + if( l == 0 || l > this->GetNumberOfHarmonics( ) ) + return; + + TComplex zp = ( *this )[ l ]; + TComplex zn = ( *this )[ -l ]; + TScalar np = std::abs( zp ); + TScalar nn = std::abs( zn ); + + // Radii + a = np + nn; + b = np - nn; + + // Rotation and phase + if( np > TScalar( 0 ) && nn > TScalar( 0 ) ) + { + zp /= np; + zn /= nn; + t = std::real( std::log( zp * zn ) / TComplex( TScalar( 0 ), TScalar( 2 ) ) ); + p = std::real( std::log( zp / zn ) / TComplex( TScalar( 0 ), TScalar( 2 * l ) ) ); + } + else + { + t = TScalar( 0 ); + p = ( np > TScalar( 0 ) )? std::arg( zp ): std::arg( zn ); + + } // fi +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +void ivq::ITK::FourierSeries< _TScalar, _VDim >:: +SetEllipse( int l, TScalar& a, TScalar& b, TScalar& t, TScalar& p ) +{ + TScalar nzp = ( a + b ) / TScalar( 2 ); + TScalar nzn = ( a - b ) / TScalar( 2 ); + TComplex et = std::polar( TScalar( 1 ), t ); + TComplex ep = std::polar( TScalar( 1 ), TScalar( l ) * p ); + ( *this )[ l ] = et * ep * nzp; + ( *this )[ -l ] = et * std::conj( ep ) * nzn; +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +TComplex ivq::ITK::FourierSeries< _TScalar, _VDim >:: +_Z( const TScalar& w, const unsigned int& n ) const +{ + static const TScalar _0 = TScalar( 0 ); + static const TScalar _1 = TScalar( 1 ); + TScalar _n = TScalar( n ); + + TComplex z( _0, _0 ); + int q = this->GetNumberOfHarmonics( ); + for( int l = -q; l <= q; ++l ) + { + TComplex v = ( *this )[ l ] * std::polar( _1, w * TScalar( l ) ); + if( n > 0 ) + v *= std::pow( TComplex( _0, TScalar( l ) ), _n ); + z += v; + + } // rof + return( z ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +void ivq::ITK::FourierSeries< _TScalar, _VDim >:: +_DFT( const std::vector< TComplex >& p, unsigned int q ) +{ + static const TScalar _2pi = TScalar( 2 ) * TScalar( vnl_math::pi ); + + this->SetNumberOfHarmonics( q ); + *this *= TScalar( 0 ); + if( p.size( ) == 0 ) + return; + + std::vector< TComplex > dft; + 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 ); + + } // 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 ) + { + 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 ); + + } // rof + */ +} + +// ------------------------------------------------------------------------- +template class ivq::ITK::FourierSeries< float, 2 >; +template class ivq::ITK::FourierSeries< float, 3 >; +template class ivq::ITK::FourierSeries< double, 2 >; +template class ivq::ITK::FourierSeries< double, 3 >; + +// eof - $RCSfile$ diff --git a/lib/ivq/ITK/FourierSeries.h b/lib/ivq/ITK/FourierSeries.h new file mode 100644 index 0000000..1353987 --- /dev/null +++ b/lib/ivq/ITK/FourierSeries.h @@ -0,0 +1,192 @@ +// ========================================================================= +// @author: Leonardo Florez-Valencia +// @email: florez-l@javeriana.edu.co +// ========================================================================= +#ifndef __ivq__ITK__FourierSeries__h__ +#define __ivq__ITK__FourierSeries__h__ + +#include +#include +#include +#include + +#include +#include +#include + +namespace ivq +{ + namespace ITK + { + /** + * 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 _TScalar, unsigned int _VDim > + class FourierSeries + : public std::deque< std::complex< _TScalar > > + { + public: + typedef FourierSeries Self; + typedef _TScalar TScalar; + itkStaticConstMacro( Dimension, unsigned int, _VDim ); + + typedef itk::Matrix< _TScalar, _VDim, _VDim > TMatrix; + typedef itk::Point< _TScalar, _VDim > TPoint; + typedef typename TPoint::VectorType TVector; + typedef std::complex< _TScalar > TComplex; + typedef std::deque< TComplex > Superclass; + + public: + FourierSeries( unsigned int q = 1 ); + template< class _TScalar2, unsigned int _VDim2 > + FourierSeries( const FourierSeries< _TScalar2, _VDim2 >& o ); + template< class _TIt > + FourierSeries( const _TIt& b, const _TIt& e, unsigned int q ); + virtual ~FourierSeries( ); + + template< class _TScalar2, unsigned int _VDim2 > + Self& operator=( const FourierSeries< _TScalar2, _VDim2 >& 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}$ + */ + TScalar GetArea( ) const; + + /** + * $\kappa(\omega)=\frac{|c'(\omega)\times c''(\omega)|}{|c'(\omega)|^{3}}$ + */ + TScalar GetCurvature( const TScalar& 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 TScalar& 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 TScalar& 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 TScalar& scale ) const; + Self& operator/=( const TScalar& scale ); + + /** + * Phase + */ + Self operator&( const TScalar& phase ) const; + Self& operator&=( const TScalar& phase ); + TScalar GetPhase( const TScalar& w = TScalar( 0 ) ) const; + + /** + * Sampling method + */ + void Sample( std::vector< TPoint >& p, const unsigned long& s ) const; + + /** + * Ellipse methods + */ + void GetEllipse( int l, TScalar& a, TScalar& b, TScalar& t, TScalar& p ) const; + void SetEllipse( int l, TScalar& a, TScalar& b, TScalar& t, TScalar& p ); + + template< class _TOtherScalar > + void SetEllipses( const std::vector< _TOtherScalar >& ellipses ); + + protected: + // Main method to compute all of series's values. + TComplex _Z( const TScalar& 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 ); + } + }; + + } // ecapseman + +} // ecapseman + +#include +#endif // __ivq__ITK__FourierSeries__h__ + +// eof - $RCSfile$ diff --git a/lib/ivq/ITK/FourierSeries.hxx b/lib/ivq/ITK/FourierSeries.hxx new file mode 100644 index 0000000..3ed11c5 --- /dev/null +++ b/lib/ivq/ITK/FourierSeries.hxx @@ -0,0 +1,101 @@ +// ========================================================================= +// @author: Leonardo Florez-Valencia +// @email: florez-l@javeriana.edu.co +// ========================================================================= +#ifndef __ivq__ITK__FourierSeries__hxx__ +#define __ivq__ITK__FourierSeries__hxx__ + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +template< class _TScalar2, unsigned int _VDim2 > +ivq::ITK::FourierSeries< _TScalar, _VDim >:: +FourierSeries( const FourierSeries< _TScalar2, _VDim2 >& o ) + : m_RealAxis( _VDim - 2 ), + m_ImagAxis( _VDim - 1 ) +{ + int q = o.GetNumberOfHarmonics( ); + this->SetNumberOfHarmonics( q ); + for( int l = -q; l <= q; ++l ) + ( *this )[ l ] = + TComplex( TScalar( std::real( o[ l ] ) ), TScalar( std::imag( o[ l ] ) ) ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +template< class _TIt > +ivq::ITK::FourierSeries< _TScalar, _VDim >:: +FourierSeries( const _TIt& b, const _TIt& e, unsigned int q ) + : m_RealAxis( _VDim - 2 ), + m_ImagAxis( _VDim - 1 ) +{ + std::vector< TComplex > aux; + for( _TIt it = b; it != e; ++it ) + aux.push_back( TComplex( TScalar( ( *it )[ 0 ] ), TScalar( ( *it )[ 1 ] ) ) ); + this->_DFT( aux, q ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +template< class _TScalar2, unsigned int _VDim2 > +typename ivq::ITK::FourierSeries< _TScalar, _VDim >:: +Self& ivq::ITK::FourierSeries< _TScalar, _VDim >:: +operator=( const FourierSeries< _TScalar2, _VDim2 >& o ) +{ + if( _VDim == _VDim2 ) + { + this->m_RealAxis = o.GetRealAxis( ); + this->m_ImagAxis = o.GetImagAxis( ); + } + else + { + this->m_RealAxis = _VDim - 2; + this->m_ImagAxis = _VDim - 1; + + } // fi + int q = o.GetNumberOfHarmonics( ); + this->SetNumberOfHarmonics( q ); + for( int l = -q; l <= q; ++l ) + ( *this )[ l ] = + TComplex( TScalar( std::real( o[ l ] ) ), TScalar( std::imag( o[ l ] ) ) ); + return( *this ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDim > +template< class _TOtherScalar > +void ivq::ITK::FourierSeries< _TScalar, _VDim >:: +SetEllipses( const std::vector< _TOtherScalar >& ellipses ) +{ + int Q = ( ellipses.size( ) - 2 ) >> 2; + Q = ( Q < 1 )? 1: Q; + this->SetNumberOfHarmonics( Q ); + auto cIt = ellipses.begin( ); + TScalar cx = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 ); + if( cIt != ellipses.end( ) ) + ++cIt; + TScalar cy = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 ); + if( cIt != ellipses.end( ) ) + ++cIt; + ( *this )[ 0 ] = TComplex( cx, cy ); + for( int l = 1; l <= Q; ++l ) + { + TScalar a = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 ); + if( cIt != ellipses.end( ) ) + ++cIt; + TScalar b = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 ); + if( cIt != ellipses.end( ) ) + ++cIt; + TScalar t = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 ); + if( cIt != ellipses.end( ) ) + ++cIt; + TScalar p = ( cIt != ellipses.end( ) )? *cIt: TScalar( 0 ); + if( cIt != ellipses.end( ) ) + ++cIt; + this->SetEllipse( l, a, b, t, p ); + + } // rof +} + +#endif // __ivq__ITK__FourierSeries__hxx__ + +// eof - $RCSfile$ -- 2.47.1