lib/cpExtensions/Algorithms/MFluxMedialness.hxx

   1 #ifndef __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__HXX__
   2 #define __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__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::MFluxMedialness< _TGradient, _TMask >::
  11 MFluxMedialness( )
  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::MFluxMedialness< _TGradient, _TMask >::
  23 ~MFluxMedialness( )
  24 {
  25 }
  26
  27 // -------------------------------------------------------------------------
  28 template< class _TGradient, class _TMask >
  29 typename cpExtensions::Algorithms::MFluxMedialness< _TGradient, _TMask >::
  30 TOutput cpExtensions::Algorithms::MFluxMedialness< _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   //const itk::Image::SpacingType& input_spacing = img->GetSpacing( );
  39
  40   double Flux1 = 0;
  41   double Flux2 = 0;
  42   double MFlux = 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;
  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         MFlux = 0;
  63         for( unsigned int I_radial = 0; I_radial < this->m_RadialSampling / 2; I_radial++ )
  64         {
  65           Flux1 = 0;
  66           Flux2 = 0;
  67
  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 ) );
  73           //dir1 *= (radius);
  74
  75           TIndex rIdx;
  76
  77           if ( img->TransformPhysicalPointToIndex( center + (dir1*radius), rIdx ) )
  78           {
  79             TVector grad_rIdx = img->GetPixel( rIdx );
  80             TVector u_i1;
  81             u_i1.SetVnlVector( ( center - ( center + dir1 ) ).GetVnlVector( ) );
  82             u_i1.Normalize( );
  83             // dot product
  84             Flux1 = grad_rIdx * u_i1;
  85           }
  86           else
  87           {
  88             //if (Self::Dimension==3)
  89             //{
  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;
  92             //}
  93             //else if (Self::Dimension==2)
  94             //{
  95             //std::cout<<"Point Edge x:"<<center[0]+dir1[0] ;
  96             //std::cout<<" y:"<<center[1]+dir1[1]<<std::endl;
  97             //}
  98           }
  99
 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 ));
 106
 107           TIndex rIdx2;
 108
 109           if ( img->TransformPhysicalPointToIndex( center + (dir2*radius), rIdx2 ) )
 110           {
 111             TVector grad_rIdx2 = img->GetPixel( rIdx2 );
 112             TVector u_i2;
 113             u_i2.SetVnlVector( ( center - ( center + dir2 ) ).GetVnlVector( ) );
 114             u_i2.Normalize( );
 115
 116             Flux2 = grad_rIdx2 * u_i2;
 117           }
 118           else
 119           {
 120             //if (Self::Dimension==3)
 121             //{
 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;
 124             //}
 125             //else if (Self::Dimension==2)
 126             //{
 127             //std::cout<<"Point Edge x:"<<center[0]+dir2[0] ;
 128             //std::cout<<" y:"<<center[1]+dir2[1]<<std::endl;
 129             //}
 130           }
 131
 132           MFlux += std::min( Flux1, Flux2 );
 133         } // rof
 134
 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;
 138
 139         MFlux *= 2;
 140         MFlux /= this->m_RadialSampling;
 141         FluxFinal.push_back(MFlux);
 142         radiusGenerated.push_back(radius);
 143
 144         radius += this->m_RadiusStep;
 145
 146       }     //elihw
 147
 148       dR= *( std::max_element( FluxFinal.begin(), FluxFinal.end() ) );
 149       optR= (dR>optR)? dR:optR;
 150
 151     } // rof
 152
 153   } // rof
 154   return( TScalar(optR) );
 155 }
 156
 157 #endif // __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__HXX__
 158
 159 // eof - $RCSfile$