yet another refactoring

[cpPlugins.git] / lib / cpExtensions / Algorithms / GulsunTekMedialness.hxx

diff --git a/lib/cpExtensions/Algorithms/GulsunTekMedialness.hxx b/lib/cpExtensions/Algorithms/GulsunTekMedialness.hxx
new file mode 100644 (file)
index 0000000..fc59929
--- /dev/null
+++ b/lib/cpExtensions/Algorithms/GulsunTekMedialness.hxx
@@ -0,0 +1,114 @@
+#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$