ENH: Expose kernel parameters in DiscreteGaussianImageFilter interface

`itk::DiscreteGaussianImageFilter` exposes parameters for generating a
Gaussian kernel to error and size specifications, but previously did
not make information on the kernel that was actually generated available
to other classes.

These changes expose protected methods to allow subclasses to create
kernels to the same specifications and also publicly exposes information
on the kernel sigma when image spacing is being accounted for.

Author: tbirdso

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
@@ -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

Modules/Filtering/Smoothing/test/itkDiscreteGaussianImageFilterTest.cxx

@@ -71,6 +71,10 @@ itkDiscreteGaussianImageFilterTest(int argc, char * argv[])
   filter->SetSigmaArray(sigma);
   ITK_TEST_SET_GET_VALUE(sigma, filter->GetSigmaArray());
 
+  // Verify spacing-adjusted variance is not available without input image
+  filter->SetUseImageSpacing(true);
+  ITK_TRY_EXPECT_EXCEPTION(filter->GetImageSpacingVarianceArray());
+
   // Set the boundary condition used by the filter
   itk::ConstantBoundaryCondition<ImageType> constantBoundaryCondition;
   filter->SetInputBoundaryCondition(&constantBoundaryCondition);
@@ -92,7 +96,7 @@ itkDiscreteGaussianImageFilterTest(int argc, char * argv[])
   filter->SetMaximumError(maximumError);
   ITK_TEST_SET_GET_VALUE(maximumError, filter->GetMaximumError());
 
-  int maximumKernelWidth = 32;
+  unsigned int maximumKernelWidth = 32;
   filter->SetMaximumKernelWidth(maximumKernelWidth);
   ITK_TEST_SET_GET_VALUE(maximumKernelWidth, filter->GetMaximumKernelWidth());
 
@@ -124,6 +128,28 @@ itkDiscreteGaussianImageFilterTest(int argc, char * argv[])
 
   ITK_TRY_EXPECT_NO_EXCEPTION(test1.Execute());
 
+  // Test variance adjustments for image spacing
+  ImageType::Pointer              inputImage = ImageType::New();
+  typename ImageType::SpacingType spacing;
+  spacing[0] = 1.5;
+  spacing[1] = 1.0;
+  spacing[2] = 0.5;
+  inputImage->SetSpacing(spacing);
+  filter->SetInput(inputImage);
+
+  filter->SetUseImageSpacing(false);
+  for (itk::IndexValueType dim = 0; dim < Dimension; ++dim)
+  {
+    ITK_TEST_EXPECT_EQUAL(sigmaValue, filter->GetVariance()[dim])
+    ITK_TEST_EXPECT_EQUAL(sigmaValue, filter->GetImageSpacingVarianceArray()[dim]);
+  }
+
+  filter->SetUseImageSpacing(true);
+  for (itk::IndexValueType dim = 0; dim < Dimension; ++dim)
+  {
+    ITK_TEST_EXPECT_EQUAL(sigmaValue, filter->GetVariance()[dim]);
+    ITK_TEST_EXPECT_EQUAL((sigmaValue / (spacing[dim] * spacing[dim])), filter->GetImageSpacingVarianceArray()[dim]);
+  }
 
   std::cout << "Test finished" << std::endl;
   return EXIT_SUCCESS;