sdiebolt · Oct 23, 2026
diff --git a/Modules/Filtering/Smoothing/include/itkSmoothingRecursiveGaussianImageFilter.h b/Modules/Filtering/Smoothing/include/itkSmoothingRecursiveGaussianImageFilter.h
@@ -37,6 +37,10 @@ namespace itk
  * filters. For multi-component images, the filter works on each
  * component independently.
  *
+ * Axes with fewer than four pixels are not smoothed. If all axes are shorter
+ * than four pixels, the input is only cast to the output pixel type.
+ * Image dimensionality and physical geometry are preserved.
+ *
  * For this filter to be able to run in-place the input and output
  * image types need to be the same and/or the same type as the
  * RealImageType.

diff --git a/Modules/Filtering/Smoothing/include/itkSmoothingRecursiveGaussianImageFilter.hxx b/Modules/Filtering/Smoothing/include/itkSmoothingRecursiveGaussianImageFilter.hxx
@@ -27,12 +27,7 @@ template <typename TInputImage, typename TOutputImage>
 SmoothingRecursiveGaussianImageFilter<TInputImage, TOutputImage>::SmoothingRecursiveGaussianImageFilter()
   : m_FirstSmoothingFilter(FirstGaussianFilterType::New())
 {
-  // NB: The first filter is the last dimension because it does not
-  // always run in-place. As this dimension provides the least amount
-  // of cache coherency, it will provide the least  amount of performance
-  // gains from running in-place. In fact, some performance tests
-  // indicate that the running in-place in the 3rd dimension the
-  // performance actually declines compared to not in-place methods.
+  // The highest eligible dimension benefits least from in-place execution, so process it first.
   m_FirstSmoothingFilter->SetOrder(GaussianOrderEnum::ZeroOrder);
   m_FirstSmoothingFilter->SetDirection(ImageDimension - 1);
   m_FirstSmoothingFilter->SetNormalizeAcrossScale(m_NormalizeAcrossScale);
@@ -49,14 +44,7 @@ SmoothingRecursiveGaussianImageFilter<TInputImage, TOutputImage>::SmoothingRecur
     m_SmoothingFilters[i]->InPlaceOn();
   }
 
-  m_SmoothingFilters[0]->SetInput(m_FirstSmoothingFilter->GetOutput());
-  for (unsigned int i = 1; i < ImageDimension - 1; ++i)
-  {
-    m_SmoothingFilters[i]->SetInput(m_SmoothingFilters[i - 1]->GetOutput());
-  }
-
   m_CastingFilter = CastingFilterType::New();
-  m_CastingFilter->SetInput(m_SmoothingFilters.back()->GetOutput());
   m_CastingFilter->InPlaceOn();
 
   this->InPlaceOff();
@@ -202,19 +190,68 @@ SmoothingRecursiveGaussianImageFilter<TInputImage, TOutputImage>::GenerateData()
 
   const typename TInputImage::ConstPointer inputImage(this->GetInput());
 
-  const typename TInputImage::RegionType region = inputImage->GetRequestedRegion();
-  const typename TInputImage::SizeType   size = region.GetSize();
-
+  const auto   size = inputImage->GetLargestPossibleRegion().GetSize();
+  unsigned int numberOfActiveDimensions = 0;
+  unsigned int firstDirection = 0;
   for (unsigned int d = 0; d < ImageDimension; ++d)
   {
-    if (size[d] < 4)
+    if (m_Sigma[d] <= 0.0)
     {
-      itkExceptionMacro(
-        "The number of pixels along dimension "
-        << d << " is less than 4. This filter requires a minimum of four pixels along the dimension to be processed.");
+      itkExceptionStringMacro("Sigma must be greater than zero.");
+    }
+    if (size[d] == 0)
+    {
+      itkExceptionMacro("The number of pixels along dimension " << d << " is zero.");
+    }
+    if (size[d] >= 4)
+    {
+      firstDirection = d;
+      ++numberOfActiveDimensions;
     }
   }
 
+  auto progress = ProgressAccumulator::New();
+  progress->SetMiniPipelineFilter(this);
+
+  m_FirstSmoothingFilter->SetInput(inputImage);
+
+  if (numberOfActiveDimensions == 0)
+  {
+    m_CastingFilter->SetInput(nullptr);
+    m_CastingFilter->GetOutput()->ReleaseData();
+    auto cast = CastImageFilter<InputImageType, OutputImageType>::New();
+    cast->SetInput(inputImage);
+    cast->SetInPlace(this->GetInPlace());
+    cast->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
+    progress->RegisterInternalFilter(cast, 1.0);
+    cast->GraftOutput(this->GetOutput());
+    cast->Update();
+    this->GraftOutput(cast->GetOutput());
+    return;
+  }
+
+  m_FirstSmoothingFilter->SetDirection(firstDirection);
+  m_FirstSmoothingFilter->SetSigma(m_Sigma[firstDirection]);
+  progress->RegisterInternalFilter(m_FirstSmoothingFilter, 1.0 / numberOfActiveDimensions);
+
+  typename RealImageType::Pointer previousOutput = m_FirstSmoothingFilter->GetOutput();
+  for (unsigned int d = 0; d < ImageDimension - 1; ++d)
+  {
+    if (d != firstDirection && size[d] >= 4)
+    {
+      auto & filter = m_SmoothingFilters[d];
+      filter->SetSigma(m_Sigma[d]);
+      filter->SetInput(previousOutput);
+      progress->RegisterInternalFilter(filter, 1.0 / numberOfActiveDimensions);
+      previousOutput = filter->GetOutput();
+    }
+    else
+    {
+      m_SmoothingFilters[d]->SetInput(nullptr);
+    }
+  }
+  m_CastingFilter->SetInput(previousOutput);
+
   // If this filter is running in-place, then set the first smoothing
   // filter to steal the bulk data, by running in-place.
   if (this->CanRunInPlace() && this->GetInPlace())
@@ -238,20 +275,6 @@ SmoothingRecursiveGaussianImageFilter<TInputImage, TOutputImage>::GenerateData()
     this->GetOutput()->ReleaseData();
   }
 
-  // Create a process accumulator for tracking the progress of this minipipeline.
-  auto progress = ProgressAccumulator::New();
-  progress->SetMiniPipelineFilter(this);
-
-  // Register the filter with the with progress accumulator using
-  // equal weight proportion.
-  for (const auto & filter : m_SmoothingFilters)
-  {
-    progress->RegisterInternalFilter(filter, 1.0 / (ImageDimension));
-  }
-
-  progress->RegisterInternalFilter(m_FirstSmoothingFilter, 1.0 / (ImageDimension));
-  m_FirstSmoothingFilter->SetInput(inputImage);
-
   // Graft our output to the internal filter to force the proper regions
   // to be generated, and the bulk data which be be from the input due
   // to the in-place option.

diff --git a/Modules/Filtering/Smoothing/test/CMakeLists.txt b/Modules/Filtering/Smoothing/test/CMakeLists.txt
@@ -399,5 +399,6 @@ set(
   itkBoxSigmaImageFilterGTest.cxx
   itkMeanImageFilterGTest.cxx
   itkMedianImageFilterGTest.cxx
+  itkSmoothingRecursiveGaussianImageFilterGTest.cxx
 )
 creategoogletestdriver(ITKSmoothing "${ITKSmoothing-Test_LIBRARIES}" "${ITKSmoothingGTests}")
