1 #ifndef __FPA__IMAGE__FUNCTORS__MFLUXMEDIALNESS__HXX__
2 #define __FPA__IMAGE__FUNCTORS__MFLUXMEDIALNESS__HXX__
5 #include <vnl/vnl_math.h>
6 #include <itkLineConstIterator.h>
8 // -------------------------------------------------------------------------
9 template< class _TGradient >
10 fpa::Image::Functors::MFluxMedialness< _TGradient >::
13 m_MinRadius( double( 0 ) ),
14 m_MaxRadius( double( 1 ) ),
15 m_RadialSampling( 4 ),
16 m_RadiusStep( double( 1 ) )
20 // -------------------------------------------------------------------------
21 template< class _TGradient >
22 fpa::Image::Functors::MFluxMedialness< _TGradient >::
27 // -------------------------------------------------------------------------
28 template< class _TGradient >
29 typename fpa::Image::Functors::MFluxMedialness< _TGradient >::
30 TOutput fpa::Image::Functors::MFluxMedialness< _TGradient >::
31 _Evaluate( const TIndex& i ) const
33 itk::Object::GlobalWarningDisplayOff( );
35 double pi2n = double( 2 ) * double( vnl_math::pi );
36 pi2n /= int( this->m_RadialSampling );
37 const _TGradient* img = this->GetInputImage( );
38 //const itk::Image::SpacingType& input_spacing = img->GetSpacing( );
44 TRCandidates FluxFinal;
45 TRCandidates radiusGenerated;
46 double dR = double( 0 );
47 double optR = double( 0 );
49 img->TransformIndexToPhysicalPoint( i, center );
52 for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ )
54 for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ )
58 radiusGenerated.clear();
59 radius = this->m_MinRadius;
60 while( radius <= this->m_MaxRadius )
63 for( unsigned int I_radial = 0; I_radial < this->m_RadialSampling / 2; I_radial++ )
68 // Direction of first profile
69 typename TPoint::VectorType dir1;
70 dir1.Fill( double( 0 ) );
71 dir1[ cx ] = std::cos( pi2n * double( I_radial ) );
72 dir1[ cy ] = std::sin( pi2n * double( I_radial ) );
77 if ( img->TransformPhysicalPointToIndex( center + (dir1*radius), rIdx ) )
79 TVector grad_rIdx = img->GetPixel( rIdx );
81 u_i1.SetVnlVector( ( center - ( center + dir1 ) ).GetVnlVector( ) );
84 Flux1 = grad_rIdx * u_i1;
88 //if (Self::Dimension==3)
90 //std::cout<<"Point Edge x:"<<center[0]+dir1[0] ;
91 //std::cout<<" y:"<<center[1]+dir1[1]<<" z:"<<center[2]+dir1[2]<<std::endl;
93 //else if (Self::Dimension==2)
95 //std::cout<<"Point Edge x:"<<center[0]+dir1[0] ;
96 //std::cout<<" y:"<<center[1]+dir1[1]<<std::endl;
100 // Direction of second profile
101 // pi2n*Iradial + 180°
102 typename TPoint::VectorType dir2;
103 dir2.Fill( double( 0 ) );
104 dir2[ cx ] = std::cos( (pi2n) * double( I_radial ) + double( vnl_math::pi ));
105 dir2[ cy ] = std::sin( (pi2n) * double( I_radial ) + double( vnl_math::pi ));
109 if ( img->TransformPhysicalPointToIndex( center + (dir2*radius), rIdx2 ) )
111 TVector grad_rIdx2 = img->GetPixel( rIdx2 );
113 u_i2.SetVnlVector( ( center - ( center + dir2 ) ).GetVnlVector( ) );
116 Flux2 = grad_rIdx2 * u_i2;
120 //if (Self::Dimension==3)
122 //std::cout<<"Point Edge x:"<<center[0]+dir2[0] ;
123 //std::cout<<" y:"<<center[1]+dir2[1]<<" z:"<<center[2]+dir2[2]<<std::endl;
125 //else if (Self::Dimension==2)
127 //std::cout<<"Point Edge x:"<<center[0]+dir2[0] ;
128 //std::cout<<" y:"<<center[1]+dir2[1]<<std::endl;
132 MFlux += std::min( Flux1, Flux2 );
135 //std::cout<<Self::Dimension<<" radius:"<<radius<<std::endl;
136 //std::cout<<"Center:"<<center[0]<<" "<<center[1]<<std::endl;
137 //std::cout<<"edge:"<<center[0]+radius*std::cos( pi2n * double( 0 ) )<<std::endl;
140 MFlux /= this->m_RadialSampling;
141 FluxFinal.push_back(MFlux);
142 radiusGenerated.push_back(radius);
144 radius += this->m_RadiusStep;
148 dR= *( std::max_element( FluxFinal.begin(), FluxFinal.end() ) );
149 optR= (dR>optR)? dR:optR;
154 return( TScalar(optR) );
157 #endif // __FPA__IMAGE__FUNCTORS__MFLUXMEDIALNESS__HXX__