diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index 30f432f4..b0a82519 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -185,3 +185,6 @@
 
 - Adam Feuer
   PIL Image import and export improvements
+
+- Geoffrey French
+  skimage.filters.rank.windowed_histogram and plot_windowed_histogram example.
\ No newline at end of file
diff --git a/doc/examples/plot_windowed_histogram.py b/doc/examples/plot_windowed_histogram.py
new file mode 100644
index 00000000..1c57721b
--- /dev/null
+++ b/doc/examples/plot_windowed_histogram.py
@@ -0,0 +1,146 @@
+from __future__ import division
+"""
+========================
+Sliding window histogram
+========================
+
+Histogram matching can be used for object detection in images [1]_.
+This example extracts a single coin from the `skimage.data.coins` image
+and uses histogram matching to attempt to locate it within the original
+image.
+
+First, a box-shaped region of the image containing the target coin is
+extracted and a histogram of its greyscale values is computed.
+
+Next, for each pixel in the test image, a histogram of the greyscale values
+in a region of the image surrounding the pixel is computed.
+`skimage.filter.rank.windowed_histogram` is used for this task, as it
+employs an efficient sliding window based algorithm that is able to compute
+these histograms quickly [2]_.
+The local histogram for the region surrounding each pixel in the image is
+compared to that of the single coin, with a similarity measure being
+computed and displayed.
+
+The histogram of the single coin is computed using `numpy.histogram` on a
+box shaped region surrounding the coin, while the sliding window histograms
+are computed using a disc shaped structural element of a slightly different
+size. This is done in aid of demonstrating that the technique still finds
+similarity in spite of these differences.
+
+To demonstrate the rotational invariance of the technique, the same
+test is performed on a version of the coins image rotated by 45 degrees.
+
+References
+----------
+.. [1] Porikli, F. "Integral Histogram: A Fast Way to Extract Histograms
+       in Cartesian Spaces" CVPR, 2005. Vol. 1. IEEE, 2005
+.. [2] S.Perreault and P.Hebert. Median filtering in constant time.
+       Trans. Image Processing, 16(9):2389-2394, 2007.
+"""
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+
+from skimage import data
+from skimage.util import img_as_ubyte
+from skimage.morphology import disk
+from skimage.filter import rank
+from skimage import transform
+
+
+matplotlib.rcParams['font.size'] = 9
+
+
+def windowed_histogram_similarity(image, selem, reference_hist, n_bins):
+    # Compute normalized windowed histogram feature vector for each pixel
+    px_histograms = rank.windowed_histogram(image, selem, n_bins=n_bins)
+
+    # Reshape coin histogram to (1,1,N) for broadcast when we want to use it in
+    # arithmetic operations with the windowed histograms from the image
+    reference_hist = reference_hist.reshape((1, 1) + reference_hist.shape)
+
+    # Compute Chi squared distance metric: sum((X-Y)^2 / (X+Y));
+    # a measure of distance between histograms
+    X = px_histograms
+    Y = reference_hist
+    num = (X-Y)*(X-Y)
+    denom = X+Y
+    frac = num / denom
+    frac[denom == 0] = 0
+    chi_sqr = np.sum(frac, axis=2) * 0.5
+
+    # Generate a similarity measure. It needs to be low when distance is high
+    # and high when distance is low; taking the reciprocal will do this.
+    # Chi squared will always be >= 0, add small value to prevent divide by 0.
+    similarity = 1 / (chi_sqr + 1.0e-4)
+
+    return similarity
+
+
+# Load the `skimage.data.coins` image
+img = img_as_ubyte(data.coins())
+
+# Quantize to 16 levels of greyscale; this way the output image will have a
+# 16-dimensional feature vector per pixel
+quantized_img = img//16
+
+# Select the coin from the 4th column, second row.
+# Co-ordinate ordering: [x1,y1,x2,y2]
+coin_coords = [184, 100, 228, 148]   # 44 x 44 region
+coin = quantized_img[coin_coords[1]:coin_coords[3],
+                     coin_coords[0]:coin_coords[2]]
+
+# Compute coin histogram and normalize
+coin_hist, _ = np.histogram(coin.flatten(), bins=16, range=(0, 16))
+coin_hist = coin_hist.astype(float) / np.sum(coin_hist)
+
+
+# Compute a disk shaped mask that will define the shape of our sliding window
+# Example coin is ~44px across, so make a disk 61px wide (2*rad+1) to be big
+# enough for other coins too.
+selem = disk(30)
+
+
+# Compute the similarity across the complete image
+similarity = windowed_histogram_similarity(quantized_img, selem, coin_hist,
+                                           coin_hist.shape[0])
+
+# Now try a rotated image
+rotated_img = img_as_ubyte(transform.rotate(img, 45.0, resize=True))
+# Quantize to 16 levels as before
+quantized_rotated_image = rotated_img//16
+# Similarity on rotated image
+rotated_similarity = windowed_histogram_similarity(quantized_rotated_image,
+                                                   selem, coin_hist,
+                                                   coin_hist.shape[0])
+
+
+# Plot it all
+fig, axes = plt.subplots(nrows=5, figsize=(6, 18))
+ax0, ax1, ax2, ax3, ax4 = axes
+
+ax0.imshow(img, cmap='gray')
+ax0.set_title('Original image')
+ax0.axis('off')
+
+ax1.imshow(quantized_img, cmap='gray')
+ax1.set_title('Quantized image')
+ax1.axis('off')
+
+ax2.imshow(coin, cmap='gray')
+ax2.set_title('Coin from 2nd row, 4th column')
+ax2.axis('off')
+
+ax3.imshow(img, cmap='gray')
+# While jet is not a great colormap, it makes the high similarity areas
+# stand out
+ax3.imshow(similarity, cmap='jet', alpha=0.5)
+ax3.set_title('Original image with overlayed similarity')
+ax3.axis('off')
+
+ax4.imshow(rotated_img, cmap='gray')
+ax4.imshow(rotated_similarity, cmap='jet', alpha=0.5)
+ax4.set_title('Rotated image with overlayed similarity')
+ax4.axis('off')
+
+plt.show()
diff --git a/doc/release/release_dev.txt b/doc/release/release_dev.txt
index cbeded27..77c8d102 100644
--- a/doc/release/release_dev.txt
+++ b/doc/release/release_dev.txt
@@ -20,6 +20,7 @@ Region Adjacency Graphs
  - Similarity RAGs (#1080)
  - Normalized Cut on RAGs (#1080)
  - RAG Drawing (#1087)
+Sliding Windowed Histogram (#1127)
 
 Improvements
 ------------
diff --git a/skimage/filter/rank/__init__.py b/skimage/filter/rank/__init__.py
index 8641b984..22b0e881 100644
--- a/skimage/filter/rank/__init__.py
+++ b/skimage/filter/rank/__init__.py
@@ -1,6 +1,6 @@
 from .generic import (autolevel, bottomhat, equalize, gradient, maximum, mean,
                       subtract_mean, median, minimum, modal, enhance_contrast,
-                      pop, threshold, tophat, noise_filter, entropy, otsu, sum)
+                      pop, threshold, tophat, noise_filter, entropy, otsu, sum, windowed_histogram)
 from ._percentile import (autolevel_percentile, gradient_percentile,
                           mean_percentile, subtract_mean_percentile,
                           enhance_contrast_percentile, percentile,
@@ -37,4 +37,5 @@ __all__ = ['autolevel',
            'noise_filter',
            'entropy',
            'otsu',
-           'percentile']
+           'percentile',
+           'windowed_histogram']
diff --git a/skimage/filter/rank/_percentile.py b/skimage/filter/rank/_percentile.py
index 01dd0b49..9b93195b 100644
--- a/skimage/filter/rank/_percentile.py
+++ b/skimage/filter/rank/_percentile.py
@@ -43,7 +43,7 @@ def _apply(func, image, selem, out, mask, shift_x, shift_y, p0, p1,
     func(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask,
          out=out, max_bin=max_bin, p0=p0, p1=p1)
 
-    return out
+    return out.reshape(out.shape[:2])
 
 
 def autolevel_percentile(image, selem, out=None, mask=None, shift_x=False,
diff --git a/skimage/filter/rank/bilateral.py b/skimage/filter/rank/bilateral.py
index d01680db..aeb318d1 100644
--- a/skimage/filter/rank/bilateral.py
+++ b/skimage/filter/rank/bilateral.py
@@ -42,7 +42,7 @@ def _apply(func, image, selem, out, mask, shift_x, shift_y, s0, s1,
     func(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask,
          out=out, max_bin=max_bin, s0=s0, s1=s1)
 
-    return out
+    return out.reshape(out.shape[:2])
 
 
 def mean_bilateral(image, selem, out=None, mask=None, shift_x=False,
@@ -158,8 +158,9 @@ def pop_bilateral(image, selem, out=None, mask=None, shift_x=False,
     return _apply(bilateral_cy._pop, image, selem, out=out,
                   mask=mask, shift_x=shift_x, shift_y=shift_y, s0=s0, s1=s1)
 
+
 def sum_bilateral(image, selem, out=None, mask=None, shift_x=False,
-                   shift_y=False, s0=10, s1=10):
+                  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
diff --git a/skimage/filter/rank/bilateral_cy.pyx b/skimage/filter/rank/bilateral_cy.pyx
index 25a81766..b4e22d3f 100644
--- a/skimage/filter/rank/bilateral_cy.pyx
+++ b/skimage/filter/rank/bilateral_cy.pyx
@@ -9,10 +9,12 @@ from libc.math cimport log
 from .core_cy cimport dtype_t, dtype_t_out, _core
 
 
-cdef inline double _kernel_mean(Py_ssize_t* histo, double pop, dtype_t g,
-                                Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                double p0, double p1,
-                                Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth,
+                              Py_ssize_t* histo,
+                              double pop, dtype_t g,
+                              Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                              double p0, double p1,
+                              Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i
     cdef Py_ssize_t bilat_pop = 0
@@ -24,17 +26,19 @@ cdef inline double _kernel_mean(Py_ssize_t* histo, double pop, dtype_t g,
                 bilat_pop += histo[i]
                 mean += histo[i] * i
         if bilat_pop:
-            return mean / bilat_pop
+            out[0] = <dtype_t_out>mean / bilat_pop
         else:
-            return 0
+            out[0] = <dtype_t_out>0
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g,
-                               Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                               double p0, double p1,
-                               Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth,
+                             Py_ssize_t* histo,
+                             double pop, dtype_t g,
+                             Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                             double p0, double p1,
+                             Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i
     cdef Py_ssize_t bilat_pop = 0
@@ -43,14 +47,17 @@ cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g,
         for i in range(max_bin):
             if (g > (i - s0)) and (g < (i + s1)):
                 bilat_pop += histo[i]
-        return bilat_pop
+        out[0] = <dtype_t_out>bilat_pop
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
-cdef inline double _kernel_sum(Py_ssize_t* histo, double pop, dtype_t g,
-                                Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                double p0, double p1,
-                                Py_ssize_t s0, Py_ssize_t s1):
+
+cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth,
+                             Py_ssize_t* histo,
+                             double pop, dtype_t g,
+                             Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                             double p0, double p1,
+                             Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i
     cdef Py_ssize_t bilat_pop = 0
@@ -62,40 +69,41 @@ cdef inline double _kernel_sum(Py_ssize_t* histo, double pop, dtype_t g,
                 bilat_pop += histo[i]
                 sum += histo[i] * i
         if bilat_pop:
-            return sum
+            out[0] = <dtype_t_out>sum
         else:
-            return 0
+            out[0] = <dtype_t_out>0
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
 def _mean(dtype_t[:, ::1] image,
           char[:, ::1] selem,
           char[:, ::1] mask,
-          dtype_t_out[:, ::1] out,
+          dtype_t_out[:, :, ::1] out,
           char shift_x, char shift_y, Py_ssize_t s0, Py_ssize_t s1,
           Py_ssize_t max_bin):
 
-    _core(_kernel_mean[dtype_t], image, selem, mask, out,
+    _core(_kernel_mean[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, s0, s1, max_bin)
 
 
 def _pop(dtype_t[:, ::1] image,
          char[:, ::1] selem,
          char[:, ::1] mask,
-         dtype_t_out[:, ::1] out,
+         dtype_t_out[:, :, ::1] out,
          char shift_x, char shift_y, Py_ssize_t s0, Py_ssize_t s1,
          Py_ssize_t max_bin):
 
-    _core(_kernel_pop[dtype_t], image, selem, mask, out,
+    _core(_kernel_pop[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, s0, s1, max_bin)
 
+
 def _sum(dtype_t[:, ::1] image,
          char[:, ::1] selem,
          char[:, ::1] mask,
-         dtype_t_out[:, ::1] out,
+         dtype_t_out[:, :, ::1] out,
          char shift_x, char shift_y, Py_ssize_t s0, Py_ssize_t s1,
          Py_ssize_t max_bin):
 
-    _core(_kernel_sum[dtype_t], image, selem, mask, out,
+    _core(_kernel_sum[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, s0, s1, max_bin)
diff --git a/skimage/filter/rank/core_cy.pxd b/skimage/filter/rank/core_cy.pxd
index 2e97e50a..9627c0d9 100644
--- a/skimage/filter/rank/core_cy.pxd
+++ b/skimage/filter/rank/core_cy.pxd
@@ -15,13 +15,13 @@ cdef dtype_t _max(dtype_t a, dtype_t b)
 cdef dtype_t _min(dtype_t a, dtype_t b)
 
 
-cdef void _core(double kernel(Py_ssize_t*, double, dtype_t,
-                              Py_ssize_t, Py_ssize_t, double,
-                              double, Py_ssize_t, Py_ssize_t),
+cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double,
+                            dtype_t, Py_ssize_t, Py_ssize_t, double,
+                            double, Py_ssize_t, Py_ssize_t),
                 dtype_t[:, ::1] image,
                 char[:, ::1] selem,
                 char[:, ::1] mask,
-                dtype_t_out[:, ::1] out,
+                dtype_t_out[:, :, ::1] out,
                 char shift_x, char shift_y,
                 double p0, double p1,
                 Py_ssize_t s0, Py_ssize_t s1,
diff --git a/skimage/filter/rank/core_cy.pyx b/skimage/filter/rank/core_cy.pyx
index 02c2c8d0..4bd21df0 100644
--- a/skimage/filter/rank/core_cy.pyx
+++ b/skimage/filter/rank/core_cy.pyx
@@ -42,13 +42,13 @@ cdef inline char is_in_mask(Py_ssize_t rows, Py_ssize_t cols,
             return 0
 
 
-cdef void _core(double kernel(Py_ssize_t*, double, dtype_t,
-                              Py_ssize_t, Py_ssize_t, double,
-                              double, Py_ssize_t, Py_ssize_t),
+cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double,
+                            dtype_t, Py_ssize_t, Py_ssize_t, double,
+                            double, Py_ssize_t, Py_ssize_t),
                 dtype_t[:, ::1] image,
                 char[:, ::1] selem,
                 char[:, ::1] mask,
-                dtype_t_out[:, ::1] out,
+                dtype_t_out[:, :, ::1] out,
                 char shift_x, char shift_y,
                 double p0, double p1,
                 Py_ssize_t s0, Py_ssize_t s1,
@@ -61,6 +61,7 @@ cdef void _core(double kernel(Py_ssize_t*, double, dtype_t,
     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 odepth = out.shape[2]
 
     cdef Py_ssize_t centre_r = <Py_ssize_t>(selem.shape[0] / 2) + shift_y
     cdef Py_ssize_t centre_c = <Py_ssize_t>(selem.shape[1] / 2) + shift_x
@@ -151,8 +152,8 @@ cdef void _core(double kernel(Py_ssize_t*, double, dtype_t,
 
     r = 0
     c = 0
-    out[r, c] = <dtype_t_out>kernel(histo, pop, image[r, c], max_bin, mid_bin,
-                                    p0, p1, s0, s1)
+    kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin,
+           p0, p1, s0, s1)
 
     # main loop
     r = 0
@@ -172,8 +173,8 @@ cdef void _core(double kernel(Py_ssize_t*, double, dtype_t,
                 if is_in_mask(rows, cols, rr, cc, mask_data):
                     histogram_decrement(histo, &pop, image[rr, cc])
 
-            out[r, c] = <dtype_t_out>kernel(histo, pop, image[r, c],
-                                            max_bin, mid_bin, p0, p1, s0, s1)
+            kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin,
+                   mid_bin, p0, p1, s0, s1)
 
         r += 1  # pass to the next row
         if r >= rows:
@@ -192,8 +193,8 @@ cdef void _core(double kernel(Py_ssize_t*, double, dtype_t,
             if is_in_mask(rows, cols, rr, cc, mask_data):
                 histogram_decrement(histo, &pop, image[rr, cc])
 
-        out[r, c] = <dtype_t_out>kernel(histo, pop, image[r, c],
-                                        max_bin, mid_bin, p0, p1, s0, s1)
+        kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin,
+               mid_bin, p0, p1, s0, s1)
 
         # ---> east to west
         for c in range(cols - 2, -1, -1):
@@ -209,8 +210,8 @@ cdef void _core(double kernel(Py_ssize_t*, double, dtype_t,
                 if is_in_mask(rows, cols, rr, cc, mask_data):
                     histogram_decrement(histo, &pop, image[rr, cc])
 
-            out[r, c] = <dtype_t_out>kernel(histo, pop, image[r, c],
-                                            max_bin, mid_bin, p0, p1, s0, s1)
+            kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin,
+                   mid_bin, p0, p1, s0, s1)
 
         r += 1  # pass to the next row
         if r >= rows:
@@ -229,8 +230,8 @@ cdef void _core(double kernel(Py_ssize_t*, double, dtype_t,
             if is_in_mask(rows, cols, rr, cc, mask_data):
                 histogram_decrement(histo, &pop, image[rr, cc])
 
-        out[r, c] = <dtype_t_out>kernel(histo, pop, image[r, c],
-                                        max_bin, mid_bin, p0, p1, s0, s1)
+        kernel(&out[r, c, 0], odepth, histo, pop, image[r, c],
+               max_bin, mid_bin, p0, p1, s0, s1)
 
     # release memory allocated by malloc
     free(se_e_r)
diff --git a/skimage/filter/rank/generic.py b/skimage/filter/rank/generic.py
index ccc1166e..a474d050 100644
--- a/skimage/filter/rank/generic.py
+++ b/skimage/filter/rank/generic.py
@@ -28,7 +28,7 @@ __all__ = ['autolevel', 'bottomhat', 'equalize', 'gradient', 'maximum', 'mean',
            'pop', 'threshold', 'tophat', 'noise_filter', 'entropy', 'otsu']
 
 
-def _handle_input(image, selem, out, mask, out_dtype=None):
+def _handle_input(image, selem, out, mask, out_dtype=None, pixel_size=1):
 
     if image.dtype not in (np.uint8, np.uint16):
         image = img_as_ubyte(image)
@@ -42,13 +42,16 @@ def _handle_input(image, selem, out, mask, out_dtype=None):
         mask = img_as_ubyte(mask)
         mask = np.ascontiguousarray(mask)
 
+    if image is out:
+        raise NotImplementedError("Cannot perform rank operation in place.")
+
     if out is None:
         if out_dtype is None:
             out_dtype = image.dtype
-        out = np.empty_like(image, dtype=out_dtype)
-
-    if image is out:
-        raise NotImplementedError("Cannot perform rank operation in place.")
+        out = np.empty(image.shape+(pixel_size,), dtype=out_dtype)
+    else:
+        if len(out.shape) == 2:
+            out = out.reshape(out.shape+(pixel_size,))
 
     is_8bit = image.dtype in (np.uint8, np.int8)
 
@@ -65,7 +68,8 @@ def _handle_input(image, selem, out, mask, out_dtype=None):
     return image, selem, out, mask, max_bin
 
 
-def _apply(func, image, selem, out, mask, shift_x, shift_y, out_dtype=None):
+def _apply_scalar_per_pixel(func, image, selem, out, mask, shift_x, shift_y,
+                            out_dtype=None):
 
     image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask,
                                                      out_dtype)
@@ -73,6 +77,19 @@ def _apply(func, image, selem, out, mask, shift_x, shift_y, out_dtype=None):
     func(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask,
          out=out, max_bin=max_bin)
 
+    return out.reshape(out.shape[:2])
+
+
+def _apply_vector_per_pixel(func, image, selem, out, mask, shift_x, shift_y,
+                            out_dtype=None, pixel_size=1):
+
+    image, selem, out, mask, max_bin = _handle_input(image, selem, out, mask,
+                                                     out_dtype,
+                                                     pixel_size=pixel_size)
+
+    func(image, selem, shift_x=shift_x, shift_y=shift_y, mask=mask,
+         out=out, max_bin=max_bin)
+
     return out
 
 
@@ -113,8 +130,9 @@ def autolevel(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._autolevel, image, selem,
-                  out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._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):
@@ -154,8 +172,9 @@ def bottomhat(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._bottomhat, image, selem,
-                  out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._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):
@@ -192,8 +211,9 @@ def equalize(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._equalize, image, selem,
-                  out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._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):
@@ -230,8 +250,9 @@ def gradient(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._gradient, image, selem,
-                  out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._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):
@@ -277,8 +298,9 @@ def maximum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._maximum, image, selem,
-                  out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._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):
@@ -315,8 +337,8 @@ def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._mean, image, selem, out=out,
-                  mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._mean, image, selem, out=out,
+                                   mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
 def subtract_mean(image, selem, out=None, mask=None, shift_x=False,
@@ -354,8 +376,9 @@ def subtract_mean(image, selem, out=None, mask=None, shift_x=False,
 
     """
 
-    return _apply(generic_cy._subtract_mean, image, selem,
-                  out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._subtract_mean, 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):
@@ -392,8 +415,9 @@ def median(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._median, image, selem,
-                  out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._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):
@@ -439,8 +463,9 @@ def minimum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._minimum, image, selem,
-                  out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._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):
@@ -479,8 +504,9 @@ def modal(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._modal, image, selem,
-                  out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._modal, image, selem,
+                                   out=out, mask=mask,
+                                   shift_x=shift_x, shift_y=shift_y)
 
 
 def enhance_contrast(image, selem, out=None, mask=None, shift_x=False,
@@ -522,8 +548,9 @@ def enhance_contrast(image, selem, out=None, mask=None, shift_x=False,
 
     """
 
-    return _apply(generic_cy._enhance_contrast, image, selem,
-                  out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._enhance_contrast, 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):
@@ -571,8 +598,9 @@ def pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._pop, image, selem, out=out,
-                  mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._pop, image, selem, out=out,
+                                   mask=mask, shift_x=shift_x,
+                                   shift_y=shift_y)
 
 
 def sum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
@@ -620,8 +648,9 @@ def sum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._sum, image, selem, out=out,
-                  mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._sum, 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):
@@ -669,8 +698,9 @@ def threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._threshold, image, selem,
-                  out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._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):
@@ -710,8 +740,9 @@ def tophat(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._tophat, image, selem,
-                  out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._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,
@@ -761,8 +792,9 @@ def noise_filter(image, selem, out=None, mask=None, shift_x=False,
     selem_cpy = selem.copy()
     selem_cpy[centre_r, centre_c] = 0
 
-    return _apply(generic_cy._noise_filter, image, selem_cpy, out=out,
-                  mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._noise_filter, 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):
@@ -806,9 +838,10 @@ def entropy(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._entropy, image, selem,
-                  out=out, mask=mask, shift_x=shift_x, shift_y=shift_y,
-                  out_dtype=np.double)
+    return _apply_scalar_per_pixel(generic_cy._entropy, image, selem,
+                                   out=out, mask=mask,
+                                   shift_x=shift_x, shift_y=shift_y,
+                                   out_dtype=np.double)
 
 
 def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
@@ -850,5 +883,59 @@ def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._otsu, image, selem, out=out,
-                  mask=mask, shift_x=shift_x, shift_y=shift_y)
+    return _apply_scalar_per_pixel(generic_cy._otsu, image, selem, out=out,
+                                   mask=mask, shift_x=shift_x,
+                                   shift_y=shift_y)
+
+
+def windowed_histogram(image, selem, out=None, mask=None,
+                       shift_x=False, shift_y=False, n_bins=None):
+    """Normalized sliding window histogram
+
+    Parameters
+    ----------
+    image : ndarray
+        Image array (uint8 array).
+    selem : 2-D array
+        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
+        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).
+    n_bins : int or None
+        The number of histogram bins. Will default to ``image.max() + 1``
+        if None is passed.
+
+    Returns
+    -------
+    out : 3-D array with float dtype of dimensions (H,W,N), where (H,W) are
+        the dimensions of the input image and N is n_bins or
+        ``image.max() + 1`` if no value is provided as a parameter.
+        Effectively, each pixel is a N-D feature vector that is the histogram.
+        The sum of the elements in the feature vector will be 1, unless no
+        pixels in the window were covered by both selem and mask, in which
+        case all elements will be 0.
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> from skimage.filter.rank import windowed_histogram
+    >>> from skimage.morphology import disk
+    >>> img = data.camera()
+    >>> hist_img = windowed_histogram(img, disk(5))
+
+    """
+
+    if n_bins is None:
+        n_bins = image.max() + 1
+
+    return _apply_vector_per_pixel(generic_cy._windowed_hist, image, selem,
+                                   out=out, mask=mask,
+                                   shift_x=shift_x, shift_y=shift_y,
+                                   out_dtype=np.double,
+                                   pixel_size=n_bins)
diff --git a/skimage/filter/rank/generic_cy.pyx b/skimage/filter/rank/generic_cy.pyx
index 1c26cf53..98c8cd8a 100644
--- a/skimage/filter/rank/generic_cy.pyx
+++ b/skimage/filter/rank/generic_cy.pyx
@@ -9,10 +9,12 @@ from libc.math cimport log
 from .core_cy cimport dtype_t, dtype_t_out, _core
 
 
-cdef inline double _kernel_autolevel(Py_ssize_t* histo, double pop, dtype_t g,
-                                     Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                     double p0, double p1,
-                                     Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_autolevel(dtype_t_out* out, Py_ssize_t odepth,
+                                   Py_ssize_t* histo,
+                                   double pop, dtype_t g,
+                                   Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                   double p0, double p1,
+                                   Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i, imin, imax, delta
 
@@ -27,17 +29,19 @@ cdef inline double _kernel_autolevel(Py_ssize_t* histo, double pop, dtype_t g,
                 break
         delta = imax - imin
         if delta > 0:
-            return <double>(max_bin - 1) * (g - imin) / delta
+            out[0] = <dtype_t_out>(max_bin - 1) * (g - imin) / delta
         else:
-            return 0
+            out[0] = <dtype_t_out>0
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_bottomhat(Py_ssize_t* histo, double pop, dtype_t g,
-                                     Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                     double p0, double p1,
-                                     Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_bottomhat(dtype_t_out* out, Py_ssize_t odepth,
+                                   Py_ssize_t* histo,
+                                   double pop, dtype_t g,
+                                   Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                   double p0, double p1,
+                                   Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i
 
@@ -45,15 +49,17 @@ cdef inline double _kernel_bottomhat(Py_ssize_t* histo, double pop, dtype_t g,
         for i in range(max_bin):
             if histo[i]:
                 break
-        return g - i
+        out[0] = <dtype_t_out>g - i
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_equalize(Py_ssize_t* histo, double pop, dtype_t g,
-                                    Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                    double p0, double p1,
-                                    Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_equalize(dtype_t_out* out, Py_ssize_t odepth,
+                                  Py_ssize_t* histo,
+                                  double pop, dtype_t g,
+                                  Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                  double p0, double p1,
+                                  Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i
     cdef Py_ssize_t sum = 0
@@ -63,15 +69,17 @@ cdef inline double _kernel_equalize(Py_ssize_t* histo, double pop, dtype_t g,
             sum += histo[i]
             if i >= g:
                 break
-        return ((max_bin - 1) * sum) / pop
+        out[0] = <dtype_t_out>(((max_bin - 1) * sum) / pop)
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_gradient(Py_ssize_t* histo, double pop, dtype_t g,
-                                    Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                    double p0, double p1,
-                                    Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_gradient(dtype_t_out* out, Py_ssize_t odepth,
+                                  Py_ssize_t* histo,
+                                  double pop, dtype_t g,
+                                  Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                  double p0, double p1,
+                                  Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i, imin, imax
 
@@ -84,65 +92,72 @@ cdef inline double _kernel_gradient(Py_ssize_t* histo, double pop, dtype_t g,
             if histo[i]:
                 imin = i
                 break
-        return imax - imin
+        out[0] = <dtype_t_out>(imax - imin)
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_maximum(Py_ssize_t* histo, double pop, dtype_t g,
-                                   Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                   double p0, double p1,
-                                   Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_maximum(dtype_t_out* out, Py_ssize_t odepth,
+                                 Py_ssize_t* histo,
+                                 double pop, dtype_t g,
+                                 Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                 double p0, double p1,
+                                 Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i
 
     if pop:
         for i in range(max_bin - 1, -1, -1):
             if histo[i]:
-                return i
+                out[0] = <dtype_t_out>i
+                return
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_mean(Py_ssize_t* histo, double pop,dtype_t g,
+cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth,
+                              Py_ssize_t* histo,
+                              double pop, dtype_t g,
+                              Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                              double p0, double p1,
+                              Py_ssize_t s0, Py_ssize_t s1):
+
+    cdef Py_ssize_t i
+    cdef Py_ssize_t mean = 0
+
+    if pop:
+        for i in range(max_bin):
+            mean += histo[i] * i
+        out[0] = <dtype_t_out>(mean / pop)
+    else:
+        out[0] = <dtype_t_out>0
+
+
+cdef inline void _kernel_subtract_mean(dtype_t_out* out, Py_ssize_t odepth,
+                                       Py_ssize_t* histo,
+                                       double pop, dtype_t g,
+                                       Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                       double p0, double p1,
+                                       Py_ssize_t s0, Py_ssize_t s1):
+
+    cdef Py_ssize_t i
+    cdef Py_ssize_t mean = 0
+
+    if pop:
+        for i in range(max_bin):
+            mean += histo[i] * i
+        out[0] = <dtype_t_out>((g - mean / pop) / 2. + 127)
+    else:
+        out[0] = <dtype_t_out>0
+
+
+cdef inline void _kernel_median(dtype_t_out* out, Py_ssize_t odepth,
+                                Py_ssize_t* histo,
+                                double pop, dtype_t g,
                                 Py_ssize_t max_bin, Py_ssize_t mid_bin,
                                 double p0, double p1,
                                 Py_ssize_t s0, Py_ssize_t s1):
 
-    cdef Py_ssize_t i
-    cdef Py_ssize_t mean = 0
-
-    if pop:
-        for i in range(max_bin):
-            mean += histo[i] * i
-        return mean / pop
-    else:
-        return 0
-
-
-cdef inline double _kernel_subtract_mean(Py_ssize_t* histo, double pop,
-                                         dtype_t g,
-                                         Py_ssize_t max_bin,
-                                         Py_ssize_t mid_bin, double p0,
-                                         double p1, Py_ssize_t s0,
-                                         Py_ssize_t s1):
-
-    cdef Py_ssize_t i
-    cdef Py_ssize_t mean = 0
-
-    if pop:
-        for i in range(max_bin):
-            mean += histo[i] * i
-        return (g - mean / pop) / 2. + 127
-    else:
-        return 0
-
-
-cdef inline double _kernel_median(Py_ssize_t* histo, double pop, dtype_t g,
-                                  Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                  double p0, double p1,
-                                  Py_ssize_t s0, Py_ssize_t s1):
-
     cdef Py_ssize_t i
     cdef double sum = pop / 2.0
 
@@ -151,30 +166,36 @@ cdef inline double _kernel_median(Py_ssize_t* histo, double pop, dtype_t g,
             if histo[i]:
                 sum -= histo[i]
                 if sum < 0:
-                    return i
+                    out[0] = <dtype_t_out>i
+                    return
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_minimum(Py_ssize_t* histo, double pop, dtype_t g,
-                                   Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                   double p0, double p1,
-                                   Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_minimum(dtype_t_out* out, Py_ssize_t odepth,
+                                 Py_ssize_t* histo,
+                                 double pop, dtype_t g,
+                                 Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                 double p0, double p1,
+                                 Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i
 
     if pop:
         for i in range(max_bin):
             if histo[i]:
-                return i
+                out[0] = <dtype_t_out>i
+                return
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_modal(Py_ssize_t* histo, double pop, dtype_t g,
-                                 Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                 double p0, double p1,
-                                 Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_modal(dtype_t_out* out, Py_ssize_t odepth,
+                               Py_ssize_t* histo,
+                               double pop, dtype_t g,
+                               Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                               double p0, double p1,
+                               Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t hmax = 0, imax = 0
 
@@ -183,17 +204,20 @@ cdef inline double _kernel_modal(Py_ssize_t* histo, double pop, dtype_t g,
             if histo[i] > hmax:
                 hmax = histo[i]
                 imax = i
-        return imax
+        out[0] = <dtype_t_out>imax
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_enhance_contrast(Py_ssize_t* histo, double pop,
-                                            dtype_t g,
-                                            Py_ssize_t max_bin,
-                                            Py_ssize_t mid_bin, double p0,
-                                            double p1, Py_ssize_t s0,
-                                            Py_ssize_t s1):
+cdef inline void _kernel_enhance_contrast(dtype_t_out* out,
+                                          Py_ssize_t odepth,
+                                          Py_ssize_t* histo,
+                                          double pop,
+                                          dtype_t g,
+                                          Py_ssize_t max_bin,
+                                          Py_ssize_t mid_bin, double p0,
+                                          double p1, Py_ssize_t s0,
+                                          Py_ssize_t s1):
 
     cdef Py_ssize_t i, imin, imax
 
@@ -207,25 +231,29 @@ cdef inline double _kernel_enhance_contrast(Py_ssize_t* histo, double pop,
                 imin = i
                 break
         if imax - g < g - imin:
-            return imax
+            out[0] = <dtype_t_out>imax
         else:
-            return imin
+            out[0] = <dtype_t_out>imin
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g,
-                               Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                               double p0, double p1,
-                               Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth,
+                             Py_ssize_t* histo,
+                             double pop, dtype_t g,
+                             Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                             double p0, double p1,
+                             Py_ssize_t s0, Py_ssize_t s1):
 
-    return pop
+    out[0] = <dtype_t_out>pop
 
 
-cdef inline double _kernel_sum(Py_ssize_t* histo, double pop,dtype_t g,
-                               Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                               double p0, double p1,
-                               Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth,
+                             Py_ssize_t* histo,
+                             double pop, dtype_t g,
+                             Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                             double p0, double p1,
+                             Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i
     cdef Py_ssize_t sum = 0
@@ -233,15 +261,17 @@ cdef inline double _kernel_sum(Py_ssize_t* histo, double pop,dtype_t g,
     if pop:
         for i in range(max_bin):
             sum += histo[i] * i
-        return sum
+        out[0] = <dtype_t_out>sum
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_threshold(Py_ssize_t* histo, double pop, dtype_t g,
-                                     Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                     double p0, double p1,
-                                     Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_threshold(dtype_t_out* out, Py_ssize_t odepth,
+                                   Py_ssize_t* histo,
+                                   double pop, dtype_t g,
+                                   Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                   double p0, double p1,
+                                   Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i
     cdef Py_ssize_t mean = 0
@@ -249,15 +279,17 @@ cdef inline double _kernel_threshold(Py_ssize_t* histo, double pop, dtype_t g,
     if pop:
         for i in range(max_bin):
             mean += histo[i] * i
-        return g > (mean / pop)
+        out[0] = <dtype_t_out>(g > (mean / pop))
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_tophat(Py_ssize_t* histo, double pop, dtype_t g,
-                                  Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                  double p0, double p1,
-                                  Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_tophat(dtype_t_out* out, Py_ssize_t odepth,
+                                Py_ssize_t* histo,
+                                double pop, dtype_t g,
+                                Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                double p0, double p1,
+                                Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i
 
@@ -265,23 +297,24 @@ cdef inline double _kernel_tophat(Py_ssize_t* histo, double pop, dtype_t g,
         for i in range(max_bin - 1, -1, -1):
             if histo[i]:
                 break
-        return i - g
+        out[0] = <dtype_t_out>(i - g)
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_noise_filter(Py_ssize_t* histo, double pop,
-                                        dtype_t g, Py_ssize_t max_bin,
-                                        Py_ssize_t mid_bin, double p0,
-                                        double p1, Py_ssize_t s0,
-                                        Py_ssize_t s1):
+cdef inline void _kernel_noise_filter(dtype_t_out* out, Py_ssize_t odepth,
+                                      Py_ssize_t* histo,
+                                      double pop, dtype_t g,
+                                      Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                      double p0, double p1,
+                                      Py_ssize_t s0, Py_ssize_t s1):
 
     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 0
+        out[0] = <dtype_t_out>0
 
     for i in range(g, -1, -1):
         if histo[i]:
@@ -291,15 +324,17 @@ cdef inline double _kernel_noise_filter(Py_ssize_t* histo, double pop,
         if histo[i]:
             break
     if i - g < min_i:
-        return i - g
+        out[0] = <dtype_t_out>(i - g)
     else:
-        return min_i
+        out[0] = <dtype_t_out>min_i
 
 
-cdef inline double _kernel_entropy(Py_ssize_t* histo, double pop, dtype_t g,
-                                   Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                   double p0, double p1,
-                                   Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_entropy(dtype_t_out* out, Py_ssize_t odepth,
+                                 Py_ssize_t* histo,
+                                 double pop, dtype_t g,
+                                 Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                 double p0, double p1,
+                                 Py_ssize_t s0, Py_ssize_t s1):
     cdef Py_ssize_t i
     cdef double e, p
 
@@ -309,15 +344,17 @@ cdef inline double _kernel_entropy(Py_ssize_t* histo, double pop, dtype_t g,
             p = histo[i] / pop
             if p > 0:
                 e -= p * log(p) / 0.6931471805599453
-        return e
+        out[0] = <dtype_t_out>e
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_otsu(Py_ssize_t* histo, double pop, dtype_t g,
-                                Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                double p0, double p1,
-                                Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_otsu(dtype_t_out* out, Py_ssize_t odepth,
+                              Py_ssize_t* histo,
+                              double pop, dtype_t g,
+                              Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                              double p0, double p1,
+                              Py_ssize_t s0, Py_ssize_t s1):
     cdef Py_ssize_t i
     cdef Py_ssize_t max_i
     cdef double P, mu1, mu2, q1, new_q1, sigma_b, max_sigma_b
@@ -329,7 +366,7 @@ cdef inline double _kernel_otsu(Py_ssize_t* histo, double pop, dtype_t g,
             mu += histo[i] * i
         mu = mu / pop
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
     # maximizing the between class variance
     max_i = 0
@@ -349,183 +386,212 @@ cdef inline double _kernel_otsu(Py_ssize_t* histo, double pop, dtype_t g,
                 max_i = i
             q1 = new_q1
 
-    return max_i
+    out[0] = <dtype_t_out>max_i
+
+
+cdef inline void _kernel_win_hist(dtype_t_out* out, Py_ssize_t odepth,
+                                  Py_ssize_t* histo,
+                                  double pop, dtype_t g,
+                                  Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                  double p0, double p1,
+                                  Py_ssize_t s0, Py_ssize_t s1):
+    cdef Py_ssize_t i
+    cdef Py_ssize_t max_i
+    cdef double scale
+    if pop:
+        scale = 1.0 / pop
+        for i in xrange(odepth):
+            out[i] = <dtype_t_out>(histo[i] * scale)
+    else:
+        for i in xrange(odepth):
+            out[i] = <dtype_t_out>0
 
 
 def _autolevel(dtype_t[:, ::1] image,
                char[:, ::1] selem,
                char[:, ::1] mask,
-               dtype_t_out[:, ::1] out,
+               dtype_t_out[:, :, ::1] out,
                char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_autolevel[dtype_t], image, selem, mask, out,
+    _core(_kernel_autolevel[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _bottomhat(dtype_t[:, ::1] image,
                char[:, ::1] selem,
                char[:, ::1] mask,
-               dtype_t_out[:, ::1] out,
+               dtype_t_out[:, :, ::1] out,
                char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_bottomhat[dtype_t], image, selem, mask, out,
+    _core(_kernel_bottomhat[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _equalize(dtype_t[:, ::1] image,
               char[:, ::1] selem,
               char[:, ::1] mask,
-              dtype_t_out[:, ::1] out,
+              dtype_t_out[:, :, ::1] out,
               char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_equalize[dtype_t], image, selem, mask, out,
+    _core(_kernel_equalize[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _gradient(dtype_t[:, ::1] image,
               char[:, ::1] selem,
               char[:, ::1] mask,
-              dtype_t_out[:, ::1] out,
+              dtype_t_out[:, :, ::1] out,
               char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_gradient[dtype_t], image, selem, mask, out,
+    _core(_kernel_gradient[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _maximum(dtype_t[:, ::1] image,
              char[:, ::1] selem,
              char[:, ::1] mask,
-             dtype_t_out[:, ::1] out,
+             dtype_t_out[:, :, ::1] out,
              char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_maximum[dtype_t], image, selem, mask, out,
+    _core(_kernel_maximum[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _mean(dtype_t[:, ::1] image,
           char[:, ::1] selem,
           char[:, ::1] mask,
-          dtype_t_out[:, ::1] out,
+          dtype_t_out[:, :, ::1] out,
           char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_mean[dtype_t], image, selem, mask, out,
+    _core(_kernel_mean[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _subtract_mean(dtype_t[:, ::1] image,
                    char[:, ::1] selem,
                    char[:, ::1] mask,
-                   dtype_t_out[:, ::1] out,
+                   dtype_t_out[:, :, ::1] out,
                    char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_subtract_mean[dtype_t], image, selem, mask,
+    _core(_kernel_subtract_mean[dtype_t_out, dtype_t], image, selem, mask,
           out, shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _median(dtype_t[:, ::1] image,
             char[:, ::1] selem,
             char[:, ::1] mask,
-            dtype_t_out[:, ::1] out,
+            dtype_t_out[:, :, ::1] out,
             char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_median[dtype_t], image, selem, mask, out,
+    _core(_kernel_median[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _minimum(dtype_t[:, ::1] image,
              char[:, ::1] selem,
              char[:, ::1] mask,
-             dtype_t_out[:, ::1] out,
+             dtype_t_out[:, :, ::1] out,
              char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_minimum[dtype_t], image, selem, mask, out,
+    _core(_kernel_minimum[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _enhance_contrast(dtype_t[:, ::1] image,
                       char[:, ::1] selem,
                       char[:, ::1] mask,
-                      dtype_t_out[:, ::1] out,
+                      dtype_t_out[:, :, ::1] out,
                       char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_enhance_contrast[dtype_t], image, selem, mask,
+    _core(_kernel_enhance_contrast[dtype_t_out, dtype_t], image, selem, mask,
           out, shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _modal(dtype_t[:, ::1] image,
            char[:, ::1] selem,
            char[:, ::1] mask,
-           dtype_t_out[:, ::1] out,
+           dtype_t_out[:, :, ::1] out,
            char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_modal[dtype_t], image, selem, mask, out,
+    _core(_kernel_modal[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _pop(dtype_t[:, ::1] image,
          char[:, ::1] selem,
          char[:, ::1] mask,
-         dtype_t_out[:, ::1] out,
+         dtype_t_out[:, :, ::1] out,
          char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_pop[dtype_t], image, selem, mask, out,
+    _core(_kernel_pop[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
+
 def _sum(dtype_t[:, ::1] image,
          char[:, ::1] selem,
          char[:, ::1] mask,
-         dtype_t_out[:, ::1] out,
+         dtype_t_out[:, :, ::1] out,
          char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_sum[dtype_t], image, selem, mask,
+    _core(_kernel_sum[dtype_t_out, dtype_t], image, selem, mask,
           out, shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _threshold(dtype_t[:, ::1] image,
                char[:, ::1] selem,
                char[:, ::1] mask,
-               dtype_t_out[:, ::1] out,
+               dtype_t_out[:, :, ::1] out,
                char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_threshold[dtype_t], image, selem, mask, out,
+    _core(_kernel_threshold[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _tophat(dtype_t[:, ::1] image,
             char[:, ::1] selem,
             char[:, ::1] mask,
-            dtype_t_out[:, ::1] out,
+            dtype_t_out[:, :, ::1] out,
             char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_tophat[dtype_t], image, selem, mask, out,
+    _core(_kernel_tophat[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _noise_filter(dtype_t[:, ::1] image,
                   char[:, ::1] selem,
                   char[:, ::1] mask,
-                  dtype_t_out[:, ::1] out,
+                  dtype_t_out[:, :, ::1] out,
                   char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_noise_filter[dtype_t], image, selem, mask, out,
+    _core(_kernel_noise_filter[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _entropy(dtype_t[:, ::1] image,
              char[:, ::1] selem,
              char[:, ::1] mask,
-             dtype_t_out[:, ::1] out,
+             dtype_t_out[:, :, ::1] out,
              char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_entropy[dtype_t], image, selem, mask, out,
+    _core(_kernel_entropy[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
 
 
 def _otsu(dtype_t[:, ::1] image,
           char[:, ::1] selem,
           char[:, ::1] mask,
-          dtype_t_out[:, ::1] out,
+          dtype_t_out[:, :, ::1] out,
           char shift_x, char shift_y, Py_ssize_t max_bin):
 
-    _core(_kernel_otsu[dtype_t], image, selem, mask, out,
+    _core(_kernel_otsu[dtype_t_out, dtype_t], image, selem, mask, out,
+          shift_x, shift_y, 0, 0, 0, 0, max_bin)
+
+
+def _windowed_hist(dtype_t[:, ::1] image,
+                   char[:, ::1] selem,
+                   char[:, ::1] mask,
+                   dtype_t_out[:, :, ::1] out,
+                   char shift_x, char shift_y, Py_ssize_t max_bin):
+
+    _core(_kernel_win_hist[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, 0, 0, 0, 0, max_bin)
diff --git a/skimage/filter/rank/percentile_cy.pyx b/skimage/filter/rank/percentile_cy.pyx
index 38d04b33..5ff584d3 100644
--- a/skimage/filter/rank/percentile_cy.pyx
+++ b/skimage/filter/rank/percentile_cy.pyx
@@ -7,10 +7,12 @@ cimport numpy as cnp
 from .core_cy cimport dtype_t, dtype_t_out, _core, _min, _max
 
 
-cdef inline double _kernel_autolevel(Py_ssize_t* histo, double pop, dtype_t g,
-                                     Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                     double p0, double p1,
-                                     Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_autolevel(dtype_t_out* out, Py_ssize_t odepth,
+                                   Py_ssize_t* histo,
+                                   double pop, dtype_t g,
+                                   Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                   double p0, double p1,
+                                   Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i, imin, imax, sum, delta
 
@@ -31,18 +33,20 @@ cdef inline double _kernel_autolevel(Py_ssize_t* histo, double pop, dtype_t g,
 
         delta = imax - imin
         if delta > 0:
-            return <double>(max_bin - 1) * (_min(_max(imin, g), imax)
-                                           - imin) / delta
+            out[0] = <dtype_t_out>((max_bin - 1) * (_min(_max(imin, g), imax)
+                                           - imin) / delta)
         else:
-            return imax - imin
+            out[0] = <dtype_t_out>(imax - imin)
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_gradient(Py_ssize_t* histo, double pop, dtype_t g,
-                                    Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                    double p0, double p1,
-                                    Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_gradient(dtype_t_out* out, Py_ssize_t odepth,
+                                  Py_ssize_t* histo,
+                                  double pop, dtype_t g,
+                                  Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                  double p0, double p1,
+                                  Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i, imin, imax, sum, delta
 
@@ -61,15 +65,17 @@ cdef inline double _kernel_gradient(Py_ssize_t* histo, double pop, dtype_t g,
                 imax = i
                 break
 
-        return imax - imin
+        out[0] = <dtype_t_out>(imax - imin)
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_mean(Py_ssize_t* histo, double pop, dtype_t g,
-                                Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                double p0, double p1,
-                                Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth,
+                              Py_ssize_t* histo,
+                              double pop, dtype_t g,
+                              Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                              double p0, double p1,
+                              Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i, sum, mean, n
 
@@ -84,16 +90,18 @@ cdef inline double _kernel_mean(Py_ssize_t* histo, double pop, dtype_t g,
                 mean += histo[i] * i
 
         if n > 0:
-            return mean / n
+            out[0] = <dtype_t_out>(mean / n)
         else:
-            return 0
+            out[0] = <dtype_t_out>0
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
-cdef inline double _kernel_sum(Py_ssize_t* histo, double pop, dtype_t g,
-                               Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                               double p0, double p1,
-                               Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth,
+                             Py_ssize_t* histo,
+                             double pop, dtype_t g,
+                             Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                             double p0, double p1,
+                             Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i, sum, sum_g, n
 
@@ -108,18 +116,19 @@ cdef inline double _kernel_sum(Py_ssize_t* histo, double pop, dtype_t g,
                 sum_g += histo[i] * i
 
         if n > 0:
-            return sum_g
+            out[0] = <dtype_t_out>(sum_g)
         else:
-            return 0
+            out[0] = <dtype_t_out>0
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
-cdef inline double _kernel_subtract_mean(Py_ssize_t* histo, double pop,
-                                         dtype_t g,
-                                         Py_ssize_t max_bin,
-                                         Py_ssize_t mid_bin, double p0,
-                                         double p1, Py_ssize_t s0,
-                                         Py_ssize_t s1):
+cdef inline void _kernel_subtract_mean(dtype_t_out* out, Py_ssize_t odepth,
+                                       Py_ssize_t* histo,
+                                       double pop, dtype_t g,
+                                       Py_ssize_t max_bin,
+                                       Py_ssize_t mid_bin, double p0,
+                                       double p1, Py_ssize_t s0,
+                                       Py_ssize_t s1):
 
     cdef Py_ssize_t i, sum, mean, n
 
@@ -133,19 +142,21 @@ cdef inline double _kernel_subtract_mean(Py_ssize_t* histo, double pop,
                 n += histo[i]
                 mean += histo[i] * i
         if n > 0:
-            return (g - (mean / n)) * .5 + mid_bin
+            out[0] = <dtype_t_out>((g - (mean / n)) * .5 + mid_bin)
         else:
-            return 0
+            out[0] = <dtype_t_out>0
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_enhance_contrast(Py_ssize_t* histo, double pop,
-                                            dtype_t g,
-                                            Py_ssize_t max_bin,
-                                            Py_ssize_t mid_bin, double p0,
-                                            double p1, Py_ssize_t s0,
-                                            Py_ssize_t s1):
+cdef inline void _kernel_enhance_contrast(dtype_t_out* out,
+                                          Py_ssize_t odepth,
+                                          Py_ssize_t* histo, double pop,
+                                          dtype_t g,
+                                          Py_ssize_t max_bin,
+                                          Py_ssize_t mid_bin, double p0,
+                                          double p1, Py_ssize_t s0,
+                                          Py_ssize_t s1):
 
     cdef Py_ssize_t i, imin, imax, sum, delta
 
@@ -164,21 +175,23 @@ cdef inline double _kernel_enhance_contrast(Py_ssize_t* histo, double pop,
                 imax = i
                 break
         if g > imax:
-            return imax
+            out[0] = <dtype_t_out>imax
         if g < imin:
-            return imin
+            out[0] = <dtype_t_out>imin
         if imax - g < g - imin:
-            return imax
+            out[0] = <dtype_t_out>imax
         else:
-            return imin
+            out[0] = <dtype_t_out>imin
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_percentile(Py_ssize_t* histo, double pop, dtype_t g,
-                                      Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                      double p0, double p1,
-                                      Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_percentile(dtype_t_out* out, Py_ssize_t odepth,
+                                    Py_ssize_t* histo,
+                                    double pop, dtype_t g,
+                                    Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                    double p0, double p1,
+                                    Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i
     cdef Py_ssize_t sum = 0
@@ -193,15 +206,17 @@ cdef inline double _kernel_percentile(Py_ssize_t* histo, double pop, dtype_t g,
                 sum += histo[i]
                 if sum > p0 * pop:
                     break
-        return i
+        out[0] = <dtype_t_out>i
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g,
-                               Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                               double p0, double p1,
-                               Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth,
+                             Py_ssize_t* histo,
+                             double pop, dtype_t g,
+                             Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                             double p0, double p1,
+                             Py_ssize_t s0, Py_ssize_t s1):
 
     cdef Py_ssize_t i, sum, n
 
@@ -212,15 +227,17 @@ cdef inline double _kernel_pop(Py_ssize_t* histo, double pop, dtype_t g,
             sum += histo[i]
             if (sum >= p0 * pop) and (sum <= p1 * pop):
                 n += histo[i]
-        return n
+        out[0] = <dtype_t_out>n
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
-cdef inline double _kernel_threshold(Py_ssize_t* histo, double pop, dtype_t g,
-                                     Py_ssize_t max_bin, Py_ssize_t mid_bin,
-                                     double p0, double p1,
-                                     Py_ssize_t s0, Py_ssize_t s1):
+cdef inline void _kernel_threshold(dtype_t_out* out, Py_ssize_t odepth,
+                                   Py_ssize_t* histo,
+                                   double pop, dtype_t g,
+                                   Py_ssize_t max_bin, Py_ssize_t mid_bin,
+                                   double p0, double p1,
+                                   Py_ssize_t s0, Py_ssize_t s1):
 
     cdef int i
     cdef Py_ssize_t sum = 0
@@ -231,103 +248,105 @@ cdef inline double _kernel_threshold(Py_ssize_t* histo, double pop, dtype_t g,
             if sum >= p0 * pop:
                 break
 
-        return (max_bin - 1) * (g >= i)
+        out[0] = <dtype_t_out>((max_bin - 1) * (g >= i))
     else:
-        return 0
+        out[0] = <dtype_t_out>0
 
 
 def _autolevel(dtype_t[:, ::1] image,
                char[:, ::1] selem,
                char[:, ::1] mask,
-               dtype_t_out[:, ::1] out,
+               dtype_t_out[:, :, ::1] out,
                char shift_x, char shift_y, double p0, double p1,
                Py_ssize_t max_bin):
 
-    _core(_kernel_autolevel[dtype_t], image, selem, mask, out,
+    _core(_kernel_autolevel[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, p0, p1, 0, 0, max_bin)
 
 
 def _gradient(dtype_t[:, ::1] image,
               char[:, ::1] selem,
               char[:, ::1] mask,
-              dtype_t_out[:, ::1] out,
+              dtype_t_out[:, :, ::1] out,
               char shift_x, char shift_y, double p0, double p1,
               Py_ssize_t max_bin):
 
-    _core(_kernel_gradient[dtype_t], image, selem, mask, out,
+    _core(_kernel_gradient[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, p0, p1, 0, 0, max_bin)
 
 
 def _mean(dtype_t[:, ::1] image,
           char[:, ::1] selem,
           char[:, ::1] mask,
-          dtype_t_out[:, ::1] out,
+          dtype_t_out[:, :, ::1] out,
           char shift_x, char shift_y, double p0, double p1,
           Py_ssize_t max_bin):
 
-    _core(_kernel_mean[dtype_t], image, selem, mask, out,
+    _core(_kernel_mean[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, p0, p1, 0, 0, max_bin)
 
+
 def _sum(dtype_t[:, ::1] image,
          char[:, ::1] selem,
          char[:, ::1] mask,
-         dtype_t_out[:, ::1] out,
+         dtype_t_out[:, :, ::1] out,
          char shift_x, char shift_y, double p0, double p1,
          Py_ssize_t max_bin):
 
-    _core(_kernel_sum[dtype_t], image, selem, mask, out,
+    _core(_kernel_sum[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, p0, p1, 0, 0, max_bin)
 
+
 def _subtract_mean(dtype_t[:, ::1] image,
                    char[:, ::1] selem,
                    char[:, ::1] mask,
-                   dtype_t_out[:, ::1] out,
+                   dtype_t_out[:, :, ::1] out,
                    char shift_x, char shift_y, double p0, double p1,
                    Py_ssize_t max_bin):
 
-    _core(_kernel_subtract_mean[dtype_t], image, selem, mask,
+    _core(_kernel_subtract_mean[dtype_t_out, dtype_t], image, selem, mask,
           out, shift_x, shift_y, p0, p1, 0, 0, max_bin)
 
 
 def _enhance_contrast(dtype_t[:, ::1] image,
                       char[:, ::1] selem,
                       char[:, ::1] mask,
-                      dtype_t_out[:, ::1] out,
+                      dtype_t_out[:, :, ::1] out,
                       char shift_x, char shift_y, double p0, double p1,
                       Py_ssize_t max_bin):
 
-    _core(_kernel_enhance_contrast[dtype_t], image, selem, mask,
+    _core(_kernel_enhance_contrast[dtype_t_out, dtype_t], image, selem, mask,
           out, shift_x, shift_y, p0, p1, 0, 0, max_bin)
 
 
 def _percentile(dtype_t[:, ::1] image,
                 char[:, ::1] selem,
                 char[:, ::1] mask,
-                dtype_t_out[:, ::1] out,
+                dtype_t_out[:, :, ::1] out,
                 char shift_x, char shift_y, double p0, double p1,
                 Py_ssize_t max_bin):
 
-    _core(_kernel_percentile[dtype_t], image, selem, mask, out,
+    _core(_kernel_percentile[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, p0, 1, 0, 0, max_bin)
 
 
 def _pop(dtype_t[:, ::1] image,
          char[:, ::1] selem,
          char[:, ::1] mask,
-         dtype_t_out[:, ::1] out,
+         dtype_t_out[:, :, ::1] out,
          char shift_x, char shift_y, double p0, double p1,
          Py_ssize_t max_bin):
 
-    _core(_kernel_pop[dtype_t], image, selem, mask, out,
+    _core(_kernel_pop[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, p0, p1, 0, 0, max_bin)
 
 
 def _threshold(dtype_t[:, ::1] image,
                char[:, ::1] selem,
                char[:, ::1] mask,
-               dtype_t_out[:, ::1] out,
+               dtype_t_out[:, :, ::1] out,
                char shift_x, char shift_y, double p0, double p1,
                Py_ssize_t max_bin):
 
-    _core(_kernel_threshold[dtype_t], image, selem, mask, out,
+    _core(_kernel_threshold[dtype_t_out, dtype_t], image, selem, mask, out,
           shift_x, shift_y, p0, 1, 0, 0, max_bin)
diff --git a/skimage/filter/rank/tests/test_rank.py b/skimage/filter/rank/tests/test_rank.py
index 5cbffd90..48d02462 100644
--- a/skimage/filter/rank/tests/test_rank.py
+++ b/skimage/filter/rank/tests/test_rank.py
@@ -542,6 +542,45 @@ def test_sum():
     rank.sum_bilateral(image=image16, selem=elem, out=out16, mask=mask,s0=1000,s1=1000)
     assert_array_equal(r, out16)
 
+def test_windowed_histogram():
+    # check the number of valid pixels in the neighborhood
+
+    image8 = 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)
+    elem = np.ones((3, 3), dtype=np.uint8)
+    outf = np.empty(image8.shape+(2,), dtype=float)
+    mask = np.ones(image8.shape, dtype=np.uint8)
+
+    # Population so we can normalize the expected output while maintaining
+    # code readability
+    pop = np.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=float)
+
+    r0 =  np.array([[3, 4, 3, 4, 3],
+                    [4, 5, 3, 5, 4],
+                    [3, 3, 0, 3, 3],
+                    [4, 5, 3, 5, 4],
+                    [3, 4, 3, 4, 3]], dtype=float) / pop
+    r1 =  np.array([[1, 2, 3, 2, 1],
+                    [2, 4, 6, 4, 2],
+                    [3, 6, 9, 6, 3],
+                    [2, 4, 6, 4, 2],
+                    [1, 2, 3, 2, 1]], dtype=float) / pop
+    rank.windowed_histogram(image=image8, selem=elem, out=outf, mask=mask)
+    assert_array_equal(r0, outf[:,:,0])
+    assert_array_equal(r1, outf[:,:,1])
+
+    # Test n_bins parameter
+    larger_output = rank.windowed_histogram(image=image8, selem=elem,
+                                            mask=mask, n_bins=5)
+    assert larger_output.shape[2] == 5
+
 
 if __name__ == "__main__":
     run_module_suite()
diff --git a/skimage/io/__init__.py b/skimage/io/__init__.py
index 144fef57..85461216 100644
--- a/skimage/io/__init__.py
+++ b/skimage/io/__init__.py
@@ -40,7 +40,11 @@ def _update_doc(doc):
     info_table = [(p, plugin_info(p).get('description', 'no description'))
                   for p in available_plugins if not p == 'test']
 
-    name_length = max([len(n) for (n, _) in info_table])
+    if len(info_table) > 0:
+        name_length = max([len(n) for (n, _) in info_table])
+    else:
+        name_length = 0
+
     description_length = WRAP_LEN - 1 - name_length
     column_lengths = [name_length, description_length]
     _format_plugin_info_table(info_table, column_lengths)