ENH: Add FFTDiscreteGaussianImageFilter

[tbirdso]

Changes:

• Implements FFTDiscreteGaussianImageFilter as a subclass of DiscreteGaussianImageFilter to perform convolution of an ND single-channel input image with a Gaussian kernel of fixed width in the Fourier domain
• Implements FFTDiscreteGaussianImageFilterEnums to select among kernel generation methods, where COMBINED_OPERATORS composes through vector multiplication the exact same ND kernel as used in DiscreteGaussianImageFilter in image format, while IMAGE_SOURCE relies on the (potentially faster) method of directly generating the kernel image through itk::GaussianImageSource
• Exposes protected methods for generating the kernel from error/width/spacing inputs in the base DiscreteGaussianImageFilter class so that FFTDiscreteGaussianImageFilter can use the same process for generating a kernel to specifications
• Adds FFTDiscreteGaussianImageFilterFactory for manually overriding DiscreteGaussianImageFilter for FFT acceleration benefits
• Updates test coverage, including a test comparing FFTDiscreteGaussianImageFilter output to DiscreteGaussianImageFilter output with the same parameters
• Includes Python wrapping

Input:

Test output screenshot:

PR Checklist

• No API changes were made (or the changes have been approved)
• No major design changes were made (or the changes have been approved)
• Added test (or behavior not changed)
• Updated API documentation (or API not changed)\
• Added Python wrapping to new files (if any) as described in ITK Software Guide Section 9.5

Refer to the ITK Software Guide for
further development details if necessary.

[thewtex]

Looking great!

[tbirdso]

@thewtex @dzenanz A few design questions for you before formally moving to the review stage.

[tbirdso]

Pasting review notes from @jhlegarreta here (great points):

i) there is an exception thrown by the class under some condition; please exercise it;
ii) the Set methods overriding the base class should be exercised, hopefully creating the necessary cases in the CMakeLists.txt so that any if/else condition in the implementation gets exercised

Will finish addressing these and then mark as ready

[tbirdso]

@thewtex @dzenanz Ping for re-review

Note that I've consolidated and exposed functionality from itk::DiscreteGaussianImageFilter for generating a kernel radius to error specifications with itk::GaussianOperator to avoid duplication in itk::FFTDiscreteGaussianImageFilter.

UseImageSpacing and MaximumError are now used for kernel generation.

[tbirdso]

For looking forward, I've gone ahead and added a factory override mechanism so that developers can opt in to having classes like itk::MultiResolutionPyramidImageFilter rely on FFTDiscreteGaussianImageFilter rather than DiscreteGaussianImageFilter for blurring. We can combine this with GPU-accelerated FFT overrides to speed up blurring with very large images and/or very large kernels.

If more discussion is required around the factory mechanism I'm fine with splitting it off in a separate PR.

[thewtex]

@tbirdso looking awesome! 🎇

Tests should also be added to verify that we get the same result with the DiscreteGaussianImageFilter and the FFTDiscreteGaussianImageFilter.

[dzenanz]

Mostly looks good, some minor suggestions.

[tbirdso]

@thewtex 's suggestion

Tests should also be added to verify that we get the same result with the DiscreteGaussianImageFilter and the FFTDiscreteGaussianImageFilter.

was well-founded, because on more testing I am getting different images for DiscreteGaussianImageFilter and FFTDiscreteGaussianImageFilter.

DiscreteGaussianImageFilterComparisonTest:

FFTDiscreteGaussianImageFilterTest:

Diff:

It looks like the main issue here is that each filter relies on a different kernel source with what appears to be fundamentally different mechanisms for Gaussian kernel generation, leading to different kernel values from the same parameters.

DiscreteGaussianImageFilter creates and applies iterative 1D itk::GaussianOperators which rely on Bessel functions for defining kernel values. FFTDiscreteGaussianImageFilter requires an ND kernel image that it attempts to generate through itk::GaussianImageSource based on the kernel radii from itk::GaussianOperators constructed to the user error specifications, however the image source relies on itk::GaussianSpatialFunction to generate kernel elements which seems to produce very different values from itk::GaussianOperator.

@thewtex @dzenanz could you take a look and weigh in on this? Two potential options that I see for moving forward:

• Dig deeper into the differences in the itk::GaussianOperator and itk::GaussianSpatialFunction procedures to understand how parameters may need to be mapped between the two for similar outputs;
• Successively generate 1D itk::GaussianOperators, copy to itk::VariableSizeMatrix objects, cross-multiply to produce an ND kernel, iteratively copy elements to a kernel image, and update image metadata to match the input image. This would increase prep time as we may need to run 3- or 4-D matrix multiplication to get "large" kernels. However I would favor this approach as it retains a dependency on itk::GaussianOperator only and time for kernel generation was previously anticipated to be small compared to convolution.

[thewtex]

Successively generate 1D itk::GaussianOperators, copy to itk::VariableSizeMatrix objects, cross-multiply to produce an ND kernel, iteratively copy elements to a kernel image, and update image metadata to match the input image. This would increase prep time as we may need to run 3- or 4-D matrix multiplication to get "large" kernels. However I would favor this approach as it retains a dependency on itk::GaussianOperator only and time for kernel generation was previously anticipated to be small compared to convolution.

👍

If it is not too much trouble and there is a performance benefit, we could keep GaussianImageSource as an option. But GaussianOperator makes sense if we are proposing a factory-enabled replacement for DiscreteGaussianImageFilter.

[tbirdso]

Thanks @thewtex . Do you know if we multiply operators to get a kernel anywhere else in ITK? I'm having trouble finding examples of how to do this efficiently.

EDIT: It looks like VariableSizeMatrix implements 2D arrays only while we are seeking a method for 3-, 4-, or N-D matrix multiplication to properly construct the kernel. Can manually implement if necessary but it would be nice to rely on a standardized method.

[thewtex]

@tbirdso nothing comes to mind offhand, but I agree -- taking a look to re-use any existing functionality makes sense, but we can also use a manual implementation.

[tbirdso]

@dzenanz @thewtex Ping for re-review.

New changes:

• Can use FFTDiscreteGaussianImageFilterEnums to specify how FFTDiscreteGaussianImageFilter should generate its kernel image for FFT convolution
• If the kernel source is set to FFTDiscreteGaussianImageFilterEnums::KernelSource::COMBINED_OPERATORS then the filter will generate an empty image and manually perform vector multiplication to compose an ND Gaussian kernel image from N 1D directional kernels
• Tests are added to exercise different kernel methods and to compare output to DiscreteGaussianImageFilter when the kernel source is COMBINED_OPERATORS

[thewtex]

Beautiful!

One minor issue noted inline.

[dzenanz]

Looks good, just a few nitpicks.

[thewtex]

LGTM!