1 // -------------------------------------------------------------------------
2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
3 // -------------------------------------------------------------------------
5 #ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__INERTIAMEDIALNESS__H__
6 #define __CPPLUGINS__EXTENSIONS__ALGORITHMS__INERTIAMEDIALNESS__H__
9 #include <itkImageFunction.h>
11 #include <cpPlugins/Extensions/Algorithms/InertiaTensorFunction.h>
12 #include <itkImageRegionConstIteratorWithIndex.h>
22 template< class I, class S = float >
23 class InertiaMedialness
24 : public itk::ImageFunction< I, S, S >
28 typedef InertiaMedialness Self;
29 typedef itk::ImageFunction< I, S, S > Superclass;
30 typedef itk::SmartPointer< Self > Pointer;
31 typedef itk::SmartPointer< const Self > ConstPointer;
33 // Types from base itk::ImageFunction
34 typedef typename Superclass::InputType TInput;
35 typedef typename Superclass::OutputType TOutput;
36 typedef typename Superclass::PointType TPoint;
37 typedef typename Superclass::ContinuousIndexType TContIndex;
38 typedef typename Superclass::IndexType TIndex;
39 typedef typename TIndex::OffsetType TOffset;
42 typedef std::map< TIndex, TOutput, typename TIndex::LexicographicCompare > TBuffer;
44 typedef InertiaTensorFunction< S, I::ImageDimension > TInertia;
48 itkTypeMacro( InertiaMedialness, itkImageFunction );
50 itkBooleanMacro( BufferResults );
51 itkGetConstMacro( BufferResults, bool );
52 itkGetConstMacro( MaxRadius, double );
54 itkSetMacro( BufferResults, bool );
55 itkSetMacro( MaxRadius, double );
58 virtual void ResetBuffer( )
60 this->m_Buffer.clear( );
63 virtual TOutput Evaluate( const TPoint& p ) const
66 this->GetInputImage( )->TransformPhysicalPointToIndex( p, i );
67 return( this->EvaluateAtIndex( i ) );
70 virtual TOutput EvaluateAtIndex( const TIndex& i ) const
72 TOutput res = TOutput( 0 );
73 bool computed = false;
74 if( this->m_BufferResults )
76 typename TBuffer::const_iterator bIt = this->m_Buffer.find( i );
77 computed = ( bIt != this->m_Buffer.end( ) );
78 res = ( computed )? bIt->second: res;
83 res = this->_Evaluate( i );
85 if( this->m_BufferResults )
86 this->m_Buffer[ i ] = res;
90 virtual TOutput EvaluateAtContinuousIndex( const TContIndex& i ) const
93 this->GetInputImage( )->TransformContinuousIndexToPhysicalPoint( i, p );
94 return( this->Evaluate( p ) );
100 m_BufferResults( false ),
101 m_MaxRadius( double( 1 ) )
103 this->m_Buffer.clear( );
106 virtual ~InertiaMedialness( )
108 this->m_Buffer.clear( );
111 virtual TOutput _Evaluate( const TIndex& idx ) const
113 const I* image = this->GetInputImage( );
115 typename I::PointType p_i;
116 image->TransformIndexToPhysicalPoint( idx, p_i );
118 typename I::PointType max_p, min_p;
119 for( unsigned int d = 0; d < I::ImageDimension; ++d )
121 max_p[ d ] = p_i[ d ] + this->m_MaxRadius;
122 min_p[ d ] = p_i[ d ] - this->m_MaxRadius;
126 image->TransformPhysicalPointToIndex( max_p, max_i );
127 image->TransformPhysicalPointToIndex( min_p, min_i );
129 typename I::RegionType in_region = image->GetRequestedRegion( );
130 TIndex in_index = in_region.GetIndex( );
131 TIndex in_last = in_index + in_region.GetSize( );
132 typename I::SizeType size;
133 for( unsigned int d = 0; d < I::ImageDimension; ++d )
135 if( min_i[ d ] < in_index[ d ] ) min_i[ d ] = in_index[ d ];
136 if( max_i[ d ] < in_index[ d ] ) max_i[ d ] = in_index[ d ];
137 if( min_i[ d ] >= in_last[ d ] ) min_i[ d ] = in_last[ d ];
138 if( max_i[ d ] >= in_last[ d ] ) max_i[ d ] = in_last[ d ];
140 size[ d ] = max_i[ d ] - min_i[ d ];
144 typename I::RegionType region;
145 region.SetIndex( min_i );
146 region.SetSize( size );
148 std::vector< typename TInertia::Pointer > inertias;
149 itk::ImageRegionConstIteratorWithIndex< I > it( image, region );
150 for( it.GoToBegin( ); !it.IsAtEnd( ); ++it )
152 TOffset off = it.GetIndex( ) - idx;
153 unsigned long l1dist = std::abs( off[ 0 ] );
154 for( unsigned int d = 1; d < I::ImageDimension; ++d )
155 l1dist = ( std::abs( off[ d ] ) > l1dist )? std::abs( off[ d ] ): l1dist;
157 typename TInertia::TPoint i_pnt;
158 image->TransformIndexToPhysicalPoint( it.GetIndex( ), i_pnt );
160 for( unsigned long l = 0; l < l1dist; ++l )
162 if( inertias.size( ) <= l )
163 inertias.push_back( TInertia::New( ) );
164 inertias[ l ]->AddMass( i_pnt.GetVectorFromOrigin( ), S( it.Get( ) ) );
167 typename TInertias::iterator inIt = inertias.find( l );
168 if( inIt == inertias.end( ) )
169 inIt = inertias.insert( std::pair< unsigned long, typename TInertia::Pointer >( l, TInertia::New( ) ) ).first;
176 if( inertias.size( ) > 0 )
179 for( unsigned int l = 0; l < inertias.size( ); ++l )
181 typename TInertia::TVector pv, r;
182 typename TInertia::TMatrix pm;
183 inertias[ l ]->GetEigenAnalysis( pm, pv, r );
185 if( l == 0 || v > res )
192 return( TOutput( 0 ) );
196 // Purposely not implemented.
197 InertiaMedialness( const Self& );
198 void operator=( const Self& );
201 mutable TBuffer m_Buffer;
202 bool m_BufferResults;
213 // TODO: #include <cpPlugins/Extensions/Algorithms/InertiaMedialness.hxx>
215 #endif // __CPPLUGINS__EXTENSIONS__ALGORITHMS__INERTIAMEDIALNESS__H__