ENH: Support singleton and short dimensions in SmoothingRecursiveGaussianImageFilter #6997

Description

SmoothingRecursiveGaussianImageFilter currently requires at least four pixels along every image dimension, even when smoothing along a particular dimension is unnecessary.

This prevents Gaussian smoothing of images containing singleton dimensions, such as 2D slices represented as 3D images (1 × Y × X). In particular, it prevents using these images with multi-resolution registration in ImageRegistrationMethodv4, unless smoothing is disabled entirely or the singleton dimension is artificially expanded.

Would it be possible to support images with dimensions smaller than four pixels, by skipping smoothing for these dimensions?

If that's something that would be interesting for ITK, I'd be up to trying my hands on a PR!

Expected behavior

• Dimensions of size <4: Skip Gaussian smoothing along these dimensions while smoothing the remaining dimensions normally.
• Preserve the original image dimensionality and physical geometry.
• Preserve existing behavior and performance for images with at least four pixels along every dimension.

Actual behavior

SmoothingRecursiveGaussianImageFilter::GenerateData() currently throws an exception if any dimension contains fewer than four pixels, regardless of whether smoothing along that dimension is necessary.

This can be reproduced with the following minimal registration example using ITK Python:

import itk
import numpy as np

# Two 3D images with a singleton Z dimension.
y, x = np.mgrid[:64, :64]

fixed = itk.image_from_array(
    np.exp(-((x - 28)**2 + (y - 32)**2) / 100)[None].astype("float32")
)
moving = itk.image_from_array(
    np.exp(-((x - 32)**2 + (y - 32)**2) / 100)[None].astype("float32")
)

ImageType = type(fixed)
RegistrationType = itk.ImageRegistrationMethodv4[ImageType, ImageType]
MetricType = itk.MeanSquaresImageToImageMetricv4[ImageType, ImageType]


def register(sigma):
    metric = MetricType.New()

    # Avoid unrelated recursive Gaussian gradient filtering.
    metric.SetUseFixedImageGradientFilter(False)
    metric.SetUseMovingImageGradientFilter(False)

    registration = RegistrationType.New(
        FixedImage=fixed,
        MovingImage=moving,
        Metric=metric,
        Optimizer=itk.RegularStepGradientDescentOptimizerv4[itk.D].New(
            LearningRate=1.0,
            MinimumStepLength=0.001,
            NumberOfIterations=25,
        ),
        InitialTransform=itk.TranslationTransform[itk.D, 3].New(),
    )
    registration.SetNumberOfLevels(1)
    registration.SetShrinkFactorsPerLevel([1])
    registration.SetSmoothingSigmasPerLevel([sigma])

    registration.Update()
    return list(registration.GetTransform().GetParameters())


for sigma in (1.0, 0.0):
    try:
        print(f"sigma={sigma}: {register(sigma)}")
    except RuntimeError as error:
        print(f"sigma={sigma}: FAILED")
        print(error)

With sigma=1.0, registration fails:

ITK ERROR: SmoothingRecursiveGaussianImageFilter:
The number of pixels along dimension 2 is less than 4.
This filter requires a minimum of four pixels
along the dimension to be processed.

With sigma=0.0, registration completes successfully.

Disabling smoothing avoids the error, but also prevents Gaussian smoothing along the non-singleton dimensions, which is undesirable for multi-resolution registration.

Additional Information

This issue isn't new at all, I found several discussions on this topic:

• ITK Discourse: SimpleITK Image registration 4.13 vs 5.1 (2020). Reports the same error during slice-to-volume registration. The discussion identifies the switch from discrete to recursive Gaussian smoothing in the registration framework and suggests expanding the singleton dimension or applying smoothing outside registration.
• ITK Discourse: Plot TRE per iteration for 2D registration (2024). Reports the same exception when converting a 2D image into a singleton-axis 3D image.
• ITK Discourse: Registering with only one slice of 3D image (2024). Reports the four-pixel restriction when registering single-slice DICOM images.
• SimpleElastix #10: Error after run elastix.Execute() (2016). Reports the same underlying recursive Gaussian restriction. Suggested workarounds include extracting a genuine 2D image or using shrinking-only image pyramids.
• ANTs #324: antsBrainExtraction.sh fails at fine-resolution (2016). Reports the recursive Gaussian four-pixel restriction when preprocessing generates images with very small dimensions.
• ANTs discussion: Slice-to-volume registration (2013). Discusses representing 2D slices as singleton-axis 3D images for slice-to-volume registration.

Interestingly, elastix's MultiResolutionGaussianSmoothingPyramidImageFilter already handles singleton dimensions:

When either a shrink factor is set to zero for a particular dimension, or when the dimension itself is 1, or both, no smoothing is performed along that dimension.

This filter uses recursive Gaussian smoothing internally but explicitly skips singleton dimensions. A similar approach in SmoothingRecursiveGaussianImageFilter could address the singleton-dimension case without requiring changes to the registration framework.