[Insight-users] recursiveGaussianImageFilter problem

[Martin Urschler]

hi

I have experienced a strange behaviour with the 
recursiveGaussianImageFilter. I'm performing a gaussian smoothing on a 
float input image using the following code:

template< typename TImageType >
typename SmallUtilityMethods<TImageType>::TImageTypePointer
SmallUtilityMethods<TImageType>::gaussSmoothVolume(
    const SmallUtilityMethods<TImageType>::TImageTypePointer& image,
    double sigma_x,
    double sigma_y,
    double sigma_z )
{
    typename SmallUtilityMethods<TImageType>::TImageTypePointer
       result = createNewImage( image );
    typename SmallUtilityMethods<TImageType>::TImageTypePointer
       temporary_result = createNewImage( image );

    typedef itk::RecursiveGaussianImageFilter<TImageType, TImageType>
       SmoothingFilterType;

    std::cout << "minimum of image: " <<
       getMinimumVoxelValue( image ) << std::endl;
    std::cout << "maximum of image: " <<
       getMaximumVoxelValue( image ) << std::endl;

    // implement the x convolution given input 'image' and store the
    // output in 'result'
    {
       typename SmoothingFilterType::Pointer filter =
          SmoothingFilterType::New();

       filter->SetDirection( 0 ); // 0 --> X direction
       filter->SetSigma( sigma_x );
       filter->SetOrder( SmoothingFilterType::ZeroOrder );
       filter->SetNormalizeAcrossScale( false );
       filter->SetInput( image );
       filter->SetReleaseDataFlag(true);
       filter->GraftOutput( result );
       //filter->Print(std::cout);

       ITK_EXCEPTION_CHECKED( "Smoothing image in direction X!",
          filter->Update(), 0 );
    }

    // implement the y convolution given input 'result' and store the
    // output in 'temporary_result'
    {
       typename SmoothingFilterType::Pointer filter =
          SmoothingFilterType::New();

       filter->SetDirection( 1 ); // 1 --> Y direction
       filter->SetSigma( sigma_y );
       filter->SetOrder( SmoothingFilterType::ZeroOrder );
       filter->SetNormalizeAcrossScale( false );
       filter->SetInput( result );
       filter->SetReleaseDataFlag(true);
       filter->GraftOutput( temporary_result );
       //filter->Print(std::cout);

       ITK_EXCEPTION_CHECKED( "Smoothing image in direction Y!",
          filter->Update(), 0 );
    }

    // implement the z convolution given input 'temporary_result' and
    // filter output in 'result'
    {
       typename SmoothingFilterType::Pointer filter =
          SmoothingFilterType::New();

       filter->SetDirection( 2 ); // 2 --> Z direction
       filter->SetSigma( sigma_z );
       filter->SetOrder( SmoothingFilterType::ZeroOrder );
       filter->SetNormalizeAcrossScale( false );
       filter->SetReleaseDataFlag(true);
       filter->SetInput( temporary_result );
       filter->GraftOutput( result );
       //filter->Print(std::cout);

       ITK_EXCEPTION_CHECKED( "Smoothing image in direction Z!",
          filter->Update(), 0 );
    }

    std::cout << "minimum of smoothed image: " <<
       getMinimumVoxelValue( result ) << std::endl;
    std::cout << "maximum of smoothed image: " <<
       getMaximumVoxelValue( result ) << std::endl;

    return result;
}





When looking at the minimum and maximum value of the images before and 
after filtering, I get the following result:

// before filtering
minimum of image: 0
maximum of image: 600

... gaussSmoothVolume ...

// after filtering
minimum of smoothed image: -1.26718
maximum of smoothed image: 560.338


All data sets are floats. I don't think it is intended that the 
intensity range of the input image should be left during filtering, 
since gaussian filters are weighted averages. Does anyone know where the 
-1.26718 might come from? Someone already stumbled upon this problem?

regards,
Martin