Add bin splitting to filters.rank.percentile_mean

skimage/filters/rank/_percentile.py

@@ -122,7 +122,7 @@ def gradient_percentile(image, footprint, out=None, mask=None, shift_x=False,
 
 def mean_percentile(image, footprint, out=None, mask=None, shift_x=False,
                     shift_y=False, p0=0, p1=1):
-    """Return local mean of an image.
+    """Return local mean of an image exlcuding outer percentiles in kernel.
 
     Only grayvalues between percentiles [p0, p1] are considered in the filter.
 
@@ -149,8 +149,35 @@ def mean_percentile(image, footprint, out=None, mask=None, shift_x=False,
     out : 2-D array (same dtype as input image)
         Output image.
 
+    Notes
+    -----
+    The algorithm is defined as follows:
+
+    1. Determine the local neighborhood and sort by value.
+    2. Drop values below percentile given by p0 and above percentile given by p1. The
+       interval is inclusive, e.g. ``p0=0`` and ``p0=0.1`` behave the same for a kernel
+       size of 10.
+    3. Calculate the arithmetic mean and drop the remainder (for backwards compability).
+
+    Note that the actual implementation differs for performance reasons.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> import skimage as ski
+    >>> image = np.array([[0, 1, 2],
+    ...                   [3, 4, 5],
+    ...                   [6, 7, 8]], dtype=np.uint8)
+    >>> fp = np.ones((3, 3), dtype=bool)
+    >>> ski.filters.rank.mean_percentile(image, footprint=fp, p0=0.25, p1=0.75)
+    array([[2, 2, 3],
+           [3, 4, 4],
+           [5, 5, 6]], dtype=uint8)
     """
-
+    if not 0. <= p0 < p1 <= 1.:
+        raise ValueError(
+            f"Percentile interval doesn't satisfy 0 <= p0 < p1 <= 1, was [{p0}, {p1}]"
+        )
     return _apply(percentile_cy._mean,
                   image, footprint, out=out, mask=mask, shift_x=shift_x,
                   shift_y=shift_y, p0=p0, p1=p1)

skimage/filters/rank/percentile_cy.pyx

@@ -5,6 +5,7 @@
 
 cimport numpy as cnp
 from .core_cy cimport dtype_t, dtype_t_out, _core, _min, _max
+from libc.math cimport floor, ceil
 cnp.import_array()
 
 cdef inline void _kernel_autolevel(dtype_t_out* out, Py_ssize_t odepth,
@@ -76,25 +77,58 @@ cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth,
                               Py_ssize_t n_bins, Py_ssize_t mid_bin,
                               cnp.float64_t p0, cnp.float64_t p1,
                               Py_ssize_t s0, Py_ssize_t s1) noexcept nogil:
+    """Return local mean of an histogram excluding optional outer percentiles.
+
+    This algorithm uses two counters -- `lower` and `inner``-- to average the
+    appropriate bins of the histogram. First, the number of pixels in each excluded bin
+    are subtracted from `lower` until the percentile is reached for which the arithmetic
+    mean is to be calculated. The same is repeated with the `inner` counter, while
+    summing the value of pixels in each bin which is later divided by the total number
+    of included pixels.
+    """
+    cdef:
+        Py_ssize_t i, lower, inner, denominator
+        cnp.uint64_t total = 0
+
+    # Counter to deplete while summing lower excluded bins in histogram
+    lower = <Py_ssize_t>ceil(pop * p0)
+    # Safely subtract 1 because border should be included in average
+    if 0 < lower:
+        lower -= 1
+    # Counter to deplete while summing inner included bins in histogram
+    inner = <Py_ssize_t>floor(pop * p1)
+    # Safely add 1 because border should be included in average
+    if inner < pop:
+        inner += 1
+    # Inner counter starts after `lower` is depleted, so adjust
+    inner -= lower
+    denominator = inner
+
+    if denominator <= 0:
+        out[0] = <dtype_t_out> 0
+        return  # Return early
+
+    i = 0
+    # Deplete counter `lower` by subtracting lower excluded bins
+    while 0 < lower:
+        lower -= histo[i]
+        i += 1
+    # If `lower` is negative, percentile border is inside bin
+    # so add as many pixel values as `lower` underflowed
+    total += -lower * (i - 1)
+    inner += lower
+    # Deplete counter `inner` by subtracting inner included bins, upper excluded bins
+    # are therefore implicitly excluded
+    while 0 < inner:
+        inner -= histo[i]
+        total += histo[i] * i
+        i += 1
+    # If `inner` is negative, percentile border is inside bin, and we added too
+    # much, so subtract as many pixel values as `inner` underflowed
+    total += inner * (i - 1)
+    # Drop remainder of mean to maintain backwards compatibility
+    out[0] = <dtype_t_out>(total / denominator)
 
-    cdef Py_ssize_t i, sum, mean, n
-
-    if pop:
-        sum = 0
-        mean = 0
-        n = 0
-        for i in range(n_bins):
-            sum += histo[i]
-            if (sum >= p0 * pop) and (sum <= p1 * pop):
-                n += histo[i]
-                mean += histo[i] * i
-
-        if n > 0:
-            out[0] = <dtype_t_out>(mean / n)
-        else:
-            out[0] = <dtype_t_out>0
-    else:
-        out[0] = <dtype_t_out>0
 
 cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth,
                              Py_ssize_t[::1] histo,

skimage/filters/rank/tests/test_rank.py

@@ -887,3 +887,54 @@ def test_input_boolean_dtype(self):
         elem = np.ones((3, 3), dtype=bool)
         with pytest.raises(ValueError):
             rank.maximum(image=image, footprint=elem)
+
+    @pytest.mark.parametrize("p0,p1", [(0.5, 0.5), (0.6, 0.4), (-1, 1), (0, 2)])
+    def test_percentile_mean_invalid(self, p0, p1):
+        image = np.ones((9, 9), dtype=np.uint8)
+        fp = np.ones((3, 3), dtype=bool)
+        with pytest.raises(ValueError, match="Percentile interval doesn't satisfy"):
+            rank.mean_percentile(image, footprint=fp, p0=p0, p1=p1)
+
+    @pytest.mark.parametrize("dtype", [np.uint8, np.uint16])
+    def test_percentile_mean_handcrafted(self, dtype):
+        image = np.array([[0, 11, 22],
+                          [33, 44, 55],
+                          [66, 77, 88]], dtype=dtype)
+        fp = np.ones((3, 3), dtype=bool)
+        result = rank.mean_percentile(image, footprint=fp, p0=0.25, p1=0.75)
+        desired = np.array([[22, 27, 33],
+                            [38, 44, 49],
+                            [55, 60, 66]], dtype=dtype)
+        np.testing.assert_equal(result, desired)
+
+    @pytest.mark.parametrize("dtype", [np.uint8, np.uint16])
+    def test_percentile_mean_pr7096(self, dtype):
+        image = np.array([[0, 0, 0],
+                          [0, 1, 1],
+                          [1, 1, 1]], dtype=dtype)
+        fp = np.ones((3, 3), dtype=bool)
+        result = rank.mean_percentile(image, footprint=fp, p0=0.25, p1=0.75)
+        desired = np.array([[0, 0, 0],
+                            [0, 0, 0],
+                            [0, 1, 1]], dtype=dtype)
+        np.testing.assert_equal(result, desired)
+
+    @pytest.mark.parametrize("p0", (0., .09, .1, .11, .3, .4))
+    @pytest.mark.parametrize("p1", (.6, .7, .89, .9, .91, 1.))
+    def test_percentile_mean_edges2(self, p0, p1):
+        image = np.arange(10, dtype=np.uint16) ** 2
+        image = image.reshape((1, 10))
+
+        lower = int(np.ceil(p0 * image.size))
+        lower = max(0, lower - 1)  # Inclusive border
+        upper = int(np.floor(p1 * image.size))
+        upper = min(image.size, upper + 1)  # Inclusive border
+
+        values = image.ravel()[lower:upper]
+        desired = np.sum(values) // len(values)
+
+        # Kernel should cover entire image for every position
+        footprint = np.ones((1, image.size * 3), dtype=bool)
+        result = rank.mean_percentile(image, footprint, p0=p0, p1=p1)
+
+        assert result[0, 0] == desired