Repository: scipy/scipy
Pull request: Add axes argument to ndimage filters
closes #10554
see also requests/discussion at: https://github.com/scipy/scipy/issues/12997#issuecomment-811268933 and https://github.com/scipy/scipy/issues/16191#issuecomment-1144781807

cc @mdhaber 

# description

This PR adds an "axes" keyword argument to most `scipy.ndimage` filters (1D filters already have an `axis` argument so this PR applies only to the others). Having an `axes` keyword-only argument allows easier filtering of a subset of the axes of an image. For example to apply a 2D footprint along 2 axes of an n-dimensional array, the user could either:
1.) expand the footprint to size one along the other axes. 
or 
2.) applied for loop(s) over non-filtered axes
However, this is a bit error prone and not always obvious to users, so I think it is helpful to provide this via an `axes` argument.

The functions for which a new `axes` argument can be introduced are:
- convolve
- correlate
- gaussian_filter
- gaussian_gradient_magnitude
- gaussian_laplace
- generic_gradient_magnitude
- generic_laplace
- laplace
- maximum_filter
- median_filter
- minimum_filter
- percentile_filter
- rank_filter
- uniform_filter
- generic_filter

For now, I reverted back to including only `gaussian_filter` in this PR. I can push other filters in a follow up.

This could potentially be extended to morphology functions (e.g. `erode`, `dilate`). 


Problem statement:
ndimage.gaussian_filter provide axis option
**Is your feature request related to a problem? Please describe.**
I have a stacked array of 2d images, which is a 3d array. I want to use ndimage.gaussian_filter for each 2d image in the array. However, current method does 3d gaussian, which I did not realize at the beginning. 

**Describe the solution you'd like**
Will the Scipy provides an option for axis and dimension in ndimage.gaussian_filter? It would be much more convenient and faster than doing it a loop. Thanks!

Hints from issue discussion:
I had the same problem, but was able to work around it by setting `sigma` to a tuple of standard deviations (instead of a single float) and fixing the standard deviation for the color channel dimension to 0: 
```python3
import numpy as np
import scipy.ndimage

#  create test image with shape (2, 2, 2), first channel is color channel dimension
img = np.stack([np.zeros((2,2)), np.ones((2,2))])

# per default this function applies blur over all axes
print(scipy.ndimage.gaussian_filter(img, sigma=1)) 
# output
#[[[0.35438556 0.35438556]
#  [0.35438556 0.35438556]]
#
# [[0.64561444 0.64561444]
#  [0.64561444 0.64561444]]]

# only apply filter to the second and third axis
print(scipy.ndimage.gaussian_filter(img, sigma=(0,1,1))) 
# output:
#[[[0. 0.]
#  [0. 0.]]
#
# [[1. 1.]
#  [1. 1.]]]
```
However, I agree that an argument to specify the axis would be useful. The default behavior of this function is surprising and quite unintuitive.
This would be a fairly straightforward feature to add at the Python level without needing to modify the underlying ndimage C code. It would basically be updating `sigma` for the users to insert 0 for any axes not to be filtered. 

I would mainly be in favor if we could provide the same API across most (all?) filtering/morphology functions. e.g. similar logic can be used in many of the other functions (e.g. expand a kernel to size 1 on any unfiltered indices for convolution functions).
Here is a related previous discussion: https://github.com/scipy/scipy/issues/12997#issuecomment-811268933
There was also a request in https://github.com/scipy/scipy/issues/16191#issuecomment-1144781807. i'd be happy to review. If the signatures are consistent enough, would you consider using a decorator to do the work (iterating over slices of the ndarray and calling the function on each slice) or would the approach you suggested be much more efficient?
> If the signatures are consistent enough, would you consider using a decorator to do the work (iterating over slices of the ndarray and calling the function on each slice) or would the approach you suggested be much more efficient?

