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 <itkGradientMagnitudeRecursiveGaussianImageFilter.h>
  18 #include <itkGradientRecursiveGaussianImageFilter.h>
  19 #include <itkHessianRecursiveGaussianImageFilter.h>
  20
  21 // -------------------------------------------------------------------------
  22 template< class I, class O >
  23 cpExtensions::Algorithms::
  24 MultiScaleGaussianImageFilter< I, O >::_Greater::
  25 _Greater( )
  26 {
  27 }
  28
  29 // -------------------------------------------------------------------------
  30 template< class I, class O >
  31 cpExtensions::Algorithms::
  32 MultiScaleGaussianImageFilter< I, O >::_Greater::
  33 ~_Greater( )
  34 {
  35 }
  36
  37 // -------------------------------------------------------------------------
  38 template< class I, class O >
  39 bool cpExtensions::Algorithms::
  40 MultiScaleGaussianImageFilter< I, O >::_Greater::
  41 operator!=( const _Greater& b ) const
  42 {
  43   return( false );
  44 }
  45
  46 // -------------------------------------------------------------------------
  47 template< class I, class O >
  48 bool cpExtensions::Algorithms::
  49 MultiScaleGaussianImageFilter< I, O >::_Greater::
  50 operator==( const _Greater& b ) const
  51 {
  52   return !( *this != b );
  53 }
  54
  55 // -------------------------------------------------------------------------
  56 template< class I, class O >
  57 typename cpExtensions::Algorithms::
  58 MultiScaleGaussianImageFilter< I, O >::_Greater::
  59 _T cpExtensions::Algorithms::
  60 MultiScaleGaussianImageFilter< I, O >::_Greater::
  61 operator()( const _T& a ) const
  62 {
  63   return( a );
  64 }
  65
  66 // -------------------------------------------------------------------------
  67 template< class I, class O >
  68 typename cpExtensions::Algorithms::
  69 MultiScaleGaussianImageFilter< I, O >::_Greater::
  70 _T cpExtensions::Algorithms::
  71 MultiScaleGaussianImageFilter< I, O >::_Greater::
  72 operator()( const _T& a, const _T& b ) const
  73 {
  74   typedef itk::NumericTraits< _T >     _TTraits;
  75   typedef typename _TTraits::ValueType _TValue;
  76   typedef vnl_vector< _TValue >        _TVector;
  77
  78   _TVector va( _TTraits::GetLength( ) );
  79   _TVector vb( _TTraits::GetLength( ) );
  80
  81   _TTraits::AssignToArray( a, va );
  82   _TTraits::AssignToArray( b, vb );
  83   return( ( vb.magnitude( ) < va.magnitude( ) )? a: b );
  84 }
  85
  86 // -------------------------------------------------------------------------
  87 template< class I, class O >
  88 void
  89 cpExtensions::Algorithms::
  90 MultiScaleGaussianImageFilter< I, O >::
  91 SetFilterToGradient( )
  92 {
  93   if(
  94     itk::NumericTraits< typename O::PixelType >::GetLength( ) ==
  95     I::ImageDimension
  96     )
  97     this->m_FilterId = Self::Gradient;
  98   else
  99     this->m_FilterId = Self::None;
 100   this->Modified( );
 101 }
 102
 103 // -------------------------------------------------------------------------
 104 template< class I, class O >
 105 void
 106 cpExtensions::Algorithms::
 107 MultiScaleGaussianImageFilter< I, O >::
 108 SetFilterToGradientMagnitude( )
 109 {
 110   if( itk::NumericTraits< typename O::PixelType >::GetLength( ) == 1 )
 111     this->m_FilterId = Self::GradientMagnitude;
 112   else
 113     this->m_FilterId = Self::None;
 114   this->Modified( );
 115 }
 116
 117 // -------------------------------------------------------------------------
 118 template< class I, class O >
 119 void
 120 cpExtensions::Algorithms::
 121 MultiScaleGaussianImageFilter< I, O >::
 122 SetFilterToHessian( )
 123 {
 124   itkExceptionMacro( << "Check for hessian definition." );
 125 }
 126
 127 // -------------------------------------------------------------------------
 128 template< class I, class O >
 129 bool
 130 cpExtensions::Algorithms::
 131 MultiScaleGaussianImageFilter< I, O >::
 132 IsGradientFilter( ) const
 133 {
 134   return( this->m_FilterId == Self::Gradient );
 135 }
 136
 137 // -------------------------------------------------------------------------
 138 template< class I, class O >
 139 bool
 140 cpExtensions::Algorithms::
 141 MultiScaleGaussianImageFilter< I, O >::
 142 IsGradientMagnitudeFilter( ) const
 143 {
 144   return( this->m_FilterId == Self::GradientMagnitude );
 145 }
 146
 147 // -------------------------------------------------------------------------
 148 template< class I, class O >
 149 bool
 150 cpExtensions::Algorithms::
 151 MultiScaleGaussianImageFilter< I, O >::
 152 IsHessianFilter( ) const
 153 {
 154   return( this->m_FilterId == Self::Hessian );
 155 }
 156
 157 // -------------------------------------------------------------------------
 158 template< class I, class O >
 159 void
 160 cpExtensions::Algorithms::
 161 MultiScaleGaussianImageFilter< I, O >::
 162 AddScale( const double& s )
 163 {
 164   if( this->m_Scales.insert( s ).second )
 165     this->Modified( );
 166 }
 167
 168 // -------------------------------------------------------------------------
 169 template< class I, class O >
 170 unsigned long
 171 cpExtensions::Algorithms::
 172 MultiScaleGaussianImageFilter< I, O >::
 173 GetNumberOfScales( ) const
 174 {
 175   return( this->m_Scales.size( ) );
 176 }
 177
 178 // -------------------------------------------------------------------------
 179 template< class I, class O >
 180 cpExtensions::Algorithms::
 181 MultiScaleGaussianImageFilter< I, O >::
 182 MultiScaleGaussianImageFilter( )
 183   : Superclass( )
 184 {
 185   this->SetFilterToGradientMagnitude( );
 186   if( !this->IsGradientMagnitudeFilter( ) )
 187   {
 188     this->SetFilterToGradient( );
 189     if( !this->IsGradientFilter( ) )
 190       this->SetFilterToHessian( );
 191
 192   } // fi
 193 }
 194
 195 // -------------------------------------------------------------------------
 196 template< class I, class O >
 197 cpExtensions::Algorithms::
 198 MultiScaleGaussianImageFilter< I, O >::
 199 ~MultiScaleGaussianImageFilter( )
 200 {
 201 }
 202
 203 // -------------------------------------------------------------------------
 204 template< class I, class O >
 205 void
 206 cpExtensions::Algorithms::
 207 MultiScaleGaussianImageFilter< I, O >::
 208 GenerateData( )
 209 {
 210   typedef itk::GradientRecursiveGaussianImageFilter< I, O >          _TGF;
 211   typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< I, O > _TGMF;
 212   typedef itk::HessianRecursiveGaussianImageFilter< I, O >           _THF;
 213
 214   if( this->IsGradientFilter( ) )
 215     this->_GenerateData< _TGF >( );
 216   else if( this->IsGradientMagnitudeFilter( ) )
 217     this->_GenerateData< _TGMF >( );
 218   else if( this->IsHessianFilter( ) )
 219     this->_GenerateData< _THF >( );
 220 }
 221
 222 // -------------------------------------------------------------------------
 223 template< class I, class O >
 224 template< class F >
 225 void
 226 cpExtensions::Algorithms::
 227 MultiScaleGaussianImageFilter< I, O >::
 228 _GenerateData( )
 229 {
 230   // Some types
 231   typedef itk::BinaryFunctorImageFilter< O, O, O, _Greater > _TMaxFilter;
 232   typedef itk::UnaryFunctorImageFilter< O, O, _Greater >     _TCopyFilter;
 233   typedef itk::ImageRegionConstIterator< O >                 _TConstIt;
 234   typedef itk::ImageRegionIterator< O >                      _TIt;
 235   typedef itk::ProgressAccumulator                           _TProgress;
 236
 237   // Some values
 238   typename I::ConstPointer input = this->GetInput( );
 239   typename O::Pointer output = this->GetOutput( );
 240   unsigned int nThreads = this->GetNumberOfThreads( );
 241   float fw = float( 1 ) / float( this->m_Scales.size( ) << 1 );
 242
 243   // Progress accumulator
 244   _TProgress::Pointer pg = _TProgress::New( );
 245   pg->SetMiniPipelineFilter( this );
 246
 247   // Copy image information
 248   output->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
 249   output->SetRequestedRegion( input->GetRequestedRegion( ) );
 250   output->SetBufferedRegion( input->GetBufferedRegion( ) );
 251   output->SetSpacing( input->GetSpacing( ) );
 252   output->SetOrigin( input->GetOrigin( ) );
 253   output->SetDirection( input->GetDirection( ) );
 254   output->Allocate( );
 255
 256   // Intermediary buffer
 257   typename O::Pointer buffer = O::New( );
 258   buffer->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
 259   buffer->SetRequestedRegion( input->GetRequestedRegion( ) );
 260   buffer->SetBufferedRegion( input->GetBufferedRegion( ) );
 261   buffer->SetSpacing( input->GetSpacing( ) );
 262   buffer->SetOrigin( input->GetOrigin( ) );
 263   buffer->SetDirection( input->GetDirection( ) );
 264   buffer->Allocate( );
 265
 266   // Perform all scales
 267   _TIt gIt( output, output->GetRequestedRegion( ) );
 268   TScalesContainer::const_iterator sIt = this->m_Scales.begin( );
 269   for( ; sIt != this->m_Scales.end( ); sIt++ )
 270   {
 271     // Single scale filter
 272     typename F::Pointer filter = F::New( );
 273     filter->SetInput( input );
 274     filter->SetNormalizeAcrossScale( true );
 275     filter->SetNumberOfThreads( nThreads );
 276     filter->SetSigma( *sIt );
 277     pg->RegisterInternalFilter( filter, fw );
 278     filter->GraftOutput( buffer );
 279     filter->Update( );
 280     buffer->Graft( filter->GetOutput( ) );
 281
 282     // Get maximum response
 283     if( sIt == this->m_Scales.begin( ) )
 284     {
 285       // Copy first result
 286       typename _TCopyFilter::Pointer copy = _TCopyFilter::New( );
 287       copy->SetInput( buffer );
 288       copy->InPlaceOff( );
 289       copy->SetNumberOfThreads( nThreads );
 290       pg->RegisterInternalFilter( copy, fw );
 291       copy->GraftOutput( output );
 292       copy->Update( );
 293       output->Graft( copy->GetOutput( ) );
 294     }
 295     else
 296     {
 297       // Maximize results
 298       typename _TMaxFilter::Pointer max = _TMaxFilter::New( );
 299       max->SetInput1( output );
 300       max->SetInput2( buffer );
 301       max->InPlaceOn( );
 302       max->SetNumberOfThreads( nThreads );
 303       pg->RegisterInternalFilter( max, fw );
 304       max->GraftOutput( output );
 305       max->Update( );
 306       output->Graft( max->GetOutput( ) );
 307
 308     } // fi
 309
 310   } // rof
 311 }
 312
 313 #endif // __CPEXTENSIONS__ALGORITHMS__MULTISCALEGAUSSIANIMAGEFILTER__HXX__
 314
 315 // eof - $RCSfile$