convert GaussianFilter kwargs to python native (#159)

* convert GaussianFilter kwargs to python native

* bump version

* add more tests

Author: gschramm

docs/changelog.rst

@@ -4,6 +4,22 @@ Changelog
 2.x
 ---
 
+2.0.2 (Jul 29, 2026)
+^^^^^^^^^^^^^^^^^^^^
+
+Bug fixes
+~~~~~~~~~
+
+- **``GaussianFilterOperator``: fully fix tensor-valued ``sigma``.** The 2.0.1
+  fix converted the input array to NumPy, but the actual trigger is a *tensor*
+  ``sigma``: ``scipy.ndimage`` builds the Gaussian kernel in ``sigma``'s array
+  namespace, so a torch / CuPy tensor ``sigma`` yields a tensor kernel that
+  ``gaussian_filter1d`` reverses with a negative-step slice (``weights[::-1]``),
+  which PyTorch does not support — failing even for a plain NumPy input array.
+  ``GaussianFilterOperator`` now coerces ``sigma`` (and other array-valued
+  kwargs) to native Python values (the input NumPy-view safeguard from 2.0.1 is
+  retained). This fully fixes downstream use via e.g. DeepInv.
+
 2.0.1 (Jul 24, 2026)
 ^^^^^^^^^^^^^^^^^^^^
 

pyproject.toml

@@ -16,7 +16,7 @@ build-backend = "hatchling.build"
 [project]
 name = "parallelproj"
 # IMPORTANT: bumping a release requires editing this one line.
-version = "2.0.2.dev0"
+version = "2.0.2"
 description = "Python tools for PET projection and reconstruction workflows."
 readme = "README.md"
 requires-python = ">=3.12"

src/parallelproj/operators.py

@@ -17,16 +17,21 @@
 import numpy as np
 import array_api_compat
 
-# GPU arrays (CuPy / PyTorch CUDA) are filtered with ``cupyx.scipy.ndimage``
-# directly (see GaussianFilterOperator); CPU arrays (NumPy, PyTorch CPU,
-# array-api-strict) are converted to a NumPy view first and then filtered with
-# ``scipy.ndimage``.  Neither path relies on scipy's array-API *delegation*
-# (where scipy would compute natively in the input's namespace): under
-# delegation ``gaussian_filter1d`` reverses the kernel with a negative-step
-# slice ``weights[::-1]`` that e.g. PyTorch does not support.  Converting to
-# NumPy ourselves makes the operator robust regardless of the scipy version or
-# the ``SCIPY_ARRAY_API`` env var (and independent of the scipy / parallelproj
-# import order).
+# GaussianFilterOperator filters GPU arrays (CuPy / PyTorch CUDA) with
+# ``cupyx.scipy.ndimage`` and CPU arrays (NumPy, PyTorch CPU, array-api-strict)
+# with ``scipy.ndimage``.  Two precautions keep it backend-robust:
+#   * filter kwargs such as ``sigma`` are coerced to native Python values
+#     (see ``_to_builtin``).  scipy builds the Gaussian kernel in ``sigma``'s
+#     array namespace, so a *tensor* ``sigma`` yields a tensor kernel that
+#     ``gaussian_filter1d`` then reverses with a negative-step slice
+#     ``weights[::-1]`` -- unsupported by e.g. PyTorch.  This is the root cause
+#     of the failure; it fires even for a plain NumPy input array.
+#   * on the CPU path the input is additionally converted to a NumPy view, so
+#     scipy never relies on its array-API *delegation* (scipy >= 1.16 /
+#     ``SCIPY_ARRAY_API``) to compute natively in the input's namespace.
+# Together the scipy call is always the canonical ``gaussian_filter(numpy,
+# float)`` form, independent of scipy version and scipy / parallelproj import
+# order.
 import scipy.ndimage as ndimage
 from array_api_compat import device, get_namespace
 
@@ -541,6 +546,30 @@ def iscomplex(self) -> bool:
         ) or self.xp.isdtype(self._values.dtype, self.xp.complex128)
 
 
+def _to_builtin(value):
+    """Coerce an array/tensor-valued filter kwarg to native Python.
+
+    scipy.ndimage builds the Gaussian kernel in the array namespace of
+    ``sigma``, so a torch / CuPy / array-api tensor ``sigma`` yields a tensor
+    kernel that ``gaussian_filter1d`` reverses with a negative-step slice
+    (``weights[::-1]``), which PyTorch does not support.  Passing native Python
+    scalars / lists avoids this for every backend.
+
+    Plain Python scalars, strings and ``None`` are returned unchanged;
+    lists / tuples are converted element-wise (preserving the container type);
+    NumPy scalars and 0-d / 1-d arrays / tensors become a Python scalar or list.
+    """
+    if isinstance(value, np.generic):
+        return value.item()
+    if value is None or isinstance(value, (bool, int, float, str)):
+        return value
+    if isinstance(value, (list, tuple)):
+        return type(value)(_to_builtin(v) for v in value)
+    # numpy / torch / cupy / array-api arrays (to_numpy_array always returns a
+    # NumPy ndarray, incl. a device->host copy for GPU arrays)
+    return to_numpy_array(value).tolist()
+
+
 class GaussianFilterOperator(LinearOperator):
     """Isotropic Gaussian smoothing operator (self-adjoint).
 
@@ -565,10 +594,15 @@ def __init__(self, in_shape: tuple[int, ...], **kwargs):
         **kwargs : dict
             passed to scipy.ndimage.gaussian_filter; most commonly ``sigma``
             (standard deviation in pixels), plus optional ``mode``, ``truncate``, etc.
+            Array/tensor-valued arguments (e.g. a torch tensor ``sigma``) are
+            coerced to native Python values so scipy builds the kernel in NumPy
+            (see :func:`_to_builtin`).
         """
         super().__init__()
         self._in_shape = in_shape
