diff --git a/skimage/filter/rank/_crank8.pyx b/skimage/filter/rank/_crank8.pyx index 84ad9204..4efd2bb0 100644 --- a/skimage/filter/rank/_crank8.pyx +++ b/skimage/filter/rank/_crank8.pyx @@ -264,6 +264,46 @@ cdef inline np.uint8_t kernel_entropy( return < np.uint8_t > e*10 +cdef inline np.uint8_t kernel_otsu( + Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0, + Py_ssize_t s1): + + cdef Py_ssize_t i + cdef Py_ssize_t max_i + cdef float P, mu1, mu2, q1,new_q1, sigma_b, max_sigma_b + cdef float mu = 0. + + # compute local mean + + if pop: + for i in range(256): + mu += histo[i] * i + mu = (mu / pop) + else: + return < np.uint8_t > (0) + + # maximizing the between class variance + max_i = 0 + q1 = histo[0]/pop + m1 = 0. + max_sigma_b = 0. + + for i in range(1,256): + P = histo[i]/pop + new_q1 = q1 + P + if new_q1>0: + mu1 = (q1*mu1 + i*P)/new_q1 + mu2 = (mu-new_q1*mu1)/(1.-new_q1) + sigma_b = new_q1*(1.-new_q1)*(mu1-mu2)**2 + if sigma_b>max_sigma_b: + max_sigma_b = sigma_b + max_i = i + q1 = new_q1 + + if g>max_i: + return < np.uint8_t > 255 + else: + return < np.uint8_t > 0 # ----------------------------------------------------------------- # python wrappers # used only internally @@ -409,3 +449,10 @@ def entropy(np.ndarray[np.uint8_t, ndim=2] image, np.ndarray[np.uint8_t, ndim=2] out=None, char shift_x=0, char shift_y=0): return _core8(kernel_entropy, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0) + +def otsu(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): + return _core8(kernel_otsu, 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 4afa546b..7e4ef405 100644 --- a/skimage/filter/rank/demo/demo_single.py +++ b/skimage/filter/rank/demo/demo_single.py @@ -10,57 +10,23 @@ from skimage.filter import denoise_bilateral if __name__ == '__main__': a8 = data.camera() a16 = data.camera().astype(np.uint16)*4 - selem = disk(10) - f8= rank.percentile_autolevel(a8,selem,p0=.0,p1=1.) - f16= rank.autolevel(a16,selem) - f16p= rank.percentile_autolevel(a16,selem,p0=.0,p1=1.) + p8 = data.page() - den = denoise_bilateral(a8,win_size=10,sigma_range=10,sigma_spatial=2)[:,:,0] - f16b= rank.bilateral_mean(a8.astype(np.uint16),disk(10),s0=10,s1=10) + selem = disk(20) - -# 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) - - selem = disk(3) - ent = rank.entropy(a16,selem) + otsu = rank.otsu(p8,selem) plt.figure() plt.subplot(1,2,1) - plt.imshow(a8) + plt.imshow(p8) + plt.colorbar() plt.subplot(1,2,2) - plt.imshow(ent) + plt.imshow(otsu) plt.colorbar() plt.show() - plt.figure() - plt.subplot(1,2,1) - plt.imshow(den) - plt.subplot(1,2,2) - plt.imshow(f16b) - plt.show() - - print f16==f16p - - plt.figure() - plt.subplot(1,3,1) - plt.imshow(f16) - plt.colorbar() - plt.subplot(1,3,2) - plt.imshow(f16p) - plt.colorbar() - plt.subplot(1,3,3) - plt.imshow(f16p-f16) - plt.colorbar() - plt.show() - - print f16 - print f16p - diff --git a/skimage/filter/rank/rank.pyx b/skimage/filter/rank/rank.pyx index ad66fc65..6a35920b 100644 --- a/skimage/filter/rank/rank.pyx +++ b/skimage/filter/rank/rank.pyx @@ -19,7 +19,7 @@ from skimage.filter.rank import _crank8, _crank16 from skimage.filter.rank.generic import find_bitdepth __all__ = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', 'mean', 'meansubstraction', 'median', 'minimum', - 'modal', 'morph_contr_enh', 'pop', 'threshold', 'tophat','noise_filter','entropy'] + 'modal', 'morph_contr_enh', 'pop', 'threshold', 'tophat','noise_filter','entropy','otsu'] def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y): @@ -655,4 +655,49 @@ def entropy(image, selem, out=None, mask=None, shift_x=False, shift_y=False): """ - return _apply(_crank8.entropy, _crank16.entropy, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) \ No newline at end of file + return _apply(_crank8.entropy, _crank16.entropy, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) + +def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False): + """Returns the image threshold using a the Otsu [otsu]_ locally . + + References + ---------- + + .. [otsu] http://en.wikipedia.org/wiki/Otsu's_method + + Parameters + ---------- + image : ndarray + Image array (uint8 array or uint16). If image is uint16, the algorithm uses max. 12bit histogram, + an exception will be raised if image has a value > 4095 + selem : ndarray + The neighborhood expressed as a 2-D array of 1's and 0's. + out : ndarray + If None, a new array will be allocated. + mask : ndarray (uint8) + 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 : uint8 array or uint16 array (same as input image) + threshold image + + + Examples + -------- + + >>> # Local entropy + >>> from skimage import data + >>> from skimage.filter.rank import otsu + >>> from skimage.morphology import disk + >>> # defining a 8- and a 16-bit test images + >>> a8 = data.camera() + >>> loc_otsu = otsu(a8,disk(5)) + + """ + + return _apply(_crank8.otsu, None, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y) \ No newline at end of file