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

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',
@@ -1486,8 +1487,11 @@ def _rank_filter(input, rank, size=None, footprint=None, output=None,
                 "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)
+        if input.ndim == 1:
+            _rank_filter_1d.rank_filter_1d(input, rank, footprint, output, mode, cval,
+                                           origin)
+        else:
+            _nd_image.rank_filter(input, rank, footprint, output, mode, cval, origins)
         if temp_needed:
             temp[...] = output
             output = temp

scipy/ndimage/_rank_filter_1d.pyx

@@ -0,0 +1,45 @@
+cimport cython
+import numpy as np
+
+ctypedef fused numeric_t:
+    cython.float
+    cython.double
+    cython.longlong
+
+cdef extern from "rank_filter_1d.cpp" nogil:
+    void rank_filter[T](T* in_arr, int rank, int arr_len, int win_len, T* out_arr,
+                        int mode, T cval, int origin)
+
+def rank_filter_1d_cpp_api(numeric_t[:] in_arr,
+                           int rank, int win_len, numeric_t[:] out_arr, int mode,
+                           numeric_t cval=0, int origin=0):
+    rank_filter(&in_arr[0], rank, in_arr.shape[0], win_len, &out_arr[0], mode, cval,
+                origin)
+    return out_arr
+
+def rank_filter_1d(x: np.ndarray, rank, footprint: np.ndarray, x_out: np.ndarray,
+                   mode: int, cval=0.0, origin=0):
+    size = footprint.size
+    rank = int(rank)
+    origin = int(origin)
+    # legacy mode handling
+    if mode == 6:
+        mode = 4
+    if mode == 5:
+        mode = 1
+    if x.dtype.name in ['float32', 'float64', 'int64']:
+        cval = x.dtype.type(cval)
+        rank_filter_1d_cpp_api(in_arr=x, rank=rank, win_len=size,out_arr=x_out,
+                               mode=mode, cval=cval, origin=origin)
+        return x_out
+    elif x.dtype.name in ['int8', 'uint8', 'int16', 'uint16', 'uint32', 'int32']:
+        x_in = x.astype(np.int64)
+        x_out_ = np.empty_like(x_in, dtype=np.int64)
+        cval = x_in.dtype.type(cval)
+        rank_filter_1d_cpp_api(in_arr=x_in, rank=rank, win_len=size, out_arr=x_out_,
+                               mode=mode, cval=cval, origin=origin)
+        np.copyto(x_out, x_out_, casting='unsafe')
+        return x_out
+    else:
+        # the original function does not support uint64 (bug) and longdouble
+        raise ValueError('Unsupported input dtype')

scipy/ndimage/meson.build

@@ -25,6 +25,16 @@ py3.extension_module('_ni_label',
   subdir: 'scipy/ndimage'
 )
 
+py3.extension_module('_rank_filter_1d',
+  cython_gen_cpp.process('_rank_filter_1d.pyx'),
+  include_directories: 'src',
+  cpp_args: cython_cpp_args,
+  dependencies: np_dep,
+  link_args: version_link_args,
+  install: true,
+  subdir: 'scipy/ndimage'
+)
+
 py3.extension_module('_ctest',
   'src/_ctest.c',
   dependencies: np_dep,

scipy/ndimage/src/rank_filter_1d.cpp

@@ -0,0 +1,223 @@
+/*
+Started working on https://ideone.com/8VVEa, I optimized by restriction of cases and proper initialization,
+also adapted for rank filter rather than the original median filter. Allowed different boundary conditons and
+Moved to C++ for polymorphism.
+*/
+
+#include <stdlib.h>
+#include <stdio.h>
+
+struct Mediator//this is used for rank keeping
+{
+   int*  pos;   //index into `heap` for each value
+   int*  heap;  //max/rank/min heap holding indexes into `data`.
+   int   N;     //allocated size.
+   int   idx;   //position in circular queue
+   int   minCt; //count of items in min heap
+   int   maxCt; //count of items in max heap
+};
+
+typedef enum {
+   NEAREST = 0,
+   WRAP = 1,
+   REFLECT = 2,
+   MIRROR = 3,
+   CONSTANT = 4,
+} Mode;
+
+/*--- Helper Functions ---*/
+
+//returns 1 if heap[i] < heap[j]
+template <typename T>
+inline int mmless(T* data, Mediator* m, int i, int j)
+{
+   return (data[m->heap[i]] < data[m->heap[j]]);
+}
+
+//swaps items i&j in heap, maintains indexes
+int mmexchange(Mediator* m, int i, int j)
+{
+   int t = m->heap[i];
+   m->heap[i] = m->heap[j];
+   m->heap[j] = t;
+   m->pos[m->heap[i]] = i;
+   m->pos[m->heap[j]] = j;
+   return 1;
+}
+
+//swaps items i & j if i < j;  returns true if swapped
+template <typename T>
+inline int mmCmpExch(T* data, Mediator* m, int i, int j)
+{
+   return (mmless(data, m,i,j) && mmexchange(m,i,j));
+}
+
+//maintains minheap property for all items below i.
+template <typename T>
+void minSortDown(T* data, Mediator* m, int i)
+{
+   for (i*=2; i <= m->minCt; i*=2)
+   {  if (i < m->minCt && mmless(data, m, i+1, i)) { ++i; }
+      if (!mmCmpExch(data, m, i, i/2)) { break; }
+   }
+}
+
+//maintains maxheap property for all items below i. (negative indexes)
+template <typename T>
+void maxSortDown(T* data, Mediator* m, int i)
+{
+   for (i*=2; i >= -m->maxCt; i*=2)
+   {  if (i > -m->maxCt && mmless(data, m, i, i-1)) { --i;}
+      if (!mmCmpExch(data, m, i/2, i)) { break; }
+   }
+}
+
+//maintains minheap property for all items above i, including the rank
+//returns true if rank changed
+template <typename T>
+inline int minSortUp(T* data, Mediator* m, int i)
+{
+   while (i>0 && mmCmpExch(data, m, i, i/2)) i/=2;
+   return (i==0);
+}
+
+//maintains maxheap property for all items above i, including rank
+//returns true if rank changed
+template <typename T>
+inline int maxSortUp(T* data, Mediator* m, int i)
+{
+   while (i<0 && mmCmpExch(data, m, i/2, i))  i/=2;
+   return (i==0);
+}
+
+/*--- Public Interface ---*/
+
+
+//creates new Mediator: to calculate `nItems` running rank.
+Mediator* MediatorNew(int nItems, int rank)
+{
+   Mediator* m =  (Mediator*)malloc(sizeof(Mediator));
+   m->pos = (int*)malloc(sizeof(int) * nItems);
+   m->heap = (int*)malloc(sizeof(int) * nItems);
+   if ((m == nullptr)||(m->pos == nullptr)||(m->heap == nullptr)){printf("out of memory\n"); exit(1);}
+   m->heap += rank; //points to rank
+   m->N = nItems;
+   m->idx = 0;
+   m->minCt = nItems - rank - 1;
+   m->maxCt = rank;
+   while (nItems--)
+   {
+      m->pos[nItems]= nItems - rank;
+      m->heap[m->pos[nItems]]=nItems;
+   }
+   return m;
+}
+
+//Inserts item, maintains rank in O(lg nItems)
+template <typename T>
+void MediatorInsert(T* data, Mediator* m, T v)
+{
+   int p = m->pos[m->idx];
+   T old = data[m->idx];
+   data[m->idx] = v;
+   m->idx++;
+   if(m->idx == m->N){m->idx = 0; }
+
+   if (p > 0) //new item is in minHeap
+   {  if (v > old) { minSortDown(data, m, p); return; }
+      if (minSortUp(data, m, p) && mmCmpExch(data, m, 0, -1)) { maxSortDown(data, m,-1); }
+   }
+   else if (p < 0) //new item is in maxheap
+   {  if ( v < old) {maxSortDown(data, m, p); return; }
+      if (maxSortUp(data, m, p) && mmCmpExch(data, m, 1, 0)) { minSortDown(data, m, 1); }
+   }
+   else //new item is at rank
+   {  if (maxSortUp(data, m, -1)) { maxSortDown(data, m, -1); }
+      if (minSortUp(data, m, 1)) { minSortDown(data, m, 1); }
+   }
+}
+
+template <typename T>
+void rank_filter(T* in_arr, int rank, int arr_len, int win_len, T* out_arr, int mode, T cval, int origin)
+{
+   int i, arr_len_thresh, lim = (win_len - 1) / 2 - origin;
+   int lim2 = arr_len - lim;
+   Mediator* m = MediatorNew(win_len, rank);
+   T* data = (T*)malloc(sizeof(T) * win_len);
+
+   switch (mode)
+   {
+      case REFLECT:
+      for (i=win_len - lim - 1; i > - 1; i--){MediatorInsert(data, m, in_arr[i]);}
+      break;
+      case CONSTANT:
+      for (i=win_len - lim; i > 0; i--){MediatorInsert(data, m, cval);}
+      break;
+      case NEAREST:
+      for (i=win_len - lim; i > 0; i--){MediatorInsert(data, m, in_arr[0]);}
+      break;
+      case MIRROR:
+      for (i=win_len - lim; i > 0; i--){MediatorInsert(data, m, in_arr[i]);}
+      break;
+      case WRAP:
+      for (i=arr_len - lim - 1 - 2 * origin; i < arr_len; i++){MediatorInsert(data, m, in_arr[i]);}
+      break;
+   }
+
+   for (i=0; i < lim; i++){MediatorInsert(data, m, in_arr[i]);}
+   for (i=lim; i < arr_len; i++)
+   {
+      MediatorInsert(data, m, in_arr[i]);
+      out_arr[i - lim] = data[m->heap[0]];
+   }
+   switch (mode)
+   {
+      case REFLECT:
+      arr_len_thresh = arr_len - 1;
+      for (i=0; i < lim; i++)
+      {
+      MediatorInsert(data, m, in_arr[arr_len_thresh - i]);
+      out_arr[lim2 + i] = data[m->heap[0]];
+      }
+      break;
+      case CONSTANT:
+      for (i=0; i < lim; i++)
+      {
+      MediatorInsert(data, m, cval);
+      out_arr[lim2 + i] = data[m->heap[0]];
+      }
+      break;
+      case NEAREST:
+      arr_len_thresh = arr_len - 1;
+      for (i=0; i < lim; i++)
+      {
+      MediatorInsert(data, m, in_arr[arr_len_thresh]);
+      out_arr[lim2 + i] = data[m->heap[0]];
+      }
+      break;
+      case MIRROR:
+      arr_len_thresh = arr_len - 2;
+      for (i=0; i < lim + 1; i++)
+      {
+      MediatorInsert(data, m, in_arr[arr_len_thresh - i]);
+      out_arr[lim2 + i] = data[m->heap[0]];
+      }
+      break;
+      case WRAP:
+      for (i=0; i < win_len; i++){
+      MediatorInsert(data, m, in_arr[i]);
+      out_arr[lim2 + i] = data[m->heap[0]];
+      }
+      break;
+   }
+
+   m->heap -= rank;
+   free(m->heap);
+   m->heap = nullptr;
+   free(m->pos);
+   m->pos = nullptr;
+   free(m);
+   m = nullptr;
+   free(data);
+   data = nullptr;
+}