-        self._kwargs = kwargs
+        # coerce array/tensor kwargs (most importantly ``sigma``) to native
+        # Python so scipy never builds the kernel in a tensor namespace.
+        self._kwargs = {k: _to_builtin(v) for k, v in kwargs.items()}
 
     @property
     def in_shape(self) -> tuple[int, ...]:
@@ -605,14 +639,12 @@ def _apply(self, x: Array) -> Array:
             return xp.asarray(xp.from_dlpack(y_cp))
 
         # CPU arrays (NumPy, PyTorch CPU, array-api-strict): filter on a NumPy
-        # view, then convert the result back to the input's namespace/device.
-        # We convert to NumPy *ourselves* instead of passing ``x`` to scipy,
-        # because modern (array-API-aware) scipy would otherwise compute
-        # natively in the input's namespace, where ``gaussian_filter1d``
-        # reverses the kernel with a negative-step slice ``weights[::-1]`` that
-        # e.g. PyTorch does not support.  On CPU the torch / array-api-strict
-        # <-> NumPy conversions are zero-copy (shared buffer); only scipy's
-        # output allocation remains, exactly as in the pure-NumPy path.
+        # view so scipy never delegates to the input's array namespace (scipy
+        # >= 1.16 / SCIPY_ARRAY_API), then convert the result back.  ``sigma``
+        # & co. were already coerced to native Python in __init__.  On CPU the
+        # torch / array-api-strict <-> NumPy conversions are zero-copy (shared
+        # buffer); only scipy's output allocation remains, as in the pure-NumPy
+        # path.
         x_np = np.asarray(to_numpy_array(x))
         result = ndimage.gaussian_filter(x_np, **self._kwargs)
         return xp.asarray(result, device=dev, dtype=x.dtype)

tests/test_gaussian_filter_torch.py

