1 // -------------------------------------------------------------------------
2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
3 // -------------------------------------------------------------------------
5 #ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
6 #define __CPPLUGINS__EXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
9 #include <vnl/vnl_math.h>
19 // -------------------------------------------------------------------------
21 cpPlugins::Extensions::Algorithms::GulsunTekMedialness< G >::
22 GulsunTekMedialness( )
24 m_MinRadius( double( 0 ) ),
25 m_MaxRadius( double( 1 ) ),
26 m_ProfileSampling( 4 ),
27 m_RadialSampling( 10 )
31 // -------------------------------------------------------------------------
33 cpPlugins::Extensions::Algorithms::GulsunTekMedialness< G >::
34 ~GulsunTekMedialness( )
38 // -------------------------------------------------------------------------
40 typename cpPlugins::Extensions::Algorithms::GulsunTekMedialness< G >::
41 TOutput cpPlugins::Extensions::Algorithms::GulsunTekMedialness< G >::
42 _Evaluate( const TIndex& i ) const
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 );
49 std::queue< TIndex > q;
50 std::set< TIndex, typename TIndex::LexicographicCompare > m;
51 std::map< TOffset, std::vector< TIndex >, typename TOffset::LexicographicCompare > profiles;
57 TIndex node = q.front( );
59 if( m.find( node ) != m.end( ) )
64 for( unsigned int d = 0; d < Self::Dimension; d++ )
66 for( int off = -1; off <= 1; off += 2 )
68 TIndex neighbor = node;
70 if( region.IsInside( neighbor ) )
72 typename G::PointType neighbor_P;
73 gradient->TransformIndexToPhysicalPoint( neighbor, neighbor_P );
74 if( double( i_P.EuclideanDistanceTo( neighbor_P ) ) > this->m_MaxRadius )
77 // Normalize offset in terms of a l1 (manhattan) distance
78 TOffset offset = neighbor - i;
79 // std::cout << "Offset: " << offset << std::endl;
83 for( unsigned int o = 0; o < Self::Dimension; ++o )
84 max_off += std::abs( offset[ o ] );
87 bool normalized = true;
88 TOffset normalized_offset;
89 for( unsigned int o = 0; o < Self::Dimension && normalized; ++o )
91 if( offset[ o ] % max_off == 0 )
92 normalized_offset[ o ] = offset[ o ] / max_off;
98 normalized_offset = offset;
99 std::cout << "offset: " << normalized_offset << std::endl;
102 profiles[ normalized_offset ].push_back( neighbor );
115 // std::cout << "HOLA: " << profiles.size( ) << std::endl;
118 return( TOutput( 0 ) );
120 const G* in = this->GetInputImage( );
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 );
127 in->TransformIndexToPhysicalPoint( i, pnt );
130 for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ )
132 for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ )
134 TProfile maxProfile( this->m_RadialSampling, double( 0 ) );
135 for( unsigned int p = 0; p < this->m_ProfileSampling; p++ )
137 double a = pi2n * double( p );
139 // Direction of this profile
141 dir.Fill( TScalar( 0 ) );
142 dir[ cx ] = TScalar( std::cos( a ) );
143 dir[ cy ] = TScalar( std::sin( a ) );
145 double maxrise = double( 0 );
146 double maxfall = double( -1 );
148 for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
150 double radius = double( r ) * rOff;
153 in->TransformPhysicalPointToIndex( pnt + ( dir * radius ), idx )
156 TVector g = in->GetPixel( idx );
157 double b = double( g.GetNorm( ) );
158 if( double( g * dir ) < double( 0 ) )
160 maxrise = ( b > maxrise )? b: maxrise;
161 if( radius >= this->m_MinRadius )
162 maxfall = ( b < maxfall )? b: maxfall;
163 profile.push_back( -b - maxrise );
166 profile.push_back( double( 0 ) );
170 for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
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;
181 for( unsigned int r = 0; r < this->m_RadialSampling; r++ )
183 double E = maxProfile[ r ] / double( this->m_RadialSampling );
184 optR = ( E > optR )? E: optR;
191 return( TScalar( optR ) );
195 #endif // __CPPLUGINS__EXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__