+++ /dev/null
-#ifndef __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
-#define __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
-
-#include <cmath>
-#include <vnl/vnl_math.h>
-#include <itkLineConstIterator.h>
-
-// -------------------------------------------------------------------------
-template< class _TGradient, class _TMask >
-cpExtensions::Algorithms::GulsunTekMedialness< _TGradient, _TMask >::
-GulsunTekMedialness( )
- : Superclass( ),
- m_MinRadius( double( 0 ) ),
- m_MaxRadius( double( 1 ) ),
- m_ProfileSampling( 4 ),
- m_RadialSampling( 10 )
-{
-}
-
-// -------------------------------------------------------------------------
-template< class _TGradient, class _TMask >
-cpExtensions::Algorithms::GulsunTekMedialness< _TGradient, _TMask >::
-~GulsunTekMedialness( )
-{
-}
-
-// -------------------------------------------------------------------------
-template< class _TGradient, class _TMask >
-typename cpExtensions::Algorithms::GulsunTekMedialness< _TGradient, _TMask >::
-TOutput cpExtensions::Algorithms::GulsunTekMedialness< _TGradient, _TMask >::
-_Evaluate( const TIndex& i ) const
-{
- itk::Object::GlobalWarningDisplayOff( );
-
- // Various values
- const _TGradient* in = this->GetInputImage( );
- double pi2n =
- double( 2 ) * double( vnl_math::pi ) /
- double( this->m_ProfileSampling );
- double rOff = this->m_MaxRadius / double( this->m_RadialSampling - 1 );
- double optR = double( 0 );
- TPoint pnt;
- in->TransformIndexToPhysicalPoint( i, pnt );
-
- // Main loop
- for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ )
- {
- for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ )
- {
- TProfile maxProfile( this->m_RadialSampling, double( 0 ) );
- for( unsigned int p = 0; p < this->m_ProfileSampling; p++ )
- {
- double a = pi2n * double( p );
-
- // Direction of this profile
- TVector dir;
- dir.Fill( TScalar( 0 ) );
- dir[ cx ] = TScalar( std::cos( a ) );
- dir[ cy ] = TScalar( std::sin( a ) );
-
- double maxrise = double( 0 );
- double maxfall = double( -1 );
- TProfile profile;
- for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
- {
- double radius = double( r ) * rOff;
- TIndex idx;
- typename TPoint::VectorType aux;
- aux.SetVnlVector( dir.GetVnlVector( ) );
- if(
- in->TransformPhysicalPointToIndex( pnt + ( aux * radius ), idx )
- )
- {
- TVector g = in->GetPixel( idx );
- double b = double( g.GetNorm( ) );
- if( double( g * dir ) < double( 0 ) )
- b *= double( -1 );
- maxrise = ( b > maxrise )? b: maxrise;
- if( radius >= this->m_MinRadius )
- maxfall = ( b < maxfall )? b: maxfall;
- profile.push_back( -b - maxrise );
- }
- else
- profile.push_back( double( 0 ) );
-
- } // rof
-
- for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
- {
- double E = profile[ r ] / -maxfall;
- E = ( E < double( 0 ) )? double( 0 ): E;
- E = ( E > double( 1 ) )? double( 1 ): E;
- maxProfile[ r ] += E;
-
- } // rof
-
- } // rof
-
- for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
- {
- double E = maxProfile[ r ] / double( this->m_RadialSampling );
- optR = ( E > optR )? E: optR;
-
- } // rof
-
- } // rof
-
- } // rof
- return( TScalar( optR ) );
-}
-
-#endif // __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
-
-// eof - $RCSfile$