lib/cpExtensions/Algorithms/MultiScaleGaussianImageFilter.hxx

   1 // -------------------------------------------------------------------------
   2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
   3 // -------------------------------------------------------------------------
   4
   5 #ifndef __cpExtensions__Algorithms__MultiScaleGaussianImageFilter__hxx__
   6 #define __cpExtensions__Algorithms__MultiScaleGaussianImageFilter__hxx__
   7
   8 #include <vnl/vnl_vector.h>
   9
  10 #include <itkImageRegionIterator.h>
  11 #include <itkImageRegionConstIterator.h>
  12 #include <itkNumericTraits.h>
  13 #include <itkProgressAccumulator.h>
  14
  15 #include <itkBinaryFunctorImageFilter.h>
  16 #include <itkUnaryFunctorImageFilter.h>
  17 #include <itkGradientRecursiveGaussianImageFilter.h>
  18
  19 // -------------------------------------------------------------------------
  20 template< class I, class O >
  21 cpExtensions::Algorithms::
  22 MultiScaleGaussianImageFilter< I, O >::_Greater::
  23 _Greater( )
  24 {
  25 }
  26
  27 // -------------------------------------------------------------------------
  28 template< class I, class O >
  29 cpExtensions::Algorithms::
  30 MultiScaleGaussianImageFilter< I, O >::_Greater::
  31 ~_Greater( )
  32 {
  33 }
  34
  35 // -------------------------------------------------------------------------
  36 template< class I, class O >
  37 bool cpExtensions::Algorithms::
  38 MultiScaleGaussianImageFilter< I, O >::_Greater::
  39 operator!=( const _Greater& b ) const
  40 {
  41   return( false );
  42 }
  43
  44 // -------------------------------------------------------------------------
  45 template< class I, class O >
  46 bool cpExtensions::Algorithms::
  47 MultiScaleGaussianImageFilter< I, O >::_Greater::
  48 operator==( const _Greater& b ) const
  49 {
  50   return !( *this != b );
  51 }
  52
  53 // -------------------------------------------------------------------------
  54 template< class I, class O >
  55 typename cpExtensions::Algorithms::
  56 MultiScaleGaussianImageFilter< I, O >::_Greater::
  57 _T cpExtensions::Algorithms::
  58 MultiScaleGaussianImageFilter< I, O >::_Greater::
  59 operator()( const _T& a ) const
  60 {
  61   return( a );
  62 }
  63
  64 // -------------------------------------------------------------------------
  65 template< class I, class O >
  66 typename cpExtensions::Algorithms::
  67 MultiScaleGaussianImageFilter< I, O >::_Greater::
  68 _T cpExtensions::Algorithms::
  69 MultiScaleGaussianImageFilter< I, O >::_Greater::
  70 operator()( const _T& a, const _T& b ) const
  71 {
  72   typedef itk::NumericTraits< _T >     _TTraits;
  73   typedef typename _TTraits::ValueType _TValue;
  74   typedef vnl_vector< _TValue >        _TVector;
  75
  76   _TVector va( _TTraits::GetLength( ) );
  77   _TVector vb( _TTraits::GetLength( ) );
  78
  79   _TTraits::AssignToArray( a, va );
  80   _TTraits::AssignToArray( b, vb );
  81   return( ( vb.magnitude( ) < va.magnitude( ) )? a: b );
  82 }
  83
  84 // -------------------------------------------------------------------------
  85 template< class I, class O >
  86 void
  87 cpExtensions::Algorithms::
  88 MultiScaleGaussianImageFilter< I, O >::
  89 AddScale( const double& s )
  90 {
  91   if( this->m_Scales.insert( s ).second )
  92     this->Modified( );
  93 }
  94
  95 // -------------------------------------------------------------------------
  96 template< class I, class O >
  97 unsigned long
  98 cpExtensions::Algorithms::
  99 MultiScaleGaussianImageFilter< I, O >::
 100 GetNumberOfScales( ) const
 101 {
 102   return( this->m_Scales.size( ) );
 103 }
 104
 105 // -------------------------------------------------------------------------
 106 template< class I, class O >
 107 cpExtensions::Algorithms::
 108 MultiScaleGaussianImageFilter< I, O >::
 109 MultiScaleGaussianImageFilter( )
 110   : Superclass( )
 111 {
 112 }
 113
 114 // -------------------------------------------------------------------------
 115 template< class I, class O >
 116 cpExtensions::Algorithms::
 117 MultiScaleGaussianImageFilter< I, O >::
 118 ~MultiScaleGaussianImageFilter( )
 119 {
 120 }
 121
 122 // -------------------------------------------------------------------------
 123 template< class I, class O >
 124 void
 125 cpExtensions::Algorithms::
 126 MultiScaleGaussianImageFilter< I, O >::
 127 GenerateData( )
 128 {
 129   typedef itk::GradientRecursiveGaussianImageFilter< I, O > _TGF;
 130   this->_GenerateData< _TGF >( );
 131 }
 132
 133 // -------------------------------------------------------------------------
 134 template< class I, class O >
 135 template< class F >
 136 void
 137 cpExtensions::Algorithms::
 138 MultiScaleGaussianImageFilter< I, O >::
 139 _GenerateData( )
 140 {
 141   // Some types
 142   typedef itk::BinaryFunctorImageFilter< O, O, O, _Greater > _TMaxFilter;
 143   typedef itk::UnaryFunctorImageFilter< O, O, _Greater >     _TCopyFilter;
 144   typedef itk::ImageRegionConstIterator< O >                 _TConstIt;
 145   typedef itk::ImageRegionIterator< O >                      _TIt;
 146   typedef itk::ProgressAccumulator                           _TProgress;
 147
 148   // Some values
 149   typename I::ConstPointer input = this->GetInput( );
 150   typename O::Pointer output = this->GetOutput( );
 151   unsigned int nThreads = this->GetNumberOfThreads( );
 152   float fw = float( 1 ) / float( this->m_Scales.size( ) << 1 );
 153
 154   // Progress accumulator
 155   _TProgress::Pointer pg = _TProgress::New( );
 156   pg->SetMiniPipelineFilter( this );
 157
 158   // Copy image information
 159   output->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
 160   output->SetRequestedRegion( input->GetRequestedRegion( ) );
 161   output->SetBufferedRegion( input->GetBufferedRegion( ) );
 162   output->SetSpacing( input->GetSpacing( ) );
 163   output->SetOrigin( input->GetOrigin( ) );
 164   output->SetDirection( input->GetDirection( ) );
 165   output->Allocate( );
 166
 167   // Intermediary buffer
 168   typename O::Pointer buffer = O::New( );
 169   buffer->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
 170   buffer->SetRequestedRegion( input->GetRequestedRegion( ) );
 171   buffer->SetBufferedRegion( input->GetBufferedRegion( ) );
 172   buffer->SetSpacing( input->GetSpacing( ) );
 173   buffer->SetOrigin( input->GetOrigin( ) );
 174   buffer->SetDirection( input->GetDirection( ) );
 175   buffer->Allocate( );
 176
 177   // Perform all scales
 178   _TIt gIt( output, output->GetRequestedRegion( ) );
 179   TScalesContainer::const_iterator sIt = this->m_Scales.begin( );
 180   for( ; sIt != this->m_Scales.end( ); sIt++ )
 181   {
 182     // Single scale filter
 183     typename F::Pointer filter = F::New( );
 184     filter->SetInput( input );
 185     filter->SetNormalizeAcrossScale( true );
 186     filter->SetNumberOfThreads( nThreads );
 187     filter->SetSigma( *sIt );
 188     pg->RegisterInternalFilter( filter, fw );
 189     filter->GraftOutput( buffer );
 190     filter->Update( );
 191     buffer->Graft( filter->GetOutput( ) );
 192
 193     // Get maximum response
 194     if( sIt == this->m_Scales.begin( ) )
 195     {
 196       // Copy first result
 197       typename _TCopyFilter::Pointer copy = _TCopyFilter::New( );
 198       copy->SetInput( buffer );
 199       copy->InPlaceOff( );
 200       copy->SetNumberOfThreads( nThreads );
 201       pg->RegisterInternalFilter( copy, fw );
 202       copy->GraftOutput( output );
 203       copy->Update( );
 204       output->Graft( copy->GetOutput( ) );
 205     }
 206     else
 207     {
 208       // Maximize results
 209       typename _TMaxFilter::Pointer max = _TMaxFilter::New( );
 210       max->SetInput1( output );
 211       max->SetInput2( buffer );
 212       max->InPlaceOn( );
 213       max->SetNumberOfThreads( nThreads );
 214       pg->RegisterInternalFilter( max, fw );
 215       max->GraftOutput( output );
 216       max->Update( );
 217       output->Graft( max->GetOutput( ) );
 218
 219     } // fi
 220
 221   } // rof
 222 }
 223
 224 #endif // __cpExtensions__Algorithms__MultiScaleGaussianImageFilter__hxx__
 225
 226 // eof - $RCSfile$