diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index 012fd53d..fd3ff235 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -123,3 +123,6 @@
 
 - Luis Pedro Coelho
   imread plugin
+
+- Steven Silvester, Karel Zuiderveld
+  Adaptive Histogram Equalization
diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py
index 1e1291d8..f57cea9c 100644
--- a/doc/examples/plot_equalize.py
+++ b/doc/examples/plot_equalize.py
@@ -18,18 +18,19 @@ that fall within the 2nd and 98th percentiles [2]_.
 
 """
 
-from skimage import data
+from skimage import data, img_as_float
 from skimage.util.dtype import dtype_range
 from skimage import exposure
 
 import matplotlib.pyplot as plt
-
 import numpy as np
 
+
 def plot_img_and_hist(img, axes, bins=256):
     """Plot an image along with its histogram and cumulative histogram.
 
     """
+    img = img_as_float(img)
     ax_img, ax_hist = axes
     ax_cdf = ax_hist.twinx()
 
@@ -38,16 +39,16 @@ def plot_img_and_hist(img, axes, bins=256):
     ax_img.set_axis_off()
 
     # Display histogram
-    ax_hist.hist(img.ravel(), bins=bins)
+    ax_hist.hist(img.ravel(), bins=bins, histtype='step', color='black')
     ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
     ax_hist.set_xlabel('Pixel intensity')
-
-    xmin, xmax = dtype_range[img.dtype.type]
-    ax_hist.set_xlim(xmin, xmax)
+    ax_hist.set_xlim(0, 1)
+    ax_hist.set_yticks([])
 
     # Display cumulative distribution
     img_cdf, bins = exposure.cumulative_distribution(img, bins)
     ax_cdf.plot(bins, img_cdf, 'r')
+    ax_cdf.set_yticks([])
 
     return ax_img, ax_hist, ax_cdf
 
@@ -61,25 +62,33 @@ p98 = np.percentile(img, 98)
 img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))
 
 # Equalization
-img_eq = exposure.equalize(img)
+img_eq = exposure.equalize_hist(img)
 
+# Adaptive Equalization
+img_adapteq = exposure.equalize_adapthist(img, clip_limit=0.03)
 
 # Display results
-f, axes = plt.subplots(2, 3, figsize=(8, 4))
+f, axes = plt.subplots(2, 4, figsize=(8, 4))
 
 ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0])
 ax_img.set_title('Low contrast image')
+
+y_min, y_max = ax_hist.get_ylim()
 ax_hist.set_ylabel('Number of pixels')
+ax_hist.set_yticks(np.linspace(0, y_max, 5))
 
 ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_rescale, axes[:, 1])
 ax_img.set_title('Contrast stretching')
 
 ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_eq, axes[:, 2])
 ax_img.set_title('Histogram equalization')
-ax_cdf.set_ylabel('Fraction of total intensity')
 
+ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_adapteq, axes[:, 3])
+ax_img.set_title('Adaptive equalization')
+
+ax_cdf.set_ylabel('Fraction of total intensity')
+ax_cdf.set_yticks(np.linspace(0, 1, 5))
 
 # prevent overlap of y-axis labels
 plt.subplots_adjust(wspace=0.4)
 plt.show()
-
diff --git a/skimage/exposure/__init__.py b/skimage/exposure/__init__.py
index 7e19d317..ae75c982 100644
--- a/skimage/exposure/__init__.py
+++ b/skimage/exposure/__init__.py
@@ -1,2 +1,3 @@
-from .exposure import histogram, equalize, cumulative_distribution
-from .exposure import rescale_intensity
+from .exposure import histogram, equalize, equalize_hist
+from .exposure import rescale_intensity, cumulative_distribution
+from ._adapthist import equalize_adapthist
diff --git a/skimage/exposure/_adapthist.py b/skimage/exposure/_adapthist.py
new file mode 100644
index 00000000..667d2cb2
--- /dev/null
+++ b/skimage/exposure/_adapthist.py
@@ -0,0 +1,325 @@
+"""
+Adapted code from "Contrast Limited Adaptive Histogram Equalization" by Karel
+Zuiderveld <karel@cv.ruu.nl>, Graphics Gems IV, Academic Press, 1994.
+
+http://tog.acm.org/resources/GraphicsGems/gems.html#gemsvi
+
+The Graphics Gems code is copyright-protected.  In other words, you cannot
+claim the text of the code as your own and resell it. Using the code is
+permitted in any program, product, or library, non-commercial or commercial.
+Giving credit is not required, though is a nice gesture.  The code comes as-is,
+and if there are any flaws or problems with any Gems code, nobody involved with
+Gems - authors, editors, publishers, or webmasters - are to be held
+responsible.  Basically, don't be a jerk, and remember that anything free
+comes with no guarantee.
+"""
+import numpy as np
+import skimage
+from skimage import color
+from skimage.exposure import rescale_intensity
+from skimage.util import view_as_blocks
+
+
+MAX_REG_X = 16  # max. # contextual regions in x-direction */
+MAX_REG_Y = 16  # max. # contextual regions in y-direction */
+NR_OF_GREY = 16384  # number of grayscale levels to use in CLAHE algorithm
+
+
+def equalize_adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01,
+                       nbins=256):
+    """Contrast Limited Adaptive Histogram Equalization.
+
+    Parameters
+    ----------
+    image : array-like
+        Input image.
+    ntiles_x : int, optional
+        Number of tile regions in the X direction.  Ranges between 2 and 16.
+    ntiles_y : int, optional
+        Number of tile regions in the Y direction.  Ranges between 2 and 16.
+    clip_limit : float: optional
+        Clipping limit, normalized between 0 and 1 (higher values give more
+        contrast).
+    nbins : int, optional
+        Number of gray bins for histogram ("dynamic range").
+
+    Returns
+    -------
+    out : ndarray
+        Equalized image.
+
+    Notes
+    -----
+    * The algorithm relies on an image whose rows and columns are even
+      multiples of the number of tiles, so the extra rows and columns are left
+      at their original values, thus  preserving the input image shape.
+    * For color images, the following steps are performed:
+       - The image is converted to LAB color space
+       - The CLAHE algorithm is run on the L channel
+       - The image is converted back to RGB space and returned
+    * For RGBA images, the original alpha channel is removed.
+
+    References
+    ----------
+    .. [1] http://tog.acm.org/resources/GraphicsGems/gems.html#gemsvi
+    .. [2] https://en.wikipedia.org/wiki/CLAHE#CLAHE
+    """
+    args = [None, ntiles_x, ntiles_y, clip_limit * nbins, nbins]
+    if image.ndim > 2:
+        lab_img = color.rgb2lab(skimage.img_as_float(image))
+        l_chan = lab_img[:, :, 0]
+        l_chan /= np.max(np.abs(l_chan))
+        l_chan = skimage.img_as_uint(l_chan)
+        args[0] = rescale_intensity(l_chan, out_range=(0, NR_OF_GREY - 1))
+        new_l = _clahe(*args).astype(float)
+        new_l = rescale_intensity(new_l, out_range=(0, 100))
+        lab_img[:new_l.shape[0], :new_l.shape[1], 0] = new_l
+        image = color.lab2rgb(lab_img)
+        image = rescale_intensity(image, out_range=(0, 1))
+    else:
+        image = skimage.img_as_uint(image)
+        args[0] = rescale_intensity(image, out_range=(0, NR_OF_GREY - 1))
+        out = _clahe(*args)
+        image[:out.shape[0], :out.shape[1]] = out
+        image = rescale_intensity(image)
+    return image
+
+
+def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
+    """Contrast Limited Adaptive Histogram Equalization.
+
+    Parameters
+    ----------
+    image : array-like
+        Input image.
+    ntiles_x : int, optional
+        Number of tile regions in the X direction.  Ranges between 2 and 16.
+    ntiles_y : int, optional
+        Number of tile regions in the Y direction.  Ranges between 2 and 16.
+    clip_limit : float, optional
+        Normalized clipping limit (higher values give more contrast).
+    nbins : int, optional
+        Number of gray bins for histogram ("dynamic range").
+
+    Returns
+    -------
+    out : ndarray
+        Equalized image.
+
+    The number of "effective" greylevels in the output image is set by `nbins`;
+    selecting a small value (eg. 128) speeds up processing and still produce
+    an output image of good quality. The output image will have the same
+    minimum and maximum value as the input image. A clip limit smaller than 1
+    results in standard (non-contrast limited) AHE.
+    """
+    ntiles_x = min(ntiles_x, MAX_REG_X)
+    ntiles_y = min(ntiles_y, MAX_REG_Y)
+    ntiles_y = max(ntiles_x, 2)
+    ntiles_x = max(ntiles_y, 2)
+
+    if clip_limit == 1.0:
+        return image  # is OK, immediately returns original image.
+
+    map_array = np.zeros((ntiles_y, ntiles_x, nbins), dtype=int)
+
+    y_res = image.shape[0] - image.shape[0] % ntiles_y
+    x_res = image.shape[1] - image.shape[1] % ntiles_x
+    image = image[: y_res, : x_res]
+
+    x_size = image.shape[1] / ntiles_x  # Actual size of contextual regions
+    y_size = image.shape[0] / ntiles_y
+    n_pixels = x_size * y_size
+
+    if clip_limit > 0.0:  # Calculate actual cliplimit
+        clip_limit = int(clip_limit * (x_size * y_size) / nbins)
+        if clip_limit < 1:
+            clip_limit = 1
+    else:
+        clip_limit = NR_OF_GREY  # Large value, do not clip (AHE)
+
+    bin_size = 1 + NR_OF_GREY / nbins
+    aLUT = np.arange(NR_OF_GREY)
+    aLUT /= bin_size
+    img_blocks = view_as_blocks(image, (y_size, x_size))
+
+    # Calculate greylevel mappings for each contextual region
+    for y in range(ntiles_y):
+        for x in range(ntiles_x):
+            sub_img = img_blocks[y, x]
+            hist = aLUT[sub_img.ravel()]
+            hist = np.bincount(hist)
+            hist = np.append(hist, np.zeros(nbins - hist.size, dtype=int))
+            hist = clip_histogram(hist, clip_limit)
+            hist = map_histogram(hist, 0, NR_OF_GREY, n_pixels)
+            map_array[y, x] = hist
+
+    # Interpolate greylevel mappings to get CLAHE image
+    ystart = 0
+    for y in range(ntiles_y + 1):
+        xstart = 0
+        if y == 0:  # special case: top row
+            ystep = y_size / 2
+            yU = 0
+            yB = 0
+        elif y == ntiles_y:  # special case: bottom row
+            ystep = y_size / 2
+            yU = ntiles_y - 1
+            yB = yU
+        else:  # default values
+            ystep = y_size
+            yU = y - 1
+            yB = yB + 1
+
+        for x in range(ntiles_x + 1):
+            if x == 0:  # special case: left column
+                xstep = x_size / 2
+                xL = 0
+                xR = 0
+            elif x == ntiles_x:  # special case: right column
+                xstep = x_size / 2
+                xL = ntiles_x - 1
+                xR = xL
+            else:  # default values
+                xstep = x_size
+                xL = x - 1
+                xR = xL + 1
+
+            mapLU = map_array[yU, xL]
+            mapRU = map_array[yU, xR]
+            mapLB = map_array[yB, xL]
+            mapRB = map_array[yB, xR]
+
+            xslice = np.arange(xstart, xstart + xstep)
+            yslice = np.arange(ystart, ystart + ystep)
+            interpolate(image, xslice, yslice,
+                        mapLU, mapRU, mapLB, mapRB, aLUT)
+
+            xstart += xstep  # set pointer on next matrix */
+
+        ystart += ystep
+
+    return image
+
+
+def clip_histogram(hist, clip_limit):
+    """Perform clipping of the histogram and redistribution of bins.
+
+    The histogram is clipped and the number of excess pixels is counted.
+    Afterwards the excess pixels are equally redistributed across the
+    whole histogram (providing the bin count is smaller than the cliplimit).
+
+    Parameters
+    ----------
+    hist : ndarray
+        Histogram array.
+    clip_limit : int
+        Maximum allowed bin count.
+
+    Returns
+    -------
+    hist : ndarray
+        Clipped histogram.
+    """
+    # calculate total number of excess pixels
+    excess_mask = hist > clip_limit
+    excess = hist[excess_mask]
+    n_excess = excess.sum() - excess.size * clip_limit
+
+    # Second part: clip histogram and redistribute excess pixels in each bin
+    bin_incr = n_excess / hist.size  # average binincrement
+    upper = clip_limit - bin_incr  # Bins larger than upper set to cliplimit
+
+    hist[excess_mask] = clip_limit
+
+    low_mask = hist < upper
+    n_excess -= hist[low_mask].size * bin_incr
+    hist[low_mask] += bin_incr
+
+    mid_mask = (hist >= upper) & (hist < clip_limit)
+    mid = hist[mid_mask]
+    n_excess -= mid.size * clip_limit - mid.sum()
+    hist[mid_mask] = clip_limit
+
+    while n_excess > 0:  # Redistribute remaining excess
+        index = 0
+        while n_excess > 0 and index < hist.size:
+            step_size = int(hist[hist < clip_limit].size / n_excess)
+            step_size = max(step_size, 1)
+            indices = np.arange(index, hist.size, step_size)
+            under = hist[indices] < clip_limit
+            hist[under] += 1
+            n_excess -= hist[under].size
+            index += 1
+
+    return hist
+
+
+def map_histogram(hist, min_val, max_val, n_pixels):
+    """Calculate the equalized lookup table (mapping).
+
+    It does so by cumulating the input histogram.
+
+    hist : ndarray
+        Clipped histogram.
+    min_val : int
+        Minimum value for mapping.
+    max_val : int
+        Maximum value for mapping.
+    n_pixels : int
+        Number of pixels in the region.
+
+    Returns
+    -------
+    out : ndarray
+       Mapped intensity LUT.
+    """
+    out = np.cumsum(hist).astype(float)
+    scale = ((float)(max_val - min_val)) / n_pixels
+    out *= scale
+    out += min_val
+    out[out > max_val] = max_val
+    return out.astype(int)
+
+
+def interpolate(image, xslice, yslice,
+                mapLU, mapRU, mapLB, mapRB, aLUT):
+    """Find the new grayscale level for a region using bilinear interpolation.
+
+    Parameters
+    ----------
+    image : ndarray
+        Full image.
+    xslice, yslice : array-like
+       Indices of the region.
+    map* : ndarray
+        Mappings of greylevels from histograms.
+    aLUT : ndarray
+        Maps grayscale levels in image to histogram levels.
+
+    Returns
+    -------
+    out : ndarray
+        Original image with the subregion replaced.
+
+    Note
+    ----
+    This function calculates the new greylevel assignments of pixels
+    within a submatrix of the image.
+    This is done by a bilinear interpolation between four different
+    mappings in order to eliminate boundary artifacts.
+    """
+    norm = xslice.size * yslice.size  # Normalization factor
+    # interpolation weight matrices
+    x_coef, y_coef = np.meshgrid(np.arange(xslice.size),
+                                 np.arange(yslice.size))
+    x_inv_coef, y_inv_coef = x_coef[:, ::-1] + 1, y_coef[::-1] + 1
+
+    view = image[yslice[0]: yslice[-1] + 1, xslice[0]: xslice[-1] + 1]
+    im_slice = aLUT[view]
+    new = ((y_inv_coef * (x_inv_coef * mapLU[im_slice]
+                          + x_coef * mapRU[im_slice])
+            + y_coef * (x_inv_coef * mapLB[im_slice]
+                        + x_coef * mapRB[im_slice]))
+           / norm)
+    view[:, :] = new
+    return image
diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index bffe660a..465b9182 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -2,6 +2,9 @@ import numpy as np
 
 from skimage import img_as_float
 from skimage.util.dtype import dtype_range
+import skimage.color as color
+from skimage.util.dtype import convert
+from skimage._shared.utils import deprecated
 
 
 __all__ = ['histogram', 'cumulative_distribution', 'equalize',
@@ -76,7 +79,12 @@ def cumulative_distribution(image, nbins=256):
     return img_cdf, bin_centers
 
 
+@deprecated('equalize_hist')
 def equalize(image, nbins=256):
+    equalize_hist(image, nbins)
+
+
+def equalize_hist(image, nbins=256):
     """Return image after histogram equalization.
 
     Parameters
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index b6ed0b12..258b4224 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -1,9 +1,10 @@
 import numpy as np
 from numpy.testing import assert_array_almost_equal as assert_close