For a common case such as an RGB or RGBA image with a small channel axis, the iteration based approach would be fine, but I think In the general case it will likely be more efficient to do things in a single call rather than iterating over any non-filtered axes so any loops get applied at the C level. I haven't thought too much yet about whether some common utility function(s) or a decorator based approach might be preferable, but I would think either is possible.

Let me try and put something together today and then I will link here for discussion.
OK. It's probably worth trying both because then the one you don't choose could be used for writing property-based tests.

Commit message: also test order argument with axes specified

Patch to review:
diff --git a/scipy/ndimage/_filters.py b/scipy/ndimage/_filters.py
index 6804ef51fa87..d6ea8d2d2571 100644
--- a/scipy/ndimage/_filters.py
+++ b/scipy/ndimage/_filters.py
@@ -278,7 +278,8 @@ def gaussian_filter1d(input, sigma, axis=-1, order=0, output=None,
 
 @_ni_docstrings.docfiller
 def gaussian_filter(input, sigma, order=0, output=None,
-                    mode="reflect", cval=0.0, truncate=4.0, *, radius=None):
+                    mode="reflect", cval=0.0, truncate=4.0, *, radius=None,
+                    axes=None):
     """Multidimensional Gaussian filter.
 
     Parameters
@@ -306,6 +307,12 @@ def gaussian_filter(input, sigma, order=0, output=None,
         for all axes. If specified, the size of the kernel along each axis
         will be ``2*radius + 1``, and `truncate` is ignored.
         Default is None.
+    axes : tuple of int or None, optional
+        If None, it is assumed all axes are to be filtered. Otherwise, only
+        the axes specified in the tuple will be considered. When `axes` is
+        specified, any tuples used for `sigma`, `order`, `mode` and/or `radius`
+        must match the length of `axes`. The ith entry in any of these tuples
+        corresponds to the ith entry in `axes.
 
     Returns
     -------
@@ -356,13 +363,15 @@ def gaussian_filter(input, sigma, order=0, output=None,
     """
     input = numpy.asarray(input)
     output = _ni_support._get_output(output, input)
