diff --git a/skimage/filter/rank/__init__.py b/skimage/filter/rank/__init__.py index 9a31a565..aad8036b 100644 --- a/skimage/filter/rank/__init__.py +++ b/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 diff --git a/skimage/filter/rank/bilateral.py b/skimage/filter/rank/bilateral.py index f1b10fec..f6e7d228 100644 --- a/skimage/filter/rank/bilateral.py +++ b/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) diff --git a/skimage/filter/rank/bilateral_cy.pyx b/skimage/filter/rank/bilateral_cy.pyx index de2b3e53..25a81766 100644 --- a/skimage/filter/rank/bilateral_cy.pyx +++ b/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) diff --git a/skimage/filter/rank/tests/test_rank.py b/skimage/filter/rank/tests/test_rank.py index b3de866a..72f9afe6 100644 --- a/skimage/filter/rank/tests/test_rank.py +++ b/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()