Complete Failure (NaN) for DiscreteGaussianImageFilter with large sigma and large images

Engineering
1

I’m having issues with a complete failure of DiscreteGaussianImageFilter as implemented in SmoothImage in ANTs as well as its internal use in the antsRegistration pipeline.

For some combinations of large resolution/large sigma, the filter produces no output (including no warnings, as what happens if the gaussian approximation is truncated to a limited with). In SmoothImage this is reasonably easy to catch, as you get an empty file, but in the case of inside antsRegistration, it completely derails the registration.

Some tests and discussion of this bug have gone on here: https://github.com/ANTsX/ANTs/issues/538

and I reported the bug to ITK Jira here: https://insightsoftwareconsortium.atlassian.net/browse/ITK-3652

As image resolutions increase we’re going to be more and more likely to run into this problem so I’m hoping I can get some help locating/solving this.

2

Hi Gabriel,

The error message noted in the bug report may identify the issue:

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.

i.e. SmoothImage may need to call SetMaximumKernelWidth.

But, for smoothing with large Gaussian, consider using these filters instead for efficiency:

HTH,
Matt

3

Hi Matt,

Thanks.

To be clear, there is no warning output for a “large enough” sigma for the blur. There is simply no output.

In addition, increasing the SetMaximumKernelWidth arbitrarily large does not expose any warning, or change the output.

4

The next step to improve itk::DiscreteGaussianImageFilter are to create a minimal pipeline that reproduces the behavior in a single .cxx with a file reader itk::DiscreteGaussianImageFilter, and a file writer, the corresponding CMakeLists.txt, and the data that causes the issue.

However, again, for the purpose of smoothing large images with large sigma, use itk::SmoothingRecursiveYvvGaussianImageFilter.

5

I have created a reproducible test with:
http://www.bic.mni.mcgill.ca/~vfonov/icbm/2009/mni_icbm152_nlin_sym_09b_nifti.zip
and
https://itk.org/Doxygen/html/Examples_2Filtering_2DiscreteGaussianImageFilter_8cxx-example.html
(only modified to do 2D-3D).

Success (warning is expected, since I have truncated the gaussian too small)

./DiscreteGaussianImageFilter mni_icbm152_nlin_sym_09b/mni_
icbm152_t1_tal_nlin_sym_09b_hires.nii smooth.nii.gz 177 32
WARNING: In /storage/software/ANTs-build/ITKv5/Modules/Core/Common/include/itkGaussianOperator.hxx, line 61
GaussianOperator (0x7ffeb651b990): 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.

WARNING: In /storage/software/ANTs-build/ITKv5/Modules/Core/Common/include/itkGaussianOperator.hxx, line 61
GaussianOperator (0x7ffeb651b990): 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.

WARNING: In /storage/software/ANTs-build/ITKv5/Modules/Core/Common/include/itkGaussianOperator.hxx, line 61
GaussianOperator (0x7ffeb651b990): 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.

WARNING: In /storage/software/ANTs-build/ITKv5/Modules/Core/Common/include/itkGaussianOperator.hxx, line 61
GaussianOperator (0x561256913710): 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.

WARNING: In /storage/software/ANTs-build/ITKv5/Modules/Core/Common/include/itkGaussianOperator.hxx, line 61
GaussianOperator (0x561256913678): 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.

WARNING: In /storage/software/ANTs-build/ITKv5/Modules/Core/Common/include/itkGaussianOperator.hxx, line 61
GaussianOperator (0x5612569135e0): 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.

Failure:

./DiscreteGaussianImageFilter mni_icbm152_nlin_sym_09b/mni_icbm152_t1_tal_nlin_sym_09b_hires.nii smooth.nii.gz 178 32
<No warning outputs and empty file>

After a bit of code annotating with itkWarningMacro, it seems that https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.hxx#L80 breaks down:

For the working case:

 ./DiscreteGaussianImageFilter mni_icbm152_nlin_sym_09b/mni_icbm152_t1_tal_nlin_sym_09b_hires.nii smooth.nii.gz 177 32
WARNING: In /storage/software/ANTs-build/ITKv5/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.hxx, line 80
DiscreteGaussianImageFilter (0x557f2da2df70): Radius 32

