[Insight-users] discreteGaussianImageFilter has problems with integer images

[Martin Urschler]

hello,

If I use the DiscreteGaussianImageFilter on e.g. a signed short itk 
image instead of a float or double itk image, then I found that the 
output image of the filter is zero.


Here is my sample code:

Int16ImageType::Pointer image =
    VolumeIOWrapper<Int16ImageType>::readITKVolume(inputFilename, "");

typedef itk::DiscreteGaussianImageFilter<Int16ImageType,
   Int16ImageType> SmoothingFilterType;
SmoothingFilterType::Pointer smoother = SmoothingFilterType::New();

SmoothingFilterType::ArrayType variances;
variances[0] = sigma_x;
variances[1] = sigma_y;
variances[2] = sigma_z;

smoother->SetVariance( variances );
smoother->SetInput( image );
ITK_EXCEPTION_CHECKED( "smoothing", smoother->Update(), EXIT_FAILURE );

Int16ImageType::Pointer smoothed_image = smoother->GetOutput();

VolumeIOWrapper<Int16ImageType>::writeITKVolume16Bit(smoothed_image,
    outputFilename, "");




I traced this error down to the itkGaussianOperator and its parent class 
the itkNeighborhoodOperator. The Gaussian operator correctly creates 
discrete gauss kernel coefficients that have floating point values 
between 0 and 1 (method "GenerateCoefficients()" in 
itkGaussianOperator.txx). However the neighborhood operator that stores 
these coefficients has signed short as pixel type due to the 
instantiation in itkDiscreteGaussianImageFilter and therefore the 
"FillCenteredDirectional(const CoefficientVector &coeff)" method of 
itkNeighborhoodOperator assigns the float coefficients to signed short, 
thereby cutting off everything behind the comma, leaving 0s for the 
operator.

this happens here:

for (data = data.Begin(); data < data.End(); ++data, ++it)
{
   *data = static_cast<TPixel>(*it);
   // *it are the float coefficients, TPixel is signed short!
}


IMHO this is a bug, since I don't understand why it shouldn't be 
possible to convolve a signed short image with a float kernel resulting 
in a signed short image. The only problem that I see are rounding errors 
due to the implicit cast from float back to signed short.

My application has to deal with large volumes (500^3 voxel), therefore I 
don't want to cast my volumes to float for smoothing!

So, what should I do in this case? Hack a little bit in the interiors of 
the GaussianImageFilter and the NeighborhoodOperator? Form a bug request?

thanks,
Martin