-    orders = _ni_support._normalize_sequence(order, input.ndim)
-    sigmas = _ni_support._normalize_sequence(sigma, input.ndim)
-    modes = _ni_support._normalize_sequence(mode, input.ndim)
-    radiuses = _ni_support._normalize_sequence(radius, input.ndim)
-    axes = list(range(input.ndim))
+
+    axes = _ni_support._check_axes(axes, input.ndim)
+    num_axes = len(axes)
+    orders = _ni_support._normalize_sequence(order, num_axes)
+    sigmas = _ni_support._normalize_sequence(sigma, num_axes)
+    modes = _ni_support._normalize_sequence(mode, num_axes)
+    radiuses = _ni_support._normalize_sequence(radius, num_axes)
     axes = [(axes[ii], sigmas[ii], orders[ii], modes[ii], radiuses[ii])
-            for ii in range(len(axes)) if sigmas[ii] > 1e-15]
+            for ii in range(num_axes) if sigmas[ii] > 1e-15]
     if len(axes) > 0:
         for axis, sigma, order, mode, radius in axes:
             gaussian_filter1d(input, sigma, axis, order, output,
diff --git a/scipy/ndimage/_ni_support.py b/scipy/ndimage/_ni_support.py
index e8f39ed540c5..dc3d6dfb3a37 100644
--- a/scipy/ndimage/_ni_support.py
+++ b/scipy/ndimage/_ni_support.py
@@ -29,6 +29,7 @@
 # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 from collections.abc import Iterable
+import operator
 import warnings
 import numpy
 
@@ -95,3 +96,23 @@ def _get_output(output, input, shape=None, complex_output=False):
     elif complex_output and output.dtype.kind != 'c':
         raise RuntimeError("output must have complex dtype")
     return output
+
+
+def _check_axes(axes, ndim):
+    if axes is None:
+        return tuple(range(ndim))
+    elif numpy.isscalar(axes):
+        axes = (operator.index(axes),)
+    elif isinstance(axes, Iterable):
+        for ax in axes:
+            axes = tuple(operator.index(ax) for ax in axes)
+            if ax < -ndim or ax > ndim - 1:
+                raise ValueError(f"specified axis: {ax} is out of range")
+        axes = tuple(ax % ndim if ax < 0 else ax for ax in axes)
+    else:
+        raise ValueError(
+            "axes must be an integer, iterable of integers or None"
+        )
+    if len(tuple(set(axes))) != len(axes):
+        raise ValueError("axes must be unique")
+    return axes
diff --git a/scipy/ndimage/tests/test_filters.py b/scipy/ndimage/tests/test_filters.py
index 1d3fe10ec4ee..3d2e60ca1965 100644
--- a/scipy/ndimage/tests/test_filters.py
+++ b/scipy/ndimage/tests/test_filters.py
@@ -1,5 +1,6 @@
 ''' Some tests for filters '''
 import functools
+import itertools
 import math
 import numpy
 
@@ -691,6 +692,100 @@ def test_gauss_memory_overlap(self):
         ndimage.gaussian_filter(input, 1.0, output=input)
         assert_array_almost_equal(output1, input)
 
+    @pytest.mark.parametrize(
+        'axes',
+        tuple(itertools.combinations(range(-3, 3), 1))
+        + tuple(itertools.combinations(range(-3, 3), 2))
+        + ((0, 1, 2),))
+    def test_gauss_axes(self, axes):
+        array = numpy.arange(6 * 8 * 12, dtype=numpy.float64).reshape(6, 8, 12)
+        sigma = 1.0
+        axes = tuple(ax % array.ndim for ax in axes)
+        if len(tuple(set(axes))) != len(axes):
+            # parametrized cases with duplicate axes raise an error
+            with pytest.raises(ValueError):
+                ndimage.gaussian_filter(array, sigma, axes=axes)
+            return
+        output = ndimage.gaussian_filter(array, sigma, axes=axes)
+
+        # result should be equivalent to sigma=0.0 on unfiltered axes
+        all_sigmas = (sigma if ax in axes else 0.0 for ax in range(array.ndim))
+        expected = ndimage.gaussian_filter(array, all_sigmas)
+        assert_allclose(output, expected)
+
+    @pytest.mark.parametrize(
+        'axes', tuple(itertools.combinations(range(-3, 3), 2))
+    )
+    def test_gauss_axes_kwargs(self, axes):
+        array = numpy.arange(6 * 8 * 12, dtype=numpy.float64).reshape(6, 8, 12)
+        sigma = (1.0, 0.5)
+        radius = (4, 2)
+        order = (0, 1)
+        mode = ('reflect', 'nearest')
+        axes = tuple(ax % array.ndim for ax in axes)
+        kwargs = dict(sigma=sigma, radius=radius, mode=mode, order=order)
+        if len(tuple(set(axes))) != len(axes):
+            # parametrized cases with duplicate axes raise an error
+            with pytest.raises(ValueError):
+                ndimage.gaussian_filter(array, axes=axes, **kwargs)
+            return
+        output = ndimage.gaussian_filter(array, axes=axes, **kwargs)
+
+        # result should be equivalent to sigma=0.0 on unfiltered axes
+        sigma_3d = [0,] * array.ndim
+        sigma_3d[axes[0]] = sigma[0]
+        sigma_3d[axes[1]] = sigma[1]
+        radius_3d = [0,] * array.ndim
+        radius_3d[axes[0]] = radius[0]
+        radius_3d[axes[1]] = radius[1]
+        mode_3d = [0,] * array.ndim
+        mode_3d[axes[0]] = mode[0]
+        mode_3d[axes[1]] = mode[1]
+        order_3d = [0,] * array.ndim
+        order_3d[axes[0]] = order[0]
+        order_3d[axes[1]] = order[1]
+        kwargs = dict(sigma=sigma_3d, radius=radius_3d, mode=mode_3d,
+                      order=order_3d)
+        expected = ndimage.gaussian_filter(array, **kwargs)
+        assert_allclose(output, expected)
+
+    @pytest.mark.parametrize(
+        'axes', [(1.5,), (0, 1, 2, 3), (3,), (-4,)]
+    )
+    def test_gauss_invalid_axes(self, axes):
+        array = numpy.arange(6 * 8 * 12, dtype=numpy.float64).reshape(6, 8, 12)
+        if any(isinstance(ax, float) for ax in axes):
+            error_class = TypeError
+        else:
+            error_class = ValueError
+        with pytest.raises(error_class):
+            ndimage.gaussian_filter(array, sigma=1.0, axes=axes)
+
+    @pytest.mark.parametrize('tuple_ndim', [1, 3])
+    def test_gauss_axes_tuple_length_mismatch(self, tuple_ndim):
+        array = numpy.arange(6 * 8 * 12, dtype=numpy.float64).reshape(6, 8, 12)
+        axes = (0, 1)
+
+        # len(radius) != len(axes)
+        radius = (3,) * tuple_ndim
+        with pytest.raises(RuntimeError):
+            ndimage.gaussian_filter(array, sigma=1.0, radius=radius, axes=axes)
+
+        # len(mode) != len(axes)
+        mode = ['constant'] * tuple_ndim
+        with pytest.raises(RuntimeError):
+            ndimage.gaussian_filter(array, sigma=1.0, mode=mode, axes=axes)
+
+        # len(sigma) != len(axes)
+        sigma = (1.0,) * tuple_ndim
+        with pytest.raises(RuntimeError):
+            ndimage.gaussian_filter(array, sigma=sigma, axes=axes)
+
+        # len(order) != len(axes)
+        order = (0,) * tuple_ndim
+        with pytest.raises(RuntimeError):
+            ndimage.gaussian_filter(array, sigma=1.0, order=order, axes=axes)
+
     @pytest.mark.parametrize('dtype', types + complex_types)
     def test_prewitt01(self, dtype):
         array = numpy.array([[3, 2, 5, 1, 4],


Code review comments:
diff_hunk: @@ -306,6 +307,12 @@ def gaussian_filter(input, sigma, order=0, output=None,
         for all axes. If specified, the size of the kernel along each axis
         will be ``2*radius + 1``, and `truncate` is ignored.
         Default is None.
+    axes : tuple of int or None, optional
+        If None, it is assumed all axes are to be filtered. Otherwise, only
path: scipy/ndimage/_filters.py
text: @user1:
Nit: the axes themselves will not be filtered.
```suggestion
        If None, `input` is filtered along all axes. Otherwise,
```

diff_hunk: @@ -306,6 +307,12 @@ def gaussian_filter(input, sigma, order=0, output=None,
         for all axes. If specified, the size of the kernel along each axis
         will be ``2*radius + 1``, and `truncate` is ignored.
         Default is None.
+    axes : tuple of int or None, optional
+        If None, it is assumed all axes are to be filtered. Otherwise, only
+        the axes specified in the tuple will be considered. When `axes` is
+        specified, any tuples used for `sigma`, `order`, `mode` and/or `radius`
+        must match the length of `axes`. The ith entry in any of these tuples
+        corresponds to the ith entry in `axes.
path: scipy/ndimage/_filters.py
text: @user1:
```suggestion
        corresponds to the ith entry in `axes`.
```

diff_hunk: @@ -306,6 +307,12 @@ def gaussian_filter(input, sigma, order=0, output=None,
         for all axes. If specified, the size of the kernel along each axis
         will be ``2*radius + 1``, and `truncate` is ignored.
         Default is None.
+    axes : tuple of int or None, optional
+        If None, it is assumed all axes are to be filtered. Otherwise, only
+        the axes specified in the tuple will be considered. When `axes` is
path: scipy/ndimage/_filters.py
text: @user1:
```suggestion
        `input` is filtered along the specified axes. When `axes` is
```