ENH: Consider requested output region in FFT convolution

Prior to this change, `itk::FFTConvolutionImageFilter` and
`itk::ConvolutionImageFilter` had the major difference where the former
would always run on the largest possible input image region whereas the
latter would run on a padded requested region. As a result spatial
convolution was much more performant when run on a small region of a
large image.

This change revises how `itk::FFTConvolutionImageFilter` processes input
and output image regions in order to perform convolution and discards
the assumption that the output requested region must always be the
largest possible region. The input image subregion is padded for the
kernel operator and then clipped of remaining image data so that forward
FFT is performed on a region exactly matching that which would be
processed under spatial convolution. The kernel image is likewise padded
to the match the input. The result of FFT multiplication and subsequent
inversion is then cropped to remove artifacts from padding, leaving only
the requested input size.

This approach has the advantage of scaling with the size of the
requested output region. FFT convolution continues to be disadvantaged
in comparison to spatial convolution where an output region is
iteratively expanded or shrunk, as under FFT convolution entirely new
complex images must be generated for any change in the output region.

Adds tests demonstrating functionality and equivalence to spatial
convolution.

Modules/Filtering/Convolution/include/itkFFTConvolutionImageFilter.h

@@ -49,6 +49,8 @@ namespace itk
  *
  * \ingroup ITKConvolution
  * \sa ConvolutionImageFilter
+ * \sa InverseDeconvolutionImageFilter
+ * \sa IterativeDeconvolutionImageFilter
  *
  */
 template <typename TInputImage,
@@ -95,6 +97,9 @@ class ITK_TEMPLATE_EXPORT FFTConvolutionImageFilter
 
   /** Internal types used by the FFT filters. */
   using InternalImageType = Image<TInternalPrecision, TInputImage::ImageDimension>;
+  using InternalRegionType = typename InternalImageType::RegionType;
+  using InternalSizeType = typename InternalImageType::SizeType;
+  using InternalIndexType = typename InternalImageType::IndexType;
   using InternalImagePointerType = typename InternalImageType::Pointer;
   using InternalComplexType = std::complex<TInternalPrecision>;
   using InternalComplexImageType = Image<InternalComplexType, TInputImage::ImageDimension>;
@@ -116,11 +121,10 @@ class ITK_TEMPLATE_EXPORT FFTConvolutionImageFilter
   using FFTFilterType = RealToHalfHermitianForwardFFTImageFilter<InternalImageType, InternalComplexImageType>;
   using IFFTFilterType = HalfHermitianToRealInverseFFTImageFilter<InternalComplexImageType, InternalImageType>;
 
-  /** FFTConvolutionImageFilter needs the entire image kernel, which in
-   * general is going to be a different size than the output requested
-   * region. As such, this filter needs to provide an implementation
-   * for GenerateInputRequestedRegion() in order to inform the
-   * pipeline execution model.
+  /** Convolution uses a spatial region equivalent to the
+   *  output region padded by the kernel radius on all sides.
+   *  The input requested region is expanded by the kernel radius
+   *  within the bounds of the input largest possible region.
    *
    * \sa ProcessObject::GenerateInputRequestedRegion()  */
   void
@@ -181,18 +185,19 @@ class ITK_TEMPLATE_EXPORT FFTConvolutionImageFilter
   void
   CropOutput(InternalImageType * paddedOutput, ProgressAccumulator * progress, float progressWeight);
 
-  /** Get the lower bound for the padding of both the kernel and input
-   * images. Assuming that the regions of the kernel and input are the
-   * same, then this lower bound can be used to move the index of the
-   * padded kernel and padded input so that they are the same. This
-   * is important to avoid exceptions in filters that operate on these
-   * images. */
-  InputSizeType
-  GetPadLowerBound() const;
+  /** Get the radius of the kernel image. Used to pad the input image
+   *  for convolution. */
+  KernelSizeType
+  GetKernelRadius() const;
 
-  /** Get the pad size. */
-  InputSizeType
-  GetPadSize() const;
+  /** Get padding around the region of interest that results from FFT
+   *  factoring requirements. FFT typically requires that image side lengths
+   *  are factorable only by a fixed set of prime numbers (often 2, 3, and 5).
+   *  After the input image is padded for the kernel width and cropped to the
+   *  region of interest the result is then padded for FFT execution. This value
+   *  is reused for kernel padding and output cropping. */
+  InternalSizeType
+  GetFFTPadSize() const;
 
   /** Get whether the X dimension has an odd size. */
   bool
@@ -202,7 +207,9 @@ class ITK_TEMPLATE_EXPORT FFTConvolutionImageFilter
   PrintSelf(std::ostream & os, Indent indent) const override;
 
 private:
-  SizeValueType m_SizeGreatestPrimeFactor;
+  SizeValueType      m_SizeGreatestPrimeFactor;
+  InternalSizeType   m_FFTPadSize{ 0 };
+  InternalRegionType m_PaddedInputRegion;
 };
 } // namespace itk