@@ -1,13 +1,16 @@
-"""Regression tests for :class:`GaussianFilterOperator` with PyTorch CPU tensors.
-
-Guards against a scipy array-API *delegation* issue: when scipy computes
-natively on a torch tensor, ``scipy.ndimage.gaussian_filter1d`` reverses the
-Gaussian kernel with a negative-step slice (``weights[::-1]``) that PyTorch does
-not support (https://github.com/pytorch/pytorch/issues/175240).  parallelproj
-must therefore filter on a NumPy view, so this never happens regardless of the
-scipy version or the ``SCIPY_ARRAY_API`` env var.
-
-These tests are torch-CPU specific and are intentionally *not* parametrized over
+"""Regression tests for :class:`GaussianFilterOperator` with PyTorch.
+
+Root cause guarded here: ``scipy.ndimage`` builds the Gaussian kernel in the
+array namespace of ``sigma``.  A *tensor* ``sigma`` therefore yields a tensor
+kernel that ``gaussian_filter1d`` reverses with a negative-step slice
+(``weights[::-1]``), which PyTorch does not support
+(https://github.com/pytorch/pytorch/issues/175240) -- this fires even for a
+plain NumPy input array.  As defence-in-depth against scipy's array-API
+delegation (scipy >= 1.16 / ``SCIPY_ARRAY_API``), the CPU input is also filtered
+via a NumPy view.  parallelproj must therefore (a) coerce ``sigma`` to native
+Python and (b) not rely on scipy's delegation.
+
+These tests are torch specific and intentionally *not* parametrized over
 ``(xp, dev)`` (so they do not import ``config.pytestmark``).
 """
 
@@ -26,8 +29,36 @@
 pytestmark = pytest.mark.skipif(not torch_available, reason="torch not installed")
 
 
-def test_gaussian_filter_operator_torch_cpu() -> None:
-    """A torch CPU tensor stays a torch CPU tensor and matches the NumPy result."""
+def test_gaussian_filter_operator_tensor_sigma() -> None:
+    """A tensor-valued ``sigma`` must not break the filter (the DeepInv bug).
+
+    Uses a plain NumPy input on purpose to show the failure is driven by the
+    ``sigma`` type, not the input array type.
+    """
+    import numpy as np
+    import torch
+    from parallelproj.operators import GaussianFilterOperator
+
+    shape = (8, 8, 4)
+    x = np.ones(shape, dtype=np.float32)
+
+    ref = GaussianFilterOperator(shape, sigma=1.5)(x)
+
+    # scalar tensor sigma and per-axis tensor sigma
+    y_scalar = GaussianFilterOperator(shape, sigma=torch.tensor(1.5))(x)
+    y_vec = GaussianFilterOperator(shape, sigma=torch.tensor([1.5, 1.5, 1.5]))(x)
+    assert np.allclose(np.asarray(y_scalar), ref, atol=1e-6)
+    assert np.allclose(np.asarray(y_vec), ref, atol=1e-6)
+
+    # and with a torch input array (result stays a torch CPU tensor)
+    xt = torch.ones(shape, dtype=torch.float32)
+    yt = GaussianFilterOperator(shape, sigma=torch.tensor(1.5))(xt)
+    assert isinstance(yt, torch.Tensor) and yt.device.type == "cpu"
+    assert np.allclose(np.asarray(yt), ref, atol=1e-6)
+
+
+def test_gaussian_filter_operator_torch_cpu_float_sigma() -> None:
+    """Basic backend check: torch CPU input + float sigma stays a torch tensor."""
     import numpy as np
     import torch
     from scipy.ndimage import gaussian_filter
@@ -36,15 +67,10 @@ def test_gaussian_filter_operator_torch_cpu() -> None:
 
     shape = (8, 8, 4)
     op = GaussianFilterOperator(shape, sigma=1.5)
+    y = op(torch.ones(shape, dtype=torch.float32))
 
-    x = torch.ones(shape, dtype=torch.float32)
-    y = op(x)
-
-    assert isinstance(y, torch.Tensor)
-    assert y.device.type == "cpu"
+    assert isinstance(y, torch.Tensor) and y.device.type == "cpu"
     assert tuple(y.shape) == shape
-    assert bool(torch.isfinite(y).all())
-
     ref = gaussian_filter(np.ones(shape, dtype=np.float32), sigma=1.5)
     assert np.allclose(np.asarray(y), ref, atol=1e-6)
 
@@ -53,10 +79,9 @@ def test_gaussian_filter_operator_torch_cpu_scipy_array_api() -> None:
     """Force scipy's array-API delegation via ``SCIPY_ARRAY_API=1``.
 
     scipy reads ``SCIPY_ARRAY_API`` at import time, so this runs in a fresh
-    interpreter with the env var set (which also keeps the rest of the suite on
-    scipy's default behavior).  Under delegation the pre-fix operator hit the
-    unsupported ``weights[::-1]`` negative-step slice on modern scipy; with the
-    NumPy-view fix it succeeds.
+    interpreter with the env var set (keeping the rest of the suite on scipy's
+    default behaviour).  Guards the CPU NumPy-view path against scipy computing
+    natively in a torch input's namespace.
     """
     code = textwrap.dedent(
         """

tests/test_operators.py

@@ -484,3 +484,40 @@ def test_resolve_namespace_device_requires_xp(xp: ModuleType, dev: str) -> None:
     rxp, _ = mop._resolve_namespace_device(None, None)
     assert rxp is mop.xp
     assert mop._resolve_namespace_device(xp, "explicit-dev") == (xp, "explicit-dev")
+
+
+def test_gaussian_filter_to_builtin(xp: ModuleType, dev: str) -> None:
+    """Cover every branch of ``operators._to_builtin`` (sigma coercion).
+
+    Ensures array/tensor-valued filter kwargs are turned into native Python so
+    scipy never builds the Gaussian kernel in a tensor namespace.
+    """
+    import numpy as realnp
+
+    to_builtin = ppo._to_builtin
+
+    # plain python scalars / str / None / bool pass through unchanged
+    assert to_builtin(1.5) == 1.5 and isinstance(to_builtin(1.5), float)
+    assert to_builtin(3) == 3
+    assert to_builtin("reflect") == "reflect"
+    assert to_builtin(None) is None
+    assert to_builtin(True) is True
+
+    # numpy scalar (np.generic) -> python scalar
+    b = to_builtin(realnp.float32(1.5))
+    assert b == 1.5 and isinstance(b, float)
+
+    # list / tuple recurse, preserving container type and coercing elements
+    assert to_builtin((1.5, 2.5)) == (1.5, 2.5)
+    assert to_builtin([realnp.float32(1.0), 2]) == [1.0, 2]
+
+    # array / tensor branch (per backend): 1-d -> list, 0-d -> scalar
+    assert to_builtin(xp.asarray([1.5, 2.0, 2.5], device=dev)) == [1.5, 2.0, 2.5]
+    assert to_builtin(xp.asarray(1.5, device=dev)) == 1.5
+
+    # end-to-end: the operator accepts an array/tensor sigma without error
+    op = ppo.GaussianFilterOperator(
+        (6, 6, 4), sigma=xp.asarray([1.0, 1.0, 1.0], device=dev)
+    )
+    y = op(xp.ones((6, 6, 4), dtype=xp.float32, device=dev))
+    assert tuple(y.shape) == (6, 6, 4)