int q = int( this->GetNumberOfHarmonics( ) );
if( q == 0 )
return( _0 );
-
+
// Get a centered copy
Self contour = *this;
contour[ 0 ] = TComplex( _0, _0 );
void ivq::ITK::FourierSeries< _TScalar, _VDim >::
_DFT( const std::vector< TComplex >& p, unsigned int q )
{
+ // Basic values and configuration
+ static const TScalar _0 = TScalar( 0 );
+ static const TScalar _1 = TScalar( 1 );
static const TScalar _2pi = TScalar( 2 ) * TScalar( vnl_math::pi );
-
this->SetNumberOfHarmonics( q );
*this *= TScalar( 0 );
if( p.size( ) == 0 )
return;
+ TScalar N = TScalar( p.size( ) );
- std::vector< TComplex > dft;
+ // Compute center
+ TComplex c( _0, _0 );
for( long m = 0; m < p.size( ); ++m )
- {
- TComplex z( TScalar( 0 ), TScalar( 0 ) );
- for( long k = 0; k < p.size( ); ++k )
- z +=
- p[ k ] *
- std::polar( TScalar( 1 ), -( _2pi * TScalar( m * k ) ) / TScalar( p.size( ) ) );
- z /= TScalar( p.size( ) );
- dft.push_back( z );
+ c += p[ m ];
+ c /= N;
- } // rof
-
- ( *this )[ 0 ] = dft[ 0 ];
- for( int l = 1; l <= q; ++l )
+ // Minimize phase
+ long minIdx = 0;
+ TScalar minImag = std::fabs( std::imag( p[ 0 ] - c ) );
+ for( long m = 1; m < p.size( ); ++m )
{
- ( *this )[ l ] = dft[ l ];
- ( *this )[ -l ] = dft[ p.size( ) - l ];
+ TScalar y = std::fabs( std::imag( p[ m ] - c ) );
+ if( y < minImag )
+ {
+ minIdx = m;
+ minImag = y;
+
+ } // fi
} // rof
- // Reset contour
- /* TODO:
- this->SetNumberOfHarmonics( q );
+ // Real DFT computation
+ std::vector< TComplex > dft;
+ dft.push_back( c );
+ TScalar _2piN = -_2pi / N;
+ for( long m = 1; m < p.size( ); ++m )
+ {
+ TComplex z( _0, _0 );
+ for( long i = 0; i < p.size( ); ++i )
+ z +=
+ p[ ( i + minIdx ) % p.size( ) ] *
+ std::polar( _1, _2piN * TScalar( m * i ) );
+ z /= N;
+ dft.push_back( z );
- // Compute number of calculations K: cut harmonics to q
- long N = long( p.size( ) );
- if( N == 0 )
- return;
- long K = ( long( q ) << 1 ) + 1;
- if( K < N - 1 )
- K = N - 1;
+ } // rof
- // Compute harmonics
- for( long l = -K; l <= K; ++l )
+ // Assign clamped coefficients
+ ( *this )[ 0 ] = dft[ 0 ];
+ for( int x = 1; x <= q; ++x )
{
- TScalar k = ( ( TScalar( l ) + TScalar( ( l < 0 )? N: 0 ) ) * _2pi ) / TScalar( N );
- ( *this )[ l ] = TComplex( TScalar( 0 ), TScalar( 0 ) );
- for( long n = 0; n < N; ++n )
- ( *this )[ l ] += p[ n ] * std::polar( TScalar( 1 ), k * TScalar( -n ) );
- ( *this )[ l ] /= TScalar( N );
+ ( *this )[ x ] = dft[ x ];
+ ( *this )[ -x ] = dft[ p.size( ) - x ];
} // rof
- */
}
// -------------------------------------------------------------------------
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+
+#include <ivq/VTK/PolyDataToFourierSeriesFilter.h>
+#include <itkQuadEdgeMesh.h>
+#include <vtkInformation.h>
+#include <vtkInformationVector.h>
+
+// -------------------------------------------------------------------------
+template< class _TFourierSeries >
+typename ivq::VTK::PolyDataToFourierSeriesFilter< _TFourierSeries >::
+Self* ivq::VTK::PolyDataToFourierSeriesFilter< _TFourierSeries >::
+New( )
+{
+ return( new Self( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TFourierSeries >
+const _TFourierSeries&
+ivq::VTK::PolyDataToFourierSeriesFilter< _TFourierSeries >::
+GetOutput( ) const
+{
+ return( this->m_FourierSeries );
+}
+
+// -------------------------------------------------------------------------
+template< class _TFourierSeries >
+void ivq::VTK::PolyDataToFourierSeriesFilter< _TFourierSeries >::
+SetNumberOfHarmonics( unsigned int q )
+{
+ if( this->m_FourierSeries.GetNumberOfHarmonics( ) != q )
+ {
+ this->m_FourierSeries.SetNumberOfHarmonics( q );
+ this->Modified( );
+
+ } // fi
+}
+
+// -------------------------------------------------------------------------
+template< class _TFourierSeries >
+unsigned int ivq::VTK::PolyDataToFourierSeriesFilter< _TFourierSeries >::
+GetNumberOfHarmonics( ) const
+{
+ return( this->m_FourierSeries.GetNumberOfHarmonics( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class _TFourierSeries >
+ivq::VTK::PolyDataToFourierSeriesFilter< _TFourierSeries >::
+PolyDataToFourierSeriesFilter( )
+ : vtkPolyDataAlgorithm( )
+{
+ this->SetNumberOfInputPorts( 1 );
+ this->SetNumberOfOutputPorts( 0 );
+ this->m_FourierSeries.SetNumberOfHarmonics( 1 );
+}
+
+// -------------------------------------------------------------------------
+template< class _TFourierSeries >
+ivq::VTK::PolyDataToFourierSeriesFilter< _TFourierSeries >::
+~PolyDataToFourierSeriesFilter( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class _TFourierSeries >
+int ivq::VTK::PolyDataToFourierSeriesFilter< _TFourierSeries >::
+RequestData(
+ vtkInformation* information,
+ vtkInformationVector** input,
+ vtkInformationVector* output
+ )
+{
+ // Get input polydata
+ vtkInformation* inInfo = input[ 0 ]->GetInformationObject( 0 );
+ vtkPolyData* cnt =
+ vtkPolyData::SafeDownCast(
+ inInfo->Get( vtkDataObject::DATA_OBJECT( ) )
+ );
+
+ // Use quad-edge to order points
+ typedef itk::QuadEdgeMesh< double, 2 > _TQEMesh;
+ typedef _TQEMesh::PointType _T2DPoint;
+ typedef std::vector< _T2DPoint > _T2DPoints;
+ _TQEMesh::Pointer mesh = _TQEMesh::New( );
+ for( unsigned int i = 0; i < cnt->GetNumberOfPoints( ); ++i )
+ {
+ _T2DPoint pnt;
+ double p[ 3 ];
+ cnt->GetPoint( i, p );
+ pnt[ 0 ] = p[ 0 ];
+ pnt[ 1 ] = p[ 1 ];
+ mesh->AddPoint( pnt );
+
+ } // rof
+
+ vtkCellArray* lines = cnt->GetLines( );
+ lines->InitTraversal( );
+ vtkIdType* ids = new vtkIdType[ 2 ];
+ vtkIdType npts;
+ while( lines->GetNextCell( npts, ids ) != 0 )
+ mesh->AddEdge( ids[ 0 ], ids[ 1 ] );
+ delete ids;
+ _TQEMesh::QEPrimal* edge = mesh->GetEdge( );
+ if( edge != NULL )
+ {
+ _T2DPoints points;
+ mesh->AddFace( edge );
+ edge = mesh->GetEdge( );
+ for(
+ _TQEMesh::QEPrimal::IteratorGeom eIt = edge->BeginGeomLnext( );
+ eIt != edge->EndGeomLnext( );
+ ++eIt
+ )
+ points.push_back( mesh->GetPoint( *eIt ) );
+ this->m_FourierSeries =
+ TFourierSeries(
+ points.begin( ), points.end( ),
+ this->m_FourierSeries.GetNumberOfHarmonics( )
+ );
+ this->m_FourierSeries.SetOrderingToCounterClockWise( );
+ }
+ else
+ this->m_FourierSeries *= ( typename _TFourierSeries::TScalar )( 0 );
+ return( 1 );
+}
+
+// -------------------------------------------------------------------------
+template< class _TFourierSeries >
+int ivq::VTK::PolyDataToFourierSeriesFilter< _TFourierSeries >::
+RequestInformation(
+ vtkInformation* information,
+ vtkInformationVector** input,
+ vtkInformationVector* output
+ )
+{
+ return( 1 );
+}
+
+// -------------------------------------------------------------------------
+#include <ivq/ITK/FourierSeries.h>
+
+template class ivq::VTK::PolyDataToFourierSeriesFilter< ivq::ITK::FourierSeries< float, 2 > >;
+template class ivq::VTK::PolyDataToFourierSeriesFilter< ivq::ITK::FourierSeries< float, 3 > >;
+template class ivq::VTK::PolyDataToFourierSeriesFilter< ivq::ITK::FourierSeries< double, 2 > >;
+template class ivq::VTK::PolyDataToFourierSeriesFilter< ivq::ITK::FourierSeries< double, 3 > >;
+
+// eof - $RCSfile$
--- /dev/null
+// =========================================================================
+// @author Leonardo Florez Valencia
+// @email florez-l@javeriana.edu.co
+// =========================================================================
+#ifndef __ivq__VTK__PolyDataToFourierSeriesFilter__h__
+#define __ivq__VTK__PolyDataToFourierSeriesFilter__h__
+
+#include <vtkPolyDataAlgorithm.h>
+
+namespace ivq
+{
+ namespace VTK
+ {
+ /**
+ */
+ template< class _TFourierSeries >
+ class PolyDataToFourierSeriesFilter
+ : public vtkPolyDataAlgorithm
+ {
+ public:
+ typedef PolyDataToFourierSeriesFilter Self;
+ typedef _TFourierSeries TFourierSeries;
+
+ public:
+ vtkTypeMacro( PolyDataToFourierSeriesFilter, vtkPolyDataAlgorithm );
+
+ public:
+ static Self* New( );
+
+ const TFourierSeries& GetOutput( ) const;
+ void SetNumberOfHarmonics( unsigned int q );
+ unsigned int GetNumberOfHarmonics( ) const;
+
+ protected:
+ PolyDataToFourierSeriesFilter( );
+ virtual ~PolyDataToFourierSeriesFilter( );
+
+ int RequestData(
+ vtkInformation* information,
+ vtkInformationVector** input,
+ vtkInformationVector* output
+ );
+ int RequestInformation(
+ vtkInformation* information,
+ vtkInformationVector** input,
+ vtkInformationVector* output
+ );
+
+ private:
+ // Purposely not implemented
+ PolyDataToFourierSeriesFilter( const Self& );
+ void operator=( const Self& );
+
+ protected:
+ TFourierSeries m_FourierSeries;
+ };
+
+ } // ecapseman
+
+} // ecapseman
+
+#endif // __ivq__VTK__PolyDataToFourierSeriesFilter__h__
+
+// eof - $RCSfile$
git://git.creatis.insa-lyon.fr / cpPlugins.git / commitdiff