-
 import skimage
 from skimage import data
 from skimage import exposure
+from skimage.color import rgb2gray
+from skimage.util.dtype import dtype_range
 
 
 # Test histogram equalization
@@ -15,7 +16,7 @@ test_img = exposure.rescale_intensity(data.camera() / 5. + 100)
 
 def test_equalize_ubyte():
     img = skimage.img_as_ubyte(test_img)
-    img_eq = exposure.equalize(img)
+    img_eq = exposure.equalize_hist(img)
 
     cdf, bin_edges = exposure.cumulative_distribution(img_eq)
     check_cdf_slope(cdf)
@@ -23,7 +24,7 @@ def test_equalize_ubyte():
 
 def test_equalize_float():
     img = skimage.img_as_float(test_img)
-    img_eq = exposure.equalize(img)
+    img_eq = exposure.equalize_hist(img)
 
     cdf, bin_edges = exposure.cumulative_distribution(img_eq)
     check_cdf_slope(cdf)
@@ -71,6 +72,99 @@ def test_rescale_out_range():
     assert_close(out, [0, 63, 127])
 
 
+# Test adaptive histogram equalization
+# ====================================
+
+def test_adapthist_scalar():
+    '''Test a scalar uint8 image
+    '''
+    img = skimage.img_as_ubyte(data.moon())
+    adapted = exposure.equalize_adapthist(img, clip_limit=0.02)
+    assert adapted.min() == 0
+    assert adapted.max() == (1 << 16) - 1
+    assert img.shape == adapted.shape
+    full_scale = skimage.exposure.rescale_intensity(skimage.img_as_uint(img))
+    assert_almost_equal = np.testing.assert_almost_equal
+    assert_almost_equal(peak_snr(full_scale, adapted), 101.231, 3)
+    assert_almost_equal(norm_brightness_err(full_scale, adapted),
+                        0.041, 3)
+    return img, adapted
+
+
+def test_adapthist_grayscale():
+    '''Test a grayscale float image
+    '''
+    img = skimage.img_as_float(data.lena())
+    img = rgb2gray(img)
+    img = np.dstack((img, img, img))
+    adapted = exposure.equalize_adapthist(img, 10, 9, clip_limit=0.01,
+                        nbins=128)
+    assert_almost_equal = np.testing.assert_almost_equal
+    assert img.shape == adapted.shape
+    assert_almost_equal(peak_snr(img, adapted), 77.584, 3)
+    assert_almost_equal(norm_brightness_err(img, adapted), 0.038, 3)
+    return data, adapted
+
+
+def test_adapthist_color():
+    '''Test a color uint16 image
+    '''
+    img = skimage.img_as_uint(data.lena())
+    adapted = exposure.equalize_adapthist(img, clip_limit=0.01)
+    assert_almost_equal = np.testing.assert_almost_equal
+    assert adapted.min() == 0
+    assert adapted.max() == 1.0
+    assert img.shape == adapted.shape
+    full_scale = skimage.exposure.rescale_intensity(img)
+    assert_almost_equal(peak_snr(full_scale, adapted), 64.717, 3)
+    assert_almost_equal(norm_brightness_err(full_scale, adapted),
+                        0.179, 3)
+    return data, adapted
+
+
+def peak_snr(img1, img2):
+    '''Peak signal to noise ratio of two images
+
+    Parameters
+    ----------
+    img1 : array-like
+    img2 : array-like
+
+    Returns
+    -------
+    peak_snr : float
+        Peak signal to noise ratio
+    '''
+    if img1.ndim == 3:
+        img1, img2 = rgb2gray(img1.copy()), rgb2gray(img2.copy())
+    img1 = skimage.img_as_float(img1)
+    img2 = skimage.img_as_float(img2)
+    mse = 1. / img1.size * np.square(img1 - img2).sum()
+    _, max_ = dtype_range[img1.dtype.type]
+    print mse, max_
+    return 20 * np.log(max_ / mse)
+
+
+def norm_brightness_err(img1, img2):
+    '''Normalized Absolute Mean Brightness Error between two images
+
+    Parameters
+    ----------
+    img1 : array-like
+    img2 : array-like
+
+    Returns
+    -------
+    norm_brightness_error : float
+        Normalized absolute mean brightness error
+    '''
+    if img1.ndim == 3:
+        img1, img2 = rgb2gray(img1), rgb2gray(img2)
+    ambe = np.abs(img1.mean() - img2.mean())
+    nbe = ambe / dtype_range[img1.dtype.type][1]
+    return nbe
+
+
 if __name__ == '__main__':
     from numpy import testing
     testing.run_module_suite()