JP-3456 - Replace use of scipy.signal.medfilt in outlier detection (#…

…8033)

Co-authored-by: Howard Bushouse <[email protected]>

CHANGES.rst

@@ -33,7 +33,15 @@ general
 
 - Add lack of python 3.12 support to project metadata [#8042]
 
-- Increase asdf maximum version to 4 [#8018]
+- Increase asdf maximum version to 4. [#8018]
+
+- Remove upper version limit for scipy. [#8033]
+
+outlier_detection
+-----------------
+
+- Remove use of ``scipy.signal.medfilt`` which is undefined for ``nan``
+  inputs. [#8033]
 
 imprint
 -------

jwst/outlier_detection/outlier_detection_ifu.py

@@ -1,17 +1,28 @@
 """Class definition for performing outlier detection on IFU data."""
 
+import logging
+
 import numpy as np
+from skimage.util import view_as_windows
+
 from stdatamodels.jwst import datamodels
-from scipy.signal import medfilt
-from .outlier_detection import OutlierDetection
 from stdatamodels.jwst.datamodels import dqflags
-import logging
+from .outlier_detection import OutlierDetection
+
 log = logging.getLogger(__name__)
 log.setLevel(logging.DEBUG)
 
 __all__ = ["OutlierDetectionIFU"]
 
 
+def medfilt(arr, kern_size):
+    # scipy.signal.medfilt (and many other median filters) have undefined behavior
+    # for nan inputs. See: https://github.com/scipy/scipy/issues/4800
+    padded = np.pad(arr, [[k // 2] for k in kern_size])
+    windows = view_as_windows(padded, kern_size, np.ones(len(kern_size), dtype='int'))
+    return np.nanmedian(windows, axis=np.arange(-len(kern_size), 0))
+
+
 class OutlierDetectionIFU(OutlierDetection):
     """Sub-class defined for performing outlier detection on IFU data.
 
@@ -220,7 +231,7 @@ def flag_outliers(self, idet, uq_det, ndet_files,
         minarr = np.nanmin(diffarr, axis=0)
 
         # Normalise the differences to a local median image to deal with ultra-bright sources
-        normarr = medfilt(minarr, kernel_size=kern_size)
+        normarr = medfilt(minarr, kern_size)
         nfloor = np.nanmedian(minarr)/3
         normarr[normarr < nfloor] = nfloor  # Ensure we never divide by a tiny number
         minarr_norm = minarr / normarr

jwst/outlier_detection/tests/test_medfilt.py

@@ -0,0 +1,42 @@
+import warnings
+
+import pytest
+import numpy as np
+import scipy.signal
+
+from jwst.outlier_detection.outlier_detection_ifu import medfilt
+
+
+
+@pytest.mark.parametrize("shape,kern_size", [
+    ([7, 7], [3, 3]),
+    ([7, 7], [3, 1]),
+    ([7, 7], [1, 3]),
+    ([7, 5], [3, 3]),
+    ([5, 7], [3, 3]),
+    ([42, 42], [7, 7]),
+    ([42, 42], [7, 5]),
+    ([42, 42], [5, 7]),
+    ([42, 7, 5], [3, 3, 3]),
+    ([5, 7, 42], [5, 5, 5]),
+])
+def test_medfilt_against_scipy(shape, kern_size):
+    arr = np.arange(np.prod(shape), dtype='uint32').reshape(shape)
+    result = medfilt(arr, kern_size)
+    expected = scipy.signal.medfilt(arr, kern_size)
+    np.testing.assert_allclose(result, expected)
+
+
+@pytest.mark.parametrize("arr,kern_size,expected", [
+    ([2, np.nan, 0], [3], [1, 1, 0]),
+    ([np.nan, np.nan, np.nan], [3], [0, np.nan, 0]),
+])
+def test_medfilt_nan(arr, kern_size, expected):
+    with warnings.catch_warnings():
+        warnings.filterwarnings(
+            "ignore",
+            message="All-NaN slice",
+            category=RuntimeWarning
+        )
+        result = medfilt(arr, kern_size)
+    np.testing.assert_allclose(result, expected)

setup.cfg

@@ -43,7 +43,8 @@ install_requires =
     poppy>=1.0.2
     pyparsing>=2.2.1
     requests>=2.22
-    scipy>=1.6.0,<1.10.0
+    scikit-image>=0.19
+    scipy>=1.6.0
     spherical-geometry>=1.2.22
     stcal>=1.4.4,<1.5.0
     stdatamodels>=1.8.3,<1.9.0
@@ -53,7 +54,7 @@ install_requires =
     tweakwcs>=0.8.3
     asdf-astropy>=0.3.0
     wiimatch>=0.2.1
-    packaging>19.0
+    packaging>20.0
     importlib-metadata>=4.11.4
     jsonschema>=4.8