From 887d3b280e7de40d45b0483767b09dcc46121f39 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Leonardo=20Fl=C3=B3rez-Valencia?= Date: Wed, 27 Sep 2017 21:41:15 -0500 Subject: [PATCH] ... --- appli/CMakeLists.txt | 4 + appli/CTArteries/CTArteries.cxx | 305 ++++++++ appli/CTArteries/CTArteries.h | 28 + appli/CTArteries/CTArteries.ui | 18 +- appli/CTArteries/FourierSeries.cxx | 728 ++++++++++++++++++ appli/CTArteries/FourierSeries.h | 187 +++++ appli/CTArteries/FourierSeries.hxx | 102 +++ .../algorithms/RandomWalkSegmentation.hxx | 2 +- 8 files changed, 1370 insertions(+), 4 deletions(-) create mode 100644 appli/CTArteries/FourierSeries.cxx create mode 100644 appli/CTArteries/FourierSeries.h create mode 100644 appli/CTArteries/FourierSeries.hxx diff --git a/appli/CMakeLists.txt b/appli/CMakeLists.txt index dbcde1a..898efa2 100644 --- a/appli/CMakeLists.txt +++ b/appli/CMakeLists.txt @@ -2,6 +2,10 @@ ## @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co) ## ========================================================================= +include_directories( + ${PROJECT_SOURCE_DIR}/lib + ${PROJECT_BINARY_DIR}/lib + ) subdirs(CTArteries) ## eof - $RCSfile$ diff --git a/appli/CTArteries/CTArteries.cxx b/appli/CTArteries/CTArteries.cxx index 2eec3c2..345b94e 100644 --- a/appli/CTArteries/CTArteries.cxx +++ b/appli/CTArteries/CTArteries.cxx @@ -11,14 +11,21 @@ #include #include +#include #include +#include +#include +#include +#include #include #include +#include #include #include +#include #include #include #include @@ -34,6 +41,72 @@ 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 ) @@ -212,6 +285,186 @@ _process( ) 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( ) @@ -242,6 +495,16 @@ _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( ); @@ -249,6 +512,23 @@ _showProcessResults( ) 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( ); } // ------------------------------------------------------------------------- @@ -303,6 +583,31 @@ _sProcess( ) 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 diff --git a/appli/CTArteries/CTArteries.h b/appli/CTArteries/CTArteries.h index 993f6ab..4875917 100644 --- a/appli/CTArteries/CTArteries.h +++ b/appli/CTArteries/CTArteries.h @@ -15,6 +15,10 @@ #include +#include "FourierSeries.h" + +class vtkPolyData; +class vtkSplineWidget; class QDialog; namespace Ui { @@ -53,6 +57,7 @@ public: 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( @@ -68,6 +73,7 @@ protected: void _showInputImage( ); void _process( ); + void _prepareQuantification( ); void _showProcessResults( ); protected slots: @@ -75,6 +81,13 @@ 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; @@ -93,8 +106,23 @@ private: 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__ diff --git a/appli/CTArteries/CTArteries.ui b/appli/CTArteries/CTArteries.ui index 5ae4c1c..5bbd884 100644 --- a/appli/CTArteries/CTArteries.ui +++ b/appli/CTArteries/CTArteries.ui @@ -16,17 +16,23 @@ - + Qt::Vertical - + Qt::Horizontal - + + + Qt::Vertical + + + + @@ -99,6 +105,12 @@
QCustomPlot.h
1 + + ivq::Qt::ImageViewerWidget + QWidget +
ivq/Qt/ImageViewerWidget.h
+ 1 + diff --git a/appli/CTArteries/FourierSeries.cxx b/appli/CTArteries/FourierSeries.cxx new file mode 100644 index 0000000..c7bf079 --- /dev/null +++ b/appli/CTArteries/FourierSeries.cxx @@ -0,0 +1,728 @@ +#include +#include +#include + +#include +#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$ diff --git a/appli/CTArteries/FourierSeries.h b/appli/CTArteries/FourierSeries.h new file mode 100644 index 0000000..9194bfa --- /dev/null +++ b/appli/CTArteries/FourierSeries.h @@ -0,0 +1,187 @@ +/* ======================================================================= + * @author: Leonardo Florez-Valencia + * @email: florez-l@javeriana.edu.co + * ======================================================================= + */ +#ifndef __FourierSeries__h__ +#define __FourierSeries__h__ + +#include +#include +#include +#include + +#include +#include +#include + +/** + * 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$ diff --git a/appli/CTArteries/FourierSeries.hxx b/appli/CTArteries/FourierSeries.hxx new file mode 100644 index 0000000..6398da8 --- /dev/null +++ b/appli/CTArteries/FourierSeries.hxx @@ -0,0 +1,102 @@ +/* ======================================================================= + * @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$ diff --git a/appli/CTArteries/algorithms/RandomWalkSegmentation.hxx b/appli/CTArteries/algorithms/RandomWalkSegmentation.hxx index 6610d17..66e1018 100644 --- a/appli/CTArteries/algorithms/RandomWalkSegmentation.hxx +++ b/appli/CTArteries/algorithms/RandomWalkSegmentation.hxx @@ -298,7 +298,7 @@ GenerateData( ) 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 ); -- 2.47.1