add:rank.sum_bilateral

skimage/filter/rank/__init__.py

@@ -5,7 +5,7 @@ from ._percentile import (autolevel_percentile, gradient_percentile,
                           mean_percentile, subtract_mean_percentile,
                           enhance_contrast_percentile, percentile,
                           pop_percentile,sum_percentile, threshold_percentile)
-from .bilateral import mean_bilateral, pop_bilateral
+from .bilateral import mean_bilateral, pop_bilateral, sum_bilateral
 
 from skimage._shared.utils import deprecated
 

skimage/filter/rank/bilateral.py

@@ -30,7 +30,7 @@ from . import bilateral_cy
 from .generic import _handle_input
 
 
-__all__ = ['mean_bilateral', 'pop_bilateral']
+__all__ = ['mean_bilateral', 'pop_bilateral', 'sum_bilateral']
 
 
 def _apply(func, image, selem, out, mask, shift_x, shift_y, s0, s1,
@@ -155,3 +155,61 @@ def pop_bilateral(image, selem, out=None, mask=None, shift_x=False,
 
     return _apply(bilateral_cy._pop, image, selem, out=out,
                   mask=mask, shift_x=shift_x, shift_y=shift_y, s0=s0, s1=s1)
+
+def sum_bilateral(image, selem, out=None, mask=None, shift_x=False,
+                   shift_y=False, s0=10, s1=10):
+    """Apply a flat kernel bilateral filter.
+
+    This is an edge-preserving and noise reducing denoising filter. It averages
+    pixels based on their spatial closeness and radiometric similarity.
+
+    Spatial closeness is measured by considering only the local pixel
+    neighborhood given by a structuring element (selem).
+
+    Radiometric similarity is defined by the greylevel interval [g-s0, g+s1]
+    where g is the current pixel greylevel. Only pixels belonging to the
+    structuring element AND having a greylevel inside this interval are
+    averaged. Return greyscale local bilateral_sum of an image.
+
+    Parameters
+    ----------
+    image : ndarray (uint8, uint16)
+        Image array.
+    selem : ndarray
+        The neighborhood expressed as a 2-D array of 1's and 0's.
+    out : ndarray (same dtype as input)
+        If None, a new array will be allocated.
+    mask : ndarray
+        Mask array that defines (>0) area of the image included in the local
+        neighborhood. If None, the complete image is used (default).
+    shift_x, shift_y : int
+        Offset added to the structuring element center point. Shift is bounded
+        to the structuring element sizes (center must be inside the given
+        structuring element).
+    s0, s1 : int
+        Define the [s0, s1] interval around the greyvalue of the center pixel
+        to be considered for computing the value.
+
+    Returns
+    -------
+    out : ndarray (same dtype as input image)
+        Output image.
+
+    See also
+    --------
+    skimage.filter.denoise_bilateral for a gaussian bilateral filter.
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> from skimage.morphology import disk
+    >>> from skimage.filter.rank import bilateral_sum
+    >>> # Load test image / cast to uint16
+    >>> ima = data.camera().astype(np.uint16)
+    >>> # bilateral filtering of cameraman image using a flat kernel
+    >>> bilat_ima = bilateral_sum(ima, disk(20), s0=10,s1=10)
+
+    """
+
+    return _apply(bilateral_cy._sum, image, selem, out=out,
+                  mask=mask, shift_x=shift_x, shift_y=shift_y, s0=s0, s1=s1)

skimage/filter/rank/bilateral_cy.pyx

@@ -47,6 +47,27 @@ cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g,
     else:
         return 0
 
+cdef inline double _kernel_sum(Py_ssize_t* histo, double pop, dtype_t g,
+                                Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                double p0, double p1,
+                                Py_ssize_t s0, Py_ssize_t s1):
+
+    cdef Py_ssize_t i
+    cdef Py_ssize_t bilat_pop = 0
+    cdef Py_ssize_t sum = 0
+
+    if pop:
+        for i in range(max_bin):
+            if (g > (i - s0)) and (g < (i + s1)):
+                bilat_pop += histo[i]
+                sum += histo[i] * i
+        if bilat_pop:
+            return sum
+        else:
+            return 0
+    else:
+        return 0
+
 
 def _mean(dtype_t[:, ::1] image,
           char[:, ::1] selem,
@@ -68,3 +89,13 @@ def _pop(dtype_t[:, ::1] image,
 
     _core(_kernel_pop[dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, s0, s1, max_bin)
+
+def _sum(dtype_t[:, ::1] image,
+         char[:, ::1] selem,
+         char[:, ::1] mask,
+         dtype_t_out[:, ::1] out,
+         char shift_x, char shift_y, Py_ssize_t s0, Py_ssize_t s1,
+         Py_ssize_t max_bin):
+
+    _core(_kernel_sum[dtype_t], image, selem, mask, out,
+          shift_x, shift_y, 0, 0, s0, s1, max_bin)

skimage/filter/rank/tests/test_rank.py

@@ -516,55 +516,30 @@ def test_sum():
     out16 = np.empty_like(image16)
     mask = np.ones(image8.shape, dtype=np.uint8)
 
+    r =  np.array([[1, 2, 3, 2, 1],
+           [2, 4, 6, 4, 2],
+           [3, 6, 9, 6, 3],
+           [2, 4, 6, 4, 2],
+           [1, 2, 3, 2, 1]], dtype=np.uint8)
     rank.sum(image=image8, selem=elem, out=out8, mask=mask)
-    r =  np.array([[1, 2, 3, 2, 1],
-           [2, 4, 6, 4, 2],
-           [3, 6, 9, 6, 3],
-           [2, 4, 6, 4, 2],
-           [1, 2, 3, 2, 1]], dtype=np.uint8)
     assert_array_equal(r, out8)
-
-    rank.sum(image=image16, selem=elem, out=out16, mask=mask)
-    r = 400* np.array([[1, 2, 3, 2, 1],
-           [2, 4, 6, 4, 2],
-           [3, 6, 9, 6, 3],
-           [2, 4, 6, 4, 2],
-           [1, 2, 3, 2, 1]], dtype=np.uint16)
-    assert_array_equal(r, out16)
-
-def test_sum_percentile():
-    # check the number of valid pixels in the neighborhood
-
-    image8 = np.array([[0, 0, 0, 0, 0],
-                       [0, 1, 1, 1, 0],
-                       [0, 1, 1, 1, 0],
-                       [0, 1, 1, 1, 0],
-                       [0, 0, 0, 0, 0]], dtype=np.uint8)
-    image16 = 400*np.array([[0, 0, 0, 0, 0],
-                       [0, 1, 1, 1, 0],
-                       [0, 1, 1, 1, 0],
-                       [0, 1, 1, 1, 0],
-                       [0, 0, 0, 0, 0]], dtype=np.uint16)
-    elem = np.ones((3, 3), dtype=np.uint8)
-    out8 = np.empty_like(image8)
-    out16 = np.empty_like(image16)
-    mask = np.ones(image16.shape, dtype=np.uint8)
-
     rank.sum_percentile(image=image8, selem=elem, out=out8, mask=mask,p0=.0,p1=1.)
-    r =  np.array([[1, 2, 3, 2, 1],
-           [2, 4, 6, 4, 2],
-           [3, 6, 9, 6, 3],
-           [2, 4, 6, 4, 2],
-           [1, 2, 3, 2, 1]], dtype=np.uint8)
+    assert_array_equal(r, out8)
+    rank.sum_bilateral(image=image8, selem=elem, out=out8, mask=mask,s0=255,s1=255)
     assert_array_equal(r, out8)
 
-    rank.sum_percentile(image=image16, selem=elem, out=out16, mask=mask,p0=.0,p1=1.)
     r = 400* np.array([[1, 2, 3, 2, 1],
            [2, 4, 6, 4, 2],
            [3, 6, 9, 6, 3],
            [2, 4, 6, 4, 2],
            [1, 2, 3, 2, 1]], dtype=np.uint16)
+    rank.sum(image=image16, selem=elem, out=out16, mask=mask)
     assert_array_equal(r, out16)
+    rank.sum_percentile(image=image16, selem=elem, out=out16, mask=mask,p0=.0,p1=1.)
+    assert_array_equal(r, out16)
+    rank.sum_bilateral(image=image16, selem=elem, out=out16, mask=mask,s0=1000,s1=1000)
+    assert_array_equal(r, out16)
+
 
 if __name__ == "__main__":
     run_module_suite()