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 <itkGradientRecursiveGaussianImageFilter.h>
19 // -------------------------------------------------------------------------
20 template< class I, class O >
21 cpExtensions::Algorithms::
22 MultiScaleGaussianImageFilter< I, O >::_Greater::
27 // -------------------------------------------------------------------------
28 template< class I, class O >
29 cpExtensions::Algorithms::
30 MultiScaleGaussianImageFilter< I, O >::_Greater::
35 // -------------------------------------------------------------------------
36 template< class I, class O >
37 bool cpExtensions::Algorithms::
38 MultiScaleGaussianImageFilter< I, O >::_Greater::
39 operator!=( const _Greater& b ) const
44 // -------------------------------------------------------------------------
45 template< class I, class O >
46 bool cpExtensions::Algorithms::
47 MultiScaleGaussianImageFilter< I, O >::_Greater::
48 operator==( const _Greater& b ) const
50 return !( *this != b );
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
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
72 typedef itk::NumericTraits< _T > _TTraits;
73 typedef typename _TTraits::ValueType _TValue;
74 typedef vnl_vector< _TValue > _TVector;
76 _TVector va( _TTraits::GetLength( ) );
77 _TVector vb( _TTraits::GetLength( ) );
79 _TTraits::AssignToArray( a, va );
80 _TTraits::AssignToArray( b, vb );
81 return( ( vb.magnitude( ) < va.magnitude( ) )? a: b );
84 // -------------------------------------------------------------------------
85 template< class I, class O >
87 cpExtensions::Algorithms::
88 MultiScaleGaussianImageFilter< I, O >::
89 AddScale( const double& s )
91 if( this->m_Scales.insert( s ).second )
95 // -------------------------------------------------------------------------
96 template< class I, class O >
98 cpExtensions::Algorithms::
99 MultiScaleGaussianImageFilter< I, O >::
100 GetNumberOfScales( ) const
102 return( this->m_Scales.size( ) );
105 // -------------------------------------------------------------------------
106 template< class I, class O >
107 cpExtensions::Algorithms::
108 MultiScaleGaussianImageFilter< I, O >::
109 MultiScaleGaussianImageFilter( )
114 // -------------------------------------------------------------------------
115 template< class I, class O >
116 cpExtensions::Algorithms::
117 MultiScaleGaussianImageFilter< I, O >::
118 ~MultiScaleGaussianImageFilter( )
122 // -------------------------------------------------------------------------
123 template< class I, class O >
125 cpExtensions::Algorithms::
126 MultiScaleGaussianImageFilter< I, O >::
129 typedef itk::GradientRecursiveGaussianImageFilter< I, O > _TGF;
130 this->_GenerateData< _TGF >( );
133 // -------------------------------------------------------------------------
134 template< class I, class O >
137 cpExtensions::Algorithms::
138 MultiScaleGaussianImageFilter< I, O >::
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;
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 );
154 // Progress accumulator
155 _TProgress::Pointer pg = _TProgress::New( );
156 pg->SetMiniPipelineFilter( this );
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( ) );
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( ) );
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++ )
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 );
191 buffer->Graft( filter->GetOutput( ) );
193 // Get maximum response
194 if( sIt == this->m_Scales.begin( ) )
197 typename _TCopyFilter::Pointer copy = _TCopyFilter::New( );
198 copy->SetInput( buffer );
200 copy->SetNumberOfThreads( nThreads );
201 pg->RegisterInternalFilter( copy, fw );
202 copy->GraftOutput( output );
204 output->Graft( copy->GetOutput( ) );
209 typename _TMaxFilter::Pointer max = _TMaxFilter::New( );
210 max->SetInput1( output );
211 max->SetInput2( buffer );
213 max->SetNumberOfThreads( nThreads );
214 pg->RegisterInternalFilter( max, fw );
215 max->GraftOutput( output );
217 output->Graft( max->GetOutput( ) );
224 #endif // __CPEXTENSIONS__ALGORITHMS__MULTISCALEGAUSSIANIMAGEFILTER__HXX__