diff --git a/skimage/filter/rank/__init__.py b/skimage/filter/rank/__init__.py index 04a4b854..9fc2c18e 100644 --- a/skimage/filter/rank/__init__.py +++ b/skimage/filter/rank/__init__.py @@ -1,11 +1,11 @@ from .generic import (autolevel, bottomhat, equalize, gradient, maximum, mean, subtract_mean, median, minimum, modal, enhance_contrast, - pop, threshold, tophat, noise_filter, entropy, otsu) + pop, threshold, tophat, noise_filter, entropy, otsu, sum) from ._percentile import (autolevel_percentile, gradient_percentile, mean_percentile, subtract_mean_percentile, enhance_contrast_percentile, percentile, - pop_percentile, threshold_percentile) -from .bilateral import mean_bilateral, pop_bilateral + pop_percentile, sum_percentile, threshold_percentile) +from .bilateral import mean_bilateral, pop_bilateral, sum_bilateral from skimage._shared.utils import deprecated @@ -51,6 +51,9 @@ __all__ = ['autolevel', 'pop', 'pop_percentile', 'pop_bilateral', + 'sum', + 'sum_bilateral', + 'sum_percentile', 'threshold', 'threshold_percentile', 'tophat', diff --git a/skimage/filter/rank/_percentile.py b/skimage/filter/rank/_percentile.py index ff3b1559..f2b5e7f4 100644 --- a/skimage/filter/rank/_percentile.py +++ b/skimage/filter/rank/_percentile.py @@ -310,6 +310,42 @@ def pop_percentile(image, selem, out=None, mask=None, shift_x=False, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1) +def sum_percentile(image, selem, out=None, mask=None, shift_x=False, + shift_y=False, p0=0, p1=1): + """Return greyscale local sum of an image. + + sum is computed on the given structuring element. Only levels between + percentiles [p0, p1] are used. Result is truncated (8bit or 16bit). + + 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). + p0, p1 : float in [0, ..., 1] + Define the [p0, p1] percentile interval to be considered for computing + the value. + + Returns + ------- + out : ndarray (same dtype as input image) + Output image. + + """ + + return _apply(percentile_cy._sum, + image, selem, out=out, mask=mask, shift_x=shift_x, + shift_y=shift_y, p0=p0, p1=p1) def threshold_percentile(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=0): diff --git a/skimage/filter/rank/bilateral.py b/skimage/filter/rank/bilateral.py index f1b10fec..4509e06b 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,62 @@ 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 + summed. Return greyscale local bilateral sum of an image. + Result is truncated (8bit or 16bit). + + 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 sum_bilateral + >>> # Load test image / cast to uint16 + >>> ima = data.camera().astype(np.uint16) + >>> # bilateral filtering of cameraman image using a flat kernel + >>> bilat_ima = sum_bilateral(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/generic.py b/skimage/filter/rank/generic.py index 34c0787b..703b66b3 100644 --- a/skimage/filter/rank/generic.py +++ b/skimage/filter/rank/generic.py @@ -508,7 +508,6 @@ def pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False): Examples -------- - >>> # Local mean >>> from skimage.morphology import square >>> import skimage.filter.rank as rank >>> ima = 255 * np.array([[0, 0, 0, 0, 0], @@ -528,6 +527,51 @@ def pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False): return _apply(generic_cy._pop, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) +def sum(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Return the sum of pixels inside the neighborhood (truncated to uint8 or uint16). + + 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). + + Returns + ------- + out : ndarray (same dtype as input image) + Output image. + + Examples + -------- + >>> from skimage.morphology import square + >>> import skimage.filter.rank as rank + >>> ima = 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) + >>> rank.sum(ima, square(3)) + 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=uint8) + + """ + + return _apply(generic_cy._sum, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + def threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """Return greyscale local threshold of an image. diff --git a/skimage/filter/rank/generic_cy.pyx b/skimage/filter/rank/generic_cy.pyx index dcf6e361..2a1a8455 100644 --- a/skimage/filter/rank/generic_cy.pyx +++ b/skimage/filter/rank/generic_cy.pyx @@ -221,6 +221,21 @@ cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g, return pop +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 sum = 0 + + if pop: + for i in range(max_bin): + sum += histo[i] * i + return sum + else: + return 0 + cdef inline double _kernel_threshold(Py_ssize_t* histo, double pop, dtype_t g, Py_ssize_t max_bin, Py_ssize_t mid_bin, @@ -455,6 +470,15 @@ def _pop(dtype_t[:, ::1] image, _core(_kernel_pop[dtype_t], image, selem, mask, out, shift_x, shift_y, 0, 0, 0, 0, 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 max_bin): + + _core(_kernel_sum[dtype_t], image, selem, mask, + out, shift_x, shift_y, 0, 0, 0, 0, max_bin) + def _threshold(dtype_t[:, ::1] image, char[:, ::1] selem, diff --git a/skimage/filter/rank/percentile_cy.pyx b/skimage/filter/rank/percentile_cy.pyx index e951a76e..38d04b33 100644 --- a/skimage/filter/rank/percentile_cy.pyx +++ b/skimage/filter/rank/percentile_cy.pyx @@ -90,6 +90,29 @@ cdef inline double _kernel_mean(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, sum, sum_g, n + + if pop: + sum = 0 + sum_g = 0 + n = 0 + for i in range(max_bin): + sum += histo[i] + if (sum >= p0 * pop) and (sum <= p1 * pop): + n += histo[i] + sum_g += histo[i] * i + + if n > 0: + return sum_g + else: + return 0 + else: + return 0 cdef inline double _kernel_subtract_mean(Py_ssize_t* histo, double pop, dtype_t g, @@ -245,6 +268,15 @@ def _mean(dtype_t[:, ::1] image, _core(_kernel_mean[dtype_t], image, selem, mask, out, shift_x, shift_y, p0, p1, 0, 0, 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, double p0, double p1, + Py_ssize_t max_bin): + + _core(_kernel_sum[dtype_t], image, selem, mask, out, + shift_x, shift_y, p0, p1, 0, 0, max_bin) def _subtract_mean(dtype_t[:, ::1] image, char[:, ::1] selem, diff --git a/skimage/filter/rank/tests/test_rank.py b/skimage/filter/rank/tests/test_rank.py index cedeb92d..72f9afe6 100644 --- a/skimage/filter/rank/tests/test_rank.py +++ b/skimage/filter/rank/tests/test_rank.py @@ -498,6 +498,48 @@ def test_percentile_median(): img_max = rank.median(img16, selem=selem) assert_array_equal(img_p0, img_max) +def test_sum(): + # 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(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) + assert_array_equal(r, out8) + rank.sum_percentile(image=image8, selem=elem, out=out8, mask=mask,p0=.0,p1=1.) + 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) + + 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()