ENH: ndimage: log(n) implementation for 1D rank filter (scipy#20543)

* ENH: rank filter for 1D cases log(n) complexity implementation

Reviewed at scipy#20543

---------

Co-authored-by: Lucas Colley <[email protected]>

Author: ggkogan

scipy/ndimage/_filters.py

@@ -38,6 +38,7 @@
 from . import _ni_support
 from . import _nd_image
 from . import _ni_docstrings
+from . import _rank_filter_1d
 
 __all__ = ['correlate1d', 'convolve1d', 'gaussian_filter1d', 'gaussian_filter',
            'prewitt', 'sobel', 'generic_laplace', 'laplace',
@@ -1553,9 +1554,30 @@ def _rank_filter(input, rank, size=None, footprint=None, output=None,
             raise RuntimeError(
                 "A sequence of modes is not supported by non-separable rank "
                 "filters")
-        mode = _ni_support._extend_mode_to_code(mode)
-        _nd_image.rank_filter(input, rank, footprint, output, mode, cval,
-                              origins)
+        mode = _ni_support._extend_mode_to_code(mode, is_filter=True)
+        if input.ndim == 1:
+            if input.dtype in (np.int64, np.float64, np.float32):
+                x = input
+                x_out = output
+            elif input.dtype == np.float16:
+                x = input.astype('float32')
+                x_out = np.empty(x.shape, dtype='float32')
+            elif np.result_type(input, np.int64) == np.int64:
+                x = input.astype('int64')
+                x_out = np.empty(x.shape, dtype='int64')
+            elif input.dtype.kind in 'biu':
+                # cast any other boolean, integer or unsigned type to int64
+                x = input.astype('int64')
+                x_out = np.empty(x.shape, dtype='int64')
+            else:
+                raise RuntimeError('Unsupported array type')
+            cval = x.dtype.type(cval)
+            _rank_filter_1d.rank_filter(x, rank, footprint.size, x_out, mode, cval,
+                                        origin)
+            if input.dtype not in (np.int64, np.float64, np.float32):
+                np.copyto(output, x_out, casting='unsafe')
+        else:
+            _nd_image.rank_filter(input, rank, footprint, output, mode, cval, origins)
         if temp_needed:
             temp[...] = output
             output = temp

scipy/ndimage/_ni_support.py

@@ -34,7 +34,7 @@
 import numpy as np
 
 
-def _extend_mode_to_code(mode):
+def _extend_mode_to_code(mode, is_filter=False):
     """Convert an extension mode to the corresponding integer code.
     """
     if mode == 'nearest':
@@ -47,8 +47,12 @@ def _extend_mode_to_code(mode):
         return 3
     elif mode == 'constant':
         return 4
+    elif mode == 'grid-wrap' and is_filter:
+        return 1
     elif mode == 'grid-wrap':
         return 5
+    elif mode == 'grid-constant' and is_filter:
+        return 4
     elif mode == 'grid-constant':
         return 6
     else:
@@ -121,7 +125,6 @@ def _check_axes(axes, ndim):
         raise ValueError("axes must be unique")
     return axes
 
-
 def _skip_if_dtype(arg):
     """'array or dtype' polymorphism.
 
@@ -138,4 +141,3 @@ def _skip_if_dtype(arg):
 
 def _skip_if_int(arg):
     return None if (arg is None or isinstance(arg, int)) else arg
-

scipy/ndimage/meson.build

@@ -25,6 +25,14 @@ py3.extension_module('_ni_label',
   subdir: 'scipy/ndimage'
 )
 
+py3.extension_module('_rank_filter_1d',
+  'src/_rank_filter_1d.cpp',
+  link_args: version_link_args,
+  install: true,
+  dependencies: np_dep,
+  subdir: 'scipy/ndimage'
+)
+
 py3.extension_module('_ctest',
   'src/_ctest.c',
   dependencies: np_dep,