1 // -------------------------------------------------------------------------
2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
3 // -------------------------------------------------------------------------
5 #ifndef __CPEXTENSIONS__ALGORITHMS__MULTISCALEGAUSSIANIMAGEFILTER__HXX__
6 #define __CPEXTENSIONS__ALGORITHMS__MULTISCALEGAUSSIANIMAGEFILTER__HXX__
8 #include <vnl/vnl_vector.h>
10 #include <itkImageRegionIterator.h>
11 #include <itkImageRegionConstIterator.h>
12 #include <itkNumericTraits.h>
13 #include <itkProgressAccumulator.h>
15 #include <itkBinaryFunctorImageFilter.h>
16 #include <itkUnaryFunctorImageFilter.h>
17 #include <itkGradientMagnitudeRecursiveGaussianImageFilter.h>
18 #include <itkGradientRecursiveGaussianImageFilter.h>
19 #include <itkHessianRecursiveGaussianImageFilter.h>
21 // -------------------------------------------------------------------------
22 template< class I, class O >
23 cpExtensions::Algorithms::
24 MultiScaleGaussianImageFilter< I, O >::_Greater::
29 // -------------------------------------------------------------------------
30 template< class I, class O >
31 cpExtensions::Algorithms::
32 MultiScaleGaussianImageFilter< I, O >::_Greater::
37 // -------------------------------------------------------------------------
38 template< class I, class O >
39 bool cpExtensions::Algorithms::
40 MultiScaleGaussianImageFilter< I, O >::_Greater::
41 operator!=( const _Greater& b ) const
46 // -------------------------------------------------------------------------
47 template< class I, class O >
48 bool cpExtensions::Algorithms::
49 MultiScaleGaussianImageFilter< I, O >::_Greater::
50 operator==( const _Greater& b ) const
52 return !( *this != b );
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
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
74 typedef itk::NumericTraits< _T > _TTraits;
75 typedef typename _TTraits::ValueType _TValue;
76 typedef vnl_vector< _TValue > _TVector;
78 _TVector va( _TTraits::GetLength( ) );
79 _TVector vb( _TTraits::GetLength( ) );
81 _TTraits::AssignToArray( a, va );
82 _TTraits::AssignToArray( b, vb );
83 return( ( vb.magnitude( ) < va.magnitude( ) )? a: b );
86 // -------------------------------------------------------------------------
87 template< class I, class O >
89 cpExtensions::Algorithms::
90 MultiScaleGaussianImageFilter< I, O >::
91 SetFilterToGradient( )
94 itk::NumericTraits< typename O::PixelType >::GetLength( ) ==
97 this->m_FilterId = Self::Gradient;
99 this->m_FilterId = Self::None;
103 // -------------------------------------------------------------------------
104 template< class I, class O >
106 cpExtensions::Algorithms::
107 MultiScaleGaussianImageFilter< I, O >::
108 SetFilterToGradientMagnitude( )
110 if( itk::NumericTraits< typename O::PixelType >::GetLength( ) == 1 )
111 this->m_FilterId = Self::GradientMagnitude;
113 this->m_FilterId = Self::None;
117 // -------------------------------------------------------------------------
118 template< class I, class O >
120 cpExtensions::Algorithms::
121 MultiScaleGaussianImageFilter< I, O >::
122 SetFilterToHessian( )
124 itkExceptionMacro( << "Check for hessian definition." );
127 // -------------------------------------------------------------------------
128 template< class I, class O >
130 cpExtensions::Algorithms::
131 MultiScaleGaussianImageFilter< I, O >::
132 IsGradientFilter( ) const
134 return( this->m_FilterId == Self::Gradient );
137 // -------------------------------------------------------------------------
138 template< class I, class O >
140 cpExtensions::Algorithms::
141 MultiScaleGaussianImageFilter< I, O >::
142 IsGradientMagnitudeFilter( ) const
144 return( this->m_FilterId == Self::GradientMagnitude );
147 // -------------------------------------------------------------------------
148 template< class I, class O >
150 cpExtensions::Algorithms::
151 MultiScaleGaussianImageFilter< I, O >::
152 IsHessianFilter( ) const
154 return( this->m_FilterId == Self::Hessian );
157 // -------------------------------------------------------------------------
158 template< class I, class O >
160 cpExtensions::Algorithms::
161 MultiScaleGaussianImageFilter< I, O >::
162 AddScale( const double& s )
164 if( this->m_Scales.insert( s ).second )
168 // -------------------------------------------------------------------------
169 template< class I, class O >
171 cpExtensions::Algorithms::
172 MultiScaleGaussianImageFilter< I, O >::
173 GetNumberOfScales( ) const
175 return( this->m_Scales.size( ) );
178 // -------------------------------------------------------------------------
179 template< class I, class O >
180 cpExtensions::Algorithms::
181 MultiScaleGaussianImageFilter< I, O >::
182 MultiScaleGaussianImageFilter( )
185 this->SetFilterToGradientMagnitude( );
186 if( !this->IsGradientMagnitudeFilter( ) )
188 this->SetFilterToGradient( );
189 if( !this->IsGradientFilter( ) )
190 this->SetFilterToHessian( );
195 // -------------------------------------------------------------------------
196 template< class I, class O >
197 cpExtensions::Algorithms::
198 MultiScaleGaussianImageFilter< I, O >::
199 ~MultiScaleGaussianImageFilter( )
203 // -------------------------------------------------------------------------
204 template< class I, class O >
206 cpExtensions::Algorithms::
207 MultiScaleGaussianImageFilter< I, O >::
210 typedef itk::GradientRecursiveGaussianImageFilter< I, O > _TGF;
211 typedef itk::GradientMagnitudeRecursiveGaussianImageFilter< I, O > _TGMF;
212 typedef itk::HessianRecursiveGaussianImageFilter< I, O > _THF;
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 >( );
222 // -------------------------------------------------------------------------
223 template< class I, class O >
226 cpExtensions::Algorithms::
227 MultiScaleGaussianImageFilter< I, O >::
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;
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 );
243 // Progress accumulator
244 _TProgress::Pointer pg = _TProgress::New( );
245 pg->SetMiniPipelineFilter( this );
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( ) );
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( ) );
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++ )
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 );
280 buffer->Graft( filter->GetOutput( ) );
282 // Get maximum response
283 if( sIt == this->m_Scales.begin( ) )
286 typename _TCopyFilter::Pointer copy = _TCopyFilter::New( );
287 copy->SetInput( buffer );
289 copy->SetNumberOfThreads( nThreads );
290 pg->RegisterInternalFilter( copy, fw );
291 copy->GraftOutput( output );
293 output->Graft( copy->GetOutput( ) );
298 typename _TMaxFilter::Pointer max = _TMaxFilter::New( );
299 max->SetInput1( output );
300 max->SetInput2( buffer );
302 max->SetNumberOfThreads( nThreads );
303 pg->RegisterInternalFilter( max, fw );
304 max->GraftOutput( output );
306 output->Graft( max->GetOutput( ) );
313 #endif // __CPEXTENSIONS__ALGORITHMS__MULTISCALEGAUSSIANIMAGEFILTER__HXX__