lib/cpPlugins/Extensions/Algorithms/InertiaMedialness.h

   1 // -------------------------------------------------------------------------
   2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
   3 // -------------------------------------------------------------------------
   4
   5 #ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__INERTIAMEDIALNESS__H__
   6 #define __CPPLUGINS__EXTENSIONS__ALGORITHMS__INERTIAMEDIALNESS__H__
   7
   8 #include <map>
   9 #include <itkImageFunction.h>
  10
  11 #include <cpPlugins/Extensions/Algorithms/InertiaTensorFunction.h>
  12 #include <itkImageRegionConstIteratorWithIndex.h>
  13
  14 namespace cpPlugins
  15 {
  16   namespace Extensions
  17   {
  18     namespace Algorithms
  19     {
  20       /**
  21        */
  22       template< class I, class S = float >
  23       class InertiaMedialness
  24         : public itk::ImageFunction< I, S, S >
  25       {
  26       public:
  27         // Standard itk types
  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;
  32
  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;
  40
  41         // Sparse buffer
  42         typedef std::map< TIndex, TOutput, typename TIndex::LexicographicCompare > TBuffer;
  43
  44         typedef InertiaTensorFunction< S, I::ImageDimension > TInertia;
  45
  46       public:
  47         itkNewMacro( Self );
  48         itkTypeMacro( InertiaMedialness, itkImageFunction );
  49
  50         itkBooleanMacro( BufferResults );
  51         itkGetConstMacro( BufferResults, bool );
  52         itkGetConstMacro( MaxRadius, double );
  53
  54         itkSetMacro( BufferResults, bool );
  55         itkSetMacro( MaxRadius, double );
  56
  57       public:
  58         virtual void ResetBuffer( )
  59           {
  60             this->m_Buffer.clear( );
  61           }
  62
  63         virtual TOutput Evaluate( const TPoint& p ) const
  64           {
  65             TIndex i;
  66             this->GetInputImage( )->TransformPhysicalPointToIndex( p, i );
  67             return( this->EvaluateAtIndex( i ) );
  68           }
  69
  70         virtual TOutput EvaluateAtIndex( const TIndex& i ) const
  71           {
  72             TOutput res = TOutput( 0 );
  73             bool computed = false;
  74             if( this->m_BufferResults )
  75             {
  76               typename TBuffer::const_iterator bIt = this->m_Buffer.find( i );
  77               computed = ( bIt != this->m_Buffer.end( ) );
  78               res = ( computed )? bIt->second: res;
  79
  80             } // fi
  81
  82             if( !computed )
  83               res = this->_Evaluate( i );
  84
  85             if( this->m_BufferResults )
  86               this->m_Buffer[ i ] = res;
  87             return( res );
  88           }
  89
  90         virtual TOutput EvaluateAtContinuousIndex( const TContIndex& i ) const
  91           {
  92             TPoint p;
  93             this->GetInputImage( )->TransformContinuousIndexToPhysicalPoint( i, p );
  94             return( this->Evaluate( p ) );
  95           }
  96
  97       protected:
  98         InertiaMedialness( )
  99           : Superclass( ),
 100             m_BufferResults( false ),
 101             m_MaxRadius( double( 1 ) )
 102           {
 103             this->m_Buffer.clear( );
 104           }
 105
 106         virtual ~InertiaMedialness( )
 107           {
 108             this->m_Buffer.clear( );
 109           }
 110
 111         virtual TOutput _Evaluate( const TIndex& idx ) const
 112           {
 113             const I* image = this->GetInputImage( );
 114
 115             typename I::PointType p_i;
 116             image->TransformIndexToPhysicalPoint( idx, p_i );
 117
 118             typename I::PointType max_p, min_p;
 119             for( unsigned int d = 0; d < I::ImageDimension; ++d )
 120             {
 121               max_p[ d ] = p_i[ d ] + this->m_MaxRadius;
 122               min_p[ d ] = p_i[ d ] - this->m_MaxRadius;
 123
 124             } // rof
 125             TIndex max_i, min_i;
 126             image->TransformPhysicalPointToIndex( max_p, max_i );
 127             image->TransformPhysicalPointToIndex( min_p, min_i );
 128
 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 )
 134             {
 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 ];
 139
 140               size[ d ] = max_i[ d ] - min_i[ d ];
 141
 142             } // rof
 143
 144             typename I::RegionType region;
 145             region.SetIndex( min_i );
 146             region.SetSize( size );
 147
 148             std::vector< typename TInertia::Pointer > inertias;
 149             itk::ImageRegionConstIteratorWithIndex< I > it( image, region );
 150             for( it.GoToBegin( ); !it.IsAtEnd( ); ++it )
 151             {
 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;
 156
 157               typename TInertia::TPoint i_pnt;
 158               image->TransformIndexToPhysicalPoint( it.GetIndex( ), i_pnt );
 159
 160               for( unsigned long l = 0; l < l1dist; ++l )
 161               {
 162                 if( inertias.size( ) <= l )
 163                   inertias.push_back( TInertia::New( ) );
 164                 inertias[ l ]->AddMass( i_pnt.GetVectorFromOrigin( ), S( it.Get( ) ) );
 165
 166                 /* TODO
 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;
 170                 */
 171
 172               } // rof
 173
 174             } // rof
 175
 176             if( inertias.size( ) > 0 )
 177             {
 178               S res = S( 0 );
 179               for( unsigned int l = 0; l < inertias.size( ); ++l )
 180               {
 181                 typename TInertia::TVector pv, r;
 182                 typename TInertia::TMatrix pm;
 183                 inertias[ l ]->GetEigenAnalysis( pm, pv, r );
 184                 S v = pv.GetNorm( );
 185                 if( l == 0 || v > res )
 186                   res = v;
 187
 188               } // rof
 189               return( res );
 190             }
 191             else
 192               return( TOutput( 0 ) );
 193           }
 194
 195       private:
 196         // Purposely not implemented.
 197         InertiaMedialness( const Self& );
 198         void operator=( const Self& );
 199
 200       protected:
 201         mutable TBuffer m_Buffer;
 202         bool m_BufferResults;
 203
 204         double m_MaxRadius;
 205       };
 206
 207     } // ecapseman
 208
 209   } // ecapseman
 210
 211 } // ecapseman
 212
 213 // TODO: #include <cpPlugins/Extensions/Algorithms/InertiaMedialness.hxx>
 214
 215 #endif // __CPPLUGINS__EXTENSIONS__ALGORITHMS__INERTIAMEDIALNESS__H__
 216
 217 // eof - $RCSfile$