diff --git a/Modules/Filtering/Smoothing/test/itkSmoothingRecursiveGaussianImageFilterGTest.cxx b/Modules/Filtering/Smoothing/test/itkSmoothingRecursiveGaussianImageFilterGTest.cxx
@@ -0,0 +1,247 @@
+/*=========================================================================
+ *
+ *  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
+ *
+ *         https://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.
+ *
+ *=========================================================================*/
+
+#include "itkSmoothingRecursiveGaussianImageFilter.h"
+#include "itkImageRegionIteratorWithIndex.h"
+#include "itkImageRegionIterator.h"
+#include "itkImageRegionConstIterator.h"
+#include "itkVectorImage.h"
+
+#include <gtest/gtest.h>
+
+namespace
+{
+using ImageType = itk::Image<float, 3>;
+using FilterType = itk::SmoothingRecursiveGaussianImageFilter<ImageType>;
+
+ImageType::Pointer
+MakeImage(const ImageType::SizeType & size)
+{
+  auto image = ImageType::New();
+  image->SetRegions(ImageType::RegionType({ { -2, 3, 5 } }, size));
+  image->SetOrigin(ImageType::PointType{ { 7.0, -11.0, 2.0 } });
+  image->SetSpacing(ImageType::SpacingType{ { 0.7, 1.3, 2.1 } });
+  ImageType::DirectionType direction;
+  direction.SetIdentity();
+  direction[0][0] = 0;
+  direction[0][1] = -1;
+  direction[1][0] = 1;
+  direction[1][1] = 0;
+  image->SetDirection(direction);
+  image->Allocate();
+  for (itk::ImageRegionIteratorWithIndex<ImageType> it(image, image->GetLargestPossibleRegion()); !it.IsAtEnd(); ++it)
+  {
+    const auto index = it.GetIndex();
+    it.Set(static_cast<float>((index[0] * index[0] + 3 * index[1] + index[2] * index[2]) % 17));
+  }
+  return image;
+}
+
+ImageType::Pointer
+SmoothEligibleAxes(const ImageType * input, const FilterType::SigmaArrayType & sigma)
+{
+  ImageType::ConstPointer previous = input;
+  ImageType::Pointer      output;
+  const auto              size = input->GetLargestPossibleRegion().GetSize();
+  unsigned int            firstDirection = 0;
+  bool                    hasActiveDirection = false;
+  for (unsigned int d = 0; d < 3; ++d)
+  {
+    if (size[d] >= 4)
+    {
+      firstDirection = d;
+      hasActiveDirection = true;
+    }
+  }
+  if (!hasActiveDirection)
+  {
+    auto copy = itk::CastImageFilter<ImageType, ImageType>::New();
+    copy->SetInput(input);
+    copy->Update();
+    output = copy->GetOutput();
+    output->DisconnectPipeline();
+    return output;
+  }
+  for (const unsigned int d : { firstDirection, 0u, 1u })
+  {
+    if (size[d] < 4 || (d == firstDirection && output.IsNotNull()))
+    {
+      continue;
+    }
+    auto recursive = itk::RecursiveGaussianImageFilter<ImageType>::New();
+    recursive->SetInput(previous);
+    recursive->SetDirection(d);
+    recursive->SetSigma(sigma[d]);
+    recursive->Update();
+    output = recursive->GetOutput();
+    output->DisconnectPipeline();
+    previous = output;
+  }
+  return output;
+}
+
+void
+ExpectSameImage(const ImageType * actual, const ImageType * expected)
+{
+  ASSERT_EQ(actual->GetLargestPossibleRegion(), expected->GetLargestPossibleRegion());
+  EXPECT_EQ(actual->GetOrigin(), expected->GetOrigin());
+  EXPECT_EQ(actual->GetSpacing(), expected->GetSpacing());
+  EXPECT_EQ(actual->GetDirection(), expected->GetDirection());
+  itk::ImageRegionConstIterator<ImageType> actualIt(actual, actual->GetLargestPossibleRegion());
+  itk::ImageRegionConstIterator<ImageType> expectedIt(expected, expected->GetLargestPossibleRegion());
+  for (; !actualIt.IsAtEnd(); ++actualIt, ++expectedIt)
+  {
+    EXPECT_FLOAT_EQ(actualIt.Get(), expectedIt.Get());
+  }
+}
+} // namespace
+
+TEST(SmoothingRecursiveGaussianImageFilter, ShortAxesAndRepeatedUpdates)
+{
+  const FilterType::SigmaArrayType sigma{ { 0.8, 1.7, 2.5 } };
+  for (const bool inPlace : { false, true })
+  {
+    auto filter = FilterType::New();
+    filter->SetSigmaArray(sigma);
+    filter->SetInPlace(inPlace);
+    for (const unsigned int x : { 1u, 2u, 3u, 4u, 9u })
+    {
+      for (const unsigned int y : { 1u, 2u, 3u, 4u, 9u })
+      {
+        for (const unsigned int z : { 1u, 2u, 3u, 4u, 9u })
+        {
+          const ImageType::SizeType size{ { x, y, z } };
+          SCOPED_TRACE(::testing::Message() << "size=" << size << ", inPlace=" << inPlace);
+          const auto image = MakeImage(size);
+          const auto expected = SmoothEligibleAxes(image, sigma);
+          filter->SetInput(image);
+          ASSERT_NO_THROW(filter->Update());
+          ExpectSameImage(filter->GetOutput(), expected);
+        }
+      }
+    }
+  }
+}
+
+TEST(SmoothingRecursiveGaussianImageFilter, AllShortAxesPreserveDoublePrecisionAndCast)
+{
+  using DoubleImageType = itk::Image<double, 3>;
+  auto image = DoubleImageType::New();
+  image->SetRegions(DoubleImageType::SizeType{ { 1, 2, 3 } });
+  image->Allocate();
+  constexpr double value = 1.0000000000000002;
+  image->FillBuffer(value);
+  for (const bool inPlace : { false, true })
+  {
+    auto filter = itk::SmoothingRecursiveGaussianImageFilter<DoubleImageType>::New();
+    filter->SetInput(image);
+    filter->SetInPlace(inPlace);
+    ASSERT_NO_THROW(filter->Update());
+    EXPECT_EQ(filter->GetOutput()->GetPixel({ { 0, 0, 0 } }), value);
+  }
+
+  using OutputImageType = itk::Image<unsigned char, 3>;
+  const auto input = MakeImage({ { 2, 3, 1 } });
+  auto       filter = itk::SmoothingRecursiveGaussianImageFilter<ImageType, OutputImageType>::New();
+  filter->SetInput(input);
+  ASSERT_NO_THROW(filter->Update());
+  itk::ImageRegionConstIterator<ImageType>       inputIt(input, input->GetLargestPossibleRegion());
+  itk::ImageRegionConstIterator<OutputImageType> outputIt(filter->GetOutput(), input->GetLargestPossibleRegion());
+  for (; !inputIt.IsAtEnd(); ++inputIt, ++outputIt)
+  {
+    EXPECT_EQ(outputIt.Get(), static_cast<unsigned char>(inputIt.Get()));
+  }
+}
+
+TEST(SmoothingRecursiveGaussianImageFilter, VectorImageShortAxes)
+{
+  using VectorImageType = itk::VectorImage<float, 3>;
+  for (const ImageType::SizeType size : { ImageType::SizeType{ { 9, 2, 1 } }, ImageType::SizeType{ { 1, 2, 3 } } })
+  {
+    const auto scalarImage = MakeImage(size);
+    auto       vectorImage = VectorImageType::New();
+    vectorImage->CopyInformation(scalarImage);
+    vectorImage->SetRegions(scalarImage->GetLargestPossibleRegion());
+    vectorImage->SetNumberOfComponentsPerPixel(2);
+    vectorImage->Allocate();
+    itk::ImageRegionConstIterator<ImageType>  scalarIt(scalarImage, scalarImage->GetLargestPossibleRegion());
+    itk::ImageRegionIterator<VectorImageType> vectorIt(vectorImage, vectorImage->GetLargestPossibleRegion());
+    for (; !scalarIt.IsAtEnd(); ++scalarIt, ++vectorIt)
+    {
+      VectorImageType::PixelType pixel(2);
+      pixel[0] = scalarIt.Get();
+      pixel[1] = 2 * scalarIt.Get();
+      vectorIt.Set(pixel);
+    }
+    auto filter = itk::SmoothingRecursiveGaussianImageFilter<VectorImageType>::New();
+    filter->SetInput(vectorImage);
+    ASSERT_NO_THROW(filter->Update());
+    const auto expected = SmoothEligibleAxes(scalarImage, FilterType::SigmaArrayType(1.0));
+    itk::ImageRegionConstIterator<ImageType>       expectedIt(expected, expected->GetLargestPossibleRegion());
+    itk::ImageRegionConstIterator<VectorImageType> outputIt(filter->GetOutput(), expected->GetLargestPossibleRegion());
+    for (; !expectedIt.IsAtEnd(); ++expectedIt, ++outputIt)
+    {
+      ASSERT_EQ(outputIt.Get().GetSize(), 2u);
+      EXPECT_FLOAT_EQ(outputIt.Get()[0], expectedIt.Get());
+      EXPECT_FLOAT_EQ(outputIt.Get()[1], 2 * expectedIt.Get());
+    }
+  }
+}
+
+TEST(SmoothingRecursiveGaussianImageFilter, IntegerPixelsWithShortFirstPipelineAxis)
+{
+  using IntegerImageType = itk::Image<unsigned char, 3>;
+  const auto                       image = MakeImage({ { 9, 9, 1 } });
+  const FilterType::SigmaArrayType sigma{ { 0.8, 1.7, 2.5 } };
+  const auto                       expected = SmoothEligibleAxes(image, sigma);
+  for (const bool inPlace : { false, true })
+  {
+    auto cast = itk::CastImageFilter<ImageType, IntegerImageType>::New();
+    cast->SetInput(image);
+    cast->Update();
+    auto filter = itk::SmoothingRecursiveGaussianImageFilter<IntegerImageType>::New();
+    filter->SetInput(cast->GetOutput());
+    filter->SetSigmaArray(sigma);
+    filter->SetInPlace(inPlace);
+    ASSERT_NO_THROW(filter->Update());
+    itk::ImageRegionConstIterator<ImageType>        expectedIt(expected, expected->GetLargestPossibleRegion());
+    itk::ImageRegionConstIterator<IntegerImageType> outputIt(filter->GetOutput(), expected->GetLargestPossibleRegion());
+    for (; !expectedIt.IsAtEnd(); ++expectedIt, ++outputIt)
+    {
+      EXPECT_EQ(outputIt.Get(), static_cast<unsigned char>(expectedIt.Get()));
+    }
+  }
+}
+
+TEST(SmoothingRecursiveGaussianImageFilter, InvalidSigmaAndEmptyAxesStillThrow)
+{
+  for (const ImageType::SizeType size : { ImageType::SizeType{ { 9, 9, 1 } }, ImageType::SizeType{ { 1, 2, 3 } } })
+  {
+    const auto image = MakeImage(size);
+    for (const double sigma : { 0.0, -1.0 })
+    {
+      auto filter = FilterType::New();
+      filter->SetInput(image);
+      filter->SetSigmaArray(FilterType::SigmaArrayType{ { 1.0, 1.0, sigma } });
+      EXPECT_THROW(filter->Update(), itk::ExceptionObject);
+    }
+  }
+  auto filter = FilterType::New();
+  filter->SetInput(MakeImage({ { 0, 9, 9 } }));
+  EXPECT_THROW(filter->Update(), itk::ExceptionObject);
+}
diff --git a/Modules/Registration/RegistrationMethodsv4/itk-module.cmake b/Modules/Registration/RegistrationMethodsv4/itk-module.cmake
@@ -10,6 +10,7 @@ itk_module(
     ITKOptimizersv4
     ITKMetricsv4
   TEST_DEPENDS
+    ITKGoogleTest
     ITKTestKernel
   DESCRIPTION "${DOCUMENTATION}"
 )