lib/cpExtensions/Algorithms/GulsunTekMedialness.hxx

   1 #ifndef __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
   2 #define __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
   3
   4 #include <cmath>
   5 #include <vnl/vnl_math.h>
   6 #include <itkLineConstIterator.h>
   7
   8 // -------------------------------------------------------------------------
   9 template< class _TGradient, class _TMask >
  10 cpExtensions::Algorithms::GulsunTekMedialness< _TGradient, _TMask >::
  11 GulsunTekMedialness( )
  12   : Superclass( ),
  13     m_MinRadius( double( 0 ) ),
  14     m_MaxRadius( double( 1 ) ),
  15     m_ProfileSampling( 4 ),
  16     m_RadialSampling( 10 )
  17 {
  18 }
  19
  20 // -------------------------------------------------------------------------
  21 template< class _TGradient, class _TMask >
  22 cpExtensions::Algorithms::GulsunTekMedialness< _TGradient, _TMask >::
  23 ~GulsunTekMedialness( )
  24 {
  25 }
  26
  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
  32 {
  33   itk::Object::GlobalWarningDisplayOff( );
  34
  35   // Various values
  36   const _TGradient* in = this->GetInputImage( );
  37   double pi2n =
  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 );
  42   TPoint pnt;
  43   in->TransformIndexToPhysicalPoint( i, pnt );
  44
  45   // Main loop
  46   for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ )
  47   {
  48     for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ )
  49     {
  50       TProfile maxProfile( this->m_RadialSampling, double( 0 ) );
  51       for( unsigned int p = 0; p < this->m_ProfileSampling; p++ )
  52       {
  53         double a = pi2n * double( p );
  54
  55         // Direction of this profile
  56         TVector dir;
  57         dir.Fill( TScalar( 0 ) );
  58         dir[ cx ] = TScalar( std::cos( a ) );
  59         dir[ cy ] = TScalar( std::sin( a ) );
  60
  61         double maxrise = double( 0 );
  62         double maxfall = double( -1 );
  63         TProfile profile;
  64         for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
  65         {
  66           double radius = double( r ) * rOff;
  67           TIndex idx;
  68           typename TPoint::VectorType aux;
  69           aux.SetVnlVector( dir.GetVnlVector( ) );
  70           if(
  71             in->TransformPhysicalPointToIndex( pnt + ( aux * radius ), idx )
  72             )
  73           {
  74             TVector g = in->GetPixel( idx );
  75             double b = double( g.GetNorm( ) );
  76             if( double( g * dir ) < double( 0 ) )
  77               b *= double( -1 );
  78             maxrise = ( b > maxrise )? b: maxrise;
  79             if( radius >= this->m_MinRadius )
  80               maxfall = ( b < maxfall )? b: maxfall;
  81             profile.push_back( -b - maxrise );
  82           }
  83           else
  84             profile.push_back( double( 0 ) );
  85
  86         } // rof
  87
  88         for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
  89         {
  90           double E = profile[ r ] / -maxfall;
  91           E = ( E < double( 0 ) )? double( 0 ): E;
  92           E = ( E > double( 1 ) )? double( 1 ): E;
  93           maxProfile[ r ] += E;
  94
  95         } // rof
  96
  97       } // rof
  98
  99       for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
 100       {
 101         double E = maxProfile[ r ] / double( this->m_RadialSampling );
 102         optR = ( E > optR )? E: optR;
 103
 104       } // rof
 105
 106     } // rof
 107
 108   } // rof
 109   return( TScalar( optR ) );
 110 }
 111
 112 #endif // __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
 113
 114 // eof - $RCSfile$