1 // -------------------------------------------------------------------------
2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
3 // -------------------------------------------------------------------------
5 #ifndef __CPEXTENSIONS__ALGORITHMS__INERTIAMEDIALNESS__H__
6 #define __CPEXTENSIONS__ALGORITHMS__INERTIAMEDIALNESS__H__
9 #include <itkImageFunction.h>
11 #include <cpExtensions/Algorithms/InertiaTensorFunction.h>
12 #include <itkImageRegionConstIteratorWithIndex.h>
14 namespace cpExtensions
20 template< class I, class S = float >
21 class InertiaMedialness
22 : public itk::ImageFunction< I, S, S >
26 typedef InertiaMedialness Self;
27 typedef itk::ImageFunction< I, S, S > Superclass;
28 typedef itk::SmartPointer< Self > Pointer;
29 typedef itk::SmartPointer< const Self > ConstPointer;
31 // Types from base itk::ImageFunction
32 typedef typename Superclass::InputType TInput;
33 typedef typename Superclass::OutputType TOutput;
34 typedef typename Superclass::PointType TPoint;
35 typedef typename Superclass::ContinuousIndexType TContIndex;
36 typedef typename Superclass::IndexType TIndex;
37 typedef typename TIndex::OffsetType TOffset;
40 typedef std::map< TIndex, TOutput, typename TIndex::LexicographicCompare > TBuffer;
42 typedef InertiaTensorFunction< S, I::ImageDimension > TInertia;
46 itkTypeMacro( InertiaMedialness, itkImageFunction );
48 itkBooleanMacro( BufferResults );
49 itkGetConstMacro( BufferResults, bool );
50 itkGetConstMacro( MaxRadius, double );
52 itkSetMacro( BufferResults, bool );
53 itkSetMacro( MaxRadius, double );
56 virtual void ResetBuffer( )
58 this->m_Buffer.clear( );
61 virtual TOutput Evaluate( const TPoint& p ) const
64 this->GetInputImage( )->TransformPhysicalPointToIndex( p, i );
65 return( this->EvaluateAtIndex( i ) );
68 virtual TOutput EvaluateAtIndex( const TIndex& i ) const
70 TOutput res = TOutput( 0 );
71 bool computed = false;
72 if( this->m_BufferResults )
74 typename TBuffer::const_iterator bIt = this->m_Buffer.find( i );
75 computed = ( bIt != this->m_Buffer.end( ) );
76 res = ( computed )? bIt->second: res;
81 res = this->_Evaluate( i );
83 if( this->m_BufferResults )
84 this->m_Buffer[ i ] = res;
88 virtual TOutput EvaluateAtContinuousIndex( const TContIndex& i ) const
91 this->GetInputImage( )->TransformContinuousIndexToPhysicalPoint( i, p );
92 return( this->Evaluate( p ) );
98 m_BufferResults( false ),
99 m_MaxRadius( double( 1 ) )
101 this->m_Buffer.clear( );
104 virtual ~InertiaMedialness( )
106 this->m_Buffer.clear( );
109 virtual TOutput _Evaluate( const TIndex& idx ) const
111 const I* image = this->GetInputImage( );
113 typename I::PointType p_i;
114 image->TransformIndexToPhysicalPoint( idx, p_i );
116 typename I::PointType max_p, min_p;
117 for( unsigned int d = 0; d < I::ImageDimension; ++d )
119 max_p[ d ] = p_i[ d ] + this->m_MaxRadius;
120 min_p[ d ] = p_i[ d ] - this->m_MaxRadius;
124 image->TransformPhysicalPointToIndex( max_p, max_i );
125 image->TransformPhysicalPointToIndex( min_p, min_i );
127 typename I::RegionType in_region = image->GetRequestedRegion( );
128 TIndex in_index = in_region.GetIndex( );
129 TIndex in_last = in_index + in_region.GetSize( );
130 typename I::SizeType size;
131 for( unsigned int d = 0; d < I::ImageDimension; ++d )
133 if( min_i[ d ] < in_index[ d ] ) min_i[ d ] = in_index[ d ];
134 if( max_i[ d ] < in_index[ d ] ) max_i[ d ] = in_index[ d ];
135 if( min_i[ d ] >= in_last[ d ] ) min_i[ d ] = in_last[ d ];
136 if( max_i[ d ] >= in_last[ d ] ) max_i[ d ] = in_last[ d ];
138 size[ d ] = max_i[ d ] - min_i[ d ];
142 typename I::RegionType region;
143 region.SetIndex( min_i );
144 region.SetSize( size );
146 std::vector< typename TInertia::Pointer > inertias;
147 itk::ImageRegionConstIteratorWithIndex< I > it( image, region );
148 for( it.GoToBegin( ); !it.IsAtEnd( ); ++it )
150 TOffset off = it.GetIndex( ) - idx;
151 unsigned long l1dist = std::abs( off[ 0 ] );
152 for( unsigned int d = 1; d < I::ImageDimension; ++d )
153 l1dist = ( std::abs( off[ d ] ) > l1dist )? std::abs( off[ d ] ): l1dist;
155 typename TInertia::TPoint i_pnt;
156 image->TransformIndexToPhysicalPoint( it.GetIndex( ), i_pnt );
158 for( unsigned long l = 0; l < l1dist; ++l )
160 if( inertias.size( ) <= l )
161 inertias.push_back( TInertia::New( ) );
162 inertias[ l ]->AddMass( i_pnt.GetVectorFromOrigin( ), S( it.Get( ) ) );
165 typename TInertias::iterator inIt = inertias.find( l );
166 if( inIt == inertias.end( ) )
167 inIt = inertias.insert( std::pair< unsigned long, typename TInertia::Pointer >( l, TInertia::New( ) ) ).first;
174 if( inertias.size( ) > 0 )
177 for( unsigned int l = 0; l < inertias.size( ); ++l )
179 typename TInertia::TVector pv, r;
180 typename TInertia::TMatrix pm;
181 inertias[ l ]->GetEigenAnalysis( pm, pv, r );
183 if( l == 0 || v > res )
190 return( TOutput( 0 ) );
194 // Purposely not implemented.
195 InertiaMedialness( const Self& );
196 void operator=( const Self& );
199 mutable TBuffer m_Buffer;
200 bool m_BufferResults;
209 #ifndef ITK_MANUAL_INSTANTIATION
210 // TODO: #include <cpExtensions/Algorithms/InertiaMedialness.hxx>
211 #endif // ITK_MANUAL_INSTANTIATION
213 #endif // __CPEXTENSIONS__ALGORITHMS__INERTIAMEDIALNESS__H__