lib/cpExtensions/Algorithms/FluxMedialness.hxx

   1 #ifndef __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__HXX__
   2 #define __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__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::FluxMedialness< _TGradient, _TMask >::
  11 FluxMedialness( )
  12   : Superclass( ),
  13     m_MinRadius( double( 0 ) ),
  14     m_MaxRadius( double( 1 ) ),
  15     m_RadialSampling( 4 ),
  16     m_RadiusStep( double( 1 ) )
  17 {
  18 }
  19
  20 // -------------------------------------------------------------------------
  21 template< class _TGradient, class _TMask >
  22 cpExtensions::Algorithms::FluxMedialness< _TGradient, _TMask >::
  23 ~FluxMedialness( )
  24 {
  25 }
  26
  27 // -------------------------------------------------------------------------
  28 template< class _TGradient, class _TMask >
  29 typename cpExtensions::Algorithms::FluxMedialness< _TGradient, _TMask >::
  30 TOutput cpExtensions::Algorithms::FluxMedialness< _TGradient, _TMask >::
  31 _Evaluate( const TIndex& i ) const
  32 {
  33   itk::Object::GlobalWarningDisplayOff( );
  34
  35   double pi2n = double( 2 ) * double( vnl_math::pi );
  36   pi2n /= int( this->m_RadialSampling );
  37   const _TGradient* img = this->GetInputImage( );
  38   // Gradient in central pixel
  39   //TVector grad_idx = img->GetPixel( i );
  40
  41   double Flux1 = 0;
  42   double Flux = 0;
  43
  44   TRCandidates FluxFinal;
  45   TRCandidates radiusGenerated;
  46   double dR = double( 0 );
  47   double optR = double( 0 );
  48   TPoint center;
  49   img->TransformIndexToPhysicalPoint( i, center );
  50   double radius = double(0);
  51
  52   for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ )
  53   {
  54     for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ )
  55     {
  56       dR = double( 0 );
  57       FluxFinal.clear();
  58       radiusGenerated.clear();
  59       radius = this->m_MinRadius;
  60       while( radius <= this->m_MaxRadius )
  61       {
  62         Flux = 0;
  63         for( unsigned int I_radial = 0; I_radial < this->m_RadialSampling ; I_radial++ )
  64         {
  65           Flux1 = 0;
  66
  67           // Direction of first profile
  68           typename TPoint::VectorType dir1;
  69           dir1.Fill( double( 0 ) );
  70           dir1[ cx ] = std::cos( pi2n * double( I_radial ) );
  71           dir1[ cy ] = std::sin( pi2n * double( I_radial ) );
  72           dir1 *= (radius);
  73           TIndex rIdx;
  74           if (img->TransformPhysicalPointToIndex( center + dir1, rIdx ))
  75           {
  76             TVector grad_rIdx = img->GetPixel( rIdx );
  77             TVector u_i1;
  78             u_i1.SetVnlVector( ( center - ( center + dir1 ) ).GetVnlVector( ) );
  79             u_i1.Normalize( );
  80             // dot product
  81             Flux1 = grad_rIdx * u_i1;
  82           }
  83
  84           Flux += Flux1;
  85
  86         } // rof
  87         //std::cout<<" radius:"<<radius<<std::endl;
  88         //std::cout<<"Center:"<<center[0]<<" "<<center[1]<<std::endl;
  89         //std::cout<<"edge:"<<center[0]+radius*std::cos( pi2n * double( 0 ) )<<std::endl;
  90         Flux /= this->m_RadialSampling;
  91         FluxFinal.push_back(Flux);
  92         radiusGenerated.push_back(radius);
  93         radius += this->m_RadiusStep;
  94
  95       }     //elihw
  96
  97       dR= *( std::max_element( FluxFinal.begin(), FluxFinal.end() ) );
  98       optR= (dR>optR)? dR:optR;
  99
 100     } // rof
 101
 102   } // rof
 103   return( TScalar(optR) );
 104 }
 105
 106 #endif // __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__HXX__
 107
 108 // eof - $RCSfile$