From: Leonardo Flórez-Valencia Date: Tue, 24 Oct 2017 17:39:19 +0000 (-0500) Subject: ... X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=7cb9afab15499837e3ef133669cd779ae44bedd6;p=cpPlugins.git ... --- diff --git a/lib/ivq/ITK/FourierSeries.cxx b/lib/ivq/ITK/FourierSeries.cxx index 04c9c2b..7fcadf0 100644 --- a/lib/ivq/ITK/FourierSeries.cxx +++ b/lib/ivq/ITK/FourierSeries.cxx @@ -487,7 +487,7 @@ GetPhase( const TScalar& w ) const int q = int( this->GetNumberOfHarmonics( ) ); if( q == 0 ) return( _0 ); - + // Get a centered copy Self contour = *this; contour[ 0 ] = TComplex( _0, _0 ); @@ -670,57 +670,61 @@ template< class _TScalar, unsigned int _VDim > 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 - */ } // ------------------------------------------------------------------------- diff --git a/lib/ivq/ITK/FourierSeries.hxx b/lib/ivq/ITK/FourierSeries.hxx index 3ed11c5..e45ece0 100644 --- a/lib/ivq/ITK/FourierSeries.hxx +++ b/lib/ivq/ITK/FourierSeries.hxx @@ -56,7 +56,10 @@ operator=( const FourierSeries< _TScalar2, _VDim2 >& o ) this->SetNumberOfHarmonics( q ); for( int l = -q; l <= q; ++l ) ( *this )[ l ] = - TComplex( TScalar( std::real( o[ l ] ) ), TScalar( std::imag( o[ l ] ) ) ); + TComplex( + TScalar( std::real( o[ l ] ) ), + TScalar( std::imag( o[ l ] ) ) + ); return( *this ); } diff --git a/lib/ivq/VTK/PolyDataToFourierSeriesFilter.cxx b/lib/ivq/VTK/PolyDataToFourierSeriesFilter.cxx new file mode 100644 index 0000000..b9d6a6a --- /dev/null +++ b/lib/ivq/VTK/PolyDataToFourierSeriesFilter.cxx @@ -0,0 +1,151 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= + +#include +#include +#include +#include + +// ------------------------------------------------------------------------- +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 + +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$ diff --git a/lib/ivq/VTK/PolyDataToFourierSeriesFilter.h b/lib/ivq/VTK/PolyDataToFourierSeriesFilter.h new file mode 100644 index 0000000..60cf5b2 --- /dev/null +++ b/lib/ivq/VTK/PolyDataToFourierSeriesFilter.h @@ -0,0 +1,64 @@ +// ========================================================================= +// @author Leonardo Florez Valencia +// @email florez-l@javeriana.edu.co +// ========================================================================= +#ifndef __ivq__VTK__PolyDataToFourierSeriesFilter__h__ +#define __ivq__VTK__PolyDataToFourierSeriesFilter__h__ + +#include + +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$