SmoothImage produces no output with "large" blur kernel and high resolution image

[gdevenyi]

I found this very odd result which I cannot yet find the origin:

#1mm Istotropic Image
> SmoothImage 3 mni_icbm152_t1_tal_nlin_sym_09c.nii 16 09c.nii 1 0
WARNING: In /opt/quarantine/ANTs/git/src/build/ITKv4-install/include/ITK-4.12/itkGaussianOperator.hxx, line 61
GaussianOperator (0x7ffc6d0f4a00): Kernel size has exceeded the specified maximum width of 32 and has been truncated to 33 elements.  You can raise the maximum width using the SetMaximumKernelWidth method.

Image:

#0.5mm isotropic image
> SmoothImage 3 mni_icbm152_t1_tal_nlin_sym_09b_hires.nii.gz 16 09b.nii 1 0
#No output

Taking a look over https://github.com/ANTsX/ANTs/blob/master/Examples/SmoothImage.cxx I don't see any place to expect an overflow or other issue, so this may be an ITK bug?

[cookpa]

I found the same behavior using the ch2.nii.gz from Mricron. SmoothImage does produce an output image, but the image is blank and there is no error message. So I would guess that the ITK DiscreteGaussianImageFilter object is returning a blank image. I suspect it has an internal maximum kernel size that it won't exceed even if the maxKernelWidth option is set.

SmoothImage 3 ch2.nii.gz 16 bigSmooth1mm.nii.gz 1 0 

generates the warning but

SmoothImage 3 ch2.nii.gz 16 bigSmooth1mm.nii.gz 1 0 0.01 64

works without any warning.

I tried upsampling ch2 to 0.5 mm and then using a kernel that would have a width of 16 voxels:

SmoothImage 3 ch2Upsample.nii.gz 8 bigSmooth05mm.nii.gz 1 0 0.01 64

and this works. I think it just can't handle a 32-voxel kernel width. I get the same blank output if I do

SmoothImage 3 ch2.nii.gz 32 bigSmooth1vox.nii.gz 1 0

on the 1mm image.

[gdevenyi]

@cookpa

Seems ITK updated their issue tracker and I no longer have access, can you open up an bug upstream?

[cookpa]

Sorry, I'm not signed up on the new system either. I think they would add you if you ask, but I expect they would want a minimal example built against the latest ITK (ANTs is often some way behind HEAD)

[gdevenyi]

Still a problem in ITKv5 HEAD

[gdevenyi]

Upstreamed to https://insightsoftwareconsortium.atlassian.net/browse/ITK-3652

[gdevenyi]

Determined to be an overflow in the bessel code used to estimate the discrete gaussian filter coefficients, as discussed https://discourse.itk.org/t/complete-failure-nan-for-discretegaussianimagefilter-with-large-sigma-and-large-images/1127

Current implementation is not a particularly good version from numerical recipes in C. Replace with version from Scipy or boost, or wait for C++ standard to include the boost version, which is coming eventually.

In the meantime, looks like ITKv5 is going to move to the RecursiveGuassian filter implementation which will eliminate this error.