1 #ifndef __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
2 #define __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
5 #include <vnl/vnl_math.h>
6 #include <itkLineConstIterator.h>
8 // -------------------------------------------------------------------------
9 template< class _TGradient, class _TMask >
10 cpExtensions::Algorithms::GulsunTekMedialness< _TGradient, _TMask >::
11 GulsunTekMedialness( )
13 m_MinRadius( double( 0 ) ),
14 m_MaxRadius( double( 1 ) ),
15 m_ProfileSampling( 4 ),
16 m_RadialSampling( 10 )
20 // -------------------------------------------------------------------------
21 template< class _TGradient, class _TMask >
22 cpExtensions::Algorithms::GulsunTekMedialness< _TGradient, _TMask >::
23 ~GulsunTekMedialness( )
27 // -------------------------------------------------------------------------
28 template< class _TGradient, class _TMask >
29 typename cpExtensions::Algorithms::GulsunTekMedialness< _TGradient, _TMask >::
30 TOutput cpExtensions::Algorithms::GulsunTekMedialness< _TGradient, _TMask >::
31 _Evaluate( const TIndex& i ) const
33 itk::Object::GlobalWarningDisplayOff( );
36 const _TGradient* in = this->GetInputImage( );
38 double( 2 ) * double( vnl_math::pi ) /
39 double( this->m_ProfileSampling );
40 double rOff = this->m_MaxRadius / double( this->m_RadialSampling - 1 );
41 double optR = double( 0 );
43 in->TransformIndexToPhysicalPoint( i, pnt );
46 for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ )
48 for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ )
50 TProfile maxProfile( this->m_RadialSampling, double( 0 ) );
51 for( unsigned int p = 0; p < this->m_ProfileSampling; p++ )
53 double a = pi2n * double( p );
55 // Direction of this profile
57 dir.Fill( TScalar( 0 ) );
58 dir[ cx ] = TScalar( std::cos( a ) );
59 dir[ cy ] = TScalar( std::sin( a ) );
61 double maxrise = double( 0 );
62 double maxfall = double( -1 );
64 for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
66 double radius = double( r ) * rOff;
68 typename TPoint::VectorType aux;
69 aux.SetVnlVector( dir.GetVnlVector( ) );
71 in->TransformPhysicalPointToIndex( pnt + ( aux * radius ), idx )
74 TVector g = in->GetPixel( idx );
75 double b = double( g.GetNorm( ) );
76 if( double( g * dir ) < double( 0 ) )
78 maxrise = ( b > maxrise )? b: maxrise;
79 if( radius >= this->m_MinRadius )
80 maxfall = ( b < maxfall )? b: maxfall;
81 profile.push_back( -b - maxrise );
84 profile.push_back( double( 0 ) );
88 for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
90 double E = profile[ r ] / -maxfall;
91 E = ( E < double( 0 ) )? double( 0 ): E;
92 E = ( E > double( 1 ) )? double( 1 ): E;
99 for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
101 double E = maxProfile[ r ] / double( this->m_RadialSampling );
102 optR = ( E > optR )? E: optR;
109 return( TScalar( optR ) );
112 #endif // __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__