add local Otsu threshold

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)

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
-
 
 

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)
+    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)