ENH: Add FFTDiscreteGaussianImageFilter

Implements `FFTDiscreteGaussianImageFilter` as a subclass of
`DiscreteGaussianImageFilter` performing convolution of an ND
single-channel input image with a Gaussian kernel of fixed width in the
Fourier domain. Includes wrappings, test and baseline image.

Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.h

@@ -298,15 +298,6 @@ class ITK_TEMPLATE_EXPORT DiscreteGaussianImageFilter : public ImageToImageFilte
   itkLegacyMacro(unsigned int GetInternalNumberOfStreamDivisions() const);
   itkLegacyMacro(void SetInternalNumberOfStreamDivisions(unsigned int));
 
-  /** DiscreteGaussianImageFilter needs a larger input requested region
-   * than the output requested region (larger by the size of the
-   * Gaussian kernel).  As such, DiscreteGaussianImageFilter needs to
-   * provide an implementation for GenerateInputRequestedRegion() in
-   * order to inform the pipeline execution model.
-   * \sa ImageToImageFilter::GenerateInputRequestedRegion() */
-  void
-  GenerateInputRequestedRegion() override;
-
 #ifdef ITK_USE_CONCEPT_CHECKING
   // Begin concept checking
 
@@ -331,6 +322,15 @@ class ITK_TEMPLATE_EXPORT DiscreteGaussianImageFilter : public ImageToImageFilte
   void
   PrintSelf(std::ostream & os, Indent indent) const override;
 