and the failure:

./DiscreteGaussianImageFilter mni_icbm152_nlin_sym_09b/mni_icbm152_t1_tal_nlin_sym_09b_hires.nii smooth.nii.gz 178 32
WARNING: In /storage/software/ANTs-build/ITKv5/Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.hxx, line 80
DiscreteGaussianImageFilter (0x56254ec58f70): Radius 1

GetRadius is from the underlying NeighborhoodOperator code wrapped by GaussianOperator, which I haven’t had a chance to take a look at yet.

That’s all I can manage tonight, I’ll follow up tomorrow.

2 Likes
6

That’s odd…

You may check how GaussianOperator::GenerateCoefficients() is behaving with those two values for the variance, 177 and 178. The coefficients determine the radius here.

1 Like
7

Thanks for the assist @phcerdan

I’ve now gotten it down to ModifiedBesselI0 and ModifiedBesselI1 returning “inf”

So, digging into those now :slight_smile:

1 Like
8

Okay, I’ve narrowed it down.

The ModifiedBesselI0 etc functions produce very large values for large variance. After their return they’re multiplied by exp(-m_variance) which then scales the values right now. Unfortunately, there’s a tipping point where we get an overflow.

It looks like given the current implementation of the bessel function, we could roll the exp(-m_variance) into the accumulator line and cancel out the overflow. Of course, the function proper wouldn’t be a modified bessel, but rather the “coeffcient generator”

1 Like
9

It looks like scipy has an implementation that’s not as susceptible to intermediate overflows, and isn’t recursive

This means we could use this implementation, and roll the exp(-m_variance) scaling inside the calculation to avoid numerical instabilities.

2 Likes
10

And finally boost has an implementation: https://www.boost.org/doc/libs/1_67_0/libs/math/doc/html/math_toolkit/bessel/mbessel.html

And a nice reference:
https://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision/

2 Likes
11

Nice research, I would suggest submitting a patch with the implementation you find more reasonable (boost looks good?), and check that the tests are passing. The ModifiedBessell functions are local to that class, so it shouldn’t interfere with any other existing code. And if it does, is probably for good.

1 Like
12

What about c++11’s math functions:

https://en.cppreference.com/w/cpp/numeric/special_math

Is there the required function there?

13

It seems they are merged into c++17 from boost. In c++11 the special math is in a technical specification (TS) which only gcc implemented… https://en.cppreference.com/w/cpp/experimental/special_math

14

I see that now, I recalled them being part of TR1, but looks like they didn’t make into c++11.

15

Thanks @phcerdan

I’m coming to the limits of my C++ knowledge here, both boost and gcc implementation here https://github.com/gcc-mirror/gcc/blob/227887ae724776730f55198d702d539bcfafe9c5/libstdc%2B%2B-v3/include/tr1/modified_bessel_func.tcc use a bit too much specal sauce for me to successfully map to my C experience.

The closest implementation I might be able to achieve is the scipy one, since its C-like.

16

That sounds good to me! The current ModifiedBessell in ITK only uses double so, we don’t need the template magic.

17

We have duplicate implementations in the GaussianOperator and the GaussianDerivativeOperator. This is a method I would expect to come from a numeric library. The VNL library we currently use seems to not have the required bessel methods.

I would suggest declaring the method in the itk::Math namespace and header.

1 Like
18

P.P.S @matt.mccormick

itk::SmoothingRecursiveYvvGaussianImageFilter is broken with the new threading updates:

terminate called after throwing an instance of 'itk::ExceptionObject'
  what():  /storage/software/ANTs-build/ITKv5/Modules/Core/Common/include/itkImageSource.hxx:280:
itk::ERROR: RecursiveLineYvvGaussianImageFilter(0x5622a8b96620): Subclass should override this method!!! If old behavior is desired invoke this->DynamicMultiThreadingOff(); before Update() is called. The best place is in class constructor.
Aborted (core dumped)
1 Like
19

Not all remote modules have been updated to new threading model. We planned to the transition to be incremental, as the need arises. As you now ask I will take a look at that module.

1 Like
20

PR created.

1 Like