From 4bc1c587ede3a5b74ef5525bd884d84914742483 Mon Sep 17 00:00:00 2001 From: Guillaume Lemaitre Date: Wed, 28 Oct 2015 20:29:48 +0100 Subject: [PATCH] Add local geometric mean filter Add the filter in cython with a specific kernnel. The implementatio is based on the log-average method. The test were added the npz. --- skimage/data/rank_filter_tests.npz | Bin 33945 -> 34762 bytes skimage/filters/rank/__init__.py | 7 ++-- skimage/filters/rank/generic.py | 42 ++++++++++++++++++++++-- skimage/filters/rank/generic_cy.pyx | 31 ++++++++++++++++- skimage/filters/rank/tests/test_rank.py | 35 ++++++++++++++++++-- 5 files changed, 107 insertions(+), 8 deletions(-) diff --git a/skimage/data/rank_filter_tests.npz b/skimage/data/rank_filter_tests.npz index 3142fcc15f05b7a5e38ac7dcbcbba15eb0a0ec94..d5169bc575222b07a04ad1cf8d538ca6a44335b1 100644 GIT binary patch delta 1916 zcmZvceM}Q~7{|||j;SMaGGQeXRAfto1;xr1VWis}HN~Yoq>NGJi%XDU2wq%LMZR*^VWcy=@iMn3zF4~2=T>H7F-{-Br z=l46Xn{s22lA%_lEzAV43p!ROIhpa8vJivm7nHvUnp<>_pdZrj6Le|mZv-9BAZh?HRa20RZfIE0MwudRstfyVHtJISj}V&r~8VAe9 z83J0v<`R)crD`_{UGgs2U%f*WD6!=i3(L$=U2BxvfqZztf%Y66=$t(o*~@SLVj`ymx#WBAAmq%_S&ie@#p|-DQ)F0k*TC6 z;$azUy*s^D9<=)t$bHanc7#D^bGS|DXwiBm;`@bbtx3;#G#^;2u&P@EgcaPMsQON{$ zXVB8sXO=q0gf4sh+{E?Z>BBLbT$V3(m356=4%pqnahK1d!JB6ZRJGfk%T(w<_K+GlFcm$e!O6f(W<+pX75 z&0gS&`OU*J5!4K8)k=#_o?EzCapsG&u@ft~1dS!Q; zTR@A4?rCjTD;7{P(_S7gOCEYAep&$l*t3#!PWo@7qHtC=N<)_68F>hmQcp(>$htf% z@sb~dx5S>wkTOM8iA>6Fy(4n8bT^FQj}jy0Ra6jMopHgz1ew^ zl)afEn^K^cBde04TgM(hBiCqYifBZIK99VEmzoo$Qbrrm8)hncCq}H9C1D!|(JPIs zUfMp3QaM8o0B8Z5ob%7bqsO9sWKpSIcDT?N(cE;bLDRW(5t=UdlTrztn3MJ0sJ#%OqpbBTXV9smF@i7&vLnMH&F1hS)gT{rT*Wnuv_HveaO17*EueFSA)<#+~V?dN_1Wi97} ziA@xI2Ni1;c>!f5NO6cj%*)aH1Qkm%d=F)XnmmWH9L!;Ms9C{WAZ`m|aX7%t`{D#+ zJ#dA&?}Eoms408BUqe}I{Qe3<>|LAm87j6U^*fX`KjSHs)sYP|r6LceF24}ws;H7* zP<0XIkD)C4Dwtl4S}1F>V1o_VvdMuB))2;o215wrXoD$)@vp%c!Z2-wisUugLPXX# znnM^`O;SvW;*;~FBqql-@j%4#fMU7Yle?N!6+rRAb*IPMr48k6-V z7Z{4^1b8zti7+E#)ZjkzbC?Dmh*~h&&=_VoS+H3OVvByWl)R$}14D6Xu3lb2CAuMx zBqBCXe%C0*q;5XBp;>_g}HBG%J8zm#o-)5M~)Kie;R(lUk%9cFDI$ zF@3e4Yz^Z_0r?M|CYQCyD}l^aOk=KxnZbu*#zj{KhRG8frNCB&ov+>qQzn3-Y@f&E zk1cXwMHd#o{d+^0fdPcMQ53E9o~++032~55s}$3^_{qDIq$Ufr^MEz}dsjCfX1WZD z#-&M<*R^UvG~Q{IQd*G8z>u4ol9`x?E#O|gb)7mnu0@QgGh?zrn>;w+Vs~-4P7Y`l zW2(%aoZqGv1@bgo$$n{=eW0) 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)