+  /** DiscreteGaussianImageFilter needs a larger input requested region
+   * than the output requested region (larger by the size of the
+   * Gaussian kernel).  As such, DiscreteGaussianImageFilter needs to
+   * provide an implementation for GenerateInputRequestedRegion() in
+   * order to inform the pipeline execution model.
+   * \sa ImageToImageFilter::GenerateInputRequestedRegion() */
+  void
+  GenerateInputRequestedRegion() override;
+
   /** Standard pipeline method. While this class does not implement a
    * ThreadedGenerateData(), its GenerateData() delegates all
    * calculations to an NeighborhoodOperatorImageFilter.  Since the

Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.h

@@ -0,0 +1,134 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkFFTDiscreteGaussianImageFilter_h
+#define itkFFTDiscreteGaussianImageFilter_h
+
+#include "itkDiscreteGaussianImageFilter.h"
+#include "itkImage.h"
+#include "itkMacro.h"
+
+namespace itk
+{
+/**
+ * \class FFTDiscreteGaussianImageFilter
+ * \brief Blurs an image by convolution with a discrete gaussian kernel
+ * in the frequency domain.
+ *
+ * This filter performs Gaussian blurring by convolution of an image
+ * and a discrete Gaussian operator (kernel) in the frequency domain
+ * by way of Fast Fourier Transform forward and inverse operations.
+ *
+ * \sa DiscreteGaussianImageFilter
+ * \sa GaussianImageSource
+ * \sa FFTConvolutionImageFilter
+ * \sa RecursiveGaussianImageFilter
+ * \sa ZeroFluxNeumannBoundaryCondition
+ *
+ * \ingroup ImageEnhancement
+ * \ingroup ImageFeatureExtraction
+ * \ingroup ITKSmoothing
+ */
+
+template <typename TInputImage, typename TOutputImage = TInputImage>
+class ITK_TEMPLATE_EXPORT FFTDiscreteGaussianImageFilter : public DiscreteGaussianImageFilter<TInputImage, TOutputImage>
+{
+public:
+  ITK_DISALLOW_COPY_AND_MOVE(FFTDiscreteGaussianImageFilter);
+
+  /** Standard class type aliases. */
+  using Self = FFTDiscreteGaussianImageFilter;
+  using Superclass = DiscreteGaussianImageFilter<TInputImage, TOutputImage>;
+  using Pointer = SmartPointer<Self>;
+  using ConstPointer = SmartPointer<const Self>;
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(FFTDiscreteGaussianImageFilter, DiscreteGaussianImageFilter);
+
+  /** Image type information. */
+  using InputImageType = typename Superclass::InputImageType;
+  using OutputImageType = typename Superclass::OutputImageType;
+
+  /** Extract some information from the image types.  Dimensionality
+   * of the two images is assumed to be the same. */
+  using OutputPixelType = typename Superclass::OutputPixelType;
+  using OutputInternalPixelType = typename Superclass::OutputInternalPixelType;
+  using InputPixelType = typename Superclass::InputPixelType;
+  using InputInternalPixelType = typename Superclass::InputInternalPixelType;
+
+  /** Pixel value type for Vector pixel types **/
+  using InputPixelValueType = typename Superclass::InputPixelValueType;
+  using OutputPixelValueType = typename Superclass::OutputPixelValueType;
+
+  /** Extract some information from the image types.  Dimensionality
+   * of the two images is assumed to be the same. */
+  static constexpr unsigned int ImageDimension = Superclass::ImageDimension;
+
+  /** Type of the pixel to use for intermediate results */
+  using RealOutputPixelType = typename Superclass::RealOutputPixelType;
+  using RealOutputImageType = typename Superclass::RealOutputImageType;
+  using RealOutputPixelValueType = typename Superclass::RealOutputPixelValueType;
+  using RealPixelType = RealOutputPixelType;
+  using RealImageType = RealOutputImageType;
+
+  /** Typedef to describe the boundary condition. */
+  using BoundaryConditionType = typename Superclass::BoundaryConditionType;
+  using InputBoundaryConditionPointerType = typename Superclass::InputBoundaryConditionPointerType;
+  using InputDefaultBoundaryConditionType = typename Superclass::InputDefaultBoundaryConditionType;
+  using RealBoundaryConditionPointerType = typename Superclass::RealBoundaryConditionPointerType;
+  using RealDefaultBoundaryConditionType = typename Superclass::RealDefaultBoundaryConditionType;
+
+  /** Typedef of double containers */
+  using ArrayType = typename Superclass::ArrayType;
+  using SigmaArrayType = typename Superclass::SigmaArrayType;
+  using ScalarRealType = typename Superclass::ScalarRealType;
+
+  /** Typedef to describe kernel parameters */
+  using typename Superclass::KernelType;
+  using typename Superclass::RadiusType;
+
+  /** Overridden accessors for unused parameters */
+
+  void
+  SetInputBoundaryCondition(const InputBoundaryConditionPointerType) override;
+
+protected:
+  FFTDiscreteGaussianImageFilter() = default;
+  ~FFTDiscreteGaussianImageFilter() override = default;
+
+  // Pad input region to kernel image size
+  void
+  GenerateInputRequestedRegion() override;
+
+  /** Standard pipeline method. While this class does not implement a
+   * ThreadedGenerateData(), its GenerateData() delegates all
+   * calculations to an NeighborhoodOperatorImageFilter.  Since the
+   * NeighborhoodOperatorImageFilter is multithreaded, this filter is
+   * multithreaded by default. */
+  void
+  GenerateData() override;
+};
+} // end namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#  include "itkFFTDiscreteGaussianImageFilter.hxx"
+#endif
+
+#endif

Modules/Filtering/Smoothing/include/itkFFTDiscreteGaussianImageFilter.hxx

