...

diff --git a/lib/ivq/ITK/ExtractNativeSlicesImageFilter.h b/lib/ivq/ITK/ExtractNativeSlicesImageFilter.h
new file mode 100644 (file)
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 <vector>
+#include <itkImage.h>
+#include <itkProcessObject.h>
+#include <itkSimpleDataObjectDecorator.h>
+#include <ivq/Config.h>
+
+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 <ivq/ITK/ExtractNativeSlicesImageFilter.hxx>
+#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 (file)
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 <itkExtractImageFilter.h>
+
+// -------------------------------------------------------------------------
+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 (file)
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 <cmath>
+#include <limits>
+#include <map>
+
+#include <vnl/vnl_math.h>
+#include <ivq/ITK/FourierSeries.h>
+
+// -------------------------------------------------------------------------
+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 (file)
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 <complex>
+#include <deque>
+#include <vector>
+#include <utility>
+
+#include <itkMatrix.h>
+#include <itkPoint.h>
+#include <itkVector.h>
+
+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 <ivq/ITK/FourierSeries.hxx>
+#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 (file)
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$