lib/cpExtensions/Algorithms/InertiaMedialness.h

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