@@ -0,0 +1,173 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkFFTDiscreteGaussianImageFilter_hxx
+#define itkFFTDiscreteGaussianImageFilter_hxx
+
+#include "itkGaussianOperator.h"
+#include "itkGaussianImageSource.h"
+#include "itkFFTConvolutionImageFilter.h"
+
+#include "itkProgressAccumulator.h"
+#include "itkImageAlgorithm.h"
+#include "itkMacro.h"
+
+namespace itk
+{
+
+template <typename TInputImage, typename TOutputImage>
+void
+FFTDiscreteGaussianImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion()
+{
+  // call the superclass' implementation of this method. this should
+  // copy the output requested region to the input requested region
+  ImageToImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion();
+
+  // Get pointer to input
+  typename Superclass::InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
+  if (inputPtr.IsNull())
+  {
+    return;
+  }
+
+  // get a copy of the input requested region (should equal the output
+  // requested region)
+  typename TInputImage::RegionType inputRequestedRegion;
+  inputRequestedRegion = inputPtr->GetRequestedRegion();
+
+  // pad the input requested region by the operator radius
+  RadiusType radius;
+  radius.Fill(0);
+  for (size_t dim = 0; dim < ImageDimension; ++dim)
+  {
+    radius[dim] = this->GetKernelRadius(dim);
+  }
+  inputRequestedRegion.PadByRadius(radius);
+
+  // crop the input requested region at the input's largest possible region
+  inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion());
+
+  inputPtr->SetRequestedRegion(inputRequestedRegion);
+}
+
+template <typename TInputImage, typename TOutputImage>
+void
+FFTDiscreteGaussianImageFilter<TInputImage, TOutputImage>::SetInputBoundaryCondition(
+  const InputBoundaryConditionPointerType)
+{
+  itkWarningMacro("FFTDiscreteGaussianImageFilter ignores InputBoundaryCondition, use RealBoundaryCondition instead");
+}
+
+template <typename TInputImage, typename TOutputImage>
+void
+FFTDiscreteGaussianImageFilter<TInputImage, TOutputImage>::GenerateData()
+{
+  TOutputImage * output = this->GetOutput();
+
+  output->SetBufferedRegion(output->GetRequestedRegion());
+  output->Allocate();
+
+  // Create an internal image to protect the input image's metdata
+  // (e.g. RequestedRegion). The StreamingImageFilter changes the
+  // requested region as part of its normal processing.
+  auto localInput = TInputImage::New();
+  localInput->Graft(this->GetInput());
+
+  // Determine the dimensionality to filter
+  unsigned int filterDimensionality = this->GetFilterDimensionality();
+  if (filterDimensionality > ImageDimension)
+  {
+    filterDimensionality = ImageDimension;
+  }
+  if (filterDimensionality == 0)
+  {
+    // no smoothing, copy input to output
+    ImageAlgorithm::Copy(localInput.GetPointer(),
+                         output,
+                         this->GetOutput()->GetRequestedRegion(),
+                         this->GetOutput()->GetRequestedRegion());
+    return;
+  }
+
+  // Create a process accumulator for tracking the progress of minipipeline
+  auto progress = ProgressAccumulator::New();
+  progress->SetMiniPipelineFilter(this);
+
+  // Create kernel image for blurring in requested dimensions
+  using GaussianImageSourceType = GaussianImageSource<RealImageType>;
+  using KernelSizeType = typename GaussianImageSourceType::SizeType;
+  using KernelMeanType = typename GaussianImageSourceType::ArrayType;
+
+  typename GaussianImageSourceType::Pointer kernelSource = GaussianImageSourceType::New();
+
+  kernelSource->SetScale(1.0);
+  kernelSource->SetNormalized(true);
+
+  kernelSource->SetSpacing(localInput->GetSpacing());
+  kernelSource->SetOrigin(localInput->GetOrigin());
+  kernelSource->SetDirection(localInput->GetDirection());
+
+  KernelSizeType kernelSize;
+  kernelSize.Fill(1);
+  for (size_t dim = 0; dim < filterDimensionality; ++dim)
+  {
+    kernelSize[dim] = static_cast<SizeValueType>(this->GetKernelRadius(dim)) * 2 + 1;
+  }
+  kernelSource->SetSize(kernelSize);
+
+  KernelMeanType mean;
+  auto           inputSpacing = localInput->GetSpacing();
+  auto           inputOrigin = localInput->GetOrigin();
+  for (size_t dim = 0; dim < ImageDimension; ++dim)
+  {
+    double radius = (kernelSize[dim] - 1) / 2;
+    mean[dim] = inputSpacing[dim] * radius + inputOrigin[dim]; // center pixel pos
+  }
+  kernelSource->SetMean(mean);
+  kernelSource->SetSigma(this->GetSigmaArray());
+
+  // Estimate kernel generation to be roughly 5% of effort
+  progress->RegisterInternalFilter(kernelSource, 0.05f);
+
+  // Perform image convolution by FFT
+  using ConvolutionImageFilterType = FFTConvolutionImageFilter<RealImageType, RealImageType, OutputImageType>;
+
+  typename ConvolutionImageFilterType::Pointer convolutionFilter = ConvolutionImageFilterType::New();
+  convolutionFilter->SetInput(this->GetInput());
+  convolutionFilter->SetKernelImage(kernelSource->GetOutput());
+  convolutionFilter->SetBoundaryCondition(this->GetRealBoundaryCondition());
+  convolutionFilter->SetNormalize(false); // Kernel is already normalized
+
+  // Estimate kernel generation to be roughly 95% of effort
+  progress->RegisterInternalFilter(convolutionFilter, 0.95f);
+
+  // Graft this filters output onto the mini-pipeline so the mini-pipeline
+  // has the correct region ivars and will write to this filters bulk data
+  // output.
+  convolutionFilter->GraftOutput(output);
+
+  // Update the last filter in the mini-pipeline
+  convolutionFilter->Update();
+
+  // Graft the last output of the mini-pipeline onto this filters output so
+  // the final output has the correct region ivars and a handle to the final
+  // bulk data
+  this->GraftOutput(output);
+}
+} // namespace itk
+
+#endif

Modules/Filtering/Smoothing/itk-module.cmake

@@ -5,8 +5,12 @@ interesting to look at the ITKAnisotropicSmoothing group of filters.")
 itk_module(ITKSmoothing
   ENABLE_SHARED
   COMPILE_DEPENDS
+    ITKConvolution
+    ITKFFT
     ITKImageFunction
+    ITKImageSources
   TEST_DEPENDS
+    ITKConvolution
     ITKTestKernel
   DESCRIPTION
     "${DOCUMENTATION}"

Modules/Filtering/Smoothing/test/Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha.sha512

@@ -0,0 +1 @@
+0eeb11ed3467d7af854c70506ea3e476fb73f57aaa1779ba3eb8abbfc30035bf4b0c11a7e47d12c65f01f1682dcac11a01e8ac20f1a0398ddb5f47f668a24a75

Modules/Filtering/Smoothing/test/CMakeLists.txt

@@ -3,6 +3,7 @@ set(ITKSmoothingTests
 itkBoxMeanImageFilterTest.cxx
 itkBoxSigmaImageFilterTest.cxx
 itkDiscreteGaussianImageFilterTest2.cxx
+itkFFTDiscreteGaussianImageFilterTest.cxx
 itkSmoothingRecursiveGaussianImageFilterTest.cxx
 itkSmoothingRecursiveGaussianImageFilterOnVectorImageTest.cxx
 itkSmoothingRecursiveGaussianImageFilterOnImageOfVectorTest.cxx
@@ -59,6 +60,15 @@ itk_add_test(NAME itkDiscreteGaussianImageFilterTest1a
       COMMAND ITKSmoothingTestDriver itkDiscreteGaussianImageFilterTest 1)
 itk_add_test(NAME itkDiscreteGaussianImageFilterTest1b
       COMMAND ITKSmoothingTestDriver itkDiscreteGaussianImageFilterTest 0)
+itk_add_test(NAME itkFFTDiscreteGaussianImageFilterTest
+      COMMAND ITKSmoothingTestDriver
+      --compare DATA{Baseline/FFTDiscreteGaussianImageFilterTestOutput.mha}
+                ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestOutput.mha
+        itkFFTDiscreteGaussianImageFilterTest
+        DATA{${ITK_DATA_ROOT}/Input/2th_cthead1.png}
+        3.5
+        35
+        ${ITK_TEST_OUTPUT_DIR}/FFTDiscreteGaussianImageFilterTestOutput.mha)
 itk_add_test(NAME itkMedianImageFilterTest
       COMMAND ITKSmoothingTestDriver itkMedianImageFilterTest)
 itk_add_test(NAME itkRecursiveGaussianImageFiltersOnTensorsTest