[Insight-users] discreteGaussianImageFilter has problems withinteger images

[Miller, James V (Research)]

How much round off error are you willing to accept?

The DiscreteGaussian uses separable convolution.  So N convolutions
are performed to smooth an ND image.  The intermediate results could
be kept as float so the second through N-1st convolution would not have
any round off error and the last convolution would bring the pixeltypes
back to short.  Or the results of all convolutions in your case could
be shorts and the round off error would acculumate as the separable 
convolutions are performed.

Do you have a preference?

Jim



-----Original Message-----
From: insight-users-bounces at itk.org
[mailto:insight-users-bounces at itk.org]On Behalf Of Martin Urschler
Sent: Thursday, April 28, 2005 8:59 AM
To: insight-users at itk.org
Subject: [Insight-users] discreteGaussianImageFilter has problems
withinteger images


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

_______________________________________________
Insight-users mailing list
Insight-users at itk.org
http://www.itk.org/mailman/listinfo/insight-users