+++ /dev/null
-#ifndef __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__HXX__
-#define __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__HXX__
-
-#include <cmath>
-#include <vnl/vnl_math.h>
-#include <itkLineConstIterator.h>
-
-// -------------------------------------------------------------------------
-template< class _TGradient, class _TMask >
-cpExtensions::Algorithms::MFluxMedialness< _TGradient, _TMask >::
-MFluxMedialness( )
- : Superclass( ),
- m_MinRadius( double( 0 ) ),
- m_MaxRadius( double( 1 ) ),
- m_RadialSampling( 4 ),
- m_RadiusStep( double( 1 ) )
-{
-}
-
-// -------------------------------------------------------------------------
-template< class _TGradient, class _TMask >
-cpExtensions::Algorithms::MFluxMedialness< _TGradient, _TMask >::
-~MFluxMedialness( )
-{
-}
-
-// -------------------------------------------------------------------------
-template< class _TGradient, class _TMask >
-typename cpExtensions::Algorithms::MFluxMedialness< _TGradient, _TMask >::
-TOutput cpExtensions::Algorithms::MFluxMedialness< _TGradient, _TMask >::
-_Evaluate( const TIndex& i ) const
-{
- itk::Object::GlobalWarningDisplayOff( );
-
- double pi2n = double( 2 ) * double( vnl_math::pi );
- pi2n /= int( this->m_RadialSampling );
- const _TGradient* img = this->GetInputImage( );
- //const itk::Image::SpacingType& input_spacing = img->GetSpacing( );
-
- double Flux1 = 0;
- double Flux2 = 0;
- double MFlux = 0;
-
- TRCandidates FluxFinal;
- TRCandidates radiusGenerated;
- double dR = double( 0 );
- double optR = double( 0 );
- TPoint center;
- img->TransformIndexToPhysicalPoint( i, center );
- double radius;
-
- for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ )
- {
- for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ )
- {
- dR = double( 0 );
- FluxFinal.clear();
- radiusGenerated.clear();
- radius = this->m_MinRadius;
- while( radius <= this->m_MaxRadius )
- {
- MFlux = 0;
- for( unsigned int I_radial = 0; I_radial < this->m_RadialSampling / 2; I_radial++ )
- {
- Flux1 = 0;
- Flux2 = 0;
-
- // Direction of first profile
- typename TPoint::VectorType dir1;
- dir1.Fill( double( 0 ) );
- dir1[ cx ] = std::cos( pi2n * double( I_radial ) );
- dir1[ cy ] = std::sin( pi2n * double( I_radial ) );
- //dir1 *= (radius);
-
- TIndex rIdx;
-
- if ( img->TransformPhysicalPointToIndex( center + (dir1*radius), rIdx ) )
- {
- TVector grad_rIdx = img->GetPixel( rIdx );
- TVector u_i1;
- u_i1.SetVnlVector( ( center - ( center + dir1 ) ).GetVnlVector( ) );
- u_i1.Normalize( );
- // dot product
- Flux1 = grad_rIdx * u_i1;
- }
- else
- {
- //if (Self::Dimension==3)
- //{
- //std::cout<<"Point Edge x:"<<center[0]+dir1[0] ;
- //std::cout<<" y:"<<center[1]+dir1[1]<<" z:"<<center[2]+dir1[2]<<std::endl;
- //}
- //else if (Self::Dimension==2)
- //{
- //std::cout<<"Point Edge x:"<<center[0]+dir1[0] ;
- //std::cout<<" y:"<<center[1]+dir1[1]<<std::endl;
- //}
- }
-
- // Direction of second profile
- // pi2n*Iradial + 180°
- typename TPoint::VectorType dir2;
- dir2.Fill( double( 0 ) );
- dir2[ cx ] = std::cos( (pi2n) * double( I_radial ) + double( vnl_math::pi ));
- dir2[ cy ] = std::sin( (pi2n) * double( I_radial ) + double( vnl_math::pi ));
-
- TIndex rIdx2;
-
- if ( img->TransformPhysicalPointToIndex( center + (dir2*radius), rIdx2 ) )
- {
- TVector grad_rIdx2 = img->GetPixel( rIdx2 );
- TVector u_i2;
- u_i2.SetVnlVector( ( center - ( center + dir2 ) ).GetVnlVector( ) );
- u_i2.Normalize( );
-
- Flux2 = grad_rIdx2 * u_i2;
- }
- else
- {
- //if (Self::Dimension==3)
- //{
- //std::cout<<"Point Edge x:"<<center[0]+dir2[0] ;
- //std::cout<<" y:"<<center[1]+dir2[1]<<" z:"<<center[2]+dir2[2]<<std::endl;
- //}
- //else if (Self::Dimension==2)
- //{
- //std::cout<<"Point Edge x:"<<center[0]+dir2[0] ;
- //std::cout<<" y:"<<center[1]+dir2[1]<<std::endl;
- //}
- }
-
- MFlux += std::min( Flux1, Flux2 );
- } // rof
-
- //std::cout<<Self::Dimension<<" radius:"<<radius<<std::endl;
- //std::cout<<"Center:"<<center[0]<<" "<<center[1]<<std::endl;
- //std::cout<<"edge:"<<center[0]+radius*std::cos( pi2n * double( 0 ) )<<std::endl;
-
- MFlux *= 2;
- MFlux /= this->m_RadialSampling;
- FluxFinal.push_back(MFlux);
- radiusGenerated.push_back(radius);
-
- radius += this->m_RadiusStep;
-
- } //elihw
-
- dR= *( std::max_element( FluxFinal.begin(), FluxFinal.end() ) );
- optR= (dR>optR)? dR:optR;
-
- } // rof
-
- } // rof
- return( TScalar(optR) );
-}
-
-#endif // __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__HXX__
-
-// eof - $RCSfile$