yet another refactoring

[cpPlugins.git] / lib / cpExtensions / Algorithms / MultiScaleGaussianImageFilter.hxx

diff --git a/lib/cpExtensions/Algorithms/MultiScaleGaussianImageFilter.hxx b/lib/cpExtensions/Algorithms/MultiScaleGaussianImageFilter.hxx
new file mode 100644 (file)
index 0000000..985cbec
--- /dev/null
+++ b/lib/cpExtensions/Algorithms/MultiScaleGaussianImageFilter.hxx
@@ -0,0 +1,226 @@
+// -------------------------------------------------------------------------
+// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
+// -------------------------------------------------------------------------
+
+#ifndef __cpExtensions__Algorithms__MultiScaleGaussianImageFilter__hxx__
+#define __cpExtensions__Algorithms__MultiScaleGaussianImageFilter__hxx__
+
+#include <vnl/vnl_vector.h>
+
+#include <itkImageRegionIterator.h>
+#include <itkImageRegionConstIterator.h>
+#include <itkNumericTraits.h>
+#include <itkProgressAccumulator.h>
+
+#include <itkBinaryFunctorImageFilter.h>
+#include <itkUnaryFunctorImageFilter.h>
+#include <itkGradientRecursiveGaussianImageFilter.h>
+
+// -------------------------------------------------------------------------
+template< class I, class O >
+cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::_Greater::
+_Greater( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class I, class O >
+cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::_Greater::
+~_Greater( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class I, class O >
+bool cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::_Greater::
+operator!=( const _Greater& b ) const
+{
+  return( false );
+}
+
+// -------------------------------------------------------------------------
+template< class I, class O >
+bool cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::_Greater::
+operator==( const _Greater& b ) const
+{
+  return !( *this != b );
+}
+
+// -------------------------------------------------------------------------
+template< class I, class O >
+typename cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::_Greater::
+_T cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::_Greater::
+operator()( const _T& a ) const
+{
+  return( a );
+}
+
+// -------------------------------------------------------------------------
+template< class I, class O >
+typename cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::_Greater::
+_T cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::_Greater::
+operator()( const _T& a, const _T& b ) const
+{
+  typedef itk::NumericTraits< _T >     _TTraits;
+  typedef typename _TTraits::ValueType _TValue;
+  typedef vnl_vector< _TValue >        _TVector;
+
+  _TVector va( _TTraits::GetLength( ) );
+  _TVector vb( _TTraits::GetLength( ) );
+
+  _TTraits::AssignToArray( a, va );
+  _TTraits::AssignToArray( b, vb );
+  return( ( vb.magnitude( ) < va.magnitude( ) )? a: b );
+}
+
+// -------------------------------------------------------------------------
+template< class I, class O >
+void
+cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::
+AddScale( const double& s )
+{
+  if( this->m_Scales.insert( s ).second )
+    this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class I, class O >
+unsigned long
+cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::
+GetNumberOfScales( ) const
+{
+  return( this->m_Scales.size( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class I, class O >
+cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::
+MultiScaleGaussianImageFilter( )
+  : Superclass( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class I, class O >
+cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::
+~MultiScaleGaussianImageFilter( )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class I, class O >
+void
+cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::
+GenerateData( )
+{
+  typedef itk::GradientRecursiveGaussianImageFilter< I, O > _TGF;
+  this->_GenerateData< _TGF >( );
+}
+
+// -------------------------------------------------------------------------
+template< class I, class O >
+template< class F >
+void
+cpExtensions::Algorithms::
+MultiScaleGaussianImageFilter< I, O >::
+_GenerateData( )
+{
+  // Some types
+  typedef itk::BinaryFunctorImageFilter< O, O, O, _Greater > _TMaxFilter;
+  typedef itk::UnaryFunctorImageFilter< O, O, _Greater >     _TCopyFilter;
+  typedef itk::ImageRegionConstIterator< O >                 _TConstIt;
+  typedef itk::ImageRegionIterator< O >                      _TIt;
+  typedef itk::ProgressAccumulator                           _TProgress;
+
+  // Some values
+  typename I::ConstPointer input = this->GetInput( );
+  typename O::Pointer output = this->GetOutput( );
+  unsigned int nThreads = this->GetNumberOfThreads( );
+  float fw = float( 1 ) / float( this->m_Scales.size( ) << 1 );
+
+  // Progress accumulator
+  _TProgress::Pointer pg = _TProgress::New( );
+  pg->SetMiniPipelineFilter( this );
+
+  // Copy image information
+  output->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
+  output->SetRequestedRegion( input->GetRequestedRegion( ) );
+  output->SetBufferedRegion( input->GetBufferedRegion( ) );
+  output->SetSpacing( input->GetSpacing( ) );
+  output->SetOrigin( input->GetOrigin( ) );
+  output->SetDirection( input->GetDirection( ) );
+  output->Allocate( );
+
+  // Intermediary buffer
+  typename O::Pointer buffer = O::New( );
+  buffer->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) );
+  buffer->SetRequestedRegion( input->GetRequestedRegion( ) );
+  buffer->SetBufferedRegion( input->GetBufferedRegion( ) );
+  buffer->SetSpacing( input->GetSpacing( ) );
+  buffer->SetOrigin( input->GetOrigin( ) );
+  buffer->SetDirection( input->GetDirection( ) );
+  buffer->Allocate( );
+
+  // Perform all scales
+  _TIt gIt( output, output->GetRequestedRegion( ) );
+  TScalesContainer::const_iterator sIt = this->m_Scales.begin( );
+  for( ; sIt != this->m_Scales.end( ); sIt++ )
+  {
+    // Single scale filter
+    typename F::Pointer filter = F::New( );
+    filter->SetInput( input );
+    filter->SetNormalizeAcrossScale( true );
+    filter->SetNumberOfThreads( nThreads );
+    filter->SetSigma( *sIt );
+    pg->RegisterInternalFilter( filter, fw );
+    filter->GraftOutput( buffer );
+    filter->Update( );
+    buffer->Graft( filter->GetOutput( ) );
+
+    // Get maximum response
+    if( sIt == this->m_Scales.begin( ) )
+    {
+      // Copy first result
+      typename _TCopyFilter::Pointer copy = _TCopyFilter::New( );
+      copy->SetInput( buffer );
+      copy->InPlaceOff( );
+      copy->SetNumberOfThreads( nThreads );
+      pg->RegisterInternalFilter( copy, fw );
+      copy->GraftOutput( output );
+      copy->Update( );
+      output->Graft( copy->GetOutput( ) );
+    }
+    else
+    {
+      // Maximize results
+      typename _TMaxFilter::Pointer max = _TMaxFilter::New( );
+      max->SetInput1( output );
+      max->SetInput2( buffer );
+      max->InPlaceOn( );
+      max->SetNumberOfThreads( nThreads );
+      pg->RegisterInternalFilter( max, fw );
+      max->GraftOutput( output );
+      max->Update( );
+      output->Graft( max->GetOutput( ) );
+
+    } // fi
+
+  } // rof
+}
+
+#endif // __cpExtensions__Algorithms__MultiScaleGaussianImageFilter__hxx__
+
+// eof - $RCSfile$