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..73147c7f 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,
diff --git a/skimage/filter/rank/bilateral_cy.pyx b/skimage/filter/rank/bilateral_cy.pyx
index 25a81766..0e232465 100644
--- a/skimage/filter/rank/bilateral_cy.pyx
+++ b/skimage/filter/rank/bilateral_cy.pyx
@@ -9,10 +9,11 @@ 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* 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 +25,18 @@ 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* 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 +45,15 @@ 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* 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 +65,40 @@ 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..7d23f12a 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[:] out, 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..3c8c9e42 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[:] out_feat, 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,
@@ -151,8 +151,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, :], histo, pop, image[r, c], max_bin, mid_bin,
+           p0, p1, s0, s1)
 
     # main loop
     r = 0
@@ -172,8 +172,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, :], 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 +192,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, :], 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 +209,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, :], 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 +229,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, :], 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..07c592af 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)
@@ -45,7 +45,10 @@ def _handle_input(image, selem, out, mask, out_dtype=None):
     if out is None:
         if out_dtype is None:
             out_dtype = image.dtype
-        out = np.empty_like(image, dtype=out_dtype)
+        out = np.empty(image.shape+(pixel_size,), dtype=out_dtype)
+    else:
+        if len(out.shape) == 2:
+            out = out.reshape(out.shape+(pixel_size,))
 
     if image is out:
         raise NotImplementedError("Cannot perform rank operation in place.")
@@ -65,7 +68,7 @@ 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 +76,17 @@ 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,7 +127,7 @@ def autolevel(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._autolevel, image, selem,
+    return _apply_scalar_per_pixel(generic_cy._autolevel, image, selem,
                   out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -154,7 +168,7 @@ def bottomhat(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._bottomhat, image, selem,
+    return _apply_scalar_per_pixel(generic_cy._bottomhat, image, selem,
                   out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -192,7 +206,7 @@ def equalize(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._equalize, image, selem,
+    return _apply_scalar_per_pixel(generic_cy._equalize, image, selem,
                   out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -230,7 +244,7 @@ def gradient(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._gradient, image, selem,
+    return _apply_scalar_per_pixel(generic_cy._gradient, image, selem,
                   out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -277,7 +291,7 @@ def maximum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._maximum, image, selem,
+    return _apply_scalar_per_pixel(generic_cy._maximum, image, selem,
                   out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -315,7 +329,7 @@ def mean(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._mean, image, selem, out=out,
+    return _apply_scalar_per_pixel(generic_cy._mean, image, selem, out=out,
                   mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -354,7 +368,7 @@ def subtract_mean(image, selem, out=None, mask=None, shift_x=False,
 
     """
 
-    return _apply(generic_cy._subtract_mean, image, selem,
+    return _apply_scalar_per_pixel(generic_cy._subtract_mean, image, selem,
                   out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -392,7 +406,7 @@ def median(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._median, image, selem,
+    return _apply_scalar_per_pixel(generic_cy._median, image, selem,
                   out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -439,7 +453,7 @@ def minimum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._minimum, image, selem,
+    return _apply_scalar_per_pixel(generic_cy._minimum, image, selem,
                   out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -479,7 +493,7 @@ def modal(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._modal, image, selem,
+    return _apply_scalar_per_pixel(generic_cy._modal, image, selem,
                   out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -522,7 +536,7 @@ def enhance_contrast(image, selem, out=None, mask=None, shift_x=False,
 
     """
 
-    return _apply(generic_cy._enhance_contrast, image, selem,
+    return _apply_scalar_per_pixel(generic_cy._enhance_contrast, image, selem,
                   out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -571,7 +585,7 @@ def pop(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._pop, image, selem, out=out,
+    return _apply_scalar_per_pixel(generic_cy._pop, image, selem, out=out,
                   mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -620,7 +634,7 @@ def sum(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._sum, image, selem, out=out,
+    return _apply_scalar_per_pixel(generic_cy._sum, image, selem, out=out,
                   mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -669,7 +683,7 @@ def threshold(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._threshold, image, selem,
+    return _apply_scalar_per_pixel(generic_cy._threshold, image, selem,
                   out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -710,7 +724,7 @@ def tophat(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._tophat, image, selem,
+    return _apply_scalar_per_pixel(generic_cy._tophat, image, selem,
                   out=out, mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -761,7 +775,7 @@ 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,
+    return _apply_scalar_per_pixel(generic_cy._noise_filter, image, selem_cpy, out=out,
                   mask=mask, shift_x=shift_x, shift_y=shift_y)
 
 
@@ -806,7 +820,7 @@ def entropy(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._entropy, image, selem,
+    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)
 
@@ -850,5 +864,48 @@ def otsu(image, selem, out=None, mask=None, shift_x=False, shift_y=False):
 
     """
 
-    return _apply(generic_cy._otsu, image, selem, out=out,
+    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):
+    """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).
+
+    Returns
+    -------
+    out : 3-D array (same dtype as input image) whose extra dimension
+        Output image.
+
+    References
+    ----------
+    .. [otsu] http://en.wikipedia.org/wiki/Otsu's_method
+
+    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))
+    >>> thresh_image = img >= local_otsu
+
+    """
+
+    return _apply_vector_per_pixel(generic_cy._windowed_hist, image, selem, out=out,
+                  mask=mask, shift_x=shift_x, shift_y=shift_y, pixel_size=image.max()+1)
diff --git a/skimage/filter/rank/generic_cy.pyx b/skimage/filter/rank/generic_cy.pyx
index 1c26cf53..88766f44 100644
--- a/skimage/filter/rank/generic_cy.pyx
+++ b/skimage/filter/rank/generic_cy.pyx
@@ -9,10 +9,11 @@ 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* 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 +28,18 @@ 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* 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 +47,16 @@ 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* 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 +66,16 @@ 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* 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 +88,68 @@ 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* 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* 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* 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* 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 +158,34 @@ 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* 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* 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 +194,19 @@ 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* 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 +220,27 @@ 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* 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* 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 +248,16 @@ 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* 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 +265,16 @@ 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* 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 +282,23 @@ 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* 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 +308,16 @@ 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* 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 +327,16 @@ 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* 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 +348,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 +368,205 @@ 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* 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
+    for i in xrange(out.shape[0]):
+        out[i] = <dtype_t_out>histo[i]
+
 
 
 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..af33546d 100644
--- a/skimage/filter/rank/percentile_cy.pyx
+++ b/skimage/filter/rank/percentile_cy.pyx
@@ -7,10 +7,11 @@ 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* 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 +32,19 @@ 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* 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 +63,16 @@ 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* 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 +87,17 @@ 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* 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 +112,18 @@ 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* 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 +137,20 @@ 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* 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 +169,22 @@ 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* 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 +199,16 @@ 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* 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 +219,16 @@ 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* 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 +239,103 @@ 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)