lib/cpExtensions/Algorithms/GulsunTekMedialness.hxx

   1 // -------------------------------------------------------------------------
   2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
   3 // -------------------------------------------------------------------------
   4
   5 #ifndef __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
   6 #define __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
   7
   8 #include <cmath>
   9 #include <vnl/vnl_math.h>
  10
  11
  12 #include <queue>
  13 #include <map>
  14 #include <set>
  15
  16
  17
  18
  19 // -------------------------------------------------------------------------
  20 template< class G >
  21 cpExtensions::Algorithms::GulsunTekMedialness< G >::
  22 GulsunTekMedialness( )
  23   : Superclass( ),
  24     m_MinRadius( double( 0 ) ),
  25     m_MaxRadius( double( 1 ) ),
  26     m_ProfileSampling( 4 ),
  27     m_RadialSampling( 10 )
  28 {
  29 }
  30
  31 // -------------------------------------------------------------------------
  32 template< class G >
  33 cpExtensions::Algorithms::GulsunTekMedialness< G >::
  34 ~GulsunTekMedialness( )
  35 {
  36 }
  37
  38 // -------------------------------------------------------------------------
  39 template< class G >
  40 typename cpExtensions::Algorithms::GulsunTekMedialness< G >::
  41 TOutput cpExtensions::Algorithms::GulsunTekMedialness< G >::
  42 _Evaluate( const TIndex& i ) const
  43 {
  44   const G* gradient = this->GetInputImage( );
  45   typename G::RegionType region = gradient->GetRequestedRegion( );
  46   typename G::PointType i_P;
  47   gradient->TransformIndexToPhysicalPoint( i, i_P );
  48
  49   std::queue< TIndex > q;
  50   std::set< TIndex, typename TIndex::LexicographicCompare > m;
  51   std::map< TOffset, std::vector< TIndex >, typename TOffset::LexicographicCompare > profiles;
  52   q.push( i );
  53
  54   while( !q.empty( ) )
  55   {
  56     // Get next node
  57     TIndex node = q.front( );
  58     q.pop( );
  59     if( m.find( node ) != m.end( ) )
  60       continue;
  61     m.insert( node );
  62
  63     // Get neighborhood
  64     for( unsigned int d = 0; d < Self::Dimension; d++ )
  65     {
  66       for( int off = -1; off <= 1; off += 2 )
  67       {
  68         TIndex neighbor = node;
  69         neighbor[ d ] += off;
  70         if( region.IsInside( neighbor ) )
  71         {
  72           typename G::PointType neighbor_P;
  73           gradient->TransformIndexToPhysicalPoint( neighbor, neighbor_P );
  74           if( double( i_P.EuclideanDistanceTo( neighbor_P ) ) > this->m_MaxRadius )
  75             continue;
  76
  77           // Normalize offset in terms of a l1 (manhattan) distance
  78           TOffset offset = neighbor - i;
  79           // std::cout << "Offset: " << offset << std::endl;
  80
  81           /*
  82           int max_off = 0;
  83           for( unsigned int o = 0; o < Self::Dimension; ++o )
  84             max_off += std::abs( offset[ o ] );
  85           if( max_off == 0 )
  86             continue;
  87           bool normalized = true;
  88           TOffset normalized_offset;
  89           for( unsigned int o = 0; o < Self::Dimension && normalized; ++o )
  90           {
  91             if( offset[ o ] % max_off == 0 )
  92               normalized_offset[ o ] = offset[ o ] / max_off;
  93             else
  94               normalized = false;
  95
  96           } // rof
  97           if( !normalized )
  98             normalized_offset = offset;
  99           std::cout << "offset: " << normalized_offset << std::endl;
 100
 101           // Update profiles
 102           profiles[ normalized_offset ].push_back( neighbor );
 103           */
 104
 105           // Update queue
 106           q.push( neighbor );
 107
 108         } // fi
 109
 110       } // rof
 111
 112     } // rof
 113
 114   } // elihw
 115   // std::cout << "HOLA: " << profiles.size( ) << std::endl;
 116
 117
 118   return( TOutput( 0 ) );
 119   /* TODO
 120      const G* in = this->GetInputImage( );
 121      double pi2n =
 122      double( 2 ) * double( vnl_math::pi ) /
 123      double( this->m_ProfileSampling );
 124      double rOff = this->m_MaxRadius / double( this->m_RadialSampling - 1 );
 125      double optR = double( 0 );
 126      TPoint pnt;
 127      in->TransformIndexToPhysicalPoint( i, pnt );
 128
 129      // Main loop
 130      for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ )
 131      {
 132      for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ )
 133      {
 134      TProfile maxProfile( this->m_RadialSampling, double( 0 ) );
 135      for( unsigned int p = 0; p < this->m_ProfileSampling; p++ )
 136      {
 137      double a = pi2n * double( p );
 138
 139      // Direction of this profile
 140      TVector dir;
 141      dir.Fill( TScalar( 0 ) );
 142      dir[ cx ] = TScalar( std::cos( a ) );
 143      dir[ cy ] = TScalar( std::sin( a ) );
 144
 145      double maxrise = double( 0 );
 146      double maxfall = double( -1 );
 147      TProfile profile;
 148      for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
 149      {
 150      double radius = double( r ) * rOff;
 151      TIndex idx;
 152      if(
 153      in->TransformPhysicalPointToIndex( pnt + ( dir * radius ), idx )
 154      )
 155      {
 156      TVector g = in->GetPixel( idx );
 157      double b = double( g.GetNorm( ) );
 158      if( double( g * dir ) < double( 0 ) )
 159      b *= double( -1 );
 160      maxrise = ( b > maxrise )? b: maxrise;
 161      if( radius >= this->m_MinRadius )
 162      maxfall = ( b < maxfall )? b: maxfall;
 163      profile.push_back( -b - maxrise );
 164      }
 165      else
 166      profile.push_back( double( 0 ) );
 167
 168      } // rof
 169
 170      for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
 171      {
 172      double E = profile[ r ] / -maxfall;
 173      E = ( E < double( 0 ) )? double( 0 ): E;
 174      E = ( E > double( 1 ) )? double( 1 ): E;
 175      maxProfile[ r ] += E;
 176
 177      } // rof
 178
 179      } // rof
 180
 181      for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
 182      {
 183      double E = maxProfile[ r ] / double( this->m_RadialSampling );
 184      optR = ( E > optR )? E: optR;
 185
 186      } // rof
 187
 188      } // rof
 189
 190      } // rof
 191      return( TScalar( optR ) );
 192   */
 193 }
 194
 195 #endif // __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
 196
 197 // eof - $RCSfile$