From a8a5e33425f73cb7e02f4054a95ec65320111a93 Mon Sep 17 00:00:00 2001 From: Olivier Debeir Date: Thu, 8 Nov 2012 10:37:55 +0100 Subject: [PATCH] fix noise_filter --- skimage/filter/rank/_crank8.pyx | 46 +++++--------------- skimage/filter/rank/demo/demo_single.py | 9 +--- skimage/filter/rank/rank.pyx | 20 ++++++++- skimage/filter/rank/tests/test_histo.py | 57 +++++++++++++++++++++++++ 4 files changed, 88 insertions(+), 44 deletions(-) create mode 100644 skimage/filter/rank/tests/test_histo.py diff --git a/skimage/filter/rank/_crank8.pyx b/skimage/filter/rank/_crank8.pyx index a0c3cc35..9ddcb12d 100644 --- a/skimage/filter/rank/_crank8.pyx +++ b/skimage/filter/rank/_crank8.pyx @@ -230,21 +230,26 @@ cdef inline np.uint8_t kernel_noise_filter( cdef Py_ssize_t i cdef Py_ssize_t min_i - for i in range(255, g, -1): + # early stop if at least one pixel of the neighborhood has the same g + if histo[g]>0: + return < np.uint8_t > 0 + + for i in range(g, -1, -1): if histo[i]: break - min_i = i-g - for i in range(0, g): + min_i = g-i + for i in range(g, 256): if histo[i]: break - if g-i < min_i: - return < np.uint8_t > (g-i) + if i-g < min_i: + return < np.uint8_t > (i-g) else: return < np.uint8_t > min_i # ----------------------------------------------------------------- # python wrappers +# used only internally # ----------------------------------------------------------------- @@ -253,8 +258,6 @@ def autolevel(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """bottom hat - """ return _core8( kernel_autolevel, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) @@ -265,8 +268,6 @@ def bottomhat(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """bottom hat - """ return _core8( kernel_bottomhat, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) @@ -277,8 +278,6 @@ def equalize(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """local egalisation of the gray level - """ return _core8( kernel_equalize, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) @@ -289,8 +288,6 @@ def gradient(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """local maximum - local minimum gray level - """ return _core8( kernel_gradient, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) @@ -301,8 +298,6 @@ def maximum(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """local maximum gray level - """ return _core8(kernel_maximum, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) @@ -311,8 +306,6 @@ def mean(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """average gray level (clipped on uint8) - """ return _core8(kernel_mean, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) @@ -321,8 +314,6 @@ def meansubstraction(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """(g - average gray level)/2+127 (clipped on uint8) - """ return _core8( kernel_meansubstraction, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) @@ -333,8 +324,6 @@ def median(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """local median - """ return _core8(kernel_median, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) @@ -343,8 +332,6 @@ def minimum(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """local minimum gray level - """ return _core8(kernel_minimum, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) @@ -353,8 +340,6 @@ def morph_contr_enh(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """morphological contrast enhancement - """ return _core8( kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) @@ -365,8 +350,6 @@ def modal(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """local mode - """ return _core8(kernel_modal, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) @@ -375,8 +358,6 @@ def pop(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """returns the number of actual pixels of the structuring element inside the mask - """ return _core8(kernel_pop, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) @@ -385,8 +366,6 @@ def threshold(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """returns 255 if gray level higher than local mean, 0 else - """ return _core8( kernel_threshold, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) @@ -397,15 +376,12 @@ def tophat(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """top hat - """ return _core8(kernel_tophat, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) + def noise_filter(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] selem, np.ndarray[np.uint8_t, ndim=2] mask=None, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): - """top hat - """ return _core8(kernel_noise_filter, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) diff --git a/skimage/filter/rank/demo/demo_single.py b/skimage/filter/rank/demo/demo_single.py index 3b7eaa95..e3fe9c1e 100644 --- a/skimage/filter/rank/demo/demo_single.py +++ b/skimage/filter/rank/demo/demo_single.py @@ -20,16 +20,11 @@ if __name__ == '__main__': f16b= rank.bilateral_mean(a8.astype(np.uint16),disk(10),s0=10,s1=10) - selem = np.ones((3,3)) - selem[1,1] = 0 - radius = 3 - selem = disk(radius) - selem[radius,radius] = 0 - print selem +# selem = np.ones((3,3),dtype=np.uint8) noise = rank.noise_filter(a8,selem) plt.imsave('noise.png',noise,cmap=plt.cm.gray) plt.imsave('cam.png',a8,cmap=plt.cm.gray) - print noise + plt.figure() plt.subplot(1,2,1) diff --git a/skimage/filter/rank/rank.pyx b/skimage/filter/rank/rank.pyx index 40ceecad..74875d1c 100644 --- a/skimage/filter/rank/rank.pyx +++ b/skimage/filter/rank/rank.pyx @@ -27,8 +27,12 @@ def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y): if mask is not None: mask = img_as_ubyte(mask) if image.dtype == np.uint8: + if func8 is None: + raise TypeError("uint8 image not supported for this filter") return func8(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, out=out) elif image.dtype == np.uint16: + if func16 is None: + raise TypeError("uint16 image not supported for this filter") bitdepth = find_bitdepth(image) if bitdepth > 11: raise ValueError("only uint16 <4096 image (12bit) supported!") @@ -581,7 +585,7 @@ def tophat(image, selem, out=None, mask=None, shift_x=False, shift_y=False): return _apply(_crank8.tophat, _crank16.tophat, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) def noise_filter(image, selem, out=None, mask=None, shift_x=False, shift_y=False): - """Return minimum absolute diffirence between a pixel and its neighborhood + """Returns the noise feature as described in [1]_ Parameters ---------- @@ -600,11 +604,23 @@ def noise_filter(image, selem, out=None, mask=None, shift_x=False, shift_y=False Offset added to the structuring element center point. Shift is bounded to the structuring element sizes (center must be inside the given structuring element). + Reference + ---------- + + .. [1] N. Hashimoto et al. Referenceless image quality evaluation for whole slide imaging. J Pathol Inform 2012;3:9. + + Returns ------- out : uint8 array or uint16 array (same as input image) The image noise . """ + # ensure that the central pixel in the structuring element is empty + centre_r = int(selem.shape[0] / 2) + shift_y + centre_c = int(selem.shape[1] / 2) + shift_x + # make a local copy + selem_cpy = selem.copy() + selem_cpy[centre_r,centre_c] = 0 - return _apply(_crank8.noise_filter, None, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) \ No newline at end of file + return _apply(_crank8.noise_filter, None, image, selem_cpy, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) \ No newline at end of file diff --git a/skimage/filter/rank/tests/test_histo.py b/skimage/filter/rank/tests/test_histo.py new file mode 100644 index 00000000..b6a7bbca --- /dev/null +++ b/skimage/filter/rank/tests/test_histo.py @@ -0,0 +1,57 @@ +import sys +print sys.path +import skimage +print skimage + +import unittest + +import numpy as np +from skimage.filter import rank + +from skimage import data +from skimage.morphology import cmorph,disk +from skimage.filter.rank import _crank8, _crank16 +from skimage.filter.rank import _crank16_percentiles + + +class TestSequenceFunctions(unittest.TestCase): + + def setUp(self): + pass + + def test_trivial_selem(self): + # check that min, max and mean returns identity if structuring element contains only central pixel + + a = np.zeros((5,5),dtype='uint8') + a[2,2] = 255 + a[2,3] = 128 + a[1,2] = 16 + elem = np.asarray([[0,0,0],[0,1,0],[0,0,0]],dtype='uint8') + f = _crank8.mean(image=a,selem = elem,shift_x=0,shift_y=0) + np.testing.assert_array_equal(a,f) + f = _crank8.minimum(image=a,selem = elem,shift_x=0,shift_y=0) + np.testing.assert_array_equal(a,f) + f = _crank8.maximum(image=a,selem = elem,shift_x=0,shift_y=0) + np.testing.assert_array_equal(a,f) + + def test_smallest_selem(self): + # check that min, max and mean returns identity if structuring element contains only central pixel + + a = np.zeros((5,5),dtype='uint8') + a[2,2] = 255 + a[2,3] = 128 + a[1,2] = 16 + elem = np.asarray([[1]],dtype='uint8') + f = _crank8.mean(image=a,selem = elem,shift_x=0,shift_y=0) + np.testing.assert_array_equal(a,f) + f = _crank8.minimum(image=a,selem = elem,shift_x=0,shift_y=0) + np.testing.assert_array_equal(a,f) + f = _crank8.maximum(image=a,selem = elem,shift_x=0,shift_y=0) + np.testing.assert_array_equal(a,f) + + + +if __name__ == '__main__': + + suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions) + unittest.TextTestRunner(verbosity=2).run(suite)