[Insight-users] Recursive Gaussian filter oscillations

[Tom Vercauteren]

Hi all,

I would like to use an accurate and fast Gaussian filter.

I have been using ITK's implementation of Deriche's filter for a while
now. It gives good results in many cases but now I get bothered by the
oscillations of the output (I have attached a code snippet showing
these oscillations).

I have been told that the accuracy of Deriche's filters could be
enhanced by considering higher orders in the IIR filtering. Are there
any plans to incorporate such high order recursive Gaussian filters in
ITK?

On a related topic, I have noticed that the
SmoothingRecursiveGaussianImageFilter was using an unnecessary
CastImageFilter. I have therefore modified this filter to avoid this
casting operation. I could contribute this piece of code to ITK if
someone is interested.

Best,
Tom Vercauteren

///////////////////////////////////////// Code snippet ///////////////////////

#include <itkImage.h>
#include <itkRecursiveGaussianImageFilter.h>
#include <itkImageRegionIteratorWithIndex.h>
#include <itkImageRegionConstIterator.h>


int main( int argc, char **argv )
{
  std::cout << "Test using a 1-D image" << std::endl;

  typedef double                        PixelType;
  typedef itk::Image< PixelType, 1 >    ImageType;

  typedef ImageType::SizeType           SizeType;
  typedef ImageType::IndexType          IndexType;
  typedef ImageType::RegionType         RegionType;
  typedef ImageType::SpacingType        SpacingType;

  typedef itk::NumericTraits< PixelType >::RealType    PixelRealType;

  SizeType size;
  size[0] = 21;

  IndexType start;
  start[0] = 0;

  RegionType region;
  region.SetIndex( start );
  region.SetSize( size );

  SpacingType spacing;
  spacing[0] = 1.0;

  ImageType::Pointer inputImage = ImageType::New();
  inputImage->SetRegions( region );
  inputImage->Allocate();
  inputImage->SetSpacing( spacing );
  inputImage->FillBuffer( itk::NumericTraits< PixelType >::Zero );

  IndexType index;
  index[0] = ( size[0] - 1 ) / 2;  // the middle pixel

  inputImage->SetPixel( index, static_cast< PixelType >( 1000.0 ) );

  typedef itk::RecursiveGaussianImageFilter<
                            ImageType, ImageType > FilterType;

  FilterType::Pointer filter = FilterType::New();

  filter->SetNormalizeAcrossScale( true );
  const double sigmaA = 2.0;
  filter->SetSigma( sigmaA );
  filter->SetZeroOrder ();

  filter->SetInput( inputImage );
  filter->Update();


  PixelType* ptrmin = std::min_element( filter->GetOutput()->GetBufferPointer(),
					filter->GetOutput()->GetBufferPointer() + size[0] );

  std::cout << "Min of filtered output is: "<<*ptrmin<<std::endl;

  PixelType* ptrmax = std::max_element( filter->GetOutput()->GetBufferPointer(),
					filter->GetOutput()->GetBufferPointer() + size[0] );

  std::cout << "Max of filtered output is: "<<*ptrmax<<std::endl;

  std::cout<<"[ ";
  for (unsigned int i=0; i<size[0]; i++)
    {
    index[0] = i;
    std::cout<<filter->GetOutput()->GetPixel( index )<<" ";
    }
  std::cout<<std::endl;
}