--- /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$