diff --git a/skimage/data/rank_filter_tests.npz b/skimage/data/rank_filter_tests.npz index 3142fcc1..d5169bc5 100644 Binary files a/skimage/data/rank_filter_tests.npz and b/skimage/data/rank_filter_tests.npz differ diff --git a/skimage/filters/rank/__init__.py b/skimage/filters/rank/__init__.py index 37b5f9cd..948adee6 100644 --- a/skimage/filters/rank/__init__.py +++ b/skimage/filters/rank/__init__.py @@ -1,7 +1,7 @@ from .generic import (autolevel, bottomhat, equalize, gradient, maximum, mean, - subtract_mean, median, minimum, modal, enhance_contrast, - pop, threshold, tophat, noise_filter, entropy, otsu, - sum, windowed_histogram) + geometric_mean, subtract_mean, median, minimum, modal, + enhance_contrast, pop, threshold, tophat, noise_filter, + entropy, otsu, sum, windowed_histogram) from ._percentile import (autolevel_percentile, gradient_percentile, mean_percentile, subtract_mean_percentile, enhance_contrast_percentile, percentile, @@ -17,6 +17,7 @@ __all__ = ['autolevel', 'gradient_percentile', 'maximum', 'mean', + 'geometric_mean', 'mean_percentile', 'mean_bilateral', 'subtract_mean', diff --git a/skimage/filters/rank/generic.py b/skimage/filters/rank/generic.py index fd718b98..8138fe21 100644 --- a/skimage/filters/rank/generic.py +++ b/skimage/filters/rank/generic.py @@ -25,8 +25,9 @@ from . import generic_cy __all__ = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', 'mean', - 'subtract_mean', 'median', 'minimum', 'modal', 'enhance_contrast', - 'pop', 'threshold', 'tophat', 'noise_filter', 'entropy', 'otsu'] + 'geometric_mean', 'subtract_mean', 'median', 'minimum', 'modal', + 'enhance_contrast', 'pop', 'threshold', 'tophat', 'noise_filter', + 'entropy', 'otsu'] def _handle_input(image, selem, out, mask, out_dtype=None, pixel_size=1): @@ -342,6 +343,43 @@ def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False): return _apply_scalar_per_pixel(generic_cy._mean, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) +def geometric_mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Return local geometric mean of an image. + + Parameters + ---------- + image : 2-D array (uint8, uint16) + Input image. + selem : 2-D array + The neighborhood expressed as a 2-D array of 1's and 0's. + out : 2-D array (same dtype as input) + If None, a new array is 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 : 2-D array (same dtype as input image) + Output image. + + Examples + -------- + >>> from skimage import data + >>> from skimage.morphology import disk + >>> from skimage.filters.rank import mean + >>> img = data.camera() + >>> avg = geomtric_mean(img, disk(5)) + + """ + + return _apply_scalar_per_pixel(generic_cy._geometric_mean, image, selem, out=out, + mask=mask, shift_x=shift_x, shift_y=shift_y) + def subtract_mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False): diff --git a/skimage/filters/rank/generic_cy.pyx b/skimage/filters/rank/generic_cy.pyx index 1356f369..9ff17dd5 100644 --- a/skimage/filters/rank/generic_cy.pyx +++ b/skimage/filters/rank/generic_cy.pyx @@ -4,7 +4,7 @@ #cython: wraparound=False cimport numpy as cnp -from libc.math cimport log +from libc.math cimport log, exp, round from .core_cy cimport dtype_t, dtype_t_out, _core @@ -133,6 +133,25 @@ cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth, out[0] = 0 +cdef inline void _kernel_geometric_mean(dtype_t_out* out, Py_ssize_t odepth, + 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) nogil: + + cdef Py_ssize_t i + cdef double mean = 0. + + if pop: + for i in range(max_bin): + if histo[i]: + mean += (histo[i] * log(i+1)) + out[0] = round(exp(mean / pop)-1) + else: + out[0] = 0 + + cdef inline void _kernel_subtract_mean(dtype_t_out* out, Py_ssize_t odepth, Py_ssize_t* histo, double pop, dtype_t g, @@ -467,6 +486,16 @@ def _mean(dtype_t[:, ::1] image, shift_x, shift_y, 0, 0, 0, 0, max_bin) +def _geometric_mean(dtype_t[:, ::1] image, + char[:, ::1] selem, + char[:, ::1] mask, + dtype_t_out[:, :, ::1] out, + signed char shift_x, signed char shift_y, Py_ssize_t max_bin): + + _core(_kernel_geometric_mean[dtype_t_out, dtype_t], image, selem, mask, out, + shift_x, shift_y, 0, 0, 0, 0, max_bin) + + def _subtract_mean(dtype_t[:, ::1] image, char[:, ::1] selem, char[:, ::1] mask, diff --git a/skimage/filters/rank/tests/test_rank.py b/skimage/filters/rank/tests/test_rank.py index 58e338f9..3fcc79ac 100644 --- a/skimage/filters/rank/tests/test_rank.py +++ b/skimage/filters/rank/tests/test_rank.py @@ -39,6 +39,8 @@ def check_all(): rank.maximum(image, selem)) assert_equal(refs["mean"], rank.mean(image, selem)) + assert_equal(refs["geometric_mean"], + rank.geometric_mean(image, selem)), assert_equal(refs["mean_percentile"], rank.mean_percentile(image, selem)) assert_equal(refs["mean_bilateral"], @@ -102,6 +104,13 @@ def test_random_sizes(): rank.mean(image=image8, selem=elem, mask=mask, out=out8, shift_x=+1, shift_y=+1) assert_equal(image8.shape, out8.shape) + + rank.geometric_mean(image=image8, selem=elem, mask=mask, out=out8, + shift_x=0, shift_y=0) + assert_equal(image8.shape, out8.shape) + rank.geometric_mean(image=image8, selem=elem, mask=mask, out=out8, + shift_x=+1, shift_y=+1) + assert_equal(image8.shape, out8.shape) image16 = np.ones((m, n), dtype=np.uint16) out16 = np.empty_like(image8, dtype=np.uint16) @@ -112,6 +121,13 @@ def test_random_sizes(): shift_x=+1, shift_y=+1) assert_equal(image16.shape, out16.shape) + rank.geometric_mean(image=image16, selem=elem, mask=mask, out=out16, + shift_x=0, shift_y=0) + assert_equal(image16.shape, out16.shape) + rank.geometric_mean(image=image16, selem=elem, mask=mask, out=out16, + shift_x=+1, shift_y=+1) + assert_equal(image16.shape, out16.shape) + rank.mean_percentile(image=image16, mask=mask, out=out16, selem=elem, shift_x=0, shift_y=0, p0=.1, p1=.9) assert_equal(image16.shape, out16.shape) @@ -292,8 +308,8 @@ def test_compare_8bit_unsigned_vs_signed(): assert_equal(image_u, img_as_ubyte(image_s)) methods = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', - 'mean', 'subtract_mean', 'median', 'minimum', 'modal', - 'enhance_contrast', 'pop', 'threshold', 'tophat'] + 'mean', 'geometric_mean', 'subtract_mean', 'median', 'minimum', + 'modal', 'enhance_contrast', 'pop', 'threshold', 'tophat'] for method in methods: func = getattr(rank, method) @@ -338,6 +354,9 @@ def test_trivial_selem8(): rank.mean(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) + rank.geometric_mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_equal(image, out) rank.minimum(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) @@ -361,6 +380,9 @@ def test_trivial_selem16(): rank.mean(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) + rank.geometric_mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_equal(image, out) rank.minimum(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) @@ -407,6 +429,9 @@ def test_smallest_selem16(): rank.mean(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) + rank.geometric_mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_equal(image, out) rank.minimum(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) @@ -432,6 +457,9 @@ def test_empty_selem(): rank.mean(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(res, out) + rank.geometric_mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_equal(res, out) rank.minimum(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(res, out) @@ -514,6 +542,9 @@ def test_selem_dtypes(): rank.mean(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out) + rank.geometric_mean(image=image, selem=elem, out=out, mask=mask, + shift_x=0, shift_y=0) + assert_equal(image, out) rank.mean_percentile(image=image, selem=elem, out=out, mask=mask, shift_x=0, shift_y=0) assert_equal(image, out)