## @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
## =========================================================================
+include_directories(
+ ${PROJECT_SOURCE_DIR}/lib
+ ${PROJECT_BINARY_DIR}/lib
+ )
subdirs(CTArteries)
## eof - $RCSfile$
#include <itkImageFileReader.h>
#include <itkImageSeriesReader.h>
+#include <itkQuadEdgeMesh.h>
#include <vtkMarchingCubes.h>
+#include <vtkMarchingSquares.h>
+#include <vtkPolyDataConnectivityFilter.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkProperty.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
+#include <vtkSplineWidget.h>
#include <QtGui>
#include <ivq/ITK/CPRImageFilter.h>
+#include <ivq/VTK/ImageActor.h>
#include <ivq/VTK/ImageViewer.h>
#include <ivq/VTK/MPRViewers.h>
#include <ivq/VTK/PolyDataActor.h>
this->m_UI->a##x, SIGNAL( triggered( ) ), this, SLOT( _s##x( ) ) \
)
+// -------------------------------------------------------------------------
+class ContourCallback
+ : public vtkCommand
+{
+public:
+ static ContourCallback* New( )
+ {
+ return( new ContourCallback( ) );
+ }
+ virtual void Execute(
+ vtkObject* caller, unsigned long eId, void* data
+ ) override
+ {
+ ivq::VTK::ImageActor* actor =
+ ivq::VTK::ImageActor::SafeDownCast( caller );
+ if( eId == vtkCommand::InteractionEvent && actor != NULL )
+ {
+ ivq::VTK::ImageViewer* viewer = this->Window->m_UI->Slice->GetViewer( );
+ if( this->Window->m_ContourActor != NULL )
+ viewer->GetRenderer( )->RemoveViewProp( this->Window->m_ContourActor->GetActor( ) );
+
+ double abnds[ 6 ];
+ actor->GetActorBounds( abnds );
+ unsigned int nSlice = actor->GetSliceNumber( );
+ for( unsigned int i = 0; i < 16; ++i )
+ {
+ double w = double( 6.28 ) * double( i ) / double( 16 );
+ CTArteries::TFourier::TPoint v = this->Window->m_Fourier[ nSlice ]( w );
+ this->Window->m_Spline->SetHandlePosition(
+ i, v[ 0 ], v[ 1 ], abnds[ 4 ]
+ );
+
+ } // rof
+ this->Window->m_Spline->SetProjectionPosition(
+ (
+ actor->GetBounds( )[ 4 ] +
+ actor->GetBounds( )[ 5 ]
+ ) / double( 2 )
+ );
+
+
+ this->Window->m_ContourActor =
+ this->Window->m_ContoursActors[ actor->GetSliceNumber( ) ];
+ viewer->GetRenderer( )->AddViewProp( this->Window->m_ContourActor->GetActor( ) );
+ viewer->Render( );
+
+ } // fi
+ }
+
+protected:
+ ContourCallback( )
+ : vtkCommand( )
+ {
+ }
+ virtual ~ContourCallback( )
+ {
+ }
+
+private:
+ // Purposely not implemented
+ ContourCallback( const ContourCallback& other );
+ ContourCallback& operator=( const ContourCallback& other );
+public:
+ CTArteries* Window;
+};
+
// -------------------------------------------------------------------------
CTArteries::
CTArteries( int argc, char* argv[], QWidget* parent )
this->m_Axis->DisconnectPipeline( );
}
+// -------------------------------------------------------------------------
+void CTArteries::
+_prepareQuantification( )
+{
+ this->m_Contours.clear( );
+ this->m_ContoursActors.clear( );
+ this->m_Fourier.clear( );
+
+ TImage::RegionType ext_region = this->m_CPRImage->GetRequestedRegion( );
+ for( long z = 0; z < ext_region.GetSize( )[ 2 ]; ++z )
+ {
+ TImage::IndexType slc_idx = ext_region.GetIndex( );
+ TImage::SizeType slc_size = ext_region.GetSize( );
+ slc_idx[ 2 ] += z;
+ slc_size[ 2 ] = 1;
+ TImage::RegionType slc_region;
+ slc_region.SetIndex( slc_idx );
+ slc_region.SetSize( slc_size );
+
+ typedef itk::RegionOfInterestImageFilter< TScalarImage, TScalarImage > _TROI;
+ _TROI::Pointer roi = _TROI::New( );
+ roi->SetInput( this->m_CPRSegmentation );
+ roi->SetRegionOfInterest( slc_region );
+
+ /* TODO
+ slc_idx[ 2 ] = 1;
+ slc_size[ 2 ] = 0;
+ slc_region.SetIndex( slc_idx );
+ slc_region.SetSize( slc_size );
+ typedef itk::Image< TBinPixel, 2 > _TBinSlice;
+ typedef itk::ExtractImageFilter< TBinImage, _TBinSlice > _TCollapse;
+ _TCollapse::Pointer collapse = _TCollapse::New( );
+ collapse->SetInput( roi->GetOutput( ) );
+ collapse->SetDirectionCollapseToIdentity( );
+ collapse->SetExtractionRegion( slc_region );
+ collapse->UpdateLargestPossibleRegion( );
+ */
+ typedef TScalarImage _TScalarSlice;
+
+ typedef itk::ImageToVTKImageFilter< _TScalarSlice > _TVTKScalarSlice;
+ _TVTKScalarSlice::Pointer vtk_slice = _TVTKScalarSlice::New( );
+ vtk_slice->SetInput( roi->GetOutput( ) );
+ vtk_slice->Update( );
+
+ double r[ 2 ];
+ vtk_slice->GetOutput( )->GetScalarRange( r );
+
+ vtkSmartPointer< vtkMarchingSquares > ms =
+ vtkSmartPointer< vtkMarchingSquares >::New( );
+ ms->SetInputData( vtk_slice->GetOutput( ) );
+ ms->SetValue( 0, 0 );
+
+ vtkSmartPointer< vtkPolyDataConnectivityFilter > conn =
+ vtkSmartPointer< vtkPolyDataConnectivityFilter >::New( );
+ conn->SetInputConnection( ms->GetOutputPort( ) );
+ conn->SetExtractionModeToClosestPointRegion( );
+ conn->SetClosestPoint( 0, 0, 0 );
+ conn->Update( );
+
+ vtkSmartPointer< vtkPolyData > cnt = vtkSmartPointer< vtkPolyData >::New( );
+ cnt->DeepCopy( conn->GetOutput( ) );
+ this->m_Contours.push_back( cnt );
+
+ typedef itk::QuadEdgeMesh< double, 2 > _TQEMesh;
+ typedef _TQEMesh::PointType _T2DPoint;
+ typedef std::vector< _T2DPoint > _T2DPoints;
+ _T2DPoints points;
+ _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;
+ mesh->AddFace( mesh->GetEdge( ) );
+ _TQEMesh::QEPrimal* edge = mesh->GetEdge( );
+ for( auto eIt = edge->BeginGeomLnext( ); eIt != edge->EndGeomLnext( ); ++eIt )
+ points.push_back( mesh->GetPoint( *eIt ) );
+
+ TFourier f( points.begin( ), points.end( ), 3 );
+ f.SetOrderingToCounterClockWise( );
+ this->m_Fourier.push_back( f );
+
+ vtkSmartPointer< ivq::VTK::PolyDataActor > actor =
+ vtkSmartPointer< ivq::VTK::PolyDataActor >::New( );
+ actor->SetInputData( cnt );
+ actor->GetMapper( )->ScalarVisibilityOff( );
+ actor->Update( );
+ actor->GetProperty( )->SetColor( 1, 0, 0 );
+ actor->GetProperty( )->SetLineWidth( 2 );
+ this->m_ContoursActors.push_back( actor );
+
+ } // rof
+
+ ivq::VTK::ImageViewer* slice_view = this->m_UI->Slice->GetViewer( );
+ vtkSmartPointer< ContourCallback > cbk =
+ vtkSmartPointer< ContourCallback >::New( );
+ cbk->Window = this;
+ slice_view->GetImageActor( )->AddObserver( vtkCommand::InteractionEvent, cbk );
+ this->m_ContourActor = NULL;
+
+ this->m_Spline = vtkSmartPointer< vtkSplineWidget >::New( );
+ this->m_Spline->SetCurrentRenderer( slice_view->GetRenderer( ) );
+ this->m_Spline->SetDefaultRenderer( slice_view->GetRenderer( ) );
+ this->m_Spline->SetInputData( slice_view->GetInput( ) );
+ this->m_Spline->SetProp3D( slice_view->GetImageActor( ) );
+ this->m_Spline->SetInteractor( slice_view->GetRenderWindow( )->GetInteractor( ) );
+ double bnds[ 6 ];
+ slice_view->GetInput( )->GetBounds( bnds );
+ this->m_Spline->PlaceWidget(
+ bnds[ 0 ], bnds[ 1 ],
+ bnds[ 2 ], bnds[ 3 ],
+ bnds[ 4 ], bnds[ 5 ]
+ );
+ this->m_Spline->ProjectToPlaneOn( );
+ this->m_Spline->SetProjectionNormalToZAxes( );
+ this->m_Spline->SetProjectionPosition(
+ (
+ slice_view->GetImageActor( )->GetBounds( )[ 4 ] +
+ slice_view->GetImageActor( )->GetBounds( )[ 5 ]
+ ) / double( 2 )
+ );
+ this->m_Spline->SetHandleSize( 0.005 );
+ this->m_Spline->ClosedOn( );
+ this->m_Spline->SetNumberOfHandles( 16 );
+ this->m_Spline->EnabledOn( );
+
+ QVector< double > xdata, dmindata, dmaxdata, areadata;
+ double data[ 3 ];
+ for( unsigned int i = 0; i < this->m_Fourier.size( ); ++i )
+ {
+ double a, b, t, p;
+ this->m_Fourier[ i ].GetEllipse( 1, a, b, t, p );
+ double area = this->m_Fourier[ i ].GetArea( );
+
+ xdata.push_back( double( i ) );
+ dmindata.push_back( b );
+ dmaxdata.push_back( a );
+ areadata.push_back( area );
+
+ } // rof
+ this->m_UI->Plot->addGraph( this->m_UI->Plot->xAxis, this->m_UI->Plot->yAxis );
+ this->m_UI->Plot->graph( 0 )->setPen( QPen( QColor( 255, 0, 0 ) ) );
+ this->m_UI->Plot->graph( 0 )->setName( "Minimum diameter (mm)" );
+ this->m_UI->Plot->graph( 0 )->setData( xdata, dmindata );
+
+ this->m_UI->Plot->addGraph( this->m_UI->Plot->xAxis, this->m_UI->Plot->yAxis );
+ this->m_UI->Plot->graph( 1 )->setPen( QPen( QColor( 0, 255, 0 ) ) );
+ this->m_UI->Plot->graph( 1 )->setName( "Maximum diameter (mm)" );
+ this->m_UI->Plot->graph( 1 )->setData( xdata, dmaxdata );
+
+ this->m_UI->Plot->addGraph( this->m_UI->Plot->xAxis, this->m_UI->Plot->yAxis2 );
+ this->m_UI->Plot->graph( 2 )->setPen( QPen( QColor( 0, 0, 255 ) ) );
+ this->m_UI->Plot->graph( 2 )->setName( "Area (mm2)" );
+ this->m_UI->Plot->graph( 2 )->setData( xdata, areadata );
+
+ this->m_UI->Plot->xAxis->setLabel( "Position" );
+ this->m_UI->Plot->yAxis->setLabel( "Diameters (mm)" );
+ this->m_UI->Plot->yAxis2->setLabel( "Areas (mm2)" );
+
+ this->m_UI->Plot->legend->setVisible( true );
+ this->m_UI->Plot->xAxis->setVisible( true );
+ this->m_UI->Plot->yAxis->setVisible( true );
+ this->m_UI->Plot->yAxis2->setVisible( true );
+
+ this->m_UI->Plot->rescaleAxes( );
+ this->m_UI->Plot->replot( );
+}
+
// -------------------------------------------------------------------------
void CTArteries::
_showProcessResults( )
this->m_Curve->AddPoint( this->m_Axis->GetPoint( i ) );
this->m_Curve->Smooth( 2 );
+ // Compute CPRs
+ this->_CPR(
+ this->m_CPRImage, this->m_VTKCPRImage,
+ this->m_Image, this->m_Curve
+ );
+ this->_CPR(
+ this->m_CPRSegmentation, this->m_VTKCPRSegmentation,
+ this->m_Segmentation, this->m_Curve
+ );
+
// Show all data
vtkRenderer* renderer =
this->m_UI->MPR->GetRendererWidget( )->GetRenderer( );
renderer->Render( );
renderer->ResetCamera( );
this->m_UI->MPR->GetRendererWidget( )->GetRenderWindow( )->Render( );
+
+ ivq::VTK::ImageViewer* cpr_view = this->m_UI->CPR->GetViewer( );
+ ivq::VTK::ImageViewer* slice_view = this->m_UI->Slice->GetViewer( );
+
+ cpr_view->SetInputData( this->m_VTKCPRImage->GetOutput( ) );
+ cpr_view->SetColorWindow( 1000 );
+ cpr_view->SetColorLevel( 100 );
+ cpr_view->SetSliceOrientationToXZ( );
+ cpr_view->Render( );
+
+ slice_view->SetInputData( this->m_VTKCPRImage->GetOutput( ) );
+ slice_view->SetColorWindow( 1000 );
+ slice_view->SetColorLevel( 100 );
+ slice_view->SetSliceOrientationToXY( );
+ slice_view->Render( );
+
+ this->_prepareQuantification( );
}
// -------------------------------------------------------------------------
this->_showProcessResults( );
}
+// -------------------------------------------------------------------------
+template< class _TImagePtr, class _TVTKImagePtr, class _TCurvePtr >
+void CTArteries::
+_CPR(
+ _TImagePtr& output, _TVTKImagePtr& vtk_output,
+ const _TImagePtr& input, const _TCurvePtr& curve
+ )
+{
+ typedef typename _TImagePtr::ObjectType _TImage;
+ typedef typename _TCurvePtr::ObjectType _TCurve;
+ typedef ivq::ITK::CPRImageFilter< _TImage, _TCurve > _TFilter;
+
+ typename _TFilter::Pointer cpr = _TFilter::New( );
+ cpr->SetInput( input );
+ cpr->SetCurve( curve );
+ cpr->SetSliceRadius( 40 );
+ cpr->Update( );
+ output = cpr->GetOutput( );
+ output->DisconnectPipeline( );
+
+ vtk_output = _TVTKImagePtr::ObjectType::New( );
+ vtk_output->SetInput( output );
+ vtk_output->Update( );
+}
+
// -------------------------------------------------------------------------
#include <QApplication>
#include <fpa/DataStructures/Image/Path.h>
+#include "FourierSeries.h"
+
+class vtkPolyData;
+class vtkSplineWidget;
class QDialog;
namespace Ui
{
typedef itk::ImageToVTKImageFilter< TScalarImage > TVTKScalarImage;
typedef fpa::DataStructures::Image::Path< Dim > TAxis;
typedef ivq::ITK::Simple3DCurve< TScalar > TCurve;
+ typedef FourierSeries< TScalar, 2 > TFourier;
public:
explicit CTArteries(
void _showInputImage( );
void _process( );
+ void _prepareQuantification( );
void _showProcessResults( );
protected slots:
void _sConfig( );
void _sProcess( );
+private:
+ template< class _TImagePtr, class _TVTKImagePtr, class _TCurvePtr >
+ void _CPR(
+ _TImagePtr& output, _TVTKImagePtr& vtk_output,
+ const _TImagePtr& input, const _TCurvePtr& curve
+ );
+
private:
Ui::CTArteries* m_UI;
TCurve::Pointer m_Curve;
vtkSmartPointer< ivq::VTK::PolyDataActor > m_Surface;
+ // CPRs
+ TImage::Pointer m_CPRImage;
+ TScalarImage::Pointer m_CPRSegmentation;
+ TVTKImage::Pointer m_VTKCPRImage;
+ TVTKScalarImage::Pointer m_VTKCPRSegmentation;
+
+ // Quantification data
+ std::vector< vtkSmartPointer< vtkPolyData > > m_Contours;
+ std::vector< vtkSmartPointer< ivq::VTK::PolyDataActor > > m_ContoursActors;
+ ivq::VTK::PolyDataActor* m_ContourActor;
+ std::vector< TFourier > m_Fourier;
+
// Widgets
vtkSmartPointer< ivq::VTK::SeedWidgetOverImageActor > m_Seeds;
+ vtkSmartPointer< vtkSplineWidget > m_Spline;
+
+ friend class ContourCallback;
};
#endif // __CTArteries__h__
<widget class="QWidget" name="centralwidget">
<layout class="QGridLayout" name="gridLayout">
<item row="0" column="0">
- <widget class="QSplitter" name="splitter_2">
+ <widget class="QSplitter" name="splitter_3">
<property name="orientation">
<enum>Qt::Vertical</enum>
</property>
<widget class="ivq::Qt::MPRViewersWidget" name="MPR" native="true"/>
- <widget class="QSplitter" name="splitter">
+ <widget class="QSplitter" name="splitter_2">
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
<widget class="QCustomPlot" name="Plot" native="true"/>
- <widget class="QWidget" name="CPR" native="true"/>
+ <widget class="QSplitter" name="splitter">
+ <property name="orientation">
+ <enum>Qt::Vertical</enum>
+ </property>
+ <widget class="ivq::Qt::ImageViewerWidget" name="Slice" native="true"/>
+ <widget class="ivq::Qt::ImageViewerWidget" name="CPR" native="true"/>
+ </widget>
</widget>
</widget>
</item>
<header>QCustomPlot.h</header>
<container>1</container>
</customwidget>
+ <customwidget>
+ <class>ivq::Qt::ImageViewerWidget</class>
+ <extends>QWidget</extends>
+ <header location="global">ivq/Qt/ImageViewerWidget.h</header>
+ <container>1</container>
+ </customwidget>
</customwidgets>
<resources/>
<connections>
--- /dev/null
+#include <cmath>
+#include <limits>
+#include <map>
+
+#include <vnl/vnl_math.h>
+#include "FourierSeries.h"
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+FourierSeries< S, D >::
+FourierSeries( unsigned int q )
+ : m_RealAxis( D - 2 ),
+ m_ImagAxis( D - 1 )
+{
+ this->resize( ( q << 1 ) + 1, TComplex( S( 0 ), S( 0 ) ) );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+FourierSeries< S, D >::
+~FourierSeries( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+const unsigned int& FourierSeries< S, D >::
+GetRealAxis( ) const
+{
+ return( this->m_RealAxis );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+const unsigned int& FourierSeries< S, D >::
+GetImagAxis( ) const
+{
+ return( this->m_ImagAxis );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void FourierSeries< S, D >::
+SetRealAxis( const unsigned int& a )
+{
+ this->m_RealAxis = a;
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void FourierSeries< S, D >::
+SetImagAxis( const unsigned int& a )
+{
+ this->m_ImagAxis = a;
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+S FourierSeries< S, D >::
+GetArea( ) const
+{
+ S a = S( 0 );
+ int q = this->GetNumberOfHarmonics( );
+ typename Self::const_iterator i = this->begin( );
+ for( int l = -q; i != this->end( ); ++i, ++l )
+ a += S( l ) * std::norm( *i );
+ return( S( vnl_math::pi ) * a );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+S FourierSeries< S, D >::
+GetCurvature( const S& w ) const
+{
+ TComplex d1 = this->_Z( w, 1 );
+ TComplex d2 = this->_Z( w, 2 );
+ S d = std::abs( d1 );
+ if( d > S( 0 ) )
+ return( std::imag( std::conj( d1 ) * d2 ) / ( d * d * d ) );
+ else
+ return( S( 0 ) );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void FourierSeries< S, D >::
+SetOrdering( bool cw )
+{
+ // Check if the area sign is different from the desired orientation
+ bool ap = ( S( 0 ) < this->GetArea( ) );
+ if( ( !ap && cw ) || ( ap && !cw ) ) // XOR
+ this->InvertOrdering( );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void FourierSeries< S, D >::
+InvertOrdering( )
+{
+ int q = this->GetNumberOfHarmonics( );
+ for( int l = 1; l <= q; l++ )
+ {
+ TComplex z = ( *this )[ l ];
+ ( *this )[ l ] = ( *this )[ -l ];
+ ( *this )[ -l ] = z;
+
+ } // rof
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void FourierSeries< S, D >::
+SetOrderingToClockWise( )
+{
+ this->SetOrdering( false );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void FourierSeries< S, D >::
+SetOrderingToCounterClockWise( )
+{
+ this->SetOrdering( true );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+int FourierSeries< S, D >::
+GetNumberOfHarmonics( ) const
+{
+ return( ( this->size( ) - 1 ) >> 1 );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void FourierSeries< S, D >::
+SetNumberOfHarmonics( const int& q )
+{
+ int diff_q = this->GetNumberOfHarmonics( ) - q;
+
+ while( diff_q != 0 )
+ {
+ if( diff_q > 0 )
+ {
+ this->pop_front( );
+ this->pop_back( );
+ diff_q--;
+ }
+ else
+ {
+ this->push_front( TComplex( S( 0 ), S( 0 ) ) );
+ this->push_back( TComplex( S( 0 ), S( 0 ) ) );
+ diff_q++;
+
+ } // fi
+
+ } // elihw
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+TComplex& FourierSeries< S, D >::
+operator[]( const int& l )
+{
+ static TComplex _zero;
+ int idx = ( ( this->size( ) - 1 ) >> 1 ) + l;
+ if( idx < 0 || idx >= this->size( ) )
+ {
+ _zero = TComplex( S( 0 ) );
+ return( _zero );
+ }
+ else
+ return( this->Superclass::operator[]( idx ) );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+const typename FourierSeries< S, D >::
+TComplex& FourierSeries< S, D >::
+operator[]( const int& l ) const
+{
+ static const TComplex _zero( S( 0 ) );
+ int idx = ( ( this->size( ) - 1 ) >> 1 ) + l;
+ if( idx < 0 || idx >= this->size( ) )
+ return( _zero );
+ else
+ return( this->Superclass::operator[]( idx ) );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+TPoint FourierSeries< S, D >::
+operator( )( const S& w ) const
+{
+ TComplex z = this->_Z( w, 0 );
+ TPoint p;
+ p.Fill( S( 0 ) );
+ p[ this->m_RealAxis ] = z.real( );
+ p[ this->m_ImagAxis ] = z.imag( );
+ return( p );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+TVector FourierSeries< S, D >::
+operator( )( const S& w, const unsigned int& n ) const
+{
+ TComplex z = this->_Z( w, n );
+ TVector v( S( 0 ) );
+ v[ this->m_RealAxis ] = z.real( );
+ v[ this->m_ImagAxis ] = z.imag( );
+ return( v );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self FourierSeries< S, D >::
+operator+( const Self& o ) const
+{
+ int q1 = o.GetNumberOfHarmonics( );
+ int q2 = this->GetNumberOfHarmonics( );
+ int q = ( q2 < q1 )? q1: q2;
+
+ Self res;
+ res.SetNumberOfHarmonics( q );
+ for( int l = -q; l <= q; ++l )
+ res[ l ] = ( *this )[ l ] + o[ l ];
+
+ return( res );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self& FourierSeries< S, D >::
+operator+=( const Self& o )
+{
+ int q1 = o.GetNumberOfHarmonics( );
+ int q2 = this->GetNumberOfHarmonics( );
+ int q = ( q2 < q1 )? q1: q2;
+
+ this->SetNumberOfHarmonics( q );
+ for( int l = -q; l <= q; ++l )
+ ( *this )[ l ] += o[ l ];
+
+ return( *this );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self FourierSeries< S, D >::
+operator-( const Self& o ) const
+{
+ int q1 = o.GetNumberOfHarmonics( );
+ int q2 = this->GetNumberOfHarmonics( );
+ int q = ( q2 < q1 )? q1: q2;
+
+ Self res;
+ res.SetNumberOfHarmonics( q );
+ for( int l = -q; l <= q; ++l )
+ res[ l ] = ( *this )[ l ] - o[ l ];
+
+ return( res );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self& FourierSeries< S, D >::
+operator-=( const Self& o )
+{
+ int q1 = o.GetNumberOfHarmonics( );
+ int q2 = this->GetNumberOfHarmonics( );
+ int q = ( q2 < q1 )? q1: q2;
+
+ this->SetNumberOfHarmonics( q );
+ for( int l = -q; l <= q; ++l )
+ ( *this )[ l ] -= o[ l ];
+
+ return( *this );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self FourierSeries< S, D >::
+operator+( const TComplex& translation ) const
+{
+ Self res = *this;
+ res[ 0 ] += translation;
+ return( res );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self& FourierSeries< S, D >::
+operator+=( const TComplex& translation )
+{
+ ( *this )[ 0 ] += translation;
+ return( *this );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self FourierSeries< S, D >::
+operator-( const TComplex& translation ) const
+{
+ Self res = *this;
+ res[ 0 ] -= translation;
+ return( res );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self& FourierSeries< S, D >::
+operator-=( const TComplex& translation )
+{
+ ( *this )[ 0 ] -= translation;
+ return( *this );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self FourierSeries< S, D >::
+operator+( const TVector& translation ) const
+{
+ return(
+ ( *this ) +
+ TComplex( translation[ this->m_RealAxis ], translation[ this->m_ImagAxis ] )
+ );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self& FourierSeries< S, D >::
+operator+=( const TVector& translation )
+{
+ ( *this ) +=
+ TComplex( translation[ this->m_RealAxis ], translation[ this->m_ImagAxis ] );
+ return( *this );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self FourierSeries< S, D >::
+operator-( const TVector& translation ) const
+{
+ return(
+ ( *this ) -
+ TComplex( translation[ this->m_RealAxis ], translation[ this->m_ImagAxis ] )
+ );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self& FourierSeries< S, D >::
+operator-=( const TVector& translation )
+{
+ ( *this ) -=
+ TComplex( translation[ this->m_RealAxis ], translation[ this->m_ImagAxis ] );
+ return( *this );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+TPoint FourierSeries< S, D >::
+GetCenter( ) const
+{
+ TComplex z = ( *this )[ 0 ];
+ TPoint p;
+ p.Fill( S( 0 ) );
+ p[ this->m_RealAxis ] = z.real( );
+ p[ this->m_ImagAxis ] = z.imag( );
+ return( p );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self FourierSeries< S, D >::
+operator*( const TComplex& scale_or_rotation ) const
+{
+ Self res;
+ res.clear( );
+ for( typename Self::const_iterator i = this->begin( ); i != this->end( ); ++i )
+ res.push_back( *i * scale_or_rotation );
+ return( res );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self& FourierSeries< S, D >::
+operator*=( const TComplex& scale_or_rotation )
+{
+ for( typename Self::iterator i = this->begin( ); i != this->end( ); ++i )
+ *i *= scale_or_rotation;
+ return( *this );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self FourierSeries< S, D >::
+operator/( const S& scale ) const
+{
+ Self res;
+ res.clear( );
+ for( typename Self::const_iterator i = this->begin( ); i != this->end( ); ++i )
+ res.push_back( *i / scale );
+ return( res );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self& FourierSeries< S, D >::
+operator/=( const S& scale )
+{
+ for( typename Self::iterator i = this->begin( ); i != this->end( ); ++i )
+ *i /= scale;
+ return( *this );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self FourierSeries< S, D >::
+operator&( const S& phase ) const
+{
+ Self res;
+ res.clear( );
+ int q = this->GetNumberOfHarmonics( );
+ typename Self::const_iterator i = this->begin( );
+ for( int l = -q; i != this->end( ); ++i, ++l )
+ res.push_back( *i * std::polar( S( 1 ), S( l ) * phase ) );
+ return( res );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+Self& FourierSeries< S, D >::
+operator&=( const S& phase )
+{
+ int q = this->GetNumberOfHarmonics( );
+ typename Self::iterator i = this->begin( );
+ for( int l = -q; i != this->end( ); ++i, ++l )
+ *i *= std::polar( S( 1 ), S( l ) * phase );
+ return( *this );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+S FourierSeries< S, D >::
+GetPhase( const S& w ) const
+{
+ // Some constant values
+ static const S _0 = S( 0 );
+ static const S _1 = S( 1 );
+ static const S _2 = S( 2 );
+ static const S _2pi = _2 * S( vnl_math::pi );
+ static const S _epsilon = S( 1e-14 );
+ static const S _tolerance = S( 1e-10 );
+ static const unsigned int _samples = 100;
+ static const unsigned int _maxIterations = 1000;
+ static const S _angleOff = _2pi / S( _samples );
+
+ // Some other values
+ int q = int( this->GetNumberOfHarmonics( ) );
+ if( q == 0 )
+ return( _0 );
+
+ // Get a centered copy
+ Self contour = *this;
+ contour[ 0 ] = TComplex( _0, _0 );
+
+ // Compute maximum coefficient norm
+ S max_norm = _0;
+ for( int l = 1; l <= q; ++l )
+ {
+ S np = std::abs( contour[ l ] );
+ S nn = std::abs( contour[ -l ] );
+ S n = ( np < nn )? nn: np;
+ if( max_norm < n )
+ max_norm = n;
+
+ } // rof
+
+ // Rotate to desired phase
+ contour *= std::polar( _1, -w );
+
+ // Try to normalize contour: if not, malformed contour!
+ if( max_norm > _0 )
+ {
+ contour /= max_norm;
+
+ // 1. Approximate initial guess by a simple dichotomy
+ std::vector< std::pair< S, S > > function;
+ S minA = std::numeric_limits< S >::max( );
+ unsigned int minIdx = -1;
+
+ // 1.1. Roughly sample phases
+ for( unsigned int s = 0; s < _samples; ++s )
+ {
+ S w = S( s ) * _angleOff;
+ S a = std::arg( contour._Z( w, 0 ) );
+ function.push_back( std::pair< S, S >( w, a ) );
+ if( a < minA )
+ {
+ minA = a;
+ minIdx = s;
+
+ } // fi
+
+ } // rof
+
+ // 1.2. Get zero cuts by zero crossing analysis and keep the farthest
+ // point in the real axis (ie. the last data in "cuts")
+ S prevA = _1;
+ std::map< S, unsigned int > cuts;
+ for( unsigned int s = 0; s < _samples; ++s )
+ {
+ unsigned int i = ( s + minIdx ) % _samples;
+ if( function[ i ].second * prevA < _0 )
+ cuts[ std::real( contour._Z( function[ i ].first, 0 ) ) ] = i;
+ prevA = function[ i ].second;
+
+ } // rof
+
+ // 1.3. Get initial guess
+ S w0 = _0;
+ if( cuts.size( ) > 0 )
+ {
+ unsigned int i = cuts.rbegin( )->second;
+ if( i == 0 )
+ i = function.size( );
+ w0 = function[ i - 1 ].first;
+
+ } // fi
+
+ // 2. Refine by Newthon-Raphson
+ bool stop = false;
+ unsigned int i = 0;
+ while( i < _maxIterations && !stop )
+ {
+ TComplex c = contour._Z( w0, 0 );
+ S d = std::imag( c * std::conj( contour._Z( w0, 1 ) ) );
+
+ // Emergency stop!!!
+ if( !( std::fabs( d ) < _epsilon ) )
+ {
+ S w1 = w0 + ( std::arg( c ) * ( std::norm( c ) / d ) );
+ if( ( std::fabs( w1 - w0 ) / std::fabs( w1 ) ) < _tolerance )
+ stop = true;
+ else
+ w0 = w1;
+ i++;
+ }
+ else
+ stop = true;
+
+ } // elihw
+ return( w0 );
+ }
+ else
+ return( _0 );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void FourierSeries< S, D >::
+Sample( std::vector< TPoint >& p, const unsigned long& s ) const
+{
+ static const S _2pi = S( 2 ) * S( vnl_math::pi );
+ S off = _2pi / S( s - 1 );
+ for( unsigned int w = 0; w < s; ++w )
+ p.push_back( ( *this )( S( w ) * off ) );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void FourierSeries< S, D >::
+GetEllipse( int l, S& a, S& b, S& t, S& p ) const
+{
+ a = b = t = p = S( 0 );
+ if( l == 0 || l > this->GetNumberOfHarmonics( ) )
+ return;
+
+ TComplex zp = ( *this )[ l ];
+ TComplex zn = ( *this )[ -l ];
+ S np = std::abs( zp );
+ S nn = std::abs( zn );
+
+ // Radii
+ a = np + nn;
+ b = np - nn;
+
+ // Rotation and phase
+ if( np > S( 0 ) && nn > S( 0 ) )
+ {
+ zp /= np;
+ zn /= nn;
+ t = std::real( std::log( zp * zn ) / TComplex( S( 0 ), S( 2 ) ) );
+ p = std::real( std::log( zp / zn ) / TComplex( S( 0 ), S( 2 * l ) ) );
+ }
+ else
+ {
+ t = S( 0 );
+ p = ( np > S( 0 ) )? std::arg( zp ): std::arg( zn );
+
+ } // fi
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void FourierSeries< S, D >::
+SetEllipse( int l, S& a, S& b, S& t, S& p )
+{
+ S nzp = ( a + b ) / S( 2 );
+ S nzn = ( a - b ) / S( 2 );
+ TComplex et = std::polar( S( 1 ), t );
+ TComplex ep = std::polar( S( 1 ), S( l ) * p );
+ ( *this )[ l ] = et * ep * nzp;
+ ( *this )[ -l ] = et * std::conj( ep ) * nzn;
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+typename FourierSeries< S, D >::
+TComplex FourierSeries< S, D >::
+_Z( const S& w, const unsigned int& n ) const
+{
+ static const S _0 = S( 0 );
+ static const S _1 = S( 1 );
+ S _n = S( n );
+
+ TComplex z( _0, _0 );
+ int q = this->GetNumberOfHarmonics( );
+ for( int l = -q; l <= q; ++l )
+ {
+ TComplex v = ( *this )[ l ] * std::polar( _1, w * S( l ) );
+ if( n > 0 )
+ v *= std::pow( TComplex( _0, S( l ) ), _n );
+ z += v;
+
+ } // rof
+ return( z );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void FourierSeries< S, D >::
+_DFT( const std::vector< TComplex >& p, unsigned int q )
+{
+ static const S _2pi = S( 2 ) * S( vnl_math::pi );
+
+ this->SetNumberOfHarmonics( q );
+ *this *= S( 0 );
+ if( p.size( ) == 0 )
+ return;
+
+ std::vector< TComplex > dft;
+ for( long m = 0; m < p.size( ); ++m )
+ {
+ TComplex z( S( 0 ), S( 0 ) );
+ for( long k = 0; k < p.size( ); ++k )
+ z +=
+ p[ k ] *
+ std::polar( S( 1 ), -( _2pi * S( m * k ) ) / S( p.size( ) ) );
+ z /= S( p.size( ) );
+ dft.push_back( z );
+
+ } // rof
+
+ ( *this )[ 0 ] = dft[ 0 ];
+ for( int l = 1; l <= q; ++l )
+ {
+ ( *this )[ l ] = dft[ l ];
+ ( *this )[ -l ] = dft[ p.size( ) - l ];
+
+ } // rof
+
+ // Reset contour
+ /* TODO:
+ this->SetNumberOfHarmonics( q );
+
+ // Compute number of calculations K: cut harmonics to q
+ long N = long( p.size( ) );
+ if( N == 0 )
+ return;
+ long K = ( long( q ) << 1 ) + 1;
+ if( K < N - 1 )
+ K = N - 1;
+
+ // Compute harmonics
+ for( long l = -K; l <= K; ++l )
+ {
+ S k = ( ( S( l ) + S( ( l < 0 )? N: 0 ) ) * _2pi ) / S( N );
+ ( *this )[ l ] = TComplex( S( 0 ), S( 0 ) );
+ for( long n = 0; n < N; ++n )
+ ( *this )[ l ] += p[ n ] * std::polar( S( 1 ), k * S( -n ) );
+ ( *this )[ l ] /= S( N );
+
+ } // rof
+ */
+}
+
+// -------------------------------------------------------------------------
+template class FourierSeries< float, 2 >;
+template class FourierSeries< float, 3 >;
+template class FourierSeries< double, 2 >;
+template class FourierSeries< double, 3 >;
+
+// eof - $RCSfile$
--- /dev/null
+/* =======================================================================
+ * @author: Leonardo Florez-Valencia
+ * @email: florez-l@javeriana.edu.co
+ * =======================================================================
+ */
+#ifndef __FourierSeries__h__
+#define __FourierSeries__h__
+
+#include <complex>
+#include <deque>
+#include <vector>
+#include <utility>
+
+#include <itkMatrix.h>
+#include <itkPoint.h>
+#include <itkVector.h>
+
+/**
+ * Planar series defined by its harmonic ellipses.
+ *
+ * The tuple [z(l),z(-l)], l != 0 defines an ellipse. The z(0)
+ * harmonic defines the series's center.
+ *
+ */
+template< class S, unsigned int D >
+class FourierSeries
+ : public std::deque< std::complex< S > >
+{
+public:
+ typedef FourierSeries Self;
+ typedef S TScalar;
+ itkStaticConstMacro( Dimension, unsigned int, D );
+
+ typedef itk::Matrix< S, D, D > TMatrix;
+ typedef itk::Point< S, D > TPoint;
+ typedef typename TPoint::VectorType TVector;
+ typedef std::complex< S > TComplex;
+ typedef std::deque< TComplex > Superclass;
+
+public:
+ FourierSeries( unsigned int q = 1 );
+ template< class S2, unsigned int D2 >
+ FourierSeries( const FourierSeries< S2, D2 >& o );
+ template< class _TIt >
+ FourierSeries( const _TIt& b, const _TIt& e, unsigned int q );
+ virtual ~FourierSeries( );
+
+ template< class S2, unsigned int D2 >
+ Self& operator=( const FourierSeries< S2, D2 >& o );
+
+ bool operator==( const Self& o ) const { return( false ); }
+ bool operator!=( const Self& o ) const { return( true ); }
+
+ const unsigned int& GetRealAxis( ) const;
+ const unsigned int& GetImagAxis( ) const;
+ void SetRealAxis( const unsigned int& a );
+ void SetImagAxis( const unsigned int& a );
+
+ /**
+ * $A=\pi\sum_{l=-q}^{q}l|z_{l}|^{2}$
+ */
+ S GetArea( ) const;
+
+ /**
+ * $\kappa(\omega)=\frac{|c'(\omega)\times c''(\omega)|}{|c'(\omega)|^{3}}$
+ */
+ S GetCurvature( const S& w ) const;
+
+ /**
+ * Sets the order of the ellipse. This order is defined as the
+ * counterclock/clock wiseness of the series. Calculations are
+ * done by computing the series's area (\see Area) and comparing
+ * its sign.
+ *
+ * \param cw counter-clock
+ */
+ void SetOrdering( bool cw );
+ void InvertOrdering( );
+ void SetOrderingToClockWise( );
+ void SetOrderingToCounterClockWise( );
+
+ /**
+ * Series manipulators.
+ */
+ int GetNumberOfHarmonics( ) const;
+ void SetNumberOfHarmonics( const int& q );
+ TComplex& operator[]( const int& l );
+ const TComplex& operator[]( const int& l ) const;
+
+ /**
+ * $c(\omega)=\sum_{l=-q}^{q}z_{l}e^{jl\omega}$
+ *
+ * where $0\le\omega<2\pi$ is the angular curvilinear parameter, q is
+ * the number of harmonics defining the series, $z_{l}\in\mathbb{C}$
+ * is the l-th harmonic and $j$ is the imaginary unity.
+ */
+ TPoint operator( )( const S& w ) const;
+
+ /**
+ * $\frac{d^{(n)}c(\omega)}{d\omega^{(n)}}=\sum_{l=-q}^{q}z_{l}e^{jl\omega}(jl)^{n}$
+ *
+ * where $n$ is the derivative number, $0\le\omega<2\pi$ is the
+ * angular curvilinear parameter, q is the number of harmonics
+ * defining the series, $z_{l}\in\mathbb{C}$ is the l-th harmonic
+ * and $j$ is the imaginary unity.
+ */
+ TVector operator( )( const S& w, const unsigned int& n ) const;
+
+ /**
+ * Coefficient-wise addition operators
+ */
+ Self operator+( const Self& o ) const;
+ Self& operator+=( const Self& o );
+ Self operator-( const Self& o ) const;
+ Self& operator-=( const Self& o );
+
+ /**
+ * Translation
+ */
+ Self operator+( const TComplex& translation ) const;
+ Self& operator+=( const TComplex& translation );
+ Self operator-( const TComplex& translation ) const;
+ Self& operator-=( const TComplex& translation );
+ Self operator+( const TVector& translation ) const;
+ Self& operator+=( const TVector& translation );
+ Self operator-( const TVector& translation ) const;
+ Self& operator-=( const TVector& translation );
+ TPoint GetCenter( ) const;
+
+ /**
+ * Scale and/or rotate
+ */
+ Self operator*( const TComplex& scale_or_rotation ) const;
+ Self& operator*=( const TComplex& scale_or_rotation );
+ Self operator/( const S& scale ) const;
+ Self& operator/=( const S& scale );
+
+ /**
+ * Phase
+ */
+ Self operator&( const S& phase ) const;
+ Self& operator&=( const S& phase );
+ S GetPhase( const S& w = S( 0 ) ) const;
+
+ /**
+ * Sampling method
+ */
+ void Sample( std::vector< TPoint >& p, const unsigned long& s ) const;
+
+ /**
+ * Ellipse methods
+ */
+ void GetEllipse( int l, S& a, S& b, S& t, S& p ) const;
+ void SetEllipse( int l, S& a, S& b, S& t, S& p );
+
+ template< class _TOtherScalar >
+ void SetEllipses( const std::vector< _TOtherScalar >& ellipses );
+
+protected:
+ // Main method to compute all of series's values.
+ TComplex _Z( const S& w, const unsigned int& n ) const;
+
+ // Fourier transform
+ void _DFT( const std::vector< TComplex >& p, unsigned int q );
+
+protected:
+ unsigned int m_RealAxis;
+ unsigned int m_ImagAxis;
+
+public:
+ friend std::ostream& operator<<( std::ostream& out, const Self& fs )
+ {
+ int q = fs.GetNumberOfHarmonics( );
+ int l = -q;
+ out << "{" << l << ":" << fs[ l ] << "}";
+ for( ++l; l <= q; l++ )
+ out << ", {" << l << ":" << fs[ l ] << "}";
+ return( out );
+ }
+};
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "FourierSeries.hxx"
+#endif // ITK_MANUAL_INSTANTIATION
+#endif // __FourierSeries__h__
+
+// eof - $RCSfile$
--- /dev/null
+/* =======================================================================
+ * @author: Leonardo Florez-Valencia
+ * @email: florez-l@javeriana.edu.co
+ * =======================================================================
+ */
+#ifndef __FourierSeries__hxx__
+#define __FourierSeries__hxx__
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+template< class S2, unsigned int D2 >
+FourierSeries< S, D >::
+FourierSeries( const FourierSeries< S2, D2 >& o )
+ : m_RealAxis( D - 2 ),
+ m_ImagAxis( D - 1 )
+{
+ int q = o.GetNumberOfHarmonics( );
+ this->SetNumberOfHarmonics( q );
+ for( int l = -q; l <= q; ++l )
+ ( *this )[ l ] =
+ TComplex( S( std::real( o[ l ] ) ), S( std::imag( o[ l ] ) ) );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+template< class _TIt >
+FourierSeries< S, D >::
+FourierSeries( const _TIt& b, const _TIt& e, unsigned int q )
+ : m_RealAxis( D - 2 ),
+ m_ImagAxis( D - 1 )
+{
+ std::vector< TComplex > aux;
+ for( _TIt it = b; it != e; ++it )
+ aux.push_back( TComplex( S( ( *it )[ 0 ] ), S( ( *it )[ 1 ] ) ) );
+ this->_DFT( aux, q );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+template< class S2, unsigned int D2 >
+typename FourierSeries< S, D >::
+Self& FourierSeries< S, D >::
+operator=( const FourierSeries< S2, D2 >& o )
+{
+ if( D == D2 )
+ {
+ this->m_RealAxis = o.GetRealAxis( );
+ this->m_ImagAxis = o.GetImagAxis( );
+ }
+ else
+ {
+ this->m_RealAxis = D - 2;
+ this->m_ImagAxis = D - 1;
+
+ } // fi
+ int q = o.GetNumberOfHarmonics( );
+ this->SetNumberOfHarmonics( q );
+ for( int l = -q; l <= q; ++l )
+ ( *this )[ l ] =
+ TComplex( S( std::real( o[ l ] ) ), S( std::imag( o[ l ] ) ) );
+ return( *this );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+template< class _TOtherScalar >
+void FourierSeries< S, D >::
+SetEllipses( const std::vector< _TOtherScalar >& ellipses )
+{
+ int Q = ( ellipses.size( ) - 2 ) >> 2;
+ Q = ( Q < 1 )? 1: Q;
+ this->SetNumberOfHarmonics( Q );
+ auto cIt = ellipses.begin( );
+ S cx = ( cIt != ellipses.end( ) )? *cIt: S( 0 );
+ if( cIt != ellipses.end( ) )
+ ++cIt;
+ S cy = ( cIt != ellipses.end( ) )? *cIt: S( 0 );
+ if( cIt != ellipses.end( ) )
+ ++cIt;
+ ( *this )[ 0 ] = TComplex( cx, cy );
+ for( int l = 1; l <= Q; ++l )
+ {
+ S a = ( cIt != ellipses.end( ) )? *cIt: S( 0 );
+ if( cIt != ellipses.end( ) )
+ ++cIt;
+ S b = ( cIt != ellipses.end( ) )? *cIt: S( 0 );
+ if( cIt != ellipses.end( ) )
+ ++cIt;
+ S t = ( cIt != ellipses.end( ) )? *cIt: S( 0 );
+ if( cIt != ellipses.end( ) )
+ ++cIt;
+ S p = ( cIt != ellipses.end( ) )? *cIt: S( 0 );
+ if( cIt != ellipses.end( ) )
+ ++cIt;
+ this->SetEllipse( l, a, b, t, p );
+
+ } // rof
+}
+
+#endif // __FourierSeries__hxx__
+
+// eof - $RCSfile$
std::cout << "3" << std::endl;
{ // begin
typename _TGaussFun::Pointer rw_fun = _TGaussFun::New( );
- rw_fun->SetBeta( this->m_Beta /*init_std / double( 2 )*/ );
+ rw_fun->SetBeta( init_std / double( 2 ) );
typename _TRandomWalker::Pointer rw = _TRandomWalker::New( );
rw->SetInputImage( input_roi );
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff