diff --git a/.gitignore b/.gitignore
index 1bd4352c..2623ea74 100644
--- a/.gitignore
+++ b/.gitignore
@@ -22,4 +22,5 @@ doc/source/auto_examples/images/plot_*.png
 doc/source/auto_examples/images/thumb
 doc/source/auto_examples/applications/
 doc/source/_static/random.js
-
+.idea/
+*.log
diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index c73b14e2..8ba62e0b 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -117,3 +117,6 @@
 
 - Petter Strandmark
   Perimeter calculation in regionprops.
+
+- Olivier Debeir
+  Rank filters (8- and 16-bits) using sliding window.
\ No newline at end of file
diff --git a/doc/__init__.py b/doc/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/doc/examples/applications/plot_rank_filters.py b/doc/examples/applications/plot_rank_filters.py
new file mode 100644
index 00000000..d52595f2
--- /dev/null
+++ b/doc/examples/applications/plot_rank_filters.py
@@ -0,0 +1,585 @@
+"""
+===============================================================
+Rank filters
+===============================================================
+
+Rank filters are non-linear filters using the local grey levels ordering to compute the filtered value.
+This ensemble of filters share a common base: the local grey-level histogram extraction computed on
+the neighborhood of a pixel (defined by a 2D structuring element).
+If the filtered value is taken as the middle value of the histogram, we get the classical median filter.
+
+Rank filters can be used for several purposes such as:
+
+* image quality enhancement
+  e.g. image smoothing, sharpening
+
+* image pre-processing
+  e.g. noise reduction, contrast enhancement
+
+* feature extraction
+  e.g. border detection, isolated point detection
+
+* post-processing
+  e.g. small object removal, object grouping, contour smoothing
+
+Some well known filters are specific cases of rank filters [1]_ e.g. morphological dilation, morphological erosion,
+median filters.
+
+The different implementation availables in ``skimage`` are compared.
+
+In this example, we will see how to filter a grey level image using some of the linear and non-linear filters
+availables in skimage.  We use the ``camera`` image from ``skimage.data``.
+
+.. [1] Pierre Soille, On morphological operators based on rank filters, Pattern Recognition 35 (2002) 527-535.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage import data
+
+ima = data.camera()
+hist = np.histogram(ima, bins=np.arange(0, 256))
+
+plt.figure(figsize=(8, 3))
+plt.subplot(121)
+plt.imshow(ima, cmap=plt.cm.gray, interpolation='nearest')
+plt.axis('off')
+plt.subplot(122)
+plt.plot(hist[1][:-1], hist[0], lw=2)
+plt.title('histogram of grey values')
+
+"""
+.. image:: PLOT2RST.current_figure
+
+Noise removal
+==============
+
+some noise is added to the image, 1% of pixels are randomly set to 255, %1% are randomly set to 0.
+The **median** filter is applied to remove the noise.
+
+.. note:: there is different implementations of median filter : ``skimage.filter.median_filter``and
+`skimage.filter.rank.median``
+"""
+
+noise = np.random.random(ima.shape)
+nima = data.camera()
+nima[noise>.99] = 255
+nima[noise<.01] = 0
+
+from skimage.filter.rank import median
+from skimage.morphology import disk
+
+fig = plt.figure(figsize=[10,7])
+
+lo = median(nima,disk(1))
+hi = median(nima,disk(5))
+ext = median(nima,disk(20))
+plt.subplot(2,2,1)
+plt.imshow(nima,cmap=plt.cm.gray,vmin=0,vmax=255)
+plt.xlabel('noised image')
+plt.subplot(2,2,2)
+plt.imshow(lo,cmap=plt.cm.gray,vmin=0,vmax=255)
+plt.xlabel('median $r=1$')
+plt.subplot(2,2,3)
+plt.imshow(hi,cmap=plt.cm.gray,vmin=0,vmax=255)
+plt.xlabel('median $r=5$')
+plt.subplot(2,2,4)
+plt.imshow(ext,cmap=plt.cm.gray,vmin=0,vmax=255)
+plt.xlabel('median $r=20$')
+
+"""
+.. image:: PLOT2RST.current_figure
+
+The added noise is efficiently removed, as the image defaults are small (1 pixel wide), a small filter radius is
+sufficient. As the radius is increasing, objects with a bigger size are filtered too such as the camera tripod.
+Median filter is commonly used for noise removal because borders are preserved.
+
+Image smoothing
+================
+
+The example hereunder shows how a local **mean** smooth the cameraman image.
+
+"""
+from skimage.filter.rank import mean
+
+fig = plt.figure(figsize=[10,7])
+
+loc_mean = mean(nima,disk(10))
+plt.subplot(1,2,1)
+plt.imshow(ima,cmap=plt.cm.gray,vmin=0,vmax=255)
+plt.xlabel('original')
+plt.subplot(1,2,2)
+plt.imshow(loc_mean,cmap=plt.cm.gray,vmin=0,vmax=255)
+plt.xlabel('local mean $r=10$')
+
+"""
+.. image:: PLOT2RST.current_figure
+
+One may be interested in smoothing an image while preserving important borders (median filters already achieved this),
+here we use the **bilateral** filter that restrict the local neighborhood to pixel having a grey level similar to the
+central one.
+
+.. note:: a different implementation is available for color images in ``skimage.filter.denoise_bilateral``.
+
+"""
+
+from skimage.filter.rank import bilateral_mean
+
+ima = data.camera()
+selem = disk(10)
+
+bilat = bilateral_mean(ima.astype(np.uint16),disk(20),s0=10,s1=10)
+
+# display results
+fig = plt.figure(figsize=[10,7])
+plt.subplot(2,2,1)
+plt.imshow(ima,cmap=plt.cm.gray)
+plt.xlabel('original')
+plt.subplot(2,2,3)
+plt.imshow(bilat,cmap=plt.cm.gray)
+plt.xlabel('bilateral mean')
+plt.subplot(2,2,2)
+plt.imshow(ima[200:350,350:450],cmap=plt.cm.gray)
+plt.subplot(2,2,4)
+plt.imshow(bilat[200:350,350:450],cmap=plt.cm.gray)
+
+"""
+.. image:: PLOT2RST.current_figure
+
+One can see that the large continuous part of the image (e.g.sky) are smoothed whereas other details are preserved.
+
+
+Contrast enhancement
+====================
+
+We compare here how the global histogram equalization is applied locally.
+
+The equalized image [2]_ has a roughly linear cumulative distribution function for each pixel neighborhood.
+The local version [3]_ of the histogram equalization emphasized every local graylevel variations.
+
+.. [2] http://en.wikipedia.org/wiki/Histogram_equalization
+.. [3] http://en.wikipedia.org/wiki/Adaptive_histogram_equalization
+
+"""
+
+from skimage import exposure
+from skimage.filter import rank
+
+ima = data.camera()
+# equalize globally and locally
+glob = exposure.equalize(ima)*255
+loc = rank.equalize(ima,disk(20))
+
+# extract histogram for each image
+hist = np.histogram(ima, bins=np.arange(0, 256))
+glob_hist = np.histogram(glob, bins=np.arange(0, 256))
+loc_hist = np.histogram(loc, bins=np.arange(0, 256))
+
+plt.figure(figsize=(10, 10))
+plt.subplot(321)
+plt.imshow(ima, cmap=plt.cm.gray, interpolation='nearest')
+plt.axis('off')
+plt.subplot(322)
+plt.plot(hist[1][:-1], hist[0], lw=2)
+plt.title('histogram of grey values')
+plt.subplot(323)
+plt.imshow(glob, cmap=plt.cm.gray, interpolation='nearest')
+plt.axis('off')
+plt.subplot(324)
+plt.plot(glob_hist[1][:-1], glob_hist[0], lw=2)
+plt.title('histogram of grey values')
+plt.subplot(325)
+plt.imshow(loc, cmap=plt.cm.gray, interpolation='nearest')
+plt.axis('off')
+plt.subplot(326)
+plt.plot(loc_hist[1][:-1], loc_hist[0], lw=2)
+plt.title('histogram of grey values')
+"""
+.. image:: PLOT2RST.current_figure
+
+an other way to maximize the number of grey level used for an image is to apply a local autoleveling,
+i.e. here a pixel grey level is proportionally remapped between local minimum and local maximum.
+
+The following example show how local autolevel enhance the camaraman picture.
+"""
+from skimage.filter.rank import autolevel
+
+ima = data.camera()
+selem = disk(10)
+
+auto = autolevel(ima.astype(np.uint16),disk(20))
+
+# display results
+fig = plt.figure(figsize=[10,7])
+plt.subplot(1,2,1)
+plt.imshow(ima,cmap=plt.cm.gray)
+plt.xlabel('original')
+plt.subplot(1,2,2)
+plt.imshow(auto,cmap=plt.cm.gray)
+plt.xlabel('local autolevel')
+"""
+.. image:: PLOT2RST.current_figure
+
+This filter is very sensitive to local outlayers, see the little white spot in the sky left part. This is due
+to a local maximum which is very high comparing to the rest of the neighborhood. One can moderate this
+using the percentile version of the autolevel filter which uses to given percentiles (one inferior, one superior)
+in place of local minimum and maximum. The example bellow illustrate how the percentile parameters influence the
+local autolevel result.
+
+"""
+from skimage.filter.rank import percentile_autolevel
+
+image = data.camera()
+
+selem = disk(20)
+loc_autolevel = autolevel(image,selem=selem)
+loc_perc_autolevel0 = percentile_autolevel(image,selem=selem,p0=.00,p1=1.0)
+loc_perc_autolevel1 = percentile_autolevel(image,selem=selem,p0=.01,p1=.99)
+loc_perc_autolevel2 = percentile_autolevel(image,selem=selem,p0=.05,p1=.95)
+loc_perc_autolevel3 = percentile_autolevel(image,selem=selem,p0=.1,p1=.9)
+
+fig, axes = plt.subplots(nrows=3, figsize=(7, 8))
+ax0, ax1, ax2 = axes
+plt.gray()
+
+ax0.imshow(np.hstack((image,loc_autolevel)))
+ax0.set_title('original / autolevel')
+
+ax1.imshow(np.hstack((loc_perc_autolevel0,loc_perc_autolevel1)),vmin=0,vmax=255)
+ax1.set_title('percentile autolevel 0%,1%')
+ax2.imshow(np.hstack((loc_perc_autolevel2,loc_perc_autolevel3)),vmin=0,vmax=255)
+ax2.set_title('percentile autolevel 5% and 10%')
+
+for ax in axes:
+    ax.axis('off')
+
+"""
+.. image:: PLOT2RST.current_figure
+
+Morphological contrast enhancement filter replaces the central pixel by local maximum
+if the original grey level value if closest to local maximum, by the minimum local otherwise.
+
+"""
+
+from skimage.filter.rank import morph_contr_enh
+
+ima = data.camera()
+
+enh = morph_contr_enh(ima,disk(5))
+
+# display results
+fig = plt.figure(figsize=[10,7])
+plt.subplot(2,2,1)
+plt.imshow(ima,cmap=plt.cm.gray)
+plt.xlabel('original')
+plt.subplot(2,2,3)
+plt.imshow(enh,cmap=plt.cm.gray)
+plt.xlabel('local morphlogical contrast enhancement')
+plt.subplot(2,2,2)
+plt.imshow(ima[200:350,350:450],cmap=plt.cm.gray)
+plt.subplot(2,2,4)
+plt.imshow(enh[200:350,350:450],cmap=plt.cm.gray)
+
+"""
+.. image:: PLOT2RST.current_figure
+
+The percentile version of the local morphological contrast enhancement, uses percentile *p0* and *p1* instead of local
+minimum and local maximum.
+
+"""
+
+from skimage.filter.rank import percentile_morph_contr_enh
+
+ima = data.camera()
+
+penh = percentile_morph_contr_enh(ima,disk(5),p0=.1,p1=.9)
+
+# display results
+fig = plt.figure(figsize=[10,7])
+plt.subplot(2,2,1)
+plt.imshow(ima,cmap=plt.cm.gray)
+plt.xlabel('original')
+plt.subplot(2,2,3)
+plt.imshow(penh,cmap=plt.cm.gray)
+plt.xlabel('local percentile morphlogical\n contrast enhancement')
+plt.subplot(2,2,2)
+plt.imshow(ima[200:350,350:450],cmap=plt.cm.gray)
+plt.subplot(2,2,4)
+plt.imshow(penh[200:350,350:450],cmap=plt.cm.gray)
+
+"""
+.. image:: PLOT2RST.current_figure
+
+Image morphology
+================
+
+Local maximum and local minimum are the base operators for grey level morphology.
+
+.. note:: ``skimage.dilate`` and ``skimage.erode`` are equivalent filters (see below for comparison).
+
+Here is an example of classical morphological grey level filters : opening, closing and morphological gradient.
+
+"""
+
+from skimage.filter.rank import maximum,minimum,gradient
+
+ima = data.camera()
+
+closing = maximum(minimum(ima,disk(5)),disk(5))
+opening = minimum(maximum(ima,disk(5)),disk(5))
+grad = gradient(ima,disk(5))
+
+# display results
+fig = plt.figure(figsize=[10,7])
+plt.subplot(2,2,1)
+plt.imshow(ima,cmap=plt.cm.gray)
+plt.xlabel('original')
+plt.subplot(2,2,2)
+plt.imshow(closing,cmap=plt.cm.gray)
+plt.xlabel('grey level closing')
+plt.subplot(2,2,3)
+plt.imshow(opening,cmap=plt.cm.gray)
+plt.xlabel('grey level opening')
+plt.subplot(2,2,4)
+plt.imshow(grad,cmap=plt.cm.gray)
+plt.xlabel('morphological gradient')
+
+"""
+.. image:: PLOT2RST.current_figure
+
+Feature extraction
+===================
+
+Local histogram can be exploited to compute local entropy, which is related to the local image complexity.
+Entropy is computed using base 2 logarithm i.e. the filter returns the minimum number of bits needed to encode local
+greylevel distribution.
+
+``skimage.rank.entropy`` returns local entropy on a given structuring element.
+The following example shows this filter applied on 8- and 16- bit images.
+
+.. note:: to better use the available image bit, the function returns 10x entropy for 8-bit images and 1000x entropy
+   for 16-bit images.
+
+"""
+
+
+
+from skimage import data
+from skimage.filter.rank import entropy
+from skimage.morphology import disk
+import numpy as np
+import matplotlib.pyplot as plt
+
+# defining a 8- and a 16-bit test images
+a8 = data.camera()
+a16 = data.camera().astype(np.uint16)*4
+
+ent8 = entropy(a8,disk(5)) # pixel value contain 10x the local entropy
+ent16 = entropy(a16,disk(5)) # pixel value contain 1000x the local entropy
+
+# display results
+plt.figure(figsize=(10, 10))
+
+plt.subplot(2,2,1)
+plt.imshow(a8, cmap=plt.cm.gray)
+plt.xlabel('8-bit image')
+plt.colorbar()
+
+plt.subplot(2,2,2)
+plt.imshow(ent8, cmap=plt.cm.jet)
+plt.xlabel('entropy*10')
+plt.colorbar()
+
+plt.subplot(2,2,3)
+plt.imshow(a16, cmap=plt.cm.gray)
+plt.xlabel('16-bit image')
+plt.colorbar()
+
+plt.subplot(2,2,4)
+plt.imshow(ent16, cmap=plt.cm.jet)
+plt.xlabel('entropy*1000')
+plt.colorbar()
+plt.show()
+
+"""
+.. image:: PLOT2RST.current_figure
+
+Implementation
+================
+
+The central part of the ``skimage.rank`` filters is build on a sliding window that update local grey level histogram.
+This approach limits the algorithm complexity to O(n) where n is the number of image pixels. The complexity is also
+limited with respect to the structuring element size.
+
+"""
+from time import time
+
+from scipy.ndimage.filters import percentile_filter
+from skimage.morphology import dilation
+from skimage.filter import median_filter
+from skimage.filter.rank import median,maximum
+
+def exec_and_timeit(func):
+    """ Decorator that returns both function results and execution time
+    (result, ms)
+    """
+    def wrapper(*arg):
+        t1 = time()
+        res = func(*arg)
+        t2 = time()
+        ms = (t2-t1)*1000.0
+        return (res,ms)
+    return wrapper
+
+
+@exec_and_timeit
+def cr_med(image,selem):
+    return median(image=image,selem = selem)
+
+@exec_and_timeit
+def cr_max(image,selem):
+    return maximum(image=image,selem = selem)
+
+@exec_and_timeit
+def cm_dil(image,selem):
+    return dilation(image=image,selem = selem)
+
+@exec_and_timeit
+def ctmf_med(image,radius):
+    return median_filter(image=image,radius=radius)
+
+@exec_and_timeit
+def ndi_med(image,n):
+    return percentile_filter(image,50,size=n*2-1)
+
+"""
+Comparison between
+
+* ``rank.maximum``
+* ``cmorph.dilate``
+
+on increasing structuring element size
+"""
+
+a = data.camera()
+
+rec = []
+e_range = range(1,20,2)
+for r in e_range:
+    elem = disk(r+1)
+    rc,ms_rc = cr_max(a,elem)
+    rcm,ms_rcm = cm_dil(a,elem)
+    rec.append((ms_rc,ms_rcm))
+
+rec = np.asarray(rec)
+
+plt.figure()
+plt.title('increasing element size')
+plt.ylabel('time (ms)')
+plt.xlabel('element radius')
+plt.plot(e_range,rec)
+plt.legend(['crank.maximum','cmorph.dilate'])
+
+"""
+and increasing image size
+
+.. image:: PLOT2RST.current_figure
+"""
+
+r = 9
+elem = disk(r+1)
+
+rec = []
+s_range = range(100,1000,100)
+for s in s_range:
+    a = (np.random.random((s,s))*256).astype('uint8')
+    (rc,ms_rc) = cr_max(a,elem)
+    (rcm,ms_rcm) = cm_dil(a,elem)
+    rec.append((ms_rc,ms_rcm))
+
+rec = np.asarray(rec)
+
+plt.figure()
+plt.title('increasing image size')
+plt.ylabel('time (ms)')
+plt.xlabel('image size')
+plt.plot(s_range,rec)
+plt.legend(['crank.maximum','cmorph.dilate'])
+
+
+"""
+.. image:: PLOT2RST.current_figure
+
+Comparison between:
+
+* ``rank.median``
+* ``ctmf.median_filter``
+* ``ndimage.percentile``
+
+on increasing structuring element size
+"""
+
+
+a = data.camera()
+
+rec = []
+e_range = range(2,30,4)
+for r in e_range:
+    elem = disk(r+1)
+    rc,ms_rc = cr_med(a,elem)
+    rctmf,ms_rctmf = ctmf_med(a,r)
+    rndi,ms_ndi = ndi_med(a,r)
+    rec.append((ms_rc,ms_rctmf,ms_ndi))
+
+rec = np.asarray(rec)
+
+plt.figure()
+plt.title('increasing element size')
+plt.plot(e_range,rec)
+plt.legend(['rank.median','ctmf.median_filter','ndimage.percentile'])
+plt.ylabel('time (ms)')
+plt.xlabel('element radius')
+"""
+.. image:: PLOT2RST.current_figure
+
+comparison of outcome of the three methods
+
+"""
+plt.figure()
+plt.imshow(np.hstack((rc,rctmf,rndi)))
+plt.xlabel('rank.median vs ctmf.median_filter vs ndimage.percentile')
+
+"""
+.. image:: PLOT2RST.current_figure
+
+and increasing image size
+
+"""
+
+r = 9
+elem = disk(r+1)
+
+rec = []
+s_range = [100,200,500,1000]
+for s in s_range:
+    a = (np.random.random((s,s))*256).astype('uint8')
+    (rc,ms_rc) = cr_med(a,elem)
+    rctmf,ms_rctmf = ctmf_med(a,r)
+    rndi,ms_ndi = ndi_med(a,r)
+    rec.append((ms_rc,ms_rctmf,ms_ndi))
+
+rec = np.asarray(rec)
+
+plt.figure()
+plt.title('increasing image size')
+plt.plot(s_range,rec)
+plt.legend(['rank.median','ctmf.median_filter','ndimage.percentile'])
+plt.ylabel('time (ms)')
+plt.xlabel('image size')
+
+"""
+.. image:: PLOT2RST.current_figure
+"""
+
+plt.show()
diff --git a/doc/examples/plot_16bitbilateral.py b/doc/examples/plot_16bitbilateral.py
new file mode 100644
index 00000000..fc30aa6b
--- /dev/null
+++ b/doc/examples/plot_16bitbilateral.py
@@ -0,0 +1,47 @@
+"""
+==============================
+Bilateral mean
+==============================
+This example compares
+
+* local mean
+* percentile mean
+* bilateral mean
+
+build on the local histogram distribution
+local mean uses all pixels belonging to the structuring element to compute average gray level,
+percentile mean uses only values between percentiles p0 and p1 (here 10% and 90%),
+whereas bilateral mean uses only pixels of the structuring element having a gray level situated inside
+g-s0 and g+s1 (here g-500 and g+500).
+The filters are applied on a 16 bit image (actual bitdepth is 12bit).
+
+Percentile and usual mean give here similar results, these filters smooth the complete image (background and details).
+Bilateral mean exhibits a high filtering rate for continuous area (i.e. background) while image higher frequencies
+remains untouched.
+
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage import data
+from skimage.morphology import disk
+import skimage.filter.rank as rank
+
+a16 = (data.coins()).astype('uint16')*16
+selem = disk(20)
+
+f1 = rank.percentile_mean(a16,selem = selem,p0=.1,p1=.9)
+f2 = rank.bilateral_mean(a16,selem = selem,s0=500,s1=500)
+f3 = rank.mean(a16,selem = selem)
+
+# display results
+fig, axes = plt.subplots(nrows=3, figsize=(15,10))
+ax0, ax1, ax2 = axes
+
+ax0.imshow(np.hstack((a16,f1)))
+ax0.set_title('percentile mean')
+ax1.imshow(np.hstack((a16,f2)))
+ax1.set_title('bilateral mean')
+ax2.imshow(np.hstack((a16,f3)))
+ax2.set_title('local mean')
+plt.show()
diff --git a/doc/examples/plot_entropy.py b/doc/examples/plot_entropy.py
new file mode 100644
index 00000000..f019d79c
--- /dev/null
+++ b/doc/examples/plot_entropy.py
@@ -0,0 +1,44 @@
+"""
+===================
+Entropy
+===================
+
+
+"""
+from skimage import data
+from skimage.filter.rank import entropy
+from skimage.morphology import disk
+import numpy as np
+import matplotlib.pyplot as plt
+
+# defining a 8- and a 16-bit test images
+a8 = data.camera()
+a16 = data.camera().astype(np.uint16)*4
+
+ent8 = entropy(a8,disk(5)) # pixel value contain 10x the local entropy
+ent16 = entropy(a16,disk(5)) # pixel value contain 1000x the local entropy
+
+# display results
+plt.figure(figsize=(10, 10))
+
+plt.subplot(2,2,1)
+plt.imshow(a8, cmap=plt.cm.gray)
+plt.xlabel('8-bit image')
+plt.colorbar()
+
+plt.subplot(2,2,2)
+plt.imshow(ent8, cmap=plt.cm.jet)
+plt.xlabel('entropy*10')
+plt.colorbar()
+
+plt.subplot(2,2,3)
+plt.imshow(a16, cmap=plt.cm.gray)
+plt.xlabel('16-bit image')
+plt.colorbar()
+
+plt.subplot(2,2,4)
+plt.imshow(ent16, cmap=plt.cm.jet)
+plt.xlabel('entropy*1000')
+plt.colorbar()
+plt.show()
+
diff --git a/doc/examples/plot_local_equalize.py b/doc/examples/plot_local_equalize.py
new file mode 100644
index 00000000..1a431f5c
--- /dev/null
+++ b/doc/examples/plot_local_equalize.py
@@ -0,0 +1,84 @@
+"""
+===============================
+Local Histogram Equalization
+===============================
+
+This examples enhances an image with low contrast, using a method called
+*local histogram equalization*, which "spreads out the most frequent intensity
+values" in an image .
+The equalized image [1]_ has a roughly linear cumulative distribution function for each pixel neighborhood.
+The local version [2]_ of the histogram equalization emphasized every local graylevel variations.
+
+.. [1] http://en.wikipedia.org/wiki/Histogram_equalization
+.. [2] http://en.wikipedia.org/wiki/Adaptive_histogram_equalization
+
+"""
+
+from skimage import data
+from skimage.util.dtype import dtype_range
+from skimage import exposure
+from skimage.morphology import disk
+
+import matplotlib.pyplot as plt
+
+import numpy as np
+from skimage.filter import rank
+
+def plot_img_and_hist(img, axes, bins=256):
+    """Plot an image along with its histogram and cumulative histogram.
+
+    """
+    ax_img, ax_hist = axes
+    ax_cdf = ax_hist.twinx()
+
+    # Display image
+    ax_img.imshow(img, cmap=plt.cm.gray)
+    ax_img.set_axis_off()
+
+    # Display histogram
+    ax_hist.hist(img.ravel(), bins=bins)
+    ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
+    ax_hist.set_xlabel('Pixel intensity')
+
+    xmin, xmax = dtype_range[img.dtype.type]
+    ax_hist.set_xlim(xmin, xmax)
+
+    # Display cumulative distribution
+    img_cdf, bins = exposure.cumulative_distribution(img, bins)
+    ax_cdf.plot(bins, img_cdf, 'r')
+
+    return ax_img, ax_hist, ax_cdf
+
+
+# Load an example image
+img = data.moon()
+
+# Contrast stretching
+p2 = np.percentile(img, 2)
+p98 = np.percentile(img, 98)
+img_rescale = exposure.equalize(img)
+
+# Equalization
+selem = disk(30)
+img_eq = rank.equalize(img,selem=selem)
+
+
+# Display results
+f, axes = plt.subplots(2, 3, figsize=(8, 4))
+
+ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0])
+ax_img.set_title('Low contrast image')
+ax_hist.set_ylabel('Number of pixels')
+
+ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_rescale, axes[:, 1])
+ax_img.set_title('Global equalise')
+
+ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_eq, axes[:, 2])
+ax_img.set_title('Local equalize')
+ax_cdf.set_ylabel('Fraction of total intensity')
+
+
+# prevent overlap of y-axis labels
+plt.subplots_adjust(wspace=0.4)
+plt.show()
+
diff --git a/doc/examples/plot_marked_watershed.py b/doc/examples/plot_marked_watershed.py
new file mode 100644
index 00000000..738a3d24
--- /dev/null
+++ b/doc/examples/plot_marked_watershed.py
@@ -0,0 +1,54 @@
+"""
+================================
+Markers for watershed transform
+================================
+
+The watershed is a classical algorithm used for **segmentation**, that
+is, for separating different objects in an image.
+
+Here a marker image is build from the region of low gradient inside the image.
+
+See Wikipedia_ for more details on the algorithm.
+
+.. _Wikipedia: http://en.wikipedia.org/wiki/Watershed_(image_processing)
+
+"""
+
+from scipy import ndimage
+import matplotlib.pyplot as plt
+from skimage.morphology import watershed,disk
+from skimage import data
+
+# original data
+from skimage.filter import rank
+
+image = data.camera()
+
+# denoise image
+denoised = rank.median(image,disk(2))
+
+# find continuous region (low gradient) --> markers
+markers = rank.gradient(denoised,disk(5))<10
+markers = ndimage.label(markers)[0]
+
+#local gradient
+gradient = rank.gradient(denoised,disk(2))
+
+# process the watershed
+labels = watershed(gradient, markers)
+
+# display results
+fig, axes = plt.subplots(ncols=4, figsize=(8, 2.7))
+ax0, ax1, ax2, ax3 = axes
+
+ax0.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
+ax1.imshow(gradient, cmap=plt.cm.spectral, interpolation='nearest')
+ax2.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest')
+ax3.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
+ax3.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest',alpha=.7)
+
+for ax in axes:
+    ax.axis('off')
+
+plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
+plt.show()
diff --git a/doc/examples/plot_watershed.py b/doc/examples/plot_watershed.py
index a1cd18cf..50857b2c 100644
--- a/doc/examples/plot_watershed.py
+++ b/doc/examples/plot_watershed.py
@@ -26,7 +26,7 @@ See Wikipedia_ for more details on the algorithm.
 """
 
 import numpy as np
-from scipy import ndimage
+
 import matplotlib.pyplot as plt
 from skimage.morphology import watershed, is_local_maximum
 
diff --git a/skimage/filter/rank/README.rst b/skimage/filter/rank/README.rst
new file mode 100644
index 00000000..aae62162
--- /dev/null
+++ b/skimage/filter/rank/README.rst
@@ -0,0 +1,35 @@
+To use this to build your Cython file use the commandline options:
+
+.. sourcecode:: text
+
+    $ python setup.py build_ext --inplace
+
+
+**To do**
+
+* add simple examples, adapt documentation on existing examples
+
+* add/check existing doc
+
+* adapting tests for each type of filter
+
+**General remarks**
+
+Basically these filters compute local histogram for each pixel. Histogram is build using a moving window in
+order to limit redundant computation. The path followed by the moving window is given hereunder
+
+ ...-----------------------\
+/--------------------------/
+\-------------------------- ...
+
+A comparison is proposed with cmorph.dilate algorithm to show how computation costs evolve with respect to image size or
+structuring element size. This implementation gives better results for large structuring elements.
+
+A local histogram is update at each pixel by introducing pixel entering the structuring element border and
+by removing those leaving it. The histogram size is 8bit (256 bins) for 8 bit images and 2 to 12 bit (up to 4096 bins)
+for 16bit image depending on the image maximum value. Image with pixels higher than 4095 raise a ValueError.
+
+The filter is applied up to the image border, the neighboorhood used is adjusted accordingly. The user may provide
+a mask image (same size as input image) where non zero value are the part of the image participating the the
+histogram computation. By default all the image is filtered.
+
diff --git a/skimage/filter/rank/__init__.py b/skimage/filter/rank/__init__.py
new file mode 100644
index 00000000..30d936db
--- /dev/null
+++ b/skimage/filter/rank/__init__.py
@@ -0,0 +1,3 @@
+from .rank import *
+from .percentile_rank import *
+from .bilateral_rank import *
diff --git a/skimage/filter/rank/_core16.pxd b/skimage/filter/rank/_core16.pxd
new file mode 100644
index 00000000..9590e277
--- /dev/null
+++ b/skimage/filter/rank/_core16.pxd
@@ -0,0 +1,18 @@
+cimport numpy as np
+
+#---------------------------------------------------------------------------
+# 16 bit core kernel receives extra information about data bitdepth
+#---------------------------------------------------------------------------
+
+# generic cdef functions
+cdef int int_max(int a, int b)
+cdef int int_min(int a, int b)
+
+cdef _core16(
+    np.uint16_t kernel(Py_ssize_t * , float, np.uint16_t, Py_ssize_t, Py_ssize_t, Py_ssize_t, float, float, Py_ssize_t, Py_ssize_t),
+    np.ndarray[np.uint16_t, ndim=2] image,
+    np.ndarray[np.uint8_t, ndim=2] selem,
+    np.ndarray[np.uint8_t, ndim=2] mask,
+    np.ndarray[np.uint16_t, ndim=2] out,
+    char shift_x, char shift_y, Py_ssize_t bitdepth,
+    float p0, float p1, Py_ssize_t s0, Py_ssize_t s1)
diff --git a/skimage/filter/rank/_core16.pyx b/skimage/filter/rank/_core16.pyx
new file mode 100644
index 00000000..4829792b
--- /dev/null
+++ b/skimage/filter/rank/_core16.pyx
@@ -0,0 +1,283 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+
+import numpy as np
+cimport numpy as np
+from libc.stdlib cimport malloc, free
+
+#---------------------------------------------------------------------------
+# 16 bit core kernel receives extra information about data bitdepth
+#---------------------------------------------------------------------------
+
+# generic cdef functions
+cdef inline int int_max(int a, int b):
+    return a if a >= b else b
+cdef inline int int_min(int a, int b):
+    return a if a <= b else b
+
+cdef inline void histogram_increment(Py_ssize_t * histo, float * pop, np.uint16_t value):
+    histo[value] += 1
+    pop[0] += 1.
+
+cdef inline void histogram_decrement(Py_ssize_t * histo, float * pop, np.uint16_t value):
+    histo[value] -= 1
+    pop[0] -= 1.
+
+cdef inline np.uint8_t is_in_mask(Py_ssize_t rows, Py_ssize_t cols, Py_ssize_t r, Py_ssize_t c, np.uint8_t * mask):
+    """ returns 1 if given(r,c) coordinate are within the image frame ([0-rows],[0-cols]) and
+        inside the given mask
+        returns 0 otherwise
+    """
+    if r < 0 or r > rows - 1 or c < 0 or c > cols - 1:
+        return 0
+    else:
+        if mask[r * cols + c]:
+            return 1
+        else:
+            return 0
+
+
+cdef inline _core16(
+    np.uint16_t kernel(Py_ssize_t *, float, np.uint16_t, Py_ssize_t, Py_ssize_t, Py_ssize_t, float, float, Py_ssize_t, Py_ssize_t),
+    np.ndarray[np.uint16_t, ndim=2] image,
+    np.ndarray[np.uint8_t, ndim=2] selem,
+    np.ndarray[np.uint8_t, ndim=2] mask,
+    np.ndarray[np.uint16_t, ndim=2] out,
+    char shift_x, char shift_y, Py_ssize_t bitdepth,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    """ Main loop, this function computes the histogram for each image point
+    - data is uint8
+    - result is uint8 casted
+    """
+
+    cdef Py_ssize_t rows = image.shape[0]
+    cdef Py_ssize_t cols = image.shape[1]
+    cdef Py_ssize_t srows = selem.shape[0]
+    cdef Py_ssize_t scols = selem.shape[1]
+
+    cdef Py_ssize_t centre_r = int(selem.shape[0] / 2) + shift_y
+    cdef Py_ssize_t centre_c = int(selem.shape[1] / 2) + shift_x
+
+    # check that structuring element center is inside the element bounding box
+    assert centre_r >= 0
+    assert centre_c >= 0
+    assert centre_r < srows
+    assert centre_c < scols
+    assert bitdepth in range(2, 13)
+
+    maxbin_list = [0, 0, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096]
+    midbin_list = [0, 0, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048]
+
+    #set maxbin and midbin
+    cdef Py_ssize_t maxbin = maxbin_list[bitdepth], midbin = midbin_list[bitdepth]
+
+    assert (image < maxbin).all()
+
+    image = np.ascontiguousarray(image)
+
+    if mask is None:
+        mask = np.ones((rows, cols), dtype=np.uint8)
+    else:
+        mask = np.ascontiguousarray(mask)
+
+    if image is out:
+        raise NotImplementedError("Cannot perform rank operation in place.")
+
+    if out is None:
+        out = np.zeros((rows, cols), dtype=np.uint16)
+    else:
+        out = np.ascontiguousarray(out)
+
+    mask = np.ascontiguousarray(mask)
+
+    # define pointers to the data
+    cdef np.uint16_t * out_data = <np.uint16_t * >out.data
+    cdef np.uint16_t * image_data = <np.uint16_t * >image.data
+    cdef np.uint8_t * mask_data = <np.uint8_t * >mask.data
+
+    # define local variable types
+    cdef Py_ssize_t r, c, rr, cc, s, value, local_max, i, even_row
+    cdef float pop                                 # number of pixels actually inside the neighborhood (float)
+
+    # allocate memory with malloc
+    cdef Py_ssize_t max_se = srows * scols
+
+    # number of element in each attack border
+    cdef Py_ssize_t num_se_n, num_se_s, num_se_e, num_se_w
+
+    # the current local histogram distribution
+    cdef Py_ssize_t * histo = <Py_ssize_t * >malloc(maxbin * sizeof(Py_ssize_t))
+
+    # these lists contain the relative pixel row and column for each of the 4 attack borders
+    # east, west, north and south
+    # e.g. se_e_r lists the rows of the east structuring element border
+
+    cdef Py_ssize_t * se_e_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_e_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_w_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_w_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_n_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_n_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_s_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_s_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+
+    # build attack and release borders
+    # by using difference along axis
+
+    t = np.hstack((selem, np.zeros((selem.shape[0], 1))))
+    t_e = np.diff(t, axis=1) == -1
+
+    t = np.hstack((np.zeros((selem.shape[0], 1)), selem))
+    t_w = np.diff(t, axis=1) == 1
+
+    t = np.vstack((selem, np.zeros((1, selem.shape[1]))))
+    t_s = np.diff(t, axis=0) == -1
+
+    t = np.vstack((np.zeros((1, selem.shape[1])), selem))
+    t_n = np.diff(t, axis=0) == 1
+
+    num_se_n = num_se_s = num_se_e = num_se_w = 0
+
+    for r in range(srows):
+        for c in range(scols):
+            if t_e[r, c]:
+                se_e_r[num_se_e] = r - centre_r
+                se_e_c[num_se_e] = c - centre_c
+                num_se_e += 1
+            if t_w[r, c]:
+                se_w_r[num_se_w] = r - centre_r
+                se_w_c[num_se_w] = c - centre_c
+                num_se_w += 1
+            if t_n[r, c]:
+                se_n_r[num_se_n] = r - centre_r
+                se_n_c[num_se_n] = c - centre_c
+                num_se_n += 1
+            if t_s[r, c]:
+                se_s_r[num_se_s] = r - centre_r
+                se_s_c[num_se_s] = c - centre_c
+                num_se_s += 1
+
+    # initial population and histogram
+    for i in range(maxbin):
+        histo[i] = 0
+
+    pop = 0
+
+    for r in range(srows):
+        for c in range(scols):
+            rr = r - centre_r
+            cc = c - centre_c
+            if selem[r, c]:
+                if is_in_mask(rows, cols, rr, cc, mask_data):
+                    histogram_increment(histo, & pop, image_data[rr * cols + cc])
+
+    r = 0
+    c = 0
+    # kernel -------------------------------------------
+    out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c],
+        bitdepth, maxbin, midbin, p0, p1, s0, s1)
+    # kernel -------------------------------------------
+
+    # main loop
+    r = 0
+    for even_row in range(0, rows, 2):
+        # ---> west to east
+        for c in range(1, cols):
+            for s in range(num_se_e):
+                rr = r + se_e_r[s]
+                cc = c + se_e_c[s]
+                if is_in_mask(rows, cols, rr, cc, mask_data):
+                    histogram_increment(histo, & pop, image_data[rr * cols + cc])
+
+            for s in range(num_se_w):
+                rr = r + se_w_r[s]
+                cc = c + se_w_c[s] - 1
+                if is_in_mask(rows, cols, rr, cc, mask_data):
+                    histogram_decrement(histo, & pop, image_data[rr * cols + cc])
+
+            # kernel -------------------------------------------
+            out_data[r * cols + c] = kernel(
+                histo, pop, image_data[r * cols + c],
+                bitdepth, maxbin, midbin, p0, p1, s0, s1)
+            # kernel -------------------------------------------
+
+        r += 1          # pass to the next row
+        if r >= rows:
+            break
+
+            # ---> north to south
+        for s in range(num_se_s):
+            rr = r + se_s_r[s]
+            cc = c + se_s_c[s]
+            if is_in_mask(rows, cols, rr, cc, mask_data):
+                histogram_increment(histo, & pop, image_data[rr * cols + cc])
+
+        for s in range(num_se_n):
+            rr = r + se_n_r[s] - 1
+            cc = c + se_n_c[s]
+            if is_in_mask(rows, cols, rr, cc, mask_data):
+                histogram_decrement(histo, & pop, image_data[rr * cols + cc])
+
+        # kernel -------------------------------------------
+        out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c],
+            bitdepth, maxbin, midbin, p0, p1, s0, s1)
+        # kernel -------------------------------------------
+
+        # ---> east to west
+        for c in range(cols - 2, -1, -1):
+            for s in range(num_se_w):
+                rr = r + se_w_r[s]
+                cc = c + se_w_c[s]
+                if is_in_mask(rows, cols, rr, cc, mask_data):
+                    histogram_increment(histo, & pop, image_data[rr * cols + cc])
+
+            for s in range(num_se_e):
+                rr = r + se_e_r[s]
+                cc = c + se_e_c[s] + 1
+                if is_in_mask(rows, cols, rr, cc, mask_data):
+                    histogram_decrement(histo, & pop, image_data[rr * cols + cc])
+
+            # kernel -------------------------------------------
+            out_data[r * cols + c] = kernel(
+                histo, pop, image_data[r * cols + c],
+                bitdepth, maxbin, midbin, p0, p1, s0, s1)
+            # kernel -------------------------------------------
+
+        r += 1           # pass to the next row
+        if r >= rows:
+            break
+
+        # ---> north to south
+        for s in range(num_se_s):
+            rr = r + se_s_r[s]
+            cc = c + se_s_c[s]
+            if is_in_mask(rows, cols, rr, cc, mask_data):
+                histogram_increment(histo, & pop, image_data[rr * cols + cc])
+
+        for s in range(num_se_n):
+            rr = r + se_n_r[s] - 1
+            cc = c + se_n_c[s]
+            if is_in_mask(rows, cols, rr, cc, mask_data):
+                histogram_decrement(histo, & pop, image_data[rr * cols + cc])
+
+        # kernel -------------------------------------------
+        out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols + c],
+            bitdepth, maxbin, midbin, p0, p1, s0, s1)
+        # kernel -------------------------------------------
+
+    # release memory allocated by malloc
+
+    free(se_e_r)
+    free(se_e_c)
+    free(se_w_r)
+    free(se_w_c)
+    free(se_n_r)
+    free(se_n_c)
+    free(se_s_r)
+    free(se_s_c)
+
+    free(histo)
+
+    return out
diff --git a/skimage/filter/rank/_core8.pxd b/skimage/filter/rank/_core8.pxd
new file mode 100644
index 00000000..9f898faa
--- /dev/null
+++ b/skimage/filter/rank/_core8.pxd
@@ -0,0 +1,17 @@
+cimport numpy as np
+
+# generic cdef functions
+cdef np.uint8_t uint8_max(np.uint8_t a, np.uint8_t b)
+cdef np.uint8_t uint8_min(np.uint8_t a, np.uint8_t b)
+
+#---------------------------------------------------------------------------
+# 8 bit core kernel receives extra information about data inferior and superior percentiles
+#---------------------------------------------------------------------------
+
+cdef _core8(
+    np.uint8_t kernel(Py_ssize_t *, float, np.uint8_t, float, float, Py_ssize_t, Py_ssize_t),
+    np.ndarray[np.uint8_t, ndim=2] image,
+    np.ndarray[np.uint8_t, ndim=2] selem,
+    np.ndarray[np.uint8_t, ndim=2] mask,
+    np.ndarray[np.uint8_t, ndim=2] out,
+    char shift_x, char shift_y, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1)
diff --git a/skimage/filter/rank/_core8.pyx b/skimage/filter/rank/_core8.pyx
new file mode 100644
index 00000000..9955a1e1
--- /dev/null
+++ b/skimage/filter/rank/_core8.pyx
@@ -0,0 +1,274 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+
+import numpy as np
+cimport numpy as np
+from libc.stdlib cimport malloc, free
+
+# generic cdef functions
+cdef inline np.uint8_t uint8_max(np.uint8_t a, np.uint8_t b):
+    return a if a >= b else b
+cdef inline np.uint8_t uint8_min(np.uint8_t a, np.uint8_t b):
+    return a if a <= b else b
+
+
+#---------------------------------------------------------------------------
+# 8 bit core kernel
+#---------------------------------------------------------------------------
+
+cdef inline void histogram_increment(Py_ssize_t * histo, float * pop, np.uint8_t value):
+    histo[value] += 1
+    pop[0] += 1.
+
+cdef inline void histogram_decrement(Py_ssize_t * histo, float * pop, np.uint8_t value):
+    histo[value] -= 1
+    pop[0] -= 1.
+
+cdef inline np.uint8_t is_in_mask(Py_ssize_t rows, Py_ssize_t cols, Py_ssize_t r, Py_ssize_t c, np.uint8_t * mask):
+    """ returns 1 if given(r,c) coordinate are within the image frame ([0-rows],[0-cols]) and
+        inside the given mask
+        returns 0 otherwise
+    """
+    if r < 0 or r > rows - 1 or c < 0 or c > cols - 1:
+        return 0
+    else:
+        if mask[r * cols + c]:
+            return 1
+        else:
+            return 0
+
+cdef inline _core8(
+    np.uint8_t kernel(Py_ssize_t * , float, np.uint8_t, float, float, Py_ssize_t, Py_ssize_t),
+    np.ndarray[np.uint8_t, ndim=2] image,
+    np.ndarray[np.uint8_t, ndim=2] selem,
+    np.ndarray[np.uint8_t, ndim=2] mask,
+    np.ndarray[np.uint8_t, ndim=2] out,
+        char shift_x, char shift_y, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    """ Main loop, this function computes the histogram for each image point
+    - data is uint8
+    - result is uint8 casted
+    """
+
+    cdef Py_ssize_t rows = image.shape[0]
+    cdef Py_ssize_t cols = image.shape[1]
+    cdef Py_ssize_t srows = selem.shape[0]
+    cdef Py_ssize_t scols = selem.shape[1]
+
+    cdef Py_ssize_t centre_r = int(selem.shape[0] / 2) + shift_y
+    cdef Py_ssize_t centre_c = int(selem.shape[1] / 2) + shift_x
+
+    # check that structuring element center is inside the element bounding box
+    assert centre_r >= 0
+    assert centre_c >= 0
+    assert centre_r < srows
+    assert centre_c < scols
+
+    image = np.ascontiguousarray(image)
+
+    if mask is None:
+        mask = np.ones((rows, cols), dtype=np.uint8)
+    else:
+        mask = np.ascontiguousarray(mask)
+
+    if image is out:
+        raise NotImplementedError("Cannot perform rank operation in place.")
+
+    if out is None:
+        out = np.zeros((rows, cols), dtype=np.uint8)
+    else:
+        out = np.ascontiguousarray(out)
+
+    mask = np.ascontiguousarray(mask)
+
+    # define pointers to the data
+
+    cdef np.uint8_t * out_data = <np.uint8_t * >out.data
+    cdef np.uint8_t * image_data = <np.uint8_t * >image.data
+    cdef np.uint8_t * mask_data = <np.uint8_t * >mask.data
+
+    # define local variable types
+    cdef Py_ssize_t r, c, rr, cc, s, value, local_max, i, even_row
+
+    # number of pixels actually inside the neighborhood (float)
+    cdef float pop
+
+    # allocate memory with malloc
+    cdef Py_ssize_t max_se = srows * scols
+
+    # number of element in each attack border
+    cdef Py_ssize_t num_se_n, num_se_s, num_se_e, num_se_w
+
+    # the current local histogram distribution
+    cdef Py_ssize_t * histo = <Py_ssize_t * >malloc(256 * sizeof(Py_ssize_t))
+
+    # these lists contain the relative pixel row and column for each of the 4 attack borders
+    # east, west, north and south
+    # e.g. se_e_r lists the rows of the east structuring element border
+
+    cdef Py_ssize_t * se_e_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_e_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_w_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_w_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_n_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_n_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_s_r = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+    cdef Py_ssize_t * se_s_c = <Py_ssize_t * >malloc(max_se * sizeof(Py_ssize_t))
+
+    # build attack and release borders
+    # by using difference along axis
+
+    t = np.hstack((selem, np.zeros((selem.shape[0], 1))))
+    t_e = np.diff(t, axis=1) == -1
+
+    t = np.hstack((np.zeros((selem.shape[0], 1)), selem))
+    t_w = np.diff(t, axis=1) == 1
+
+    t = np.vstack((selem, np.zeros((1, selem.shape[1]))))
+    t_s = np.diff(t, axis=0) == -1
+
+    t = np.vstack((np.zeros((1, selem.shape[1])), selem))
+    t_n = np.diff(t, axis=0) == 1
+
+    num_se_n = num_se_s = num_se_e = num_se_w = 0
+
+    for r in range(srows):
+        for c in range(scols):
+            if t_e[r, c]:
+                se_e_r[num_se_e] = r - centre_r
+                se_e_c[num_se_e] = c - centre_c
+                num_se_e += 1
+            if t_w[r, c]:
+                se_w_r[num_se_w] = r - centre_r
+                se_w_c[num_se_w] = c - centre_c
+                num_se_w += 1
+            if t_n[r, c]:
+                se_n_r[num_se_n] = r - centre_r
+                se_n_c[num_se_n] = c - centre_c
+                num_se_n += 1
+            if t_s[r, c]:
+                se_s_r[num_se_s] = r - centre_r
+                se_s_c[num_se_s] = c - centre_c
+                num_se_s += 1
+
+    # initial population and histogram (kernel is centered on the first row and column)
+    for i in range(256):
+        histo[i] = 0
+
+    pop = 0
+
+    for r in range(srows):
+        for c in range(scols):
+            rr = r - centre_r
+            cc = c - centre_c
+            if selem[r, c]:
+                if is_in_mask(rows, cols, rr, cc, mask_data):
+                    histogram_increment(histo, & pop, image_data[rr * cols + cc])
+
+    r = 0
+    c = 0
+    # kernel --------------------------------------------------------------------
+    out_data[r * cols + c] = kernel(histo, pop, image_data[r * cols +
+                                    c], p0, p1, s0, s1)
+    # kernel --------------------------------------------------------------------
+
+    # main loop
+    r = 0
+    for even_row in range(0, rows, 2):
+        # ---> west to east
+        for c in range(1, cols):
+            for s in range(num_se_e):
+                rr = r + se_e_r[s]
+                cc = c + se_e_c[s]
+                if is_in_mask(rows, cols, rr, cc, mask_data):
+                    histogram_increment(histo, & pop, image_data[rr * cols + cc])
+
+            for s in range(num_se_w):
+                rr = r + se_w_r[s]
+                cc = c + se_w_c[s] - 1
+                if is_in_mask(rows, cols, rr, cc, mask_data):
+                    histogram_decrement(histo, & pop, image_data[rr * cols + cc])
+
+            # kernel --------------------------------------------------------------------
+            out_data[r * cols + c] = kernel(
+                histo, pop, image_data[r * cols + c], p0, p1, s0, s1)
+            # kernel --------------------------------------------------------------------
+
+        r += 1          # pass to the next row
+        if r >= rows:
+            break
+
+            # ---> north to south
+        for s in range(num_se_s):
+            rr = r + se_s_r[s]
+            cc = c + se_s_c[s]
+            if is_in_mask(rows, cols, rr, cc, mask_data):
+                histogram_increment(histo, & pop, image_data[rr * cols + cc])
+
+        for s in range(num_se_n):
+            rr = r + se_n_r[s] - 1
+            cc = c + se_n_c[s]
+            if is_in_mask(rows, cols, rr, cc, mask_data):
+                histogram_decrement(histo, & pop, image_data[rr * cols + cc])
+
+        # kernel --------------------------------------------------------------------
+        out_data[r * cols + c] = kernel(histo, pop, image_data[r *
+                                        cols + c], p0, p1, s0, s1)
+        # kernel --------------------------------------------------------------------
+
+        # ---> east to west
+        for c in range(cols - 2, -1, -1):
+            for s in range(num_se_w):
+                rr = r + se_w_r[s]
+                cc = c + se_w_c[s]
+                if is_in_mask(rows, cols, rr, cc, mask_data):
+                    histogram_increment(histo, & pop, image_data[rr * cols + cc])
+
+            for s in range(num_se_e):
+                rr = r + se_e_r[s]
+                cc = c + se_e_c[s] + 1
+                if is_in_mask(rows, cols, rr, cc, mask_data):
+                    histogram_decrement(histo, & pop, image_data[rr * cols + cc])
+
+            # kernel --------------------------------------------------------------------
+            out_data[r * cols + c] = kernel(
+                histo, pop, image_data[r * cols + c], p0, p1, s0, s1)
+            # kernel --------------------------------------------------------------------
+
+        r += 1           # pass to the next row
+        if r >= rows:
+            break
+
+        # ---> north to south
+        for s in range(num_se_s):
+            rr = r + se_s_r[s]
+            cc = c + se_s_c[s]
+            if is_in_mask(rows, cols, rr, cc, mask_data):
+                histogram_increment(histo, & pop, image_data[rr * cols + cc])
+
+        for s in range(num_se_n):
+            rr = r + se_n_r[s] - 1
+            cc = c + se_n_c[s]
+            if is_in_mask(rows, cols, rr, cc, mask_data):
+                histogram_decrement(histo, & pop, image_data[rr * cols + cc])
+
+        # kernel --------------------------------------------------------------------
+        out_data[r * cols + c] = kernel(histo, pop, image_data[r *
+                                        cols + c], p0, p1, s0, s1)
+        # kernel --------------------------------------------------------------------
+
+    # release memory allocated by malloc
+
+    free(se_e_r)
+    free(se_e_c)
+    free(se_w_r)
+    free(se_w_c)
+    free(se_n_r)
+    free(se_n_c)
+    free(se_s_r)
+    free(se_s_c)
+
+    free(histo)
+
+    return out
diff --git a/skimage/filter/rank/_crank16.pyx b/skimage/filter/rank/_crank16.pyx
new file mode 100644
index 00000000..73e8e0bd
--- /dev/null
+++ b/skimage/filter/rank/_crank16.pyx
@@ -0,0 +1,364 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+
+import numpy as np
+cimport numpy as np
+from libc.math cimport log2
+
+# import main loop
+from skimage.filter.rank._core16 cimport _core16
+
+# -----------------------------------------------------------------
+# kernels uint16 take extra parameter for defining the bitdepth
+# -----------------------------------------------------------------
+
+cdef inline np.uint16_t kernel_autolevel(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i, imin, imax, delta
+
+    if pop:
+        for i in range(maxbin - 1, -1, -1):
+            if histo[i]:
+                imax = i
+                break
+        for i in range(maxbin):
+            if histo[i]:
+                imin = i
+                break
+    delta = imax - imin
+    if delta > 0:
+        return < np.uint16_t > (1. * (maxbin - 1) * (g - imin) / delta)
+    else:
+        return < np.uint16_t > (imax - imin)
+
+cdef inline np.uint16_t kernel_bottomhat(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i
+
+    for i in range(maxbin):
+        if histo[i]:
+            break
+
+    return < np.uint16_t > (g - i)
+
+
+cdef inline np.uint16_t kernel_equalize(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i
+    cdef float sum = 0.
+
+    if pop:
+        for i in range(maxbin):
+            sum += histo[i]
+            if i >= g:
+                break
+
+        return < np.uint16_t > (((maxbin - 1) * sum) / pop)
+    else:
+        return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_gradient(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i, imin, imax
+
+    if pop:
+        for i in range(maxbin - 1, -1, -1):
+            if histo[i]:
+                imax = i
+                break
+        for i in range(maxbin):
+            if histo[i]:
+                imin = i
+                break
+        return < np.uint16_t > (imax - imin)
+    else:
+        return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_maximum(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i
+
+    if pop:
+        for i in range(maxbin - 1, -1, -1):
+            if histo[i]:
+                return < np.uint16_t > (i)
+
+    return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_mean(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i
+    cdef float mean = 0.
+
+    if pop:
+        for i in range(maxbin):
+            mean += histo[i] * i
+        return < np.uint16_t > (mean / pop)
+    else:
+        return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_meansubstraction(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i
+    cdef float mean = 0.
+
+    if pop:
+        for i in range(maxbin):
+            mean += histo[i] * i
+        return < np.uint16_t > ((g - mean / pop) / 2. + (midbin - 1))
+    else:
+        return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_median(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i
+    cdef float sum = pop / 2.0
+
+    if pop:
+        for i in range(maxbin):
+            if histo[i]:
+                sum -= histo[i]
+                if sum < 0:
+                    return < np.uint16_t > (i)
+
+    return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_minimum(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i
+
+    if pop:
+        for i in range(maxbin):
+            if histo[i]:
+                return < np.uint16_t > (i)
+
+    return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_modal(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t hmax = 0, imax = 0
+
+    if pop:
+        for i in range(maxbin):
+            if histo[i] > hmax:
+                hmax = histo[i]
+                imax = i
+        return < np.uint16_t > (imax)
+
+    return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_morph_contr_enh(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i, imin, imax
+
+    if pop:
+        for i in range(maxbin - 1, -1, -1):
+            if histo[i]:
+                imax = i
+                break
+        for i in range(maxbin):
+            if histo[i]:
+                imin = i
+                break
+        if imax - g < g - imin:
+            return < np.uint16_t > (imax)
+        else:
+            return < np.uint16_t > (imin)
+    else:
+        return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_pop(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    return < np.uint16_t > (pop)
+
+cdef inline np.uint16_t kernel_threshold(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i
+    cdef float mean = 0.
+
+    if pop:
+        for i in range(maxbin):
+            mean += histo[i] * i
+        return < np.uint16_t > (g > (mean / pop))
+    else:
+        return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_tophat(
+    Py_ssize_t * histo, float pop, np.uint16_t g,
+    Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i
+
+    for i in range(maxbin - 1, -1, -1):
+        if histo[i]:
+            break
+
+    return < np.uint16_t > (i - g)
+
+
+cdef inline np.uint16_t kernel_entropy(
+        Py_ssize_t * histo, float pop, np.uint16_t g,
+        Py_ssize_t bitdepth, Py_ssize_t maxbin, Py_ssize_t midbin,
+        float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i
+    cdef float e,p
+
+    e = 0.
+
+    for i in range(maxbin):
+        p = histo[i]/pop
+        if p>0:
+            e -= p*log2(p)
+
+    return < np.uint16_t > e*1000
+
+# -----------------------------------------------------------------
+# python wrappers
+# -----------------------------------------------------------------
+
+
+def autolevel(np.ndarray[np.uint16_t, ndim=2] image,
+              np.ndarray[np.uint8_t, ndim=2] selem,
+              np.ndarray[np.uint8_t, ndim=2] mask=None,
+              np.ndarray[np.uint16_t, ndim=2] out=None,
+              char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_autolevel, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def bottomhat(np.ndarray[np.uint16_t, ndim=2] image,
+              np.ndarray[np.uint8_t, ndim=2] selem,
+              np.ndarray[np.uint8_t, ndim=2] mask=None,
+              np.ndarray[np.uint16_t, ndim=2] out=None,
+              char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_bottomhat, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def equalize(np.ndarray[np.uint16_t, ndim=2] image,
+             np.ndarray[np.uint8_t, ndim=2] selem,
+             np.ndarray[np.uint8_t, ndim=2] mask=None,
+             np.ndarray[np.uint16_t, ndim=2] out=None,
+             char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_equalize, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def gradient(np.ndarray[np.uint16_t, ndim=2] image,
+             np.ndarray[np.uint8_t, ndim=2] selem,
+             np.ndarray[np.uint8_t, ndim=2] mask=None,
+             np.ndarray[np.uint16_t, ndim=2] out=None,
+             char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_gradient, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def maximum(np.ndarray[np.uint16_t, ndim=2] image,
+            np.ndarray[np.uint8_t, ndim=2] selem,
+            np.ndarray[np.uint8_t, ndim=2] mask=None,
+            np.ndarray[np.uint16_t, ndim=2] out=None,
+            char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_maximum, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def mean(np.ndarray[np.uint16_t, ndim=2] image,
+         np.ndarray[np.uint8_t, ndim=2] selem,
+         np.ndarray[np.uint8_t, ndim=2] mask=None,
+         np.ndarray[np.uint16_t, ndim=2] out=None,
+         char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_mean, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def meansubstraction(np.ndarray[np.uint16_t, ndim=2] image,
+                     np.ndarray[np.uint8_t, ndim=2] selem,
+                     np.ndarray[np.uint8_t, ndim=2] mask=None,
+                     np.ndarray[np.uint16_t, ndim=2] out=None,
+                     char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_meansubstraction, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def median(np.ndarray[np.uint16_t, ndim=2] image,
+           np.ndarray[np.uint8_t, ndim=2] selem,
+           np.ndarray[np.uint8_t, ndim=2] mask=None,
+           np.ndarray[np.uint16_t, ndim=2] out=None,
+           char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_median, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def minimum(np.ndarray[np.uint16_t, ndim=2] image,
+            np.ndarray[np.uint8_t, ndim=2] selem,
+            np.ndarray[np.uint8_t, ndim=2] mask=None,
+            np.ndarray[np.uint16_t, ndim=2] out=None,
+            char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_minimum, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def morph_contr_enh(np.ndarray[np.uint16_t, ndim=2] image,
+                    np.ndarray[np.uint8_t, ndim=2] selem,
+                    np.ndarray[np.uint8_t, ndim=2] mask=None,
+                    np.ndarray[np.uint16_t, ndim=2] out=None,
+                    char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def modal(np.ndarray[np.uint16_t, ndim=2] image,
+          np.ndarray[np.uint8_t, ndim=2] selem,
+          np.ndarray[np.uint8_t, ndim=2] mask=None,
+          np.ndarray[np.uint16_t, ndim=2] out=None,
+          char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_modal, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def pop(np.ndarray[np.uint16_t, ndim=2] image,
+        np.ndarray[np.uint8_t, ndim=2] selem,
+        np.ndarray[np.uint8_t, ndim=2] mask=None,
+        np.ndarray[np.uint16_t, ndim=2] out=None,
+        char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_pop, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def threshold(np.ndarray[np.uint16_t, ndim=2] image,
+              np.ndarray[np.uint8_t, ndim=2] selem,
+              np.ndarray[np.uint8_t, ndim=2] mask=None,
+              np.ndarray[np.uint16_t, ndim=2] out=None,
+              char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_threshold, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def tophat(np.ndarray[np.uint16_t, ndim=2] image,
+           np.ndarray[np.uint8_t, ndim=2] selem,
+           np.ndarray[np.uint8_t, ndim=2] mask=None,
+           np.ndarray[np.uint16_t, ndim=2] out=None,
+           char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_tophat, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+def entropy(np.ndarray[np.uint16_t, ndim=2] image,
+           np.ndarray[np.uint8_t, ndim=2] selem,
+           np.ndarray[np.uint8_t, ndim=2] mask=None,
+           np.ndarray[np.uint16_t, ndim=2] out=None,
+           char shift_x=0, char shift_y=0, Py_ssize_t bitdepth=8):
+    return _core16(kernel_entropy, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
diff --git a/skimage/filter/rank/_crank16_bilateral.pyx b/skimage/filter/rank/_crank16_bilateral.pyx
new file mode 100644
index 00000000..c013b779
--- /dev/null
+++ b/skimage/filter/rank/_crank16_bilateral.pyx
@@ -0,0 +1,72 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+
+import numpy as np
+cimport numpy as np
+
+# import main loop
+from skimage.filter.rank._core16 cimport _core16
+
+# -----------------------------------------------------------------
+# kernels uint16 take extra parameter for defining the bitdepth
+# -----------------------------------------------------------------
+
+
+cdef inline np.uint16_t kernel_mean(
+    Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
+        Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+
+    cdef int i, bilat_pop = 0
+    cdef float mean = 0.
+
+    if pop:
+        for i in range(maxbin):
+            if (g > (i - s0)) and (g < (i + s1)):
+                bilat_pop += histo[i]
+                mean += histo[i] * i
+        if bilat_pop:
+            return < np.uint16_t > (mean / bilat_pop)
+        else:
+            return < np.uint16_t > (0)
+    else:
+        return < np.uint16_t > (0)
+
+
+cdef inline np.uint16_t kernel_pop(
+    Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth, Py_ssize_t maxbin,
+        Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+
+    cdef int i, bilat_pop = 0
+
+    if pop:
+        for i in range(maxbin):
+            if (g > (i - s0)) and (g < (i + s1)):
+                bilat_pop += histo[i]
+        return < np.uint16_t > (bilat_pop)
+    else:
+        return < np.uint16_t > (0)
+
+
+# -----------------------------------------------------------------
+# python wrappers
+# -----------------------------------------------------------------
+def mean(np.ndarray[np.uint16_t, ndim=2] image,
+         np.ndarray[np.uint8_t, ndim=2] selem,
+         np.ndarray[np.uint8_t, ndim=2] mask=None,
+         np.ndarray[np.uint16_t, ndim=2] out=None,
+         char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
+    """average gray level (clipped on uint8)
+    """
+    return _core16(kernel_mean, image, selem, mask, out, shift_x, shift_y, bitdepth, 0., 0., s0, s1)
+
+
+def pop(np.ndarray[np.uint16_t, ndim=2] image,
+        np.ndarray[np.uint8_t, ndim=2] selem,
+        np.ndarray[np.uint8_t, ndim=2] mask=None,
+        np.ndarray[np.uint16_t, ndim=2] out=None,
+        char shift_x=0, char shift_y=0, int bitdepth=8, int s0=1, int s1=1):
+    """returns the number of actual pixels of the structuring element inside the mask
+    """
+    return _core16(kernel_pop, image, selem, mask, out, shift_x, shift_y, bitdepth, .0, .0, s0, s1)
diff --git a/skimage/filter/rank/_crank16_percentiles.pyx b/skimage/filter/rank/_crank16_percentiles.pyx
new file mode 100644
index 00000000..0d37b77c
--- /dev/null
+++ b/skimage/filter/rank/_crank16_percentiles.pyx
@@ -0,0 +1,297 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+
+import numpy as np
+cimport numpy as np
+
+# import main loop
+from skimage.filter.rank._core16 cimport _core16, int_min, int_max
+
+# -----------------------------------------------------------------
+# kernels uint16 (SOFT version using percentiles)
+# -----------------------------------------------------------------
+
+cdef inline np.uint16_t kernel_autolevel(
+    Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
+        Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+
+    cdef int i, imin, imax, sum, delta
+
+    if pop:
+        sum = 0
+        p1 = 1.0 - p1
+        for i in range(maxbin):
+            sum += histo[i]
+            if sum > p0 * pop:
+                imin = i
+                break
+        sum = 0
+        for i in range(maxbin - 1, -1, -1):
+            sum += histo[i]
+            if sum > p1 * pop:
+                imax = i
+                break
+
+        delta = imax - imin
+        if delta > 0:
+            return < np.uint16_t > (1.0 * (maxbin - 1) * (int_min(int_max(imin, g), imax) - imin) / delta)
+        else:
+            return < np.uint16_t > (imax - imin)
+    else:
+        return < np.uint16_t > (0)
+
+
+cdef inline np.uint16_t kernel_gradient(
+    Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
+        Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+
+    cdef int i, imin, imax, sum, delta
+
+    if pop:
+        sum = 0
+        p1 = 1.0 - p1
+        for i in range(maxbin):
+            sum += histo[i]
+            if sum >= p0 * pop:
+                imin = i
+                break
+        sum = 0
+        for i in range((maxbin - 1), -1, -1):
+            sum += histo[i]
+            if sum >= p1 * pop:
+                imax = i
+                break
+
+        return < np.uint16_t > (imax - imin)
+    else:
+        return < np.uint16_t > (0)
+
+
+cdef inline np.uint16_t kernel_mean(
+    Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
+        Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+
+    cdef int i, sum, mean, n
+
+    if pop:
+        sum = 0
+        mean = 0
+        n = 0
+        for i in range(maxbin):
+            sum += histo[i]
+            if (sum >= p0 * pop) and (sum <= p1 * pop):
+                n += histo[i]
+                mean += histo[i] * i
+
+        if n > 0:
+            return < np.uint16_t > (1.0 * mean / n)
+        else:
+            return < np.uint16_t > (0)
+    else:
+        return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_mean_substraction(
+    Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
+        Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+
+    cdef int i, sum, mean, n
+
+    if pop:
+        sum = 0
+        mean = 0
+        n = 0
+        for i in range(maxbin):
+            sum += histo[i]
+            if (sum >= p0 * pop) and (sum <= p1 * pop):
+                n += histo[i]
+                mean += histo[i] * i
+        if n > 0:
+            return < np.uint16_t > ((g - (mean / n)) * .5 + midbin)
+        else:
+            return < np.uint16_t > (0)
+    else:
+        return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_morph_contr_enh(
+    Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
+        Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+
+    cdef int i, imin, imax, sum, delta
+
+    if pop:
+        sum = 0
+        p1 = 1.0 - p1
+        for i in range(maxbin):
+            sum += histo[i]
+            if sum > p0 * pop:
+                imin = i
+                break
+        sum = 0
+        for i in range((maxbin - 1), -1, -1):
+            sum += histo[i]
+            if sum > p1 * pop:
+                imax = i
+                break
+        if g > imax:
+            return < np.uint16_t > imax
+        if g < imin:
+            return < np.uint16_t > imin
+        if imax - g < g - imin:
+            return < np.uint16_t > imax
+        else:
+            return < np.uint16_t > imin
+    else:
+        return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_percentile(
+    Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
+        Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+
+    cdef int i
+    cdef float sum = 0.
+
+    if pop:
+        for i in range(maxbin):
+            sum += histo[i]
+            if sum >= p0 * pop:
+                break
+
+        return < np.uint16_t > (i)
+    else:
+        return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_pop(
+    Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
+        Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+
+    cdef int i, sum, n
+
+    if pop:
+        sum = 0
+        n = 0
+        for i in range(maxbin):
+            sum += histo[i]
+            if (sum >= p0 * pop) and (sum <= p1 * pop):
+                n += histo[i]
+        return < np.uint16_t > (n)
+    else:
+        return < np.uint16_t > (0)
+
+cdef inline np.uint16_t kernel_threshold(
+    Py_ssize_t * histo, float pop, np.uint16_t g, Py_ssize_t bitdepth,
+        Py_ssize_t maxbin, Py_ssize_t midbin, float p0, float p1, Py_ssize_t s0, Py_ssize_t s1):
+
+    cdef int i
+    cdef float sum = 0.
+
+    if pop:
+        for i in range(maxbin):
+            sum += histo[i]
+            if sum >= p0 * pop:
+                break
+
+        return < np.uint16_t > ((maxbin - 1) * (g >= i))
+    else:
+        return < np.uint16_t > (0)
+
+# -----------------------------------------------------------------
+# python wrappers
+# -----------------------------------------------------------------
+
+
+def autolevel(np.ndarray[np.uint16_t, ndim=2] image,
+              np.ndarray[np.uint8_t, ndim=2] selem,
+              np.ndarray[np.uint8_t, ndim=2] mask=None,
+              np.ndarray[np.uint16_t, ndim=2] out=None,
+              char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
+    """bottom hat
+    """
+    return _core16(
+        kernel_autolevel, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1,
+        < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def gradient(np.ndarray[np.uint16_t, ndim=2] image,
+             np.ndarray[np.uint8_t, ndim=2] selem,
+             np.ndarray[np.uint8_t, ndim=2] mask=None,
+             np.ndarray[np.uint16_t, ndim=2] out=None,
+             char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
+    """return p0,p1 percentile gradient
+    """
+    return _core16(
+        kernel_gradient, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def mean(np.ndarray[np.uint16_t, ndim=2] image,
+         np.ndarray[np.uint8_t, ndim=2] selem,
+         np.ndarray[np.uint8_t, ndim=2] mask=None,
+         np.ndarray[np.uint16_t, ndim=2] out=None,
+         char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
+    """return mean between [p0 and p1] percentiles
+    """
+    return _core16(
+        kernel_mean, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def mean_substraction(np.ndarray[np.uint16_t, ndim=2] image,
+                      np.ndarray[np.uint8_t, ndim=2] selem,
+                      np.ndarray[np.uint8_t, ndim=2] mask=None,
+                      np.ndarray[np.uint16_t, ndim=2] out=None,
+                      char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
+    """return original - mean between [p0 and p1] percentiles *.5 +127
+    """
+    return _core16(
+        kernel_mean_substraction, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1,
+        < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def morph_contr_enh(np.ndarray[np.uint16_t, ndim=2] image,
+                    np.ndarray[np.uint8_t, ndim=2] selem,
+                    np.ndarray[np.uint8_t, ndim=2] mask=None,
+                    np.ndarray[np.uint16_t, ndim=2] out=None,
+                    char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
+    """reforce contrast using percentiles
+    """
+    return _core16(
+        kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1,
+        < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def percentile(np.ndarray[np.uint16_t, ndim=2] image,
+               np.ndarray[np.uint8_t, ndim=2] selem,
+               np.ndarray[np.uint8_t, ndim=2] mask=None,
+               np.ndarray[np.uint16_t, ndim=2] out=None,
+               char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
+    """return p0 percentile
+    """
+    return _core16(
+        kernel_percentile, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1,
+        < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def pop(np.ndarray[np.uint16_t, ndim=2] image,
+        np.ndarray[np.uint8_t, ndim=2] selem,
+        np.ndarray[np.uint8_t, ndim=2] mask=None,
+        np.ndarray[np.uint16_t, ndim=2] out=None,
+        char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
+    """return nb of pixels between [p0 and p1]
+    """
+    return _core16(
+        kernel_pop, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1,
+        < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def threshold(np.ndarray[np.uint16_t, ndim=2] image,
+              np.ndarray[np.uint8_t, ndim=2] selem,
+              np.ndarray[np.uint8_t, ndim=2] mask=None,
+              np.ndarray[np.uint16_t, ndim=2] out=None,
+              char shift_x=0, char shift_y=0, int bitdepth=8, float p0=0., float p1=0.):
+    """return (maxbin-1) if g > percentile p0
+    """
+    return _core16(
+        kernel_threshold, image, selem, mask, out, shift_x, shift_y, bitdepth, p0, p1,
+        < Py_ssize_t > 0, < Py_ssize_t > 0)
diff --git a/skimage/filter/rank/_crank8.pyx b/skimage/filter/rank/_crank8.pyx
new file mode 100644
index 00000000..4efd2bb0
--- /dev/null
+++ b/skimage/filter/rank/_crank8.pyx
@@ -0,0 +1,458 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+
+import numpy as np
+cimport numpy as np
+
+from libc.math cimport log2
+
+# import main loop
+from skimage.filter.rank._core8 cimport _core8
+
+# -----------------------------------------------------------------
+# kernels uint8
+# -----------------------------------------------------------------
+
+cdef inline np.uint8_t kernel_autolevel(
+    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, imin, imax, delta
+
+    if pop:
+        for i in range(255, -1, -1):
+            if histo[i]:
+                imax = i
+                break
+        for i in range(256):
+            if histo[i]:
+                imin = i
+                break
+        delta = imax - imin
+        if delta > 0:
+            return < np.uint8_t > (255. * (g - imin) / delta)
+        else:
+            return < np.uint8_t > (imax - imin)
+    else:
+        return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_bottomhat(
+    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
+
+    for i in range(256):
+        if histo[i]:
+            break
+
+    return < np.uint8_t > (g - i)
+
+
+cdef inline np.uint8_t kernel_equalize(
+    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 float sum = 0.
+
+    if pop:
+        for i in range(256):
+            sum += histo[i]
+            if i >= g:
+                break
+
+        return < np.uint8_t > ((255 * sum) / pop)
+    else:
+        return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_gradient(
+    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, imin, imax
+
+    if pop:
+        for i in range(255, -1, -1):
+            if histo[i]:
+                imax = i
+                break
+        for i in range(256):
+            if histo[i]:
+                imin = i
+                break
+        return < np.uint8_t > (imax - imin)
+    else:
+        return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_maximum(
+    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
+
+    if pop:
+        for i in range(255, -1, -1):
+            if histo[i]:
+                return < np.uint8_t > (i)
+
+    return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_mean(
+    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 float mean = 0.
+
+    if pop:
+        for i in range(256):
+            mean += histo[i] * i
+        return < np.uint8_t > (mean / pop)
+    else:
+        return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_meansubstraction(
+    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 float mean = 0.
+
+    if pop:
+        for i in range(256):
+            mean += histo[i] * i
+        return < np.uint8_t > ((g - mean / pop) / 2. + 127)
+    else:
+        return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_median(
+    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 float sum = pop / 2.0
+
+    if pop:
+        for i in range(256):
+            if histo[i]:
+                sum -= histo[i]
+                if sum < 0:
+                    return < np.uint8_t > (i)
+
+    return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_minimum(
+    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
+
+    if pop:
+        for i in range(256):
+            if histo[i]:
+                return < np.uint8_t > (i)
+
+    return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_modal(
+        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 hmax = 0, imax = 0
+
+    if pop:
+        for i in range(256):
+            if histo[i] > hmax:
+                hmax = histo[i]
+                imax = i
+        return < np.uint8_t > (imax)
+
+    return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_morph_contr_enh(
+    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, imin, imax
+
+    if pop:
+        for i in range(255, -1, -1):
+            if histo[i]:
+                imax = i
+                break
+        for i in range(256):
+            if histo[i]:
+                imin = i
+                break
+        if imax - g < g - imin:
+            return < np.uint8_t > (imax)
+        else:
+            return < np.uint8_t > (imin)
+    else:
+        return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_pop(
+    Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
+        Py_ssize_t s1):
+
+    return < np.uint8_t > (pop)
+
+cdef inline np.uint8_t kernel_threshold(
+    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 float mean = 0.
+
+    if pop:
+        for i in range(256):
+            mean += histo[i] * i
+        return < np.uint8_t > (g > (mean / pop))
+    else:
+        return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_tophat(
+    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
+
+    for i in range(255, -1, -1):
+        if histo[i]:
+            break
+
+    return < np.uint8_t > (i - g)
+
+cdef inline np.uint8_t kernel_noise_filter(
+        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 min_i
+
+    # 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 = g-i
+    for i in range(g, 256):
+        if histo[i]:
+            break
+    if i-g < min_i:
+        return < np.uint8_t > (i-g)
+    else:
+        return < np.uint8_t > min_i
+
+cdef inline np.uint8_t kernel_entropy(
+        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 float e,p
+
+    e = 0.
+
+    for i in range(256):
+        p = histo[i]/pop
+        if p>0:
+            e -= p*log2(p)
+
+    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
+# -----------------------------------------------------------------
+
+
+def autolevel(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_autolevel, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def bottomhat(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_bottomhat, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def equalize(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_equalize, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def gradient(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_gradient, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def maximum(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_maximum, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def mean(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_mean, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def meansubstraction(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_meansubstraction, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def median(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_median, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def minimum(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_minimum, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def morph_contr_enh(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_morph_contr_enh, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def modal(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_modal, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def pop(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_pop, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def threshold(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_threshold, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def tophat(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_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):
+    return _core8(kernel_noise_filter, image, selem, mask, out, shift_x, shift_y, .0, .0, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+def entropy(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_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/_crank8_percentiles.pyx b/skimage/filter/rank/_crank8_percentiles.pyx
new file mode 100644
index 00000000..618a6452
--- /dev/null
+++ b/skimage/filter/rank/_crank8_percentiles.pyx
@@ -0,0 +1,286 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+
+import numpy as np
+cimport numpy as np
+
+# import main loop
+from skimage.filter.rank._core8 cimport _core8, uint8_max, uint8_min
+
+# -----------------------------------------------------------------
+# kernels uint8 (SOFT version using percentiles)
+# -----------------------------------------------------------------
+
+cdef inline np.uint8_t kernel_autolevel(
+    Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
+        Py_ssize_t s1):
+    cdef int i, imin, imax, sum, delta
+
+    if pop:
+        sum = 0
+        p1 = 1.0 - p1
+        imin = 0
+        imax = 255
+
+        for i in range(256):
+            sum += histo[i]
+            if sum > (p0 * pop):
+                imin = i
+                break
+        sum = 0
+        for i in range(255, -1, -1):
+            sum += histo[i]
+            if sum > (p1 * pop):
+                imax = i
+                break
+        delta = imax - imin
+        if delta > 0:
+            return < np.uint8_t > (255 * (uint8_min(uint8_max(imin, g), imax) - imin) / delta)
+        else:
+            return < np.uint8_t > (imax - imin)
+    else:
+        return < np.uint8_t > (128)
+
+
+cdef inline np.uint8_t kernel_gradient(
+    Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
+        Py_ssize_t s1):
+    cdef int i, imin, imax, sum, delta
+
+    if pop:
+        sum = 0
+        p1 = 1.0 - p1
+        for i in range(256):
+            sum += histo[i]
+            if sum >= p0 * pop:
+                imin = i
+                break
+        sum = 0
+        for i in range(255, -1, -1):
+            sum += histo[i]
+            if sum >= p1 * pop:
+                imax = i
+                break
+
+        return < np.uint8_t > (imax - imin)
+    else:
+        return < np.uint8_t > (0)
+
+
+cdef inline np.uint8_t kernel_mean(
+    Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
+        Py_ssize_t s1):
+    cdef int i, sum, mean, n
+
+    if pop:
+        sum = 0
+        mean = 0
+        n = 0
+        for i in range(256):
+            sum += histo[i]
+            if (sum >= p0 * pop) and (sum <= p1 * pop):
+                n += histo[i]
+                mean += histo[i] * i
+        if n > 0:
+            return < np.uint8_t > (1.0 * mean / n)
+        else:
+            return < np.uint8_t > (0)
+    else:
+        return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_mean_substraction(
+    Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1,
+        Py_ssize_t s0, Py_ssize_t s1):
+    cdef int i, sum, mean, n
+
+    if pop:
+        sum = 0
+        mean = 0
+        n = 0
+        for i in range(256):
+            sum += histo[i]
+            if (sum >= p0 * pop) and (sum <= p1 * pop):
+                n += histo[i]
+                mean += histo[i] * i
+        if n > 0:
+            return < np.uint8_t > ((g - (mean / n)) * .5 + 127)
+        else:
+            return < np.uint8_t > (0)
+    else:
+        return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_morph_contr_enh(
+    Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1,
+        Py_ssize_t s0, Py_ssize_t s1):
+    cdef int i, imin, imax, sum, delta
+
+    if pop:
+        sum = 0
+        p1 = 1.0 - p1
+        for i in range(256):
+            sum += histo[i]
+            if sum >= p0 * pop:
+                imin = i
+                break
+        sum = 0
+        for i in range(255, -1, -1):
+            sum += histo[i]
+            if sum >= p1 * pop:
+                imax = i
+                break
+        if g > imax:
+            return < np.uint8_t > imax
+        if g < imin:
+            return < np.uint8_t > imin
+        if imax - g < g - imin:
+            return < np.uint8_t > imax
+        else:
+            return < np.uint8_t > imin
+    else:
+        return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_percentile(
+    Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1,
+        Py_ssize_t s0, Py_ssize_t s1):
+    cdef int i
+    cdef float sum = 0.
+
+    if pop:
+        for i in range(256):
+            sum += histo[i]
+            if sum >= p0 * pop:
+                break
+
+        return < np.uint8_t > (i)
+    else:
+        return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_pop(
+    Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1,
+        Py_ssize_t s0, Py_ssize_t s1):
+    cdef int i, sum, n
+
+    if pop:
+        sum = 0
+        n = 0
+        for i in range(256):
+            sum += histo[i]
+            if (sum >= p0 * pop) and (sum <= p1 * pop):
+                n += histo[i]
+        return < np.uint8_t > (n)
+    else:
+        return < np.uint8_t > (0)
+
+cdef inline np.uint8_t kernel_threshold(
+    Py_ssize_t * histo, float pop, np.uint8_t g, float p0, float p1, Py_ssize_t s0,
+        Py_ssize_t s1):
+    cdef int i
+    cdef float sum = 0.
+
+    if pop:
+        for i in range(256):
+            sum += histo[i]
+            if sum >= p0 * pop:
+                break
+
+        return < np.uint8_t > (255 * (g >= i))
+    else:
+        return < np.uint8_t > (0)
+
+# -----------------------------------------------------------------
+# python wrappers
+# -----------------------------------------------------------------
+
+
+def autolevel(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, float p0=0., float p1=0.):
+    """autolevel
+    """
+    return _core8(
+        kernel_autolevel, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def gradient(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, float p0=0., float p1=0.):
+    """return p0,p1 percentile gradient
+    """
+    return _core8(
+        kernel_gradient, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def mean(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, float p0=0., float p1=0.):
+    """return mean between [p0 and p1] percentiles
+    """
+    return _core8(kernel_mean, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def mean_substraction(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, float p0=0., float p1=0.):
+    """return original - mean between [p0 and p1] percentiles *.5 +127
+    """
+    return _core8(
+        kernel_mean_substraction, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def morph_contr_enh(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, float p0=0., float p1=0.):
+    """reforce contrast using percentiles
+    """
+    return _core8(
+        kernel_morph_contr_enh, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def percentile(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, float p0=0., float p1=0.):
+    """return p0 percentile
+    """
+    return _core8(
+        kernel_percentile, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
+
+
+def pop(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, float p0=0., float p1=0.):
+    """return nb of pixels between [p0 and p1]
+    """
+    return _core8(kernel_pop, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0, < Py_ssize_t > 0)
+
+
+def threshold(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, float p0=0., float p1=0.):
+    """return 255 if g > percentile p0
+    """
+    return _core8(
+        kernel_threshold, image, selem, mask, out, shift_x, shift_y, p0, p1, < Py_ssize_t > 0,
+        < Py_ssize_t > 0)
diff --git a/skimage/filter/rank/bilateral_rank.pyx b/skimage/filter/rank/bilateral_rank.pyx
new file mode 100644
index 00000000..f181735b
--- /dev/null
+++ b/skimage/filter/rank/bilateral_rank.pyx
@@ -0,0 +1,165 @@
+"""bilateral_rank.py - approximate bilateral rankfilter for local (custom kernel) mean
+
+The local histogram is computed using a sliding window similar to the method described in
+
+Reference: Huang, T. ,Yang, G. ;  Tang, G.. "A fast two-dimensional median filtering algorithm",
+IEEE Transactions on Acoustics, Speech and Signal Processing, Feb 1979. Volume: 27 , Issue: 1, Page(s): 13 - 18.
+
+input image can be 8 bit or 16 bit with a value < 4096 (i.e. 12 bit),
+8 bit images are casted in 16 bit
+the number of histogram bins is determined from the maximum value present in the image
+
+The pixel neighborhood is defined by:
+
+* the given structuring element
+
+* an interval [g-s0,g+s1] in gray level around g the processed pixel gray level
+
+The kernel is flat (i.e. each pixel belonging to the neighborhood contributes equally)
+
+result image is 16 bit with respect to the input image
+
+"""
+
+from skimage import img_as_ubyte
+
+import numpy as np
+from skimage.filter.rank import _crank16_bilateral
+
+from skimage.filter.rank.generic import find_bitdepth
+
+
+__all__ = ['bilateral_mean', 'bilateral_pop']
+
+
+def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y, s0, s1):
+    selem = img_as_ubyte(selem)
+    if mask is not None:
+        mask = img_as_ubyte(mask)
+    if image.dtype == np.uint8:
+        image = image.astype(np.uint16)
+    elif image.dtype == np.uint16:
+        pass
+    else:
+        raise TypeError("only uint8 and uint16 image supported!")
+    bitdepth = find_bitdepth(image)
+    if bitdepth > 11:
+        raise ValueError("only uint16 <4096 image (12bit) supported!")
+    return func16(
+        image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, bitdepth=bitdepth + 1, out=out,
+        s0=s0, s1=s1)
+
+
+def bilateral_mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False, s0=10, s1=10):
+    """Apply a flat kernel bilateral filter.
+
+    This is an edge-preserving and noise reducing denoising filter. It averages
+    pixels based on their spatial closeness and radiometric similarity.
+
+    Spatial closeness is measured by considering only the local pixel neighborhood given by a
+    structuring element (selem).
+
+    Radiometric similarity is defined by the gray level interval [g-s0,g+s1] where g is the current pixel gray level.
+    Only pixels belonging to the structuring element AND having a gray level inside this interval are averaged.
+    Return greyscale local bilateral_mean of an image.
+
+    Parameters
+    ----------
+    image : ndarray
+        Image array (uint8 array or uint16). If image is uint16, as 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).
+    s0, s1 : int
+        define the [s0,s1] interval to be considered for computing the value.
+
+    Returns
+    -------
+    out : uint16 array (uint8 image are casted to uint16)
+        The result of the local bilateral mean.
+
+    See also
+    --------
+    skimage.filter.denoise_bilateral() for a gaussian bilateral filter.
+
+    Notes
+    -----
+
+    * input image can be 8 bit or 16 bit with a value < 4096 (i.e. 12 bit)
+
+    * 8 bit images are casted in 16 bit
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> from skimage.morphology import disk
+    >>> from skimage.filter.rank import bilateral_mean
+    >>> # Load test image
+    >>> ima = data.camera()
+    >>> # bilateral filtering of cameraman image using a flat kernel
+    >>> bilat_ima = bilateral_mean(ima, disk(20), s0=10,s1=10)
+    """
+
+    return _apply(
+        None, _crank16_bilateral.mean, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y,
+        s0=s0, s1=s1)
+
+
+def bilateral_pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False, s0=10, s1=10):
+    """Return the number (population) of pixels actually inside the bilateral neighborhood,
+    i.e. being inside the structuring element AND having a gray level inside the interval [g-s0,g+s1].
+
+    Parameters
+    ----------
+    image : ndarray
+        Image array (uint8 array or uint16). If image is uint16, as 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).
+    s0, s1 : int
+        define the [s0,s1] interval to be considered for computing the value.
+
+    Returns
+    -------
+    out : uint16 array (uint8 image are casted to uint16)
+        the local number of pixels inside the bilateral neighborhood
+
+    Examples
+    --------
+    >>> # Local mean
+    >>> from skimage.morphology import square
+    >>> import skimage.filter.rank as rank
+    >>> ima8 = 255*np.array([[0, 0, 0, 0, 0],
+    ...                           [0, 1, 1, 1, 0],
+    ...                           [0, 1, 1, 1, 0],
+    ...                           [0, 1, 1, 1, 0],
+    ...                           [0, 0, 0, 0, 0]], dtype=np.uint8)
+    >>> rank.bilateral_pop(ima8, square(3), s0=10,s1=10)
+    array([[3, 4, 3, 4, 3],
+           [4, 4, 6, 4, 4],
+           [3, 6, 9, 6, 3],
+           [4, 4, 6, 4, 4],
+           [3, 4, 3, 4, 3]], dtype=uint16)
+
+    """
+
+    return _apply(
+        None, _crank16_bilateral.pop, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y,
+        s0=s0, s1=s1)
+
diff --git a/skimage/filter/rank/demo/demo_all.py b/skimage/filter/rank/demo/demo_all.py
new file mode 100644
index 00000000..038c749b
--- /dev/null
+++ b/skimage/filter/rank/demo/demo_all.py
@@ -0,0 +1,77 @@
+import matplotlib.pyplot as plt
+from pprint import pprint
+
+from skimage import data
+from skimage.morphology.selem import disk
+import skimage.filter.rank as rank
+
+def plot_all():
+    a8 = data.camera()
+    a16 = a8.astype('uint16')*16
+    selem = disk(5)
+
+    name_list = sorted([n for n in dir(rank) if n[0] is  not '_'])
+    print name_list
+
+    for n in name_list:
+        if n.rfind('bilateral')==0:
+            print n
+            method = eval('rank.%s'%n)
+            if type(method) == type(rank.maximum):
+                print method
+                f8 = method(a8,selem = selem,s0=10,s1=10)
+                f16 = method(a16,selem = selem,s0=10,s1=10)
+                plt.figure()
+                plt.subplot(2,2,1)
+                plt.imshow(a8)
+                plt.colorbar()
+                plt.subplot(2,2,2)
+                plt.imshow(f8)
+                plt.colorbar()
+                plt.subplot(2,2,3)
+                plt.imshow(f16)
+                plt.colorbar()
+                plt.title(method)
+    for n in name_list:
+        if n.rfind('percentile')==0:
+            print n
+            method = eval('rank.%s'%n)
+            if type(method) == type(rank.maximum):
+                print method
+                f8 = method(a8,selem = selem,p0=.1,p1=.9)
+                f16 = method(a16,selem = selem,p0=.1,p1=.9)
+                plt.figure()
+                plt.subplot(2,2,1)
+                plt.imshow(a8)
+                plt.colorbar()
+                plt.subplot(2,2,2)
+                plt.imshow(f8)
+                plt.colorbar()
+                plt.subplot(2,2,3)
+                plt.imshow(f16)
+                plt.colorbar()
+                plt.title(method)
+    for n in name_list:
+        if n.find('percentile')==-1 and n.find('bilateral')==-1:
+            print n
+            method = eval('rank.%s'%n)
+            if type(method) == type(rank.maximum):
+                print method
+                f8 = method(a8,selem = selem)
+                f16 = method(a16,selem = selem)
+                plt.figure()
+                plt.subplot(2,2,1)
+                plt.imshow(a8)
+                plt.colorbar()
+                plt.subplot(2,2,2)
+                plt.imshow(f8)
+                plt.colorbar()
+                plt.subplot(2,2,3)
+                plt.imshow(f16)
+                plt.colorbar()
+                plt.title(method)
+    plt.show()
+
+if __name__ == '__main__':
+    plot_all()
+    pprint(dir(rank))
\ No newline at end of file
diff --git a/skimage/filter/rank/demo/demo_single.py b/skimage/filter/rank/demo/demo_single.py
new file mode 100644
index 00000000..7e4ef405
--- /dev/null
+++ b/skimage/filter/rank/demo/demo_single.py
@@ -0,0 +1,32 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage import data
+from skimage.morphology.selem import disk
+import skimage.filter.rank as rank
+
+from skimage.filter import denoise_bilateral
+
+if __name__ == '__main__':
+    a8 = data.camera()
+    a16 = data.camera().astype(np.uint16)*4
+
+    p8 = data.page()
+
+    selem = disk(20)
+
+    otsu = rank.otsu(p8,selem)
+
+
+    plt.figure()
+    plt.subplot(1,2,1)
+    plt.imshow(p8)
+    plt.colorbar()
+    plt.subplot(1,2,2)
+    plt.imshow(otsu)
+    plt.colorbar()
+    plt.show()
+
+
+
+
diff --git a/skimage/filter/rank/demo/iko_pan_Ja1.tif b/skimage/filter/rank/demo/iko_pan_Ja1.tif
new file mode 100644
index 00000000..47201695
Binary files /dev/null and b/skimage/filter/rank/demo/iko_pan_Ja1.tif differ
diff --git a/skimage/filter/rank/demo/test_morph_contr_enh.py b/skimage/filter/rank/demo/test_morph_contr_enh.py
new file mode 100644
index 00000000..f2f0f7c9
--- /dev/null
+++ b/skimage/filter/rank/demo/test_morph_contr_enh.py
@@ -0,0 +1,28 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import gdal
+
+from skimage.morphology import disk
+import skimage.filter.rank as rank
+
+filename = 'iko_pan_Ja1.tif'
+im16 = gdal.Open(filename).ReadAsArray().astype(np.uint16)
+
+plt.figure()
+plt.imshow(im16,cmap=plt.cm.gray)
+plt.colorbar()
+
+f0 = rank.median(im16,disk(1))
+f1 = rank.bilateral_mean(im16,disk(20),s0=200,s1=200)
+f2 = rank.equalize(f1,disk(10))
+f3 = rank.bottomhat(f1,disk(1))
+
+plt.figure()
+plt.imshow(f2,cmap=plt.cm.gray,interpolation='nearest')
+plt.colorbar()
+
+plt.show()
+
+
+
+
diff --git a/skimage/filter/rank/demo/test_rank.py b/skimage/filter/rank/demo/test_rank.py
new file mode 100644
index 00000000..09cfdcd5
--- /dev/null
+++ b/skimage/filter/rank/demo/test_rank.py
@@ -0,0 +1,42 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from skimage import data
+from skimage.morphology.selem import disk
+import skimage.filter.rank as rank
+
+print dir(rank)
+
+print rank.mean
+print rank.percentile_mean
+print rank.bilateral_mean
+
+a8 = data.camera()
+a16 = a8.astype('uint16')*16
+selem = disk(10)
+
+f8 = rank.mean(a8,selem)
+f16 = rank.mean(a16,selem)
+
+plt.figure()
+plt.imshow(np.hstack((a8,f8)))
+plt.colorbar()
+plt.figure()
+plt.imshow(np.hstack((a16,f16)))
+plt.colorbar()
+
+f8 = rank.percentile_mean(a8,selem,p0=.1,p1=.9)
+f16 = rank.percentile_mean(a16,selem,p0=.1,p1=.9)
+
+plt.figure()
+plt.imshow(np.hstack((a8,f8)))
+plt.colorbar()
+plt.figure()
+plt.imshow(np.hstack((a16,f16)))
+plt.colorbar()
+
+plt.show()
+
+
+
+
diff --git a/skimage/filter/rank/generic.py b/skimage/filter/rank/generic.py
new file mode 100644
index 00000000..94fc3130
--- /dev/null
+++ b/skimage/filter/rank/generic.py
@@ -0,0 +1,11 @@
+import numpy as np
+
+
+def find_bitdepth(image):
+    """returns the max bith depth of a uint16 image
+    """
+    umax = np.max(image)
+    if umax > 2:
+        return int(np.log2(umax))
+    else:
+        return 1
diff --git a/skimage/filter/rank/percentile_rank.pyx b/skimage/filter/rank/percentile_rank.pyx
new file mode 100644
index 00000000..1bd89eb8
--- /dev/null
+++ b/skimage/filter/rank/percentile_rank.pyx
@@ -0,0 +1,337 @@
+"""percentile_rank.py - inferior and superior ranks, provided by the user, are passed to the kernel function
+to provide a softer version of the rank filters. E.g. percentile_autolevel will stretch image levels between
+percentile [p0,p1] instead of using [min,max]. It means that isolate bright or dark pixels will not produce halos.
+
+The local histogram is computed using a sliding window similar to the method described in
+
+Reference: Huang, T. ,Yang, G. ;  Tang, G.. "A fast two-dimensional median filtering algorithm",
+IEEE Transactions on Acoustics, Speech and Signal Processing, Feb 1979. Volume: 27 , Issue: 1, Page(s): 13 - 18.
+
+input image can be 8 bit or 16 bit with a value < 4096 (i.e. 12 bit),
+for 16 bit input images, the number of histogram bins is determined from the maximum value present in the image
+
+result image is 8 or 16 bit with respect to the input image
+
+"""
+
+from skimage import img_as_ubyte
+import numpy as np
+
+from skimage.filter.rank.generic import find_bitdepth
+from skimage.filter.rank import _crank16_percentiles, _crank8_percentiles
+
+__all__ = ['percentile_autolevel', 'percentile_gradient',
+           'percentile_mean', 'percentile_mean_substraction',
+           'percentile_morph_contr_enh', 'percentile', 'percentile_pop', 'percentile_threshold']
+
+
+def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y, p0, p1):
+    selem = img_as_ubyte(selem)
+    if mask is not None:
+        mask = img_as_ubyte(mask)
+    if image.dtype == np.uint8:
+        return func8(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, out=out, p0=p0, p1=p1)
+    elif image.dtype == np.uint16:
+        bitdepth = find_bitdepth(image)
+        if bitdepth > 11:
+            raise ValueError("only uint16 <4096 image (12bit) supported!")
+        return func16(
+            image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, bitdepth=bitdepth + 1, out=out,
+            p0=p0, p1=p1)
+    else:
+        raise TypeError("only uint8 and uint16 image supported!")
+
+
+def percentile_autolevel(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
+    """Return greyscale local autolevel of an image.
+
+    Autolevel is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
+
+    Parameters
+    ----------
+    image : ndarray
+        Image array (uint8 array or uint16). If image is uint16, as 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).
+    p0, p1 : float in [0.,...,1.]
+        define the [p0,p1] percentile interval to be considered for computing the value.
+
+    Returns
+    -------
+    local autolevel : uint8 array or uint16 array depending on input image
+        The result of the local autolevel.
+
+
+    """
+
+    return _apply(
+        _crank8_percentiles.autolevel, _crank16_percentiles.autolevel, image, selem, out=out, mask=mask,
+        shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
+
+
+def percentile_gradient(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
+    """Return greyscale local percentile_gradient of an image.
+
+    percentile_gradient is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
+
+    Parameters
+    ----------
+    image : ndarray
+        Image array (uint8 array or uint16). If image is uint16, as 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).
+    p0, p1 : float in [0.,...,1.]
+        define the [p0,p1] percentile interval to be considered for computing the value.
+
+    Returns
+    -------
+    local percentile_gradient : uint8 array or uint16 array depending on input image
+        The result of the local percentile_gradient.
+
+
+
+    """
+
+    return _apply(
+        _crank8_percentiles.gradient, _crank16_percentiles.gradient, image, selem, out=out, mask=mask,
+        shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
+
+
+def percentile_mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
+    """Return greyscale local mean of an image.
+
+    Mean is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
+
+    Parameters
+    ----------
+    image : ndarray
+        Image array (uint8 array or uint16). If image is uint16, as 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).
+    p0, p1 : float in [0.,...,1.]
+        define the [p0,p1] percentile interval to be considered for computing the value.
+
+    Returns
+    -------
+    local mean : uint8 array or uint16 array depending on input image
+        The result of the local mean.
+
+
+    """
+
+    return _apply(
+        _crank8_percentiles.mean, _crank16_percentiles.mean, image, selem, out=out, mask=mask,
+        shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
+
+
+def percentile_mean_substraction(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
+    """Return greyscale local mean_substraction of an image.
+
+    mean_substraction is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
+
+    Parameters
+    ----------
+    image : ndarray
+        Image array (uint8 array or uint16). If image is uint16, as 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).
+    p0, p1 : float in [0.,...,1.]
+        define the [p0,p1] percentile interval to be considered for computing the value.
+
+    Returns
+    -------
+    local mean_substraction : uint8 array or uint16 array depending on input image
+        The result of the local mean_substraction.
+
+
+
+    """
+
+    return _apply(
+        _crank8_percentiles.mean_substraction, _crank16_percentiles.mean_substraction, image, selem, out=out,
+        mask=mask, shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
+
+
+def percentile_morph_contr_enh(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
+    """Return greyscale local morph_contr_enh of an image.
+
+    morph_contr_enh is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
+
+    Parameters
+    ----------
+    image : ndarray
+        Image array (uint8 array or uint16). If image is uint16, as 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).
+    p0, p1 : float in [0.,...,1.]
+        define the [p0,p1] percentile interval to be considered for computing the value.
+
+    Returns
+    -------
+    local morph_contr_enh : uint8 array or uint16 array depending on input image
+        The result of the local morph_contr_enh.
+
+
+
+    """
+
+    return _apply(
+        _crank8_percentiles.morph_contr_enh, _crank16_percentiles.morph_contr_enh, image, selem, out=out,
+        mask=mask, shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
+
+
+def percentile(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
+    """Return greyscale local percentile of an image.
+
+    percentile is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
+
+    Parameters
+    ----------
+    image : ndarray
+        Image array (uint8 array or uint16). If image is uint16, as 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).
+    p0, p1 : float in [0.,...,1.]
+        define the [p0,p1] percentile interval to be considered for computing the value.
+
+    Returns
+    -------
+    local percentile : uint8 array or uint16 array depending on input image
+        The result of the local percentile.
+
+
+
+    """
+
+    return _apply(
+        _crank8_percentiles.percentile, _crank16_percentiles.percentile, image, selem, out=out, mask=mask,
+        shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
+
+
+def percentile_pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
+    """Return greyscale local pop of an image.
+
+    pop is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
+
+    Parameters
+    ----------
+    image : ndarray
+        Image array (uint8 array or uint16). If image is uint16, as 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).
+    p0, p1 : float in [0.,...,1.]
+        define the [p0,p1] percentile interval to be considered for computing the value.
+
+    Returns
+    -------
+    local pop : uint8 array or uint16 array depending on input image
+        The result of the local pop.
+
+
+
+    """
+
+    return _apply(
+        _crank8_percentiles.pop, _crank16_percentiles.pop, image, selem, out=out, mask=mask, shift_x=shift_x,
+        shift_y=shift_y, p0=p0, p1=p1)
+
+
+def percentile_threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False, p0=.0, p1=1.):
+    """Return greyscale local threshold of an image.
+
+    threshold is computed on the given structuring element. Only levels between percentiles [p0,p1] ,are used.
+
+    Parameters
+    ----------
+    image : ndarray
+        Image array (uint8 array or uint16). If image is uint16, as 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).
+    p0, p1 : float in [0.,...,1.]
+        define the [p0,p1] percentile interval to be considered for computing the value.
+
+    Returns
+    -------
+    local threshold : uint8 array or uint16 array depending on input image
+        The result of the local threshold.
+
+
+    """
+
+    return _apply(
+        _crank8_percentiles.threshold, _crank16_percentiles.threshold, image, selem, out=out, mask=mask,
+        shift_x=shift_x, shift_y=shift_y, p0=p0, p1=p1)
+
+
diff --git a/skimage/filter/rank/rank.pyx b/skimage/filter/rank/rank.pyx
new file mode 100644
index 00000000..6a35920b
--- /dev/null
+++ b/skimage/filter/rank/rank.pyx
@@ -0,0 +1,703 @@
+"""rank.py - rankfilter for local (custom kernel) maximum, minimum, median, mean, auto-level, equalization, etc
+
+The local histogram is computed using a sliding window similar to the method described in
+
+Reference: Huang, T. ,Yang, G. ;  Tang, G.. "A fast two-dimensional median filtering algorithm",
+IEEE Transactions on Acoustics, Speech and Signal Processing, Feb 1979. Volume: 27 , Issue: 1, Page(s): 13 - 18.
+
+input image can be 8 bit or 16 bit with a value < 4096 (i.e. 12 bit),
+for 16 bit input images, the number of histogram bins is determined from the maximum value present in the image
+
+result image is 8 or 16 bit with respect to the input image
+
+"""
+
+from skimage import img_as_ubyte
+import numpy as np
+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','otsu']
+
+
+def _apply(func8, func16, image, selem, out, mask, shift_x, shift_y):
+    selem = img_as_ubyte(selem)
+    if mask is not None:
+        mask = img_as_ubyte(mask)
+    if image.dtype == np.uint8:
+        if func8 is None:
+            raise TypeError("not implemented for uint8 image")
+        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("not implemented for uint16 image")
+        bitdepth = find_bitdepth(image)
+        if bitdepth > 11:
+            raise ValueError("only uint16 <4096 image (12bit) supported!")
+        return func16(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask, bitdepth=bitdepth + 1, out=out)
+    else:
+        raise TypeError("only uint8 and uint16 image supported!")
+
+
+def autolevel(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Autolevel image using local histogram.
+
+    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)
+        The result of the local autolevel.
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> from skimage.morphology import disk
+    >>> from skimage.filter.rank import autolevel
+    >>> # Load test image
+    >>> ima = data.camera()
+    >>> # Stretch image contrast locally
+    >>> auto = autolevel(ima, disk(20))
+
+    """
+
+    return _apply(
+        _crank8.autolevel, _crank16.autolevel, image, selem, out=out, mask=mask, shift_x=shift_x,
+        shift_y=shift_y)
+
+
+def bottomhat(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Returns greyscale local bottomhat of an image.
+
+    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
+    -------
+    local bottomhat : uint8 array or uint16 array depending on input image
+        The result of the local bottomhat.
+
+
+    """
+
+    return _apply(
+        _crank8.bottomhat, _crank16.bottomhat, image, selem, out=out, mask=mask, shift_x=shift_x,
+        shift_y=shift_y)
+
+
+def equalize(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Equalize image using local histogram.
+
+    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)
+        The result of the local equalize.
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> from skimage.morphology import disk
+    >>> from skimage.filter.rank import equalize
+    >>> # Load test image
+    >>> ima = data.camera()
+    >>> # Local equalization
+    >>> equ = equalize(ima, disk(20))
+    """
+
+    return _apply(
+        _crank8.equalize, _crank16.equalize, image, selem, out=out, mask=mask, shift_x=shift_x,
+        shift_y=shift_y)
+
+
+def gradient(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Return greyscale local gradient of an image (i.e. local maximum - local minimum).
+
+
+    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)
+        The local gradient.
+
+    """
+
+    return _apply(
+        _crank8.gradient, _crank16.gradient, image, selem, out=out, mask=mask, shift_x=shift_x,
+        shift_y=shift_y)
+
+
+def maximum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Return greyscale local maximum of an image.
+
+
+    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)
+        The local maximum.
+
+    See also
+    --------
+    skimage.morphology.dilation()
+
+    Note
+    ----
+    * input image can be 8 bit or 16 bit with a value < 4096 (i.e. 12 bit)
+
+    * the lower algorithm complexity makes the rank.maximum() more efficient for larger images and structuring elements
+
+    """
+
+    return _apply(_crank8.maximum, _crank16.maximum, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+
+
+def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Return greyscale local mean of an image.
+
+    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)
+        The local mean.
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> from skimage.morphology import disk
+    >>> from skimage.filter.rank import mean
+    >>> # Load test image
+    >>> ima = data.camera()
+    >>> # Local mean
+    >>> avg = mean(ima, disk(20))
+    """
+
+    return _apply(_crank8.mean, _crank16.mean, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+
+
+def meansubstraction(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Return image substracted from its local mean.
+
+    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)
+        The result of the local meansubstraction.
+
+
+
+    """
+
+    return _apply(
+        _crank8.meansubstraction, _crank16.meansubstraction, image, selem, out=out, mask=mask,
+        shift_x=shift_x, shift_y=shift_y)
+
+
+def median(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Return greyscale local median of an image.
+
+
+    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)
+        The local median.
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> from skimage.morphology import disk
+    >>> from skimage.filter.rank import median
+    >>> # Load test image
+    >>> ima = data.camera()
+    >>> # Local mean
+    >>> avg = median(ima, disk(20))
+    """
+
+    return _apply(_crank8.median, _crank16.median, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+
+
+def minimum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Return greyscale local minimum of an image.
+
+    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)
+        The local minimum.
+
+    See also
+    --------
+    skimage.morphology.erosion()
+
+    Note
+    ----
+    * input image can be 8 bit or 16 bit with a value < 4096 (i.e. 12 bit)
+
+    * the lower algorithm complexity makes the rank.minimum() more efficient for larger images and structuring elements
+
+    """
+
+    return _apply(_crank8.minimum, _crank16.minimum, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+
+
+def modal(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Return greyscale local mode of an image.
+
+    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)
+        The local modal.
+
+
+    """
+
+    return _apply(_crank8.modal, _crank16.modal, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+
+
+def morph_contr_enh(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Enhance an image replacing each pixel by the local maximum if pixel graylevel is closest to maximimum
+     than local minimum OR local minimum otherwise.
+
+    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)
+        The result of the local morph_contr_enh.
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> from skimage.morphology import disk
+    >>> from skimage.filter.rank import morph_contr_enh
+    >>> # Load test image
+    >>> ima = data.camera()
+    >>> # Local mean
+    >>> avg = morph_contr_enh(ima, disk(20))
+    """
+
+    return _apply(
+        _crank8.morph_contr_enh, _crank16.morph_contr_enh, image, selem, out=out, mask=mask, shift_x=shift_x,
+        shift_y=shift_y)
+
+
+def pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Return the number (population) of pixels actually inside the neighborhood.
+
+    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)
+        The number of pixels belonging to the neighborhood.
+
+    Examples
+    --------
+    >>> # Local mean
+    >>> from skimage.morphology import square
+    >>> import skimage.filter.rank as rank
+    >>> ima = 255*np.array([[0, 0, 0, 0, 0],
+    ...                           [0, 1, 1, 1, 0],
+    ...                           [0, 1, 1, 1, 0],
+    ...                           [0, 1, 1, 1, 0],
+    ...                           [0, 0, 0, 0, 0]], dtype=np.uint8)
+    >>> rank.pop(ima, square(3))
+    array([[4, 6, 6, 6, 4],
+           [6, 9, 9, 9, 6],
+           [6, 9, 9, 9, 6],
+           [6, 9, 9, 9, 6],
+           [4, 6, 6, 6, 4]], dtype=uint8)
+
+
+    """
+
+    return _apply(_crank8.pop, _crank16.pop, image, selem, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+
+
+def threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Return greyscale local threshold of an image.
+
+    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)
+        The result of the local threshold.
+
+    Examples
+    --------
+    >>> # Local threshold
+    >>> from skimage.morphology import square
+    >>> from skimage.filter.rank import threshold
+    >>> ima = 255*np.array([[0, 0, 0, 0, 0],
+    ...                           [0, 1, 1, 1, 0],
+    ...                           [0, 1, 1, 1, 0],
+    ...                           [0, 1, 1, 1, 0],
+    ...                           [0, 0, 0, 0, 0]], dtype=np.uint8)
+    >>> threshold(ima, square(3))
+    array([[0, 0, 0, 0, 0],
+           [0, 1, 1, 1, 0],
+           [0, 1, 0, 1, 0],
+           [0, 1, 1, 1, 0],
+           [0, 0, 0, 0, 0]], dtype=uint8)
+
+
+    """
+
+    return _apply(
+        _crank8.threshold, _crank16.threshold, image, selem, out=out, mask=mask, shift_x=shift_x,
+        shift_y=shift_y)
+
+
+def tophat(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Return greyscale local tophat of an image.
+
+    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)
+        The image tophat.
+
+    """
+
+    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):
+    """Returns the noise feature as described in [Hashimoto12]_
+
+    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. Central element is removed during the filtering.
+    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).
+
+    Reference
+    ----------
+
+    .. [Hashimoto12] 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_cpy, out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+
+def entropy(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
+    """Returns the entropy [wiki_entropy]_ computed locally. Entropy is computed using base 2 logarithm i.e.
+    the filter returns the minimum number of bits needed to encode local greylevel distribution.
+
+    References
+    ----------
+    .. [wiki_entropy] http://en.wikipedia.org/wiki/Entropy_(information_theory)
+
+    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)
+        entropy x10 (uint8 images) and entropy x1000 (uint16 images)
+
+
+    Examples
+    --------
+
+    >>> # Local entropy
+    >>> from skimage import data
+    >>> from skimage.filter.rank import entropy
+    >>> from skimage.morphology import disk
+    >>> # defining a 8- and a 16-bit test images
+    >>> a8 = data.camera()
+    >>> a16 = data.camera().astype(np.uint16)*4
+    >>> ent8 = entropy(a8,disk(5)) # pixel value contain 10x the local entropy
+    >>> ent16 = entropy(a16,disk(5)) # pixel value contain 1000x the local entropy
+
+    """
+
+    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
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)
diff --git a/skimage/filter/rank/tests/test_suite.py b/skimage/filter/rank/tests/test_suite.py
new file mode 100644
index 00000000..66b6537c
--- /dev/null
+++ b/skimage/filter/rank/tests/test_suite.py
@@ -0,0 +1,191 @@
+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_random_sizes(self):
+        # make sure the size is not a problem
+
+        niter = 10
+        elem = np.asarray([[1,1,1],[1,1,1],[1,1,1]],dtype='uint8')
+        for m,n in np.random.random_integers(1,100,size=(10,2)):
+            a8 = np.ones((m,n),dtype='uint8')
+            r = _crank8.mean(image=a8,selem = elem,shift_x=0,shift_y=0)
+            self.assertTrue(a8.shape == r.shape)
+            r = _crank8.mean(image=a8,selem = elem,shift_x=+1,shift_y=+1)
+            self.assertTrue(a8.shape == r.shape)
+
+        for m,n in np.random.random_integers(1,100,size=(10,2)):
+            a16 = np.ones((m,n),dtype='uint16')
+            r = _crank16.mean(image=a16,selem = elem,shift_x=0,shift_y=0)
+            self.assertTrue(a16.shape == r.shape)
+            r = _crank16.mean(image=a16,selem = elem,shift_x=+1,shift_y=+1)
+            self.assertTrue(a16.shape == r.shape)
+
+        for m,n in np.random.random_integers(1,100,size=(10,2)):
+            a16 = np.ones((m,n),dtype='uint16')
+            r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9)
+            self.assertTrue(a16.shape == r.shape)
+            r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=+1,shift_y=+1,p0=.1,p1=.9)
+            self.assertTrue(a16.shape == r.shape)
+
+    def test_compare_with_cmorph_dilate(self):
+        #compare the result of maximum filter with dilate
+
+        a = (np.random.random((500,500))*256).astype('uint8')
+
+        for r in range(1,20,1):
+            elem = np.ones((r,r),dtype='uint8')
+            #        elem = (np.random.random((r,r))>.5).astype('uint8')
+            rc = _crank8.maximum(image=a,selem = elem)
+            cm = cmorph.dilate(image=a,selem = elem)
+            self.assertTrue((rc==cm).all())
+
+    def test_compare_with_cmorph_erode(self):
+        #compare the result of maximum filter with erode
+
+        a = (np.random.random((500,500))*256).astype('uint8')
+
+        for r in range(1,20,1):
+            elem = np.ones((r,r),dtype='uint8')
+            #        elem = (np.random.random((r,r))>.5).astype('uint8')
+            rc = _crank8.minimum(image=a,selem = elem)
+            cm = cmorph.erode(image=a,selem = elem)
+            self.assertTrue((rc==cm).all())
+
+    def test_bitdepth(self):
+        # test the different bit depth for rank16
+
+        elem = np.ones((3,3),dtype='uint8')
+        a16 = np.ones((100,100),dtype='uint16')*255
+        r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=8)
+        a16 = np.ones((100,100),dtype='uint16')*255*2
+        r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=9)
+        a16 = np.ones((100,100),dtype='uint16')*255*4
+        r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=10)
+        a16 = np.ones((100,100),dtype='uint16')*255*8
+        r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=11)
+        a16 = np.ones((100,100),dtype='uint16')*255*16
+        r = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=12)
+
+    def test_population(self):
+        # check the number of valid pixels in the neighborhood
+
+        a = np.zeros((5,5),dtype='uint8')
+        elem = np.ones((3,3),dtype='uint8')
+        p = _crank8.pop(image=a,selem = elem)
+        r = np.asarray([[4, 6, 6, 6, 4],
+            [6, 9, 9, 9, 6],
+            [6, 9, 9, 9, 6],
+            [6, 9, 9, 9, 6],
+            [4, 6, 6, 6, 4]])
+        np.testing.assert_array_equal(r,p)
+
+    def test_structuring_element(self):
+        # check the output for a custom structuring element
+
+        a = np.zeros((6,6),dtype='uint8')
+        a[2,2] = 255
+        elem = np.asarray([[1,1,0],[1,1,1],[0,0,1]],dtype='uint8')
+        f = _crank8.maximum(image=a,selem = elem,shift_x=1,shift_y=1)
+        r = np.asarray([[  0,   0,   0,   0,   0,   0],
+            [  0,   0,   0,   0,   0,   0],
+            [  0,   0, 255,   0,   0,   0],
+            [  0,   0, 255, 255, 255,   0],
+            [  0,   0,   0, 255, 255,   0],
+            [  0,   0,   0,   0,   0,   0]])
+        np.testing.assert_array_equal(r,f)
+
+
+    @unittest.expectedFailure
+    def test_fail_on_bitdepth(self):
+        # should fail because data bitdepth is too high for the function
+
+        a16 = np.ones((100,100),dtype='uint16')*255
+        elem = np.ones((3,3),dtype='uint8')
+        f = _crank16_percentiles.mean(image=a16,selem = elem,shift_x=0,shift_y=0,p0=.1,p1=.9,bitdepth=4)
+
+    def test_output(self):
+        #check rank function with external OUT output array
+
+        selem = disk(20)
+        a = (np.random.random((500,500))*256).astype('uint8')
+        out = np.zeros_like(a)
+        f1 = rank.mean(a,selem,out=out)
+        f2 = rank.mean(a,selem)
+        np.testing.assert_array_equal(f1,f2)
+        np.testing.assert_array_equal(out,f2)
+
+    @unittest.expectedFailure
+    def test_inplace_output(self):
+        #rank filters are not supposed to filter inplace
+
+        selem = disk(20)
+        a = (np.random.random((500,500))*256).astype('uint8')
+        out = a
+        f = rank.mean(a,selem,out=out)
+        np.testing.assert_array_equal(f,out)
+
+
+    def test_compare_autolevels(self):
+        # compare autolevel and percentile autolevel with p0=0.0 and p1=1.0
+        # should returns the same arrays
+
+        image = data.camera()
+
+        selem = disk(20)
+        loc_autolevel = rank.autolevel(image,selem=selem)
+        loc_perc_autolevel = rank.percentile_autolevel(image,selem=selem,p0=.0,p1=1.)
+
+        assert (loc_autolevel==loc_perc_autolevel).all()
+
+    def test_compare_autolevels_16bit(self):
+        # compare autolevel(16bit) and percentile autolevel(16bit) with p0=0.0 and p1=1.0
+        # should returns the same arrays
+
+        image = data.camera().astype(np.uint16)*4
+
+        selem = disk(20)
+        loc_autolevel = rank.autolevel(image,selem=selem)
+        loc_perc_autolevel = rank.percentile_autolevel(image,selem=selem,p0=.0,p1=1.)
+
+        assert (loc_autolevel==loc_perc_autolevel).all()
+
+    def test_compare_8bit_vs_16bit(self):
+        # filters applied on 8bit image ore 16bit image (having only real 8bit of dynamic)
+        # should be identical
+
+        i8 = data.camera()
+        i16 = i8.astype(np.uint16)
+        assert (i8==i16).all()
+
+        methods = ['autolevel','bottomhat','equalize','gradient','maximum','mean'
+            ,'meansubstraction','median','minimum','modal','morph_contr_enh','pop','threshold', 'tophat']
+
+        for method in methods:
+            func = eval('rank.%s'%method)
+            f8 = func(i8,disk(3))
+            f16 = func(i16,disk(3))
+            assert (f8==f16).all()
+
+
+if __name__ == '__main__':
+
+    suite = unittest.TestLoader().loadTestsFromTestCase(TestSequenceFunctions)
+    unittest.TextTestRunner(verbosity=2).run(suite)
diff --git a/skimage/filter/setup.py b/skimage/filter/setup.py
index b996055b..56c1e9e5 100644
--- a/skimage/filter/setup.py
+++ b/skimage/filter/setup.py
@@ -14,11 +14,47 @@ def configuration(parent_package='', top_path=None):
 
     cython(['_ctmf.pyx'], working_path=base_path)
     cython(['_denoise.pyx'], working_path=base_path)
+    cython(['rank/_core8.pyx'], working_path=base_path)
+    cython(['rank/_core16.pyx'], working_path=base_path)
+    cython(['rank/_crank8.pyx'], working_path=base_path)
+    cython(['rank/_crank8_percentiles.pyx'], working_path=base_path)
+    cython(['rank/_crank16.pyx'], working_path=base_path)
+    cython(['rank/_crank16_percentiles.pyx'], working_path=base_path)
+    cython(['rank/_crank16_bilateral.pyx'], working_path=base_path)
+    cython(['rank/rank.pyx'], working_path=base_path)
+    cython(['rank/percentile_rank.pyx'], working_path=base_path)
+    cython(['rank/bilateral_rank.pyx'], working_path=base_path)
 
     config.add_extension('_ctmf', sources=['_ctmf.c'],
-                         include_dirs=[get_numpy_include_dirs()])
+        include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_denoise', sources=['_denoise.c'],
-                         include_dirs=[get_numpy_include_dirs(), '../_shared'])
+        include_dirs=[get_numpy_include_dirs(), '../_shared'])
+    config.add_extension('rank/_core8', sources=['rank/_core8.c'],
+        include_dirs=[get_numpy_include_dirs()])
+    config.add_extension('rank/_core16', sources=['rank/_core16.c'],
+        include_dirs=[get_numpy_include_dirs()])
+    config.add_extension('rank/_crank8', sources=['rank/_crank8.c'],
+        include_dirs=[get_numpy_include_dirs()])
+    config.add_extension(
+        'rank/_crank8_percentiles', sources=['rank/_crank8_percentiles.c'],
+        include_dirs=[get_numpy_include_dirs()])
+    config.add_extension('rank/_crank16', sources=['rank/_crank16.c'],
+        include_dirs=[get_numpy_include_dirs()])
+    config.add_extension(
+        'rank/_crank16_percentiles', sources=['rank/_crank16_percentiles.c'],
+        include_dirs=[get_numpy_include_dirs()])
+    config.add_extension(
+        'rank/_crank16_bilateral', sources=['rank/_crank16_bilateral.c'],
+        include_dirs=[get_numpy_include_dirs()])
+    config.add_extension(
+        'rank/rank', sources=['rank/rank.c'],
+        include_dirs=[get_numpy_include_dirs()])
+    config.add_extension(
+        'rank/percentile_rank', sources=['rank/percentile_rank.c'],
+        include_dirs=[get_numpy_include_dirs()])
+    config.add_extension(
+        'rank/bilateral_rank', sources=['rank/bilateral_rank.c'],
+        include_dirs=[get_numpy_include_dirs()])
 
     return config