ENH: Add FFTDiscreteGaussianImageFilter

Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.h

@@ -18,6 +18,7 @@
 #ifndef itkDiscreteGaussianImageFilter_h
 #define itkDiscreteGaussianImageFilter_h
 
+#include "itkGaussianOperator.h"
 #include "itkImageToImageFilter.h"
 #include "itkImage.h"
 #include "itkZeroFluxNeumannBoundaryCondition.h"
@@ -113,6 +114,10 @@ class ITK_TEMPLATE_EXPORT DiscreteGaussianImageFilter : public ImageToImageFilte
   using SigmaArrayType = ArrayType;
   using ScalarRealType = double;
 
+  /** Type of kernel to be used in blurring */
+  using KernelType = GaussianOperator<RealOutputPixelValueType, ImageDimension>;
+  using RadiusType = typename KernelType::RadiusType;
+
   /** The variance for the discrete Gaussian kernel.  Sets the variance
    * independently for each dimension, but
    * see also SetVariance(const double v). The default is 0.0 in each
@@ -130,8 +135,8 @@ class ITK_TEMPLATE_EXPORT DiscreteGaussianImageFilter : public ImageToImageFilte
 
   /** Set the kernel to be no wider than MaximumKernelWidth pixels,
    *  even if MaximumError demands it. The default is 32 pixels. */
-  itkGetConstMacro(MaximumKernelWidth, int);
-  itkSetMacro(MaximumKernelWidth, int);
+  itkGetConstMacro(MaximumKernelWidth, unsigned int);
+  itkSetMacro(MaximumKernelWidth, unsigned int);
 
   /** Set the number of dimensions to smooth. Defaults to the image
    * dimension. Can be set to less than ImageDimension, smoothing all
@@ -293,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
 
@@ -326,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
@@ -334,6 +339,18 @@ class ITK_TEMPLATE_EXPORT DiscreteGaussianImageFilter : public ImageToImageFilte
   void
   GenerateData() override;
 
+  /** Build a directional kernel to match user specifications */
+  void
+  GenerateKernel(const unsigned int dimension, KernelType & oper) const;
+
+  /** Get the radius of the generated directional kernel */
+  unsigned int
+  GetKernelRadius(const unsigned int dimension) const;
+
+  /** Get the variance, optionally adjusted for pixel spacing */
+  ArrayType
+  GetKernelVarianceArray() const;
+
 private:
   /** The variance of the gaussian blurring kernel in each dimensional
     direction. */
@@ -346,7 +363,7 @@ class ITK_TEMPLATE_EXPORT DiscreteGaussianImageFilter : public ImageToImageFilte
 
   /** Maximum allowed kernel width for any dimension of the discrete Gaussian
       approximation */
-  int m_MaximumKernelWidth;
+  unsigned int m_MaximumKernelWidth;
 
   /** Number of dimensions to process. Default is all dimensions */
   unsigned int m_FilterDimensionality;

Modules/Filtering/Smoothing/include/itkDiscreteGaussianImageFilter.hxx

@@ -26,6 +26,59 @@
 
 namespace itk
 {
+template <typename TInputImage, typename TOutputImage>
+void
+DiscreteGaussianImageFilter<TInputImage, TOutputImage>::GenerateKernel(const unsigned int dimension,
+                                                                       KernelType &       oper) const
+{
+  // Determine the size of the operator in this dimension.  Note that the
+  // Gaussian is built as a 1D operator in each of the specified directions.
+  oper.SetDirection(dimension);
+  oper.SetMaximumError(m_MaximumError[dimension]);
+  oper.SetMaximumKernelWidth(m_MaximumKernelWidth);
+
+  oper.SetVariance(this->GetKernelVarianceArray()[dimension]);
+
+  oper.CreateDirectional();
+}
+
+template <typename TInputImage, typename TOutputImage>
+unsigned int
+DiscreteGaussianImageFilter<TInputImage, TOutputImage>::GetKernelRadius(const unsigned int dimension) const
+{
+  KernelType oper;
+  this->GenerateKernel(dimension, oper);
+  return oper.GetRadius(dimension);
+}
+
+template <typename TInputImage, typename TOutputImage>
+typename DiscreteGaussianImageFilter<TInputImage, TOutputImage>::ArrayType
+DiscreteGaussianImageFilter<TInputImage, TOutputImage>::GetKernelVarianceArray() const
+{
+  if (m_UseImageSpacing)
+  {
+    if (this->GetInput() == nullptr)
+    {
+      itkExceptionMacro("UseImageSpacing is ON but no input image was provided");
+    }
+
+    ArrayType adjustedVariance;
+    // Adjusted variance = var / (spacing ^ 2)
+    for (unsigned int dim = 0; dim < ImageDimension; ++dim)
+    {
+      // convert the variance from physical units to pixels
+      double s = this->GetInput()->GetSpacing()[dim];
+      s = s * s;
+      adjustedVariance[dim] = m_Variance[dim] / s;
+    }
+    return adjustedVariance;
+  }
+  else
+  {
+    return this->GetVariance();
+  }
+}
+
 template <typename TInputImage, typename TOutputImage>
 void
 DiscreteGaussianImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion()
@@ -42,39 +95,18 @@ DiscreteGaussianImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRe
     return;
   }
 
-  // Build an operator so that we can determine the kernel size
-  GaussianOperator<OutputPixelValueType, ImageDimension> oper;
-
-  typename TInputImage::SizeType radius;
-
+  // Determine the kernel size in each direction
+  RadiusType radius;
   for (unsigned int i = 0; i < TInputImage::ImageDimension; ++i)
   {
-    // Determine the size of the operator in this dimension.  Note that the
-    // Gaussian is built as a 1D operator in each of the specified directions.
-    oper.SetDirection(i);
-    if (m_UseImageSpacing == true)
+    if (i < m_FilterDimensionality)
     {
-      if (this->GetInput()->GetSpacing()[i] == 0.0)
-      {
-        itkExceptionMacro(<< "Pixel spacing cannot be zero");
-      }
-      else
-      {
-        // convert the variance from physical units to pixels
-        double s = this->GetInput()->GetSpacing()[i];
-        s = s * s;
-        oper.SetVariance(m_Variance[i] / s);
-      }
+      radius[i] = GetKernelRadius(i);
     }
     else
     {
-      oper.SetVariance(m_Variance[i]);
+      radius[i] = 0;
     }
-    oper.SetMaximumError(m_MaximumError[i]);
-    oper.SetMaximumKernelWidth(m_MaximumKernelWidth);
-    oper.CreateDirectional();
-
-    radius[i] = oper.GetRadius(i);
   }
 
   // get a copy of the input requested region (should equal the output
@@ -86,26 +118,9 @@ DiscreteGaussianImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRe
   inputRequestedRegion.PadByRadius(radius);
 
   // crop the input requested region at the input's largest possible region
-  if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()))
-  {
-    inputPtr->SetRequestedRegion(inputRequestedRegion);
-    return;
-  }
-  else
-  {
-    // Couldn't crop the region (requested region is outside the largest
-    // possible region).  Throw an exception.
-
-    // store what we tried to request (prior to trying to crop)
-    inputPtr->SetRequestedRegion(inputRequestedRegion);
-
-    // build an exception
-    InvalidRequestedRegionError e(__FILE__, __LINE__);
-    e.SetLocation(ITK_LOCATION);
-    e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
-    e.SetDataObject(inputPtr);
-    throw e;
-  }
+  inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion());
+
+  inputPtr->SetRequestedRegion(inputRequestedRegion);
 }
 
 template <typename TInputImage, typename TOutputImage>
@@ -160,9 +175,8 @@ DiscreteGaussianImageFilter<TInputImage, TOutputImage>::GenerateData()
   using SingleFilterPointer = typename SingleFilterType::Pointer;
 
   // Create a series of operators
-  using OperatorType = GaussianOperator<RealOutputPixelValueType, ImageDimension>;
 
-  std::vector<OperatorType> oper;
+  std::vector<KernelType> oper;
   oper.resize(filterDimensionality);
 
   // Create a process accumulator for tracking the progress of minipipeline
@@ -177,30 +191,7 @@ DiscreteGaussianImageFilter<TInputImage, TOutputImage>::GenerateData()
     // the largest dimension will be split slice wise for streaming
     unsigned int reverse_i = filterDimensionality - i - 1;
 
-    // Set up the operator for this dimension
-    oper[reverse_i].SetDirection(i);
-    if (m_UseImageSpacing == true)
-    {
-      if (localInput->GetSpacing()[i] == 0.0)
-      {
-        itkExceptionMacro(<< "Pixel spacing cannot be zero");
-      }
-      else
-      {
-        // convert the variance from physical units to pixels
-        double s = localInput->GetSpacing()[i];
-        s = s * s;
-        oper[reverse_i].SetVariance(m_Variance[i] / s);
-      }
-    }
-    else
-    {
-      oper[reverse_i].SetVariance(m_Variance[i]);
-    }
-
-    oper[reverse_i].SetMaximumKernelWidth(m_MaximumKernelWidth);
-    oper[reverse_i].SetMaximumError(m_MaximumError[i]);
-    oper[reverse_i].CreateDirectional();
+    this->GenerateKernel(i, oper[reverse_i]);
   }
 
   // Create a chain of filters