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 // TODO: #include <cpExtensions/Algorithms/InertiaMedialness.hxx>
 210
 211 #endif // __CPEXTENSIONS__ALGORITHMS__INERTIAMEDIALNESS__H__
 212
 213 // eof - $RCSfile$