From 34af272c54cfe69c1f8cd2a71be30a419f51773a Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sun, 5 Aug 2012 23:07:24 -0500
Subject: [PATCH 01/31] Initial commit of adapthist modifications.

---
 doc/examples/plot_equalize.py           |   8 +-
 skimage/exposure/_adapthist.pyx         |  40 +++
 skimage/exposure/clahe.h                | 329 ++++++++++++++++++++++++
 skimage/exposure/exposure.py            |  58 ++++-
 skimage/exposure/setup.py               |  31 +++
 skimage/exposure/tests/test_exposure.py |  71 +++++
 skimage/setup.py                        |   1 +
 7 files changed, 536 insertions(+), 2 deletions(-)
 create mode 100644 skimage/exposure/_adapthist.pyx
 create mode 100644 skimage/exposure/clahe.h
 create mode 100644 skimage/exposure/setup.py

diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py
index 1e1291d8..2cb9cbea 100644
--- a/doc/examples/plot_equalize.py
+++ b/doc/examples/plot_equalize.py
@@ -63,9 +63,12 @@ img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))
 # Equalization
 img_eq = exposure.equalize(img)
 
+# Adaptive Equalization
+img_adapteq = exposure.adapthist(img)
+
 
 # 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')
@@ -78,6 +81,9 @@ 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_eq, axes[:, 3])
+ax_img.set_title('Adaptive equalization')
+ax_cdf.set_ylabel('Fraction of total intensity')
 
 # prevent overlap of y-axis labels
 plt.subplots_adjust(wspace=0.4)
diff --git a/skimage/exposure/_adapthist.pyx b/skimage/exposure/_adapthist.pyx
new file mode 100644
index 00000000..16e84336
--- /dev/null
+++ b/skimage/exposure/_adapthist.pyx
@@ -0,0 +1,40 @@
+import numpy as np
+cimport numpy as np
+
+cdef extern from "clahe.h":
+    ctypedef unsigned int kz_pixel_t
+    int CLAHE(kz_pixel_t * pImage, unsigned int uiXRes, unsigned int uiYRes,
+        kz_pixel_t Min,
+        kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,
+        unsigned int uiNrBins, float fCliplimit)
+
+
+def _adapthist(np.ndarray[kz_pixel_t, ndim=2] image, min, max, nx, ny, nbins,
+                 clip_limit=0):
+    '''
+    image -  to the input/output image
+    min - Minimum greyvalue of input image (also becomes min of output image)
+    max - Maximum greyvalue of input image (also becomes max of output image)
+    nx - Number of tile regions in the X direction (min 2, max uiMAX_REG_X)
+    ny - Number of tile regions in the Y direction (min 2, max uiMAX_REG_Y)
+    nbins - Number of greybins for histogram ("dynamic range")
+    clip_limit - Normalized cliplimit (higher values give more contrast)
+
+    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.
+    '''
+    # we need the image to be divisible by nx and ny
+    uiYRes = image.shape[0] - image.shape[0] % ny
+    uiXRes = image.shape[1] - image.shape[1] % nx
+
+    cdef np.ndarray[kz_pixel_t, ndim = 1] flattened
+    flattened = image[:uiYRes, :uiXRes].ravel()
+
+    ret = CLAHE(< kz_pixel_t * > flattened.data, uiXRes, uiYRes, min,
+                max, nx, ny, nbins, clip_limit * nbins)
+    if ret < 0:
+        print 'clahe error: ', ret
+    return flattened.reshape((uiYRes, uiXRes))
diff --git a/skimage/exposure/clahe.h b/skimage/exposure/clahe.h
new file mode 100644
index 00000000..8b115aa1
--- /dev/null
+++ b/skimage/exposure/clahe.h
@@ -0,0 +1,329 @@
+/*
+ANSI C code from the article
+ * "Contrast Limited Adaptive Histogram Equalization"
+ * by Karel Zuiderveld, karel@cv.ruu.nl
+ * in "Graphics Gems IV", Academic Press, 1994
+=============
+
+ACM Software License Agreement
+
+All software, both binary and source published by the Association for Computing Machinery (hereafter, Software) is copyrighted by the Association (hereafter, ACM) and ownership of all right, title and interest in and to the Software remains with ACM. By using or copying the Software, User agrees to abide by the terms of this Agreement.
+Noncommercial Use
+
+The ACM grants to you (hereafter, User) a royalty-free, nonexclusive right to execute, copy, modify and distribute both the binary and source code solely for academic, research and other similar noncommercial uses, subject to the following conditions:
+
+    User acknowledges that the Software is still in the development stage and that it is being supplied "as is," without any support services from ACM. Neither ACM nor the author makes any representations or warranties, express or implied, including, without limitation, any representations or warranties of the merchantability or fitness for any particular purpose, or that the application of the software, will not infringe on any patents or other proprietary rights of others.
+    ACM shall not be held liable for direct, indirect, incidental or consequential damages arising from any claim by User or any third party with respect to uses allowed under this Agreement, or from any use of the Software.
+    User agrees to fully indemnify and hold harmless ACM and/or the author(s) of the original work from and against any and all claims, demands, suits, losses, damages, costs and expenses arising out of the User's use of the Software, including, without limitation, arising out of the User's modification of the Software.
+    User may modify the Software and distribute that modified work to third parties provided that: (a) if posted separately, it clearly acknowledges that it contains material copyrighted by ACM (b) no charge is associated with such copies, (c) User agrees to notify ACM and the Author(s) of the distribution, and (d) User clearly notifies secondary users that such modified work is not the original Software.
+    User agrees that ACM, the authors of the original work and others may enjoy a royalty-free, non-exclusive license to use, copy, modify and redistribute these modifications to the Software made by the User and distributed to third parties as a derivative work under this agreement.
+    This agreement will terminate immediately upon User's breach of, or non-compliance with, any of its terms. User may be held liable for any copyright infringement or the infringement of any other proprietary rights in the Software that is caused or facilitated by the User's failure to abide by the terms of this agreement.
+    This agreement will be construed and enforced in accordance with the law of the state of New York applicable to contracts performed entirely within the State. The parties irrevocably consent to the exclusive jurisdiction of the state or federal courts located in the City of New York for all disputes concerning this agreement. 
+
+Commercial Use
+
+Any User wishing to make a commercial use of the Software must contact ACM at permissions@acm.org to arrange an appropriate license. Commercial use includes (1) integrating or incorporating all or part of the source code into a product for sale or license by, or on behalf of, User to third parties, or (2) distribution of the binary or source code to third parties for use with a commercial product sold or licensed by, or on behalf of, User.
+
+Revised 6/98
+*/
+
+#ifdef BYTE_IMAGE
+typedef unsigned char kz_pixel_t;	 /* for 8 bit-per-pixel images */
+#define uiNR_OF_GREY (256)
+#else
+typedef unsigned short kz_pixel_t;	 /* for 12 bit-per-pixel images (default) */
+# define uiNR_OF_GREY (4096)
+#endif
+
+/*
+ * ANSI C code from the article
+ * "Contrast Limited Adaptive Histogram Equalization"
+ * by Karel Zuiderveld, karel@cv.ruu.nl
+ * in "Graphics Gems IV", Academic Press, 1994
+ *
+ *
+ *  These functions implement Contrast Limited Adaptive Histogram Equalization.
+ *  The main routine (CLAHE) expects an input image that is stored contiguously in
+ *  memory;  the CLAHE output image overwrites the original input image and has the
+ *  same minimum and maximum values (which must be provided by the user).
+ *  This implementation assumes that the X- and Y image resolutions are an integer
+ *  multiple of the X- and Y sizes of the contextual regions. A check on various other
+ *  error conditions is performed.
+ *
+ *  #define the symbol BYTE_IMAGE to make this implementation suitable for
+ *  8-bit images. The maximum number of contextual regions can be redefined
+ *  by changing uiMAX_REG_X and/or uiMAX_REG_Y; the use of more than 256
+ *  contextual regions is not recommended.
+ *
+ *  The code is ANSI-C and is also C++ compliant.
+ *
+ *  Author: Karel Zuiderveld, Computer Vision Research Group,
+ *	     Utrecht, The Netherlands (karel@cv.ruu.nl)
+ */
+
+/*********************** Local prototypes ************************/
+static void ClipHistogram (unsigned long*, unsigned int, unsigned long);
+static void MakeHistogram (kz_pixel_t*, unsigned int, unsigned int, unsigned int,
+		unsigned long*, unsigned int, kz_pixel_t*);
+static void MapHistogram (unsigned long*, kz_pixel_t, kz_pixel_t,
+	       unsigned int, unsigned long);
+static void MakeLut (kz_pixel_t*, kz_pixel_t, kz_pixel_t, unsigned int);
+static void Interpolate (kz_pixel_t*, int, unsigned long*, unsigned long*,
+	unsigned long*, unsigned long*, unsigned int, unsigned int, kz_pixel_t*);
+
+/**************	 Start of actual code **************/
+#include <stdlib.h>			 /* To get prototypes of malloc() and free() */
+
+
+const unsigned int uiMAX_REG_X = 16;	  /* max. # contextual regions in x-direction */
+const unsigned int uiMAX_REG_Y = 16;	  /* max. # contextual regions in y-direction */
+
+
+
+/************************** main function CLAHE ******************/
+int CLAHE (kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes,
+	 kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,
+	      unsigned int uiNrBins, float fCliplimit)
+/*   pImage - Pointer to the input/output image
+ *   uiXRes - Image resolution in the X direction
+ *   uiYRes - Image resolution in the Y direction
+ *   Min - Minimum greyvalue of input image (also becomes minimum of output image)
+ *   Max - Maximum greyvalue of input image (also becomes maximum of output image)
+ *   uiNrX - Number of contextial regions in the X direction (min 2, max uiMAX_REG_X)
+ *   uiNrY - Number of contextial regions in the Y direction (min 2, max uiMAX_REG_Y)
+ *   uiNrBins - Number of greybins for histogram ("dynamic range")
+ *   float fCliplimit - Normalized cliplimit (higher values give more contrast)
+ * The number of "effective" greylevels in the output image is set by uiNrBins; 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.
+ */
+{
+    unsigned int uiX, uiY;		  /* counters */
+    unsigned int uiXSize, uiYSize, uiSubX, uiSubY; /* size of context. reg. and subimages */
+    unsigned int uiXL, uiXR, uiYU, uiYB;  /* auxiliary variables interpolation routine */
+    unsigned long ulClipLimit, ulNrPixels;/* clip limit and region pixel count */
+    kz_pixel_t* pImPointer;		   /* pointer to image */
+    kz_pixel_t aLUT[uiNR_OF_GREY];	    /* lookup table used for scaling of input image */
+    unsigned long* pulHist, *pulMapArray; /* pointer to histogram and mappings*/
+    unsigned long* pulLU, *pulLB, *pulRU, *pulRB; /* auxiliary pointers interpolation */
+    
+    if (uiNrX > uiMAX_REG_X) return -1;	   /* # of regions x-direction too large */
+    if (uiNrY > uiMAX_REG_Y) return -2;	   /* # of regions y-direction too large */
+    if (uiXRes % uiNrX) return -3;	  /* x-resolution no multiple of uiNrX */
+    if (uiYRes % uiNrY) return -4;	  /* y-resolution no multiple of uiNrY */
+    if (Max >= uiNR_OF_GREY) return -5;	   /* maximum too large */
+    if (Min >= Max) return -6;		  /* minimum equal or larger than maximum */
+    if (uiNrX < 2 || uiNrY < 2) return -7;/* at least 4 contextual regions required */
+    if (fCliplimit == 1.0) return 0;	  /* is OK, immediately returns original image. */
+    if (uiNrBins == 0) uiNrBins = 128;	  /* default value when not specified */
+
+    pulMapArray=(unsigned long *)malloc(sizeof(unsigned long)*uiNrX*uiNrY*uiNrBins);
+    if (pulMapArray == 0) return -8;	  /* Not enough memory! (try reducing uiNrBins) */
+
+    uiXSize = uiXRes/uiNrX; uiYSize = uiYRes/uiNrY;  /* Actual size of contextual regions */
+    ulNrPixels = (unsigned long)uiXSize * (unsigned long)uiYSize;
+
+    if(fCliplimit > 0.0) {		  /* Calculate actual cliplimit	 */
+       ulClipLimit = (unsigned long) (fCliplimit * (uiXSize * uiYSize) / uiNrBins);
+       ulClipLimit = (ulClipLimit < 1UL) ? 1UL : ulClipLimit;
+    }
+    else ulClipLimit = 1UL<<14;		  /* Large value, do not clip (AHE) */
+    MakeLut(aLUT, Min, Max, uiNrBins);	  /* Make lookup table for mapping of greyvalues */
+    /* Calculate greylevel mappings for each contextual region */
+    for (uiY = 0, pImPointer = pImage; uiY < uiNrY; uiY++) {
+	for (uiX = 0; uiX < uiNrX; uiX++, pImPointer += uiXSize) {
+	    pulHist = &pulMapArray[uiNrBins * (uiY * uiNrX + uiX)];
+	    MakeHistogram(pImPointer,uiXRes,uiXSize,uiYSize,pulHist,uiNrBins,aLUT);
+	    ClipHistogram(pulHist, uiNrBins, ulClipLimit);
+	    MapHistogram(pulHist, Min, Max, uiNrBins, ulNrPixels);
+	}
+	pImPointer += (uiYSize - 1) * uiXRes;		  /* skip lines, set pointer */
+    }
+    /* Interpolate greylevel mappings to get CLAHE image */
+    for (pImPointer = pImage, uiY = 0; uiY <= uiNrY; uiY++) {
+	if (uiY == 0) {					  /* special case: top row */
+	    uiSubY = uiYSize >> 1;  uiYU = 0; uiYB = 0;
+	}
+	else {
+	    if (uiY == uiNrY) {				  /* special case: bottom row */
+		uiSubY = uiYSize >> 1;	uiYU = uiNrY-1;	 uiYB = uiYU;
+	    }
+	    else {					  /* default values */
+		uiSubY = uiYSize; uiYU = uiY - 1; uiYB = uiYU + 1;
+	    }
+	}
+	for (uiX = 0; uiX <= uiNrX; uiX++) {
+	    if (uiX == 0) {				  /* special case: left column */
+		uiSubX = uiXSize >> 1; uiXL = 0; uiXR = 0;
+	    }
+	    else {
+		if (uiX == uiNrX) {			  /* special case: right column */
+		    uiSubX = uiXSize >> 1;  uiXL = uiNrX - 1; uiXR = uiXL;
+		}
+		else {					  /* default values */
+		    uiSubX = uiXSize; uiXL = uiX - 1; uiXR = uiXL + 1;
+		}
+	    }
+
+	    pulLU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXL)];
+	    pulRU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXR)];
+	    pulLB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXL)];
+	    pulRB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXR)];
+	    Interpolate(pImPointer,uiXRes,pulLU,pulRU,pulLB,pulRB,uiSubX,uiSubY,aLUT);
+	    pImPointer += uiSubX;			  /* set pointer on next matrix */
+	}
+	pImPointer += (uiSubY - 1) * uiXRes;
+    }
+    free(pulMapArray);					  /* free space for histograms */
+    return 0;						  /* return status OK */
+}
+void ClipHistogram (unsigned long* pulHistogram, unsigned int
+		     uiNrGreylevels, unsigned long ulClipLimit)
+/* This function performs 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).
+ */
+{
+    unsigned long* pulBinPointer, *pulEndPointer, *pulHisto;
+    unsigned long ulNrExcess, ulUpper, ulBinIncr, ulStepSize, i;
+    long lBinExcess;
+
+    ulNrExcess = 0;  pulBinPointer = pulHistogram;
+    for (i = 0; i < uiNrGreylevels; i++) { /* calculate total number of excess pixels */
+	lBinExcess = (long) pulBinPointer[i] - (long) ulClipLimit;
+	if (lBinExcess > 0) ulNrExcess += lBinExcess;	  /* excess in current bin */
+    };
+
+    /* Second part: clip histogram and redistribute excess pixels in each bin */
+    ulBinIncr = ulNrExcess / uiNrGreylevels;		  /* average binincrement */
+    ulUpper =  ulClipLimit - ulBinIncr;	 /* Bins larger than ulUpper set to cliplimit */
+
+    for (i = 0; i < uiNrGreylevels; i++) {
+      if (pulHistogram[i] > ulClipLimit) pulHistogram[i] = ulClipLimit; /* clip bin */
+      else {
+	  if (pulHistogram[i] > ulUpper) {		/* high bin count */
+	      ulNrExcess -= pulHistogram[i] - ulUpper; pulHistogram[i]=ulClipLimit;
+	  }
+	  else {					/* low bin count */
+	      ulNrExcess -= ulBinIncr; pulHistogram[i] += ulBinIncr;
+	  }
+       }
+    }
+
+    while (ulNrExcess) {   /* Redistribute remaining excess  */
+	pulEndPointer = &pulHistogram[uiNrGreylevels]; pulHisto = pulHistogram;
+
+	while (ulNrExcess && pulHisto < pulEndPointer) {
+	    ulStepSize = uiNrGreylevels / ulNrExcess;
+	    if (ulStepSize < 1) ulStepSize = 1;		  /* stepsize at least 1 */
+	    for (pulBinPointer=pulHisto; pulBinPointer < pulEndPointer && ulNrExcess;
+		 pulBinPointer += ulStepSize) {
+		if (*pulBinPointer < ulClipLimit) {
+		    (*pulBinPointer)++;	 ulNrExcess--;	  /* reduce excess */
+		}
+	    }
+	    pulHisto++;		  /* restart redistributing on other bin location */
+	}
+    }
+}
+void MakeHistogram (kz_pixel_t* pImage, unsigned int uiXRes,
+		unsigned int uiSizeX, unsigned int uiSizeY,
+		unsigned long* pulHistogram,
+		unsigned int uiNrGreylevels, kz_pixel_t* pLookupTable)
+/* This function classifies the greylevels present in the array image into
+ * a greylevel histogram. The pLookupTable specifies the relationship
+ * between the greyvalue of the pixel (typically between 0 and 4095) and
+ * the corresponding bin in the histogram (usually containing only 128 bins).
+ */
+{
+    kz_pixel_t* pImagePointer;
+    unsigned int i;
+
+    for (i = 0; i < uiNrGreylevels; i++) pulHistogram[i] = 0L; /* clear histogram */
+
+    for (i = 0; i < uiSizeY; i++) {
+	pImagePointer = &pImage[uiSizeX];
+	while (pImage < pImagePointer) pulHistogram[pLookupTable[*pImage++]]++;
+	pImagePointer += uiXRes;
+	pImage = &pImagePointer[-uiSizeX];
+    }
+}
+
+void MapHistogram (unsigned long* pulHistogram, kz_pixel_t Min, kz_pixel_t Max,
+	       unsigned int uiNrGreylevels, unsigned long ulNrOfPixels)
+/* This function calculates the equalized lookup table (mapping) by
+ * cumulating the input histogram. Note: lookup table is rescaled in range [Min..Max].
+ */
+{
+    unsigned int i;  unsigned long ulSum = 0;
+    const float fScale = ((float)(Max - Min)) / ulNrOfPixels;
+    const unsigned long ulMin = (unsigned long) Min;
+
+    for (i = 0; i < uiNrGreylevels; i++) {
+	ulSum += pulHistogram[i]; pulHistogram[i]=(unsigned long)(ulMin+ulSum*fScale);
+	if (pulHistogram[i] > Max) pulHistogram[i] = Max;
+    }
+}
+
+void MakeLut (kz_pixel_t * pLUT, kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrBins)
+/* To speed up histogram clipping, the input image [Min,Max] is scaled down to
+ * [0,uiNrBins-1]. This function calculates the LUT.
+ */
+{
+    int i;
+    const kz_pixel_t BinSize = (kz_pixel_t) (1 + (Max - Min) / uiNrBins);
+
+    for (i = Min; i <= Max; i++)  pLUT[i] = (i - Min) / BinSize;
+}
+
+void Interpolate (kz_pixel_t * pImage, int uiXRes, unsigned long * pulMapLU,
+     unsigned long * pulMapRU, unsigned long * pulMapLB,  unsigned long * pulMapRB,
+     unsigned int uiXSize, unsigned int uiYSize, kz_pixel_t * pLUT)
+/* pImage      - pointer to input/output image
+ * uiXRes      - resolution of image in x-direction
+ * pulMap*     - mappings of greylevels from histograms
+ * uiXSize     - uiXSize of image submatrix
+ * uiYSize     - uiYSize of image submatrix
+ * pLUT	       - lookup table containing mapping greyvalues to bins
+ * This function calculates the new greylevel assignments of pixels within a submatrix
+ * of the image with size uiXSize and uiYSize. This is done by a bilinear interpolation
+ * between four different mappings in order to eliminate boundary artifacts.
+ * It uses a division; since division is often an expensive operation, I added code to
+ * perform a logical shift instead when feasible.
+ */
+{
+    const unsigned int uiIncr = uiXRes-uiXSize; /* Pointer increment after processing row */
+    kz_pixel_t GreyValue; unsigned int uiNum = uiXSize*uiYSize; /* Normalization factor */
+
+    unsigned int uiXCoef, uiYCoef, uiXInvCoef, uiYInvCoef, uiShift = 0;
+
+    if (uiNum & (uiNum - 1))   /* If uiNum is not a power of two, use division */
+    for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
+	 uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {
+	for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
+	     uiXCoef++, uiXInvCoef--) {
+	    GreyValue = pLUT[*pImage];		   /* get histogram bin value */
+	    *pImage++ = (kz_pixel_t ) ((uiYInvCoef * (uiXInvCoef*pulMapLU[GreyValue]
+				      + uiXCoef * pulMapRU[GreyValue])
+				+ uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
+				      + uiXCoef * pulMapRB[GreyValue])) / uiNum);
+	}
+    }
+    else {			   /* avoid the division and use a right shift instead */
+	while (uiNum >>= 1) uiShift++;		   /* Calculate 2log of uiNum */
+	for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
+	     uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {
+	     for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
+	       uiXCoef++, uiXInvCoef--) {
+	       GreyValue = pLUT[*pImage];	  /* get histogram bin value */
+	       *pImage++ = (kz_pixel_t)((uiYInvCoef* (uiXInvCoef * pulMapLU[GreyValue]
+				      + uiXCoef * pulMapRU[GreyValue])
+				+ uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
+				      + uiXCoef * pulMapRB[GreyValue])) >> uiShift);
+	    }
+	}
+    }
+}
+
diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index 5a722308..71869915 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -2,10 +2,13 @@ import numpy as np
 
 import skimage
 from skimage.util.dtype import dtype_range
+from skimage.color import rgb2gray
+from skimage.util.dtype import convert
 
+from _adapthist import _adapthist
 
 __all__ = ['histogram', 'cumulative_distribution', 'equalize',
-           'rescale_intensity']
+           'rescale_intensity', 'adapthist']
 
 
 def histogram(image, nbins=256):
@@ -183,3 +186,56 @@ def rescale_intensity(image, in_range=None, out_range=None):
 
     image = (image - imin) / float(imax - imin)
     return dtype(image * (omax - omin) + omin)
+
+
+def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
+    '''Contrast Limited Adaptive Histogram Equalization
+
+    Parameters
+    ----------
+    image : array-like
+        original image
+    nx : int, optional
+        Tile regions in the X direction (2, 16)
+    ny : int, optional
+        Tile regions in the Y direction (2, 16)
+    clip_limit : float: optional
+        Normalized cliplimit (higher values give more contrast)
+    nbins : int, optional
+        Greybins for histogram ("dynamic range")
+    out_range : str, optional
+        Range of the output image data.
+           - 'original' - Use original image limits
+           - 'full' - Use full range of image data type
+
+    Returns
+    -------
+    out - np.ndarray :
+        equalized image - may be a different shape than the original
+    '''
+    in_type = image.dtype.type
+    if out_range == 'full':
+        out_range = None
+    else:
+        out_range = (image.min(), image.max())
+    # must be converted to 12 bit for CLAHE
+    image = skimage.img_as_uint(image)
+    MAX_VAL = 2 ** 12 - 1
+    image = rescale_intensity(image, out_range=(0, MAX_VAL))
+    # handle color images - CLAHE accepts scalar images only
+    args = [image.copy(), 0, MAX_VAL, nx, ny, nbins, clip_limit]
+    if image.ndim == 3:
+        image = image[:, :, :3]
+        for channel in range(3):
+            args[0] = image[:, :, channel]
+            out = _adapthist(*args)
+            image[:out.shape[0], :out.shape[1], channel] = out
+    else:
+        out = _adapthist(*args)
+        image[:out.shape[0], :out.shape[1]] = out
+    # restore to desired output type and output limits
+    image = rescale_intensity(image)
+    if in_type != np.uint16:
+        image = convert(image, in_type)
+    image = rescale_intensity(image, out_range=out_range)
+    return image
diff --git a/skimage/exposure/setup.py b/skimage/exposure/setup.py
new file mode 100644
index 00000000..43bab23f
--- /dev/null
+++ b/skimage/exposure/setup.py
@@ -0,0 +1,31 @@
+#!/usr/bin/env python
+
+import os
+from skimage._build import cython
+
+base_path = os.path.abspath(os.path.dirname(__file__))
+
+
+def configuration(parent_package='', top_path=None):
+    from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
+
+    config = Configuration('exposure', parent_package, top_path)
+    config.add_data_dir('tests')
+
+    cython(['_adapthist.pyx'], working_path=base_path)
+
+    config.add_extension('_adapthist', sources=['_adapthist.c'],
+                         include_dirs=[get_numpy_include_dirs()])
+
+    return config
+
+if __name__ == '__main__':
+    from numpy.distutils.core import setup
+    setup(maintainer='scikits-image developers',
+          author='scikits-image developers',
+          maintainer_email='scikits-image@googlegroups.com',
+          description='Exposure',
+          url='https://github.com/scikits-image/scikits-image',
+          license='SciPy License (BSD Style)',
+          **(configuration(top_path='').todict())
+          )
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index b6ed0b12..fbad207b 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -4,6 +4,7 @@ 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
 
 
 # Test histogram equalization
@@ -71,6 +72,76 @@ def test_rescale_out_range():
     assert_close(out, [0, 63, 127])
 
 
+# Test rescale intensity
+# ======================
+
+def test_adapthist_ubyte():
+    '''Test a scalar uint8 image
+    '''
+    img = skimage.img_as_ubyte(data.moon())
+    adapted = exposure.adapthist(img, clip_limit=0.02)
+    assert adapted.min() == 0
+    assert adapted.max() == 255
+    assert img.shape == adapted.shape
+    assert peak_snr(img, adapted) > 22
+    assert norm_brightness_err(img, adapted) < 0.05
+    return img, adapted
+
+
+def test_adapthist_float():
+    '''Test an RGB float image
+    '''
+    img = skimage.img_as_float(data.lena())
+    adapted = exposure.adapthist(img, nx=10, ny=9, clip_limit=0.01,
+                        nbins=128, out_range='original')
+    assert_almost_equal(adapted.min() , img.min())
+    assert_almost_equal(adapted.min(), img.min())
+    assert img.shape == adapted.shape
+    assert peak_snr(img, adapted) > 136
+    assert norm_brightness_err(img, adapted) < 0.02
+    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())
+    mse = 1. / img1.size * np.square(img1 - img2).sum()
+    _, max_ = dtype_range[img1.dtype.type]
+    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
+        Normalize 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()
diff --git a/skimage/setup.py b/skimage/setup.py
index 0afc5bc9..80fdc0cb 100644
--- a/skimage/setup.py
+++ b/skimage/setup.py
@@ -17,6 +17,7 @@ def configuration(parent_package='', top_path=None):
     config.add_subpackage('morphology')
     config.add_subpackage('transform')
     config.add_subpackage('util')
+    config.add_subpackage('exposure')
 
     def add_test_directories(arg, dirname, fnames):
         if dirname.split(os.path.sep)[-1] == 'tests':

From 819332baa3f624637b52a26c89788b0bf4bb5f74 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sun, 5 Aug 2012 23:27:24 -0500
Subject: [PATCH 02/31] Tests pass and example works on build.

---
 doc/examples/plot_equalize.py           | 4 ++--
 skimage/exposure/__init__.py            | 2 +-
 skimage/exposure/tests/test_exposure.py | 2 ++
 3 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py
index 2cb9cbea..8b06890a 100644
--- a/doc/examples/plot_equalize.py
+++ b/doc/examples/plot_equalize.py
@@ -64,7 +64,7 @@ img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))
 img_eq = exposure.equalize(img)
 
 # Adaptive Equalization
-img_adapteq = exposure.adapthist(img)
+img_adapteq = exposure.adapthist(img, clip_limit=0.03)
 
 
 # Display results
@@ -81,7 +81,7 @@ 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_eq, axes[:, 3])
+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')
 
diff --git a/skimage/exposure/__init__.py b/skimage/exposure/__init__.py
index 7e19d317..a5e6c25a 100644
--- a/skimage/exposure/__init__.py
+++ b/skimage/exposure/__init__.py
@@ -1,2 +1,2 @@
-from .exposure import histogram, equalize, cumulative_distribution
+from .exposure import histogram, equalize, cumulative_distribution, adapthist
 from .exposure import rescale_intensity
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index fbad207b..963758f3 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -5,6 +5,7 @@ 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
@@ -94,6 +95,7 @@ def test_adapthist_float():
     img = skimage.img_as_float(data.lena())
     adapted = exposure.adapthist(img, nx=10, ny=9, clip_limit=0.01,
                         nbins=128, out_range='original')
+    assert_almost_equal = np.testing.assert_almost_equal
     assert_almost_equal(adapted.min() , img.min())
     assert_almost_equal(adapted.min(), img.min())
     assert img.shape == adapted.shape

From 045f6279b441045abf7aeec5c32cda68fb2fec92 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sun, 5 Aug 2012 23:36:39 -0500
Subject: [PATCH 03/31] Added Notes to the docstring.

---
 skimage/exposure/exposure.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index 71869915..df98da1f 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -212,6 +212,14 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
     -------
     out - np.ndarray :
         equalized image - may be a different shape than the original
+        
+    Notes
+    -----
+    * The underlying 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 RGB or RGBA images, the algorithm is run on each channel.
+    * For RGBA images, the original alpha channel is removed.
     '''
     in_type = image.dtype.type
     if out_range == 'full':

From 7dc17334a2431fc617f1467cfa47b274c3344159 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sun, 5 Aug 2012 23:42:44 -0500
Subject: [PATCH 04/31] Pep8 fixes.

---
 doc/examples/plot_equalize.py           | 2 +-
 skimage/exposure/exposure.py            | 2 +-
 skimage/exposure/tests/test_exposure.py | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py
index 8b06890a..e09b5a07 100644
--- a/doc/examples/plot_equalize.py
+++ b/doc/examples/plot_equalize.py
@@ -26,6 +26,7 @@ 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.
 
@@ -88,4 +89,3 @@ ax_cdf.set_ylabel('Fraction of total intensity')
 # prevent overlap of y-axis labels
 plt.subplots_adjust(wspace=0.4)
 plt.show()
-
diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index df98da1f..fbc2384e 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -212,7 +212,7 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
     -------
     out - np.ndarray :
         equalized image - may be a different shape than the original
-        
+
     Notes
     -----
     * The underlying algorithm relies on an image whose rows and columns are even multiples of
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 963758f3..9d70969a 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -96,7 +96,7 @@ def test_adapthist_float():
     adapted = exposure.adapthist(img, nx=10, ny=9, clip_limit=0.01,
                         nbins=128, out_range='original')
     assert_almost_equal = np.testing.assert_almost_equal
-    assert_almost_equal(adapted.min() , img.min())
+    assert_almost_equal(adapted.min(), img.min())
     assert_almost_equal(adapted.min(), img.min())
     assert img.shape == adapted.shape
     assert peak_snr(img, adapted) > 136

From 0e957dab3913f1095b11de58cdab081c94300e77 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Mon, 6 Aug 2012 19:06:57 -0500
Subject: [PATCH 05/31] Replaced generic ACM license with the Graphics Gems IV
 license.

---
 skimage/exposure/clahe.h | 31 ++++++++++++-------------------
 1 file changed, 12 insertions(+), 19 deletions(-)

diff --git a/skimage/exposure/clahe.h b/skimage/exposure/clahe.h
index 8b115aa1..a2721b98 100644
--- a/skimage/exposure/clahe.h
+++ b/skimage/exposure/clahe.h
@@ -5,26 +5,19 @@ ANSI C code from the article
  * in "Graphics Gems IV", Academic Press, 1994
 =============
 
-ACM Software License Agreement
+http://tog.acm.org/resources/GraphicsGems/
 
-All software, both binary and source published by the Association for Computing Machinery (hereafter, Software) is copyrighted by the Association (hereafter, ACM) and ownership of all right, title and interest in and to the Software remains with ACM. By using or copying the Software, User agrees to abide by the terms of this Agreement.
-Noncommercial Use
+EULA: 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.
 
-The ACM grants to you (hereafter, User) a royalty-free, nonexclusive right to execute, copy, modify and distribute both the binary and source code solely for academic, research and other similar noncommercial uses, subject to the following conditions:
-
-    User acknowledges that the Software is still in the development stage and that it is being supplied "as is," without any support services from ACM. Neither ACM nor the author makes any representations or warranties, express or implied, including, without limitation, any representations or warranties of the merchantability or fitness for any particular purpose, or that the application of the software, will not infringe on any patents or other proprietary rights of others.
-    ACM shall not be held liable for direct, indirect, incidental or consequential damages arising from any claim by User or any third party with respect to uses allowed under this Agreement, or from any use of the Software.
-    User agrees to fully indemnify and hold harmless ACM and/or the author(s) of the original work from and against any and all claims, demands, suits, losses, damages, costs and expenses arising out of the User's use of the Software, including, without limitation, arising out of the User's modification of the Software.
-    User may modify the Software and distribute that modified work to third parties provided that: (a) if posted separately, it clearly acknowledges that it contains material copyrighted by ACM (b) no charge is associated with such copies, (c) User agrees to notify ACM and the Author(s) of the distribution, and (d) User clearly notifies secondary users that such modified work is not the original Software.
-    User agrees that ACM, the authors of the original work and others may enjoy a royalty-free, non-exclusive license to use, copy, modify and redistribute these modifications to the Software made by the User and distributed to third parties as a derivative work under this agreement.
-    This agreement will terminate immediately upon User's breach of, or non-compliance with, any of its terms. User may be held liable for any copyright infringement or the infringement of any other proprietary rights in the Software that is caused or facilitated by the User's failure to abide by the terms of this agreement.
-    This agreement will be construed and enforced in accordance with the law of the state of New York applicable to contracts performed entirely within the State. The parties irrevocably consent to the exclusive jurisdiction of the state or federal courts located in the City of New York for all disputes concerning this agreement. 
-
-Commercial Use
-
-Any User wishing to make a commercial use of the Software must contact ACM at permissions@acm.org to arrange an appropriate license. Commercial use includes (1) integrating or incorporating all or part of the source code into a product for sale or license by, or on behalf of, User to third parties, or (2) distribution of the binary or source code to third parties for use with a commercial product sold or licensed by, or on behalf of, User.
-
-Revised 6/98
 */
 
 #ifdef BYTE_IMAGE
@@ -107,7 +100,7 @@ int CLAHE (kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes,
     kz_pixel_t aLUT[uiNR_OF_GREY];	    /* lookup table used for scaling of input image */
     unsigned long* pulHist, *pulMapArray; /* pointer to histogram and mappings*/
     unsigned long* pulLU, *pulLB, *pulRU, *pulRB; /* auxiliary pointers interpolation */
-    
+
     if (uiNrX > uiMAX_REG_X) return -1;	   /* # of regions x-direction too large */
     if (uiNrY > uiMAX_REG_Y) return -2;	   /* # of regions y-direction too large */
     if (uiXRes % uiNrX) return -3;	  /* x-resolution no multiple of uiNrX */

From 55a3f22ba293bec479d70079f4b28a147659d568 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sat, 11 Aug 2012 15:48:18 -0500
Subject: [PATCH 06/31] Improved support for color images using lab transform.

---
 skimage/exposure/exposure.py            | 60 ++++++++++++++++++++-----
 skimage/exposure/tests/test_exposure.py | 38 ++++++++++++----
 2 files changed, 77 insertions(+), 21 deletions(-)

diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index fbc2384e..47c98584 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -2,7 +2,7 @@ import numpy as np
 
 import skimage
 from skimage.util.dtype import dtype_range
-from skimage.color import rgb2gray
+import skimage.color as color
 from skimage.util.dtype import convert
 
 from _adapthist import _adapthist
@@ -210,7 +210,7 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
 
     Returns
     -------
-    out - np.ndarray :
+    out : np.ndarray
         equalized image - may be a different shape than the original
 
     Notes
@@ -218,8 +218,17 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
     * The underlying 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 RGB or RGBA images, the algorithm is run on each channel.
+    * For grayscale images, CLAHE is performed on one channel, and a grayscale is returned
+    * 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/
+    .. [2] https://en.wikipedia.org/wiki/CLAHE#CLAHE
     '''
     in_type = image.dtype.type
     if out_range == 'full':
@@ -227,23 +236,50 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
     else:
         out_range = (image.min(), image.max())
     # must be converted to 12 bit for CLAHE
-    image = skimage.img_as_uint(image)
+    int_image = skimage.img_as_uint(image)
     MAX_VAL = 2 ** 12 - 1
-    image = rescale_intensity(image, out_range=(0, MAX_VAL))
+    int_image = rescale_intensity(int_image, out_range=(0, MAX_VAL))
     # handle color images - CLAHE accepts scalar images only
-    args = [image.copy(), 0, MAX_VAL, nx, ny, nbins, clip_limit]
+    args = [int_image.copy(), 0, MAX_VAL, nx, ny, nbins, clip_limit]
     if image.ndim == 3:
-        image = image[:, :, :3]
-        for channel in range(3):
-            args[0] = image[:, :, channel]
+        # check for grayscale
+        if (np.allclose(image[:, :, 0], image[:, :, 1]) and
+            np.allclose(image[:, :, 2], image[:, :, 3])):
+            args[0] = image[:, :, 0]
             out = _adapthist(*args)
-            image[:out.shape[0], :out.shape[1], channel] = out
+            image = int_image[:, :, :3]
+            for channel in range(3):
+                image[:out.shape[0], :out.shape[1], channel] = out
+        # for color images, convert to LAB space for processing
+        else:
+            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, MAX_VAL))
+            new_L = _adapthist(*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:
         out = _adapthist(*args)
+        image = int_image
         image[:out.shape[0], :out.shape[1]] = out
     # restore to desired output type and output limits
     image = rescale_intensity(image)
-    if in_type != np.uint16:
-        image = convert(image, in_type)
+    image = convert(image, in_type)
     image = rescale_intensity(image, out_range=out_range)
     return image
+
+if __name__ == '__main__':
+    from skimage import data
+    import matplotlib.pyplot as plt
+    img = skimage.img_as_uint(data.lena())
+    adapted = adapthist(img, nx=10, ny=9, clip_limit=0.01,
+                            nbins=128, out_range='original')
+    plt.imshow(img)
+    plt.figure(); plt.imshow(skimage.img_as_ubyte(adapted))
+    plt.figure(); plt.imshow(color.lab2rgb(color.rgb2lab(img)))
+    plt.show()
+    print 'Done'
\ No newline at end of file
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 9d70969a..8ce7c8c6 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -76,7 +76,7 @@ def test_rescale_out_range():
 # Test rescale intensity
 # ======================
 
-def test_adapthist_ubyte():
+def test_adapthist_scalar():
     '''Test a scalar uint8 image
     '''
     img = skimage.img_as_ubyte(data.moon())
@@ -84,23 +84,43 @@ def test_adapthist_ubyte():
     assert adapted.min() == 0
     assert adapted.max() == 255
     assert img.shape == adapted.shape
-    assert peak_snr(img, adapted) > 22
-    assert norm_brightness_err(img, adapted) < 0.05
+    full_scale = skimage.exposure.rescale_intensity(img)
+    assert_almost_equal = np.testing.assert_almost_equal
+    assert_almost_equal(peak_snr(full_scale, adapted), 22.19920073)
+    assert_almost_equal(norm_brightness_err(full_scale, adapted),
+                        0.04161278)
     return img, adapted
 
 
-def test_adapthist_float():
-    '''Test an RGB float image
+def test_adapthist_grayscale():
+    '''Test a grayscale float image
     '''
     img = skimage.img_as_float(data.lena())
+    img = rgb2gray(img)
     adapted = exposure.adapthist(img, nx=10, ny=9, clip_limit=0.01,
                         nbins=128, out_range='original')
     assert_almost_equal = np.testing.assert_almost_equal
     assert_almost_equal(adapted.min(), img.min())
-    assert_almost_equal(adapted.min(), img.min())
+    assert_almost_equal(adapted.max(), img.max())
     assert img.shape == adapted.shape
-    assert peak_snr(img, adapted) > 136
-    assert norm_brightness_err(img, adapted) < 0.02
+    assert_almost_equal(peak_snr(img, adapted), 131.4962063)
+    assert_almost_equal(norm_brightness_err(img, adapted), 0.0208805)
+    return data, adapted
+
+
+def test_adapthist_color():
+    '''Test a color uint16 image
+    '''
+    img = skimage.img_as_uint(data.lena())
+    adapted = exposure.adapthist(img, clip_limit=0.01)
+    assert_almost_equal = np.testing.assert_almost_equal
+    assert adapted.min() == 0
+    assert adapted.max() == 65535
+    assert img.shape == adapted.shape
+    full_scale = skimage.exposure.rescale_intensity(img)
+    assert_almost_equal(peak_snr(full_scale, adapted), 64.29546231)
+    assert_almost_equal(norm_brightness_err(full_scale, adapted),
+                        0.181473754)
     return data, adapted
 
 
@@ -135,7 +155,7 @@ def norm_brightness_err(img1, img2):
     Returns
     -------
     norm_brightness_error : float
-        Normalize absolute mean brightness error
+        Normalized absolute mean brightness error
     '''
     if img1.ndim == 3:
         img1, img2 = rgb2gray(img1), rgb2gray(img2)

From d9926d8de4fb148dc76a38648757f5149df75316 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sat, 11 Aug 2012 21:28:08 -0500
Subject: [PATCH 07/31] Passed all tests with full coverage, PEP8 and Pylint
 performed.

---
 skimage/exposure/exposure.py            | 44 ++++++++++---------------
 skimage/exposure/tests/test_exposure.py |  1 +
 2 files changed, 18 insertions(+), 27 deletions(-)

diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index 47c98584..32231b00 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -188,16 +188,17 @@ def rescale_intensity(image, in_range=None, out_range=None):
     return dtype(image * (omax - omin) + omin)
 
 
-def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
+def adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, nbins=256,
+               out_range='full'):
     '''Contrast Limited Adaptive Histogram Equalization
 
     Parameters
     ----------
     image : array-like
         original image
-    nx : int, optional
+    ntiles_x : int, optional
         Tile regions in the X direction (2, 16)
-    ny : int, optional
+    ntiles_y : int, optional
         Tile regions in the Y direction (2, 16)
     clip_limit : float: optional
         Normalized cliplimit (higher values give more contrast)
@@ -237,15 +238,16 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
         out_range = (image.min(), image.max())
     # must be converted to 12 bit for CLAHE
     int_image = skimage.img_as_uint(image)
-    MAX_VAL = 2 ** 12 - 1
-    int_image = rescale_intensity(int_image, out_range=(0, MAX_VAL))
+    max_val = 2 ** 12 - 1
+    int_image = rescale_intensity(int_image, out_range=(0, max_val))
     # handle color images - CLAHE accepts scalar images only
-    args = [int_image.copy(), 0, MAX_VAL, nx, ny, nbins, clip_limit]
+    args = [int_image.copy(), 0, max_val, ntiles_x, ntiles_y, nbins,
+            clip_limit]
     if image.ndim == 3:
         # check for grayscale
         if (np.allclose(image[:, :, 0], image[:, :, 1]) and
-            np.allclose(image[:, :, 2], image[:, :, 3])):
-            args[0] = image[:, :, 0]
+            np.allclose(image[:, :, 1], image[:, :, 2])):
+            args[0] = int_image[:, :, 0]
             out = _adapthist(*args)
             image = int_image[:, :, :3]
             for channel in range(3):
@@ -253,13 +255,13 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
         # for color images, convert to LAB space for processing
         else:
             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, MAX_VAL))
-            new_L = _adapthist(*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
+            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, max_val))
+            new_l = _adapthist(*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:
@@ -271,15 +273,3 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
     image = convert(image, in_type)
     image = rescale_intensity(image, out_range=out_range)
     return image
-
-if __name__ == '__main__':
-    from skimage import data
-    import matplotlib.pyplot as plt
-    img = skimage.img_as_uint(data.lena())
-    adapted = adapthist(img, nx=10, ny=9, clip_limit=0.01,
-                            nbins=128, out_range='original')
-    plt.imshow(img)
-    plt.figure(); plt.imshow(skimage.img_as_ubyte(adapted))
-    plt.figure(); plt.imshow(color.lab2rgb(color.rgb2lab(img)))
-    plt.show()
-    print 'Done'
\ No newline at end of file
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 8ce7c8c6..2e64e741 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -97,6 +97,7 @@ def test_adapthist_grayscale():
     '''
     img = skimage.img_as_float(data.lena())
     img = rgb2gray(img)
+    img = np.dstack((img, img, img))
     adapted = exposure.adapthist(img, nx=10, ny=9, clip_limit=0.01,
                         nbins=128, out_range='original')
     assert_almost_equal = np.testing.assert_almost_equal

From c778de04ef5b56827cc105222e1681adc4eddc57 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sun, 21 Oct 2012 13:25:18 -0500
Subject: [PATCH 08/31] Translated clahe C code to python, added to list of
 contributors.

---
 CONTRIBUTORS.txt                        |   3 +
 skimage/exposure/__init__.py            |   3 +-
 skimage/exposure/_adapthist.py          | 337 ++++++++++++++++++++++++
 skimage/exposure/_adapthist.pyx         |  40 ---
 skimage/exposure/clahe.h                | 322 ----------------------
 skimage/exposure/exposure.py            |  92 +------
 skimage/exposure/setup.py               |  31 ---
 skimage/exposure/tests/test_exposure.py |  28 +-
 8 files changed, 357 insertions(+), 499 deletions(-)
 create mode 100644 skimage/exposure/_adapthist.py
 delete mode 100644 skimage/exposure/_adapthist.pyx
 delete mode 100644 skimage/exposure/clahe.h
 delete mode 100644 skimage/exposure/setup.py

diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index baaf064b..9b4dadae 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -107,3 +107,6 @@
   Polygon, circle and ellipse drawing functions
   Adaptive thresholding
   Implementation of Matlab's `regionprops`
+
+- Steven Silvester, Karel Zuiderveld
+  Adaptive Histogram Equalization
\ No newline at end of file
diff --git a/skimage/exposure/__init__.py b/skimage/exposure/__init__.py
index a5e6c25a..f8a4f046 100644
--- a/skimage/exposure/__init__.py
+++ b/skimage/exposure/__init__.py
@@ -1,2 +1,3 @@
-from .exposure import histogram, equalize, cumulative_distribution, adapthist
+from .exposure import histogram, equalize, cumulative_distribution
 from .exposure import rescale_intensity
+from ._adaphist import adapthist
diff --git a/skimage/exposure/_adapthist.py b/skimage/exposure/_adapthist.py
new file mode 100644
index 00000000..f35a5da5
--- /dev/null
+++ b/skimage/exposure/_adapthist.py
@@ -0,0 +1,337 @@
+'''
+Adapted code from the article
+ * "Contrast Limited Adaptive Histogram Equalization"
+ * by Karel Zuiderveld, karel@cv.ruu.nl
+ * in "Graphics Gems IV", Academic Press, 1994
+=============
+
+http://tog.acm.org/resources/GraphicsGems/
+
+EULA: 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.
+ *
+ *  Author: Karel Zuiderveld, Computer Vision Research Group,
+ *	     Utrecht, The Netherlands (karel@cv.ruu.nl)
+'''
+import numpy as np
+import skimage
+from skimage import color
+from skimage.exposure import rescale_intensity
+
+
+MAX_REG_X = 16  # max. # contextual regions in x-direction */
+MAX_REG_Y = 16  # max. # contextual regions in y-direction */
+NR_OF_GREY = 1 << 14  # number of grayscale levels to use in CLAHE algorithm
+
+
+def adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, nbins=256):
+    '''Contrast Limited Adaptive Histogram Equalization
+
+    Parameters
+    ----------
+    image : array-like
+        original image
+    ntiles_x : int, optional
+        Tile regions in the X direction (2, 16)
+    ntiles_y : int, optional
+        Tile regions in the Y direction (2, 16)
+    clip_limit : float: optional
+        Normalized cliplimit (higher values give more contrast)
+    nbins : int, optional
+        Greybins for histogram ("dynamic range")
+
+    Returns
+    -------
+    out : np.ndarray
+        equalized image - grayscale images are uint16, color images are float
+
+    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 grayscale images, CLAHE is performed on one channel,
+      and a grayscale is returned
+    * 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/
+    .. [2] https://en.wikipedia.org/wiki/CLAHE#CLAHE
+    '''
+    # convert to uint image
+    int_image = skimage.img_as_uint(image)
+    int_image = rescale_intensity(int_image, out_range=(0, NR_OF_GREY - 1))
+    # handle color images - CLAHE accepts scalar images only
+    args = [int_image.copy(), ntiles_x, ntiles_y, clip_limit * nbins, nbins]
+    if image.ndim == 3:
+        # check for grayscale
+        if (np.allclose(image[:, :, 0], image[:, :, 1]) and
+            np.allclose(image[:, :, 1], image[:, :, 2])):
+            args[0] = int_image[:, :, 0]
+            out = _clahe(*args)
+            image = int_image[:, :, :3]
+            for channel in range(3):
+                image[:out.shape[0], :out.shape[1], channel] = out
+        # for color images, convert to LAB space for processing
+        else:
+            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:
+        out = _clahe(*args)
+        image = int_image
+        image[:out.shape[0], :out.shape[1]] = out
+    return rescale_intensity(image)
+
+
+def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
+    '''Contrast Limited Adaptive Histogram Equalization
+
+    Parameters
+    ----------
+    image : array-like
+        original image
+    ntiles_x : int, optional
+        Tile regions in the X direction (2, 16)
+    ntiles_y : int, optional
+        Tile regions in the Y direction (2, 16)
+    clip_limit : float: optional
+        Normalized cliplimit (higher values give more contrast)
+    nbins : int, optional
+        Greybins for histogram ("dynamic range")
+
+    Returns
+    -------
+    out : np.ndarray
+
+    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_x * ntiles_y, 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
+    # Calculate greylevel mappings for each contextual region
+    ystart = 0
+    for y in range(ntiles_y):
+        xstart = 0
+        for x in range(ntiles_x):
+            sub_img = image[ystart: ystart + y_size,
+                            xstart: xstart + x_size]
+            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 * ntiles_x + x] = hist
+            xstart += x_size
+        ystart += y_size
+    # 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 * ntiles_x + xL]
+            mapRU = map_array[yU * ntiles_x + xR]
+            mapLB = map_array[yB * ntiles_x + xL]
+            mapRB = map_array[yB * ntiles_x + xR]
+            interpolate(image, xstart, xstep, ystart, ystep,
+                        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 : np.ndarray
+        histogram array
+    clip_limit : int
+        maximum allowed bin count
+
+    Returns
+    -------
+    hist : np.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):
+    '''Calculates the equalized lookup table (mapping)
+
+    It does so by cumulating the input histogram.
+
+    hist : np.ndarray
+        clipped histogram
+    min_val : int
+        min value for mapping
+    max_val : int
+        max value for mapping
+    n_pixels : int
+        number of pixels in the region
+
+    Returns
+    -------
+    out : np.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, xstart, xstep, ystart, ystep,
+                mapLU, mapRU, mapLB, mapRB, aLUT):
+    '''Find the new grayscale level for a region using bilinear interpolation
+
+    Parameters
+    ----------
+    image : np.ndarray
+        full image
+    xstart, xstop : int
+       indices of xslice
+    ystart, ystop : int
+        indices of yslice
+    map* : np.ndarray
+        mappings of greylevels from histograms
+    aLUT : np.ndarray
+        maps grayscale levels in image to histogram levels
+
+    Returns
+    -------
+    out : np.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 = xstep * ystep  # Normalization factor
+
+    # interpolation weight matrices
+    x_coef, y_coef = np.meshgrid(np.arange(xstep),
+                                 np.arange(ystep))
+    x_inv_coef, y_inv_coef = x_coef[:, ::-1] + 1, y_coef[::-1] + 1
+
+    im_slice = image[ystart: ystart + ystep, xstart: xstart + xstep]
+    im_slice = aLUT[im_slice]
+    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)
+    image[ystart: ystart + ystep, xstart: xstart + xstep] = new
diff --git a/skimage/exposure/_adapthist.pyx b/skimage/exposure/_adapthist.pyx
deleted file mode 100644
index 16e84336..00000000
--- a/skimage/exposure/_adapthist.pyx
+++ /dev/null
@@ -1,40 +0,0 @@
-import numpy as np
-cimport numpy as np
-
-cdef extern from "clahe.h":
-    ctypedef unsigned int kz_pixel_t
-    int CLAHE(kz_pixel_t * pImage, unsigned int uiXRes, unsigned int uiYRes,
-        kz_pixel_t Min,
-        kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,
-        unsigned int uiNrBins, float fCliplimit)
-
-
-def _adapthist(np.ndarray[kz_pixel_t, ndim=2] image, min, max, nx, ny, nbins,
-                 clip_limit=0):
-    '''
-    image -  to the input/output image
-    min - Minimum greyvalue of input image (also becomes min of output image)
-    max - Maximum greyvalue of input image (also becomes max of output image)
-    nx - Number of tile regions in the X direction (min 2, max uiMAX_REG_X)
-    ny - Number of tile regions in the Y direction (min 2, max uiMAX_REG_Y)
-    nbins - Number of greybins for histogram ("dynamic range")
-    clip_limit - Normalized cliplimit (higher values give more contrast)
-
-    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.
-    '''
-    # we need the image to be divisible by nx and ny
-    uiYRes = image.shape[0] - image.shape[0] % ny
-    uiXRes = image.shape[1] - image.shape[1] % nx
-
-    cdef np.ndarray[kz_pixel_t, ndim = 1] flattened
-    flattened = image[:uiYRes, :uiXRes].ravel()
-
-    ret = CLAHE(< kz_pixel_t * > flattened.data, uiXRes, uiYRes, min,
-                max, nx, ny, nbins, clip_limit * nbins)
-    if ret < 0:
-        print 'clahe error: ', ret
-    return flattened.reshape((uiYRes, uiXRes))
diff --git a/skimage/exposure/clahe.h b/skimage/exposure/clahe.h
deleted file mode 100644
index a2721b98..00000000
--- a/skimage/exposure/clahe.h
+++ /dev/null
@@ -1,322 +0,0 @@
-/*
-ANSI C code from the article
- * "Contrast Limited Adaptive Histogram Equalization"
- * by Karel Zuiderveld, karel@cv.ruu.nl
- * in "Graphics Gems IV", Academic Press, 1994
-=============
-
-http://tog.acm.org/resources/GraphicsGems/
-
-EULA: 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.
-
-*/
-
-#ifdef BYTE_IMAGE
-typedef unsigned char kz_pixel_t;	 /* for 8 bit-per-pixel images */
-#define uiNR_OF_GREY (256)
-#else
-typedef unsigned short kz_pixel_t;	 /* for 12 bit-per-pixel images (default) */
-# define uiNR_OF_GREY (4096)
-#endif
-
-/*
- * ANSI C code from the article
- * "Contrast Limited Adaptive Histogram Equalization"
- * by Karel Zuiderveld, karel@cv.ruu.nl
- * in "Graphics Gems IV", Academic Press, 1994
- *
- *
- *  These functions implement Contrast Limited Adaptive Histogram Equalization.
- *  The main routine (CLAHE) expects an input image that is stored contiguously in
- *  memory;  the CLAHE output image overwrites the original input image and has the
- *  same minimum and maximum values (which must be provided by the user).
- *  This implementation assumes that the X- and Y image resolutions are an integer
- *  multiple of the X- and Y sizes of the contextual regions. A check on various other
- *  error conditions is performed.
- *
- *  #define the symbol BYTE_IMAGE to make this implementation suitable for
- *  8-bit images. The maximum number of contextual regions can be redefined
- *  by changing uiMAX_REG_X and/or uiMAX_REG_Y; the use of more than 256
- *  contextual regions is not recommended.
- *
- *  The code is ANSI-C and is also C++ compliant.
- *
- *  Author: Karel Zuiderveld, Computer Vision Research Group,
- *	     Utrecht, The Netherlands (karel@cv.ruu.nl)
- */
-
-/*********************** Local prototypes ************************/
-static void ClipHistogram (unsigned long*, unsigned int, unsigned long);
-static void MakeHistogram (kz_pixel_t*, unsigned int, unsigned int, unsigned int,
-		unsigned long*, unsigned int, kz_pixel_t*);
-static void MapHistogram (unsigned long*, kz_pixel_t, kz_pixel_t,
-	       unsigned int, unsigned long);
-static void MakeLut (kz_pixel_t*, kz_pixel_t, kz_pixel_t, unsigned int);
-static void Interpolate (kz_pixel_t*, int, unsigned long*, unsigned long*,
-	unsigned long*, unsigned long*, unsigned int, unsigned int, kz_pixel_t*);
-
-/**************	 Start of actual code **************/
-#include <stdlib.h>			 /* To get prototypes of malloc() and free() */
-
-
-const unsigned int uiMAX_REG_X = 16;	  /* max. # contextual regions in x-direction */
-const unsigned int uiMAX_REG_Y = 16;	  /* max. # contextual regions in y-direction */
-
-
-
-/************************** main function CLAHE ******************/
-int CLAHE (kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes,
-	 kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,
-	      unsigned int uiNrBins, float fCliplimit)
-/*   pImage - Pointer to the input/output image
- *   uiXRes - Image resolution in the X direction
- *   uiYRes - Image resolution in the Y direction
- *   Min - Minimum greyvalue of input image (also becomes minimum of output image)
- *   Max - Maximum greyvalue of input image (also becomes maximum of output image)
- *   uiNrX - Number of contextial regions in the X direction (min 2, max uiMAX_REG_X)
- *   uiNrY - Number of contextial regions in the Y direction (min 2, max uiMAX_REG_Y)
- *   uiNrBins - Number of greybins for histogram ("dynamic range")
- *   float fCliplimit - Normalized cliplimit (higher values give more contrast)
- * The number of "effective" greylevels in the output image is set by uiNrBins; 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.
- */
-{
-    unsigned int uiX, uiY;		  /* counters */
-    unsigned int uiXSize, uiYSize, uiSubX, uiSubY; /* size of context. reg. and subimages */
-    unsigned int uiXL, uiXR, uiYU, uiYB;  /* auxiliary variables interpolation routine */
-    unsigned long ulClipLimit, ulNrPixels;/* clip limit and region pixel count */
-    kz_pixel_t* pImPointer;		   /* pointer to image */
-    kz_pixel_t aLUT[uiNR_OF_GREY];	    /* lookup table used for scaling of input image */
-    unsigned long* pulHist, *pulMapArray; /* pointer to histogram and mappings*/
-    unsigned long* pulLU, *pulLB, *pulRU, *pulRB; /* auxiliary pointers interpolation */
-
-    if (uiNrX > uiMAX_REG_X) return -1;	   /* # of regions x-direction too large */
-    if (uiNrY > uiMAX_REG_Y) return -2;	   /* # of regions y-direction too large */
-    if (uiXRes % uiNrX) return -3;	  /* x-resolution no multiple of uiNrX */
-    if (uiYRes % uiNrY) return -4;	  /* y-resolution no multiple of uiNrY */
-    if (Max >= uiNR_OF_GREY) return -5;	   /* maximum too large */
-    if (Min >= Max) return -6;		  /* minimum equal or larger than maximum */
-    if (uiNrX < 2 || uiNrY < 2) return -7;/* at least 4 contextual regions required */
-    if (fCliplimit == 1.0) return 0;	  /* is OK, immediately returns original image. */
-    if (uiNrBins == 0) uiNrBins = 128;	  /* default value when not specified */
-
-    pulMapArray=(unsigned long *)malloc(sizeof(unsigned long)*uiNrX*uiNrY*uiNrBins);
-    if (pulMapArray == 0) return -8;	  /* Not enough memory! (try reducing uiNrBins) */
-
-    uiXSize = uiXRes/uiNrX; uiYSize = uiYRes/uiNrY;  /* Actual size of contextual regions */
-    ulNrPixels = (unsigned long)uiXSize * (unsigned long)uiYSize;
-
-    if(fCliplimit > 0.0) {		  /* Calculate actual cliplimit	 */
-       ulClipLimit = (unsigned long) (fCliplimit * (uiXSize * uiYSize) / uiNrBins);
-       ulClipLimit = (ulClipLimit < 1UL) ? 1UL : ulClipLimit;
-    }
-    else ulClipLimit = 1UL<<14;		  /* Large value, do not clip (AHE) */
-    MakeLut(aLUT, Min, Max, uiNrBins);	  /* Make lookup table for mapping of greyvalues */
-    /* Calculate greylevel mappings for each contextual region */
-    for (uiY = 0, pImPointer = pImage; uiY < uiNrY; uiY++) {
-	for (uiX = 0; uiX < uiNrX; uiX++, pImPointer += uiXSize) {
-	    pulHist = &pulMapArray[uiNrBins * (uiY * uiNrX + uiX)];
-	    MakeHistogram(pImPointer,uiXRes,uiXSize,uiYSize,pulHist,uiNrBins,aLUT);
-	    ClipHistogram(pulHist, uiNrBins, ulClipLimit);
-	    MapHistogram(pulHist, Min, Max, uiNrBins, ulNrPixels);
-	}
-	pImPointer += (uiYSize - 1) * uiXRes;		  /* skip lines, set pointer */
-    }
-    /* Interpolate greylevel mappings to get CLAHE image */
-    for (pImPointer = pImage, uiY = 0; uiY <= uiNrY; uiY++) {
-	if (uiY == 0) {					  /* special case: top row */
-	    uiSubY = uiYSize >> 1;  uiYU = 0; uiYB = 0;
-	}
-	else {
-	    if (uiY == uiNrY) {				  /* special case: bottom row */
-		uiSubY = uiYSize >> 1;	uiYU = uiNrY-1;	 uiYB = uiYU;
-	    }
-	    else {					  /* default values */
-		uiSubY = uiYSize; uiYU = uiY - 1; uiYB = uiYU + 1;
-	    }
-	}
-	for (uiX = 0; uiX <= uiNrX; uiX++) {
-	    if (uiX == 0) {				  /* special case: left column */
-		uiSubX = uiXSize >> 1; uiXL = 0; uiXR = 0;
-	    }
-	    else {
-		if (uiX == uiNrX) {			  /* special case: right column */
-		    uiSubX = uiXSize >> 1;  uiXL = uiNrX - 1; uiXR = uiXL;
-		}
-		else {					  /* default values */
-		    uiSubX = uiXSize; uiXL = uiX - 1; uiXR = uiXL + 1;
-		}
-	    }
-
-	    pulLU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXL)];
-	    pulRU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXR)];
-	    pulLB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXL)];
-	    pulRB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXR)];
-	    Interpolate(pImPointer,uiXRes,pulLU,pulRU,pulLB,pulRB,uiSubX,uiSubY,aLUT);
-	    pImPointer += uiSubX;			  /* set pointer on next matrix */
-	}
-	pImPointer += (uiSubY - 1) * uiXRes;
-    }
-    free(pulMapArray);					  /* free space for histograms */
-    return 0;						  /* return status OK */
-}
-void ClipHistogram (unsigned long* pulHistogram, unsigned int
-		     uiNrGreylevels, unsigned long ulClipLimit)
-/* This function performs 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).
- */
-{
-    unsigned long* pulBinPointer, *pulEndPointer, *pulHisto;
-    unsigned long ulNrExcess, ulUpper, ulBinIncr, ulStepSize, i;
-    long lBinExcess;
-
-    ulNrExcess = 0;  pulBinPointer = pulHistogram;
-    for (i = 0; i < uiNrGreylevels; i++) { /* calculate total number of excess pixels */
-	lBinExcess = (long) pulBinPointer[i] - (long) ulClipLimit;
-	if (lBinExcess > 0) ulNrExcess += lBinExcess;	  /* excess in current bin */
-    };
-
-    /* Second part: clip histogram and redistribute excess pixels in each bin */
-    ulBinIncr = ulNrExcess / uiNrGreylevels;		  /* average binincrement */
-    ulUpper =  ulClipLimit - ulBinIncr;	 /* Bins larger than ulUpper set to cliplimit */
-
-    for (i = 0; i < uiNrGreylevels; i++) {
-      if (pulHistogram[i] > ulClipLimit) pulHistogram[i] = ulClipLimit; /* clip bin */
-      else {
-	  if (pulHistogram[i] > ulUpper) {		/* high bin count */
-	      ulNrExcess -= pulHistogram[i] - ulUpper; pulHistogram[i]=ulClipLimit;
-	  }
-	  else {					/* low bin count */
-	      ulNrExcess -= ulBinIncr; pulHistogram[i] += ulBinIncr;
-	  }
-       }
-    }
-
-    while (ulNrExcess) {   /* Redistribute remaining excess  */
-	pulEndPointer = &pulHistogram[uiNrGreylevels]; pulHisto = pulHistogram;
-
-	while (ulNrExcess && pulHisto < pulEndPointer) {
-	    ulStepSize = uiNrGreylevels / ulNrExcess;
-	    if (ulStepSize < 1) ulStepSize = 1;		  /* stepsize at least 1 */
-	    for (pulBinPointer=pulHisto; pulBinPointer < pulEndPointer && ulNrExcess;
-		 pulBinPointer += ulStepSize) {
-		if (*pulBinPointer < ulClipLimit) {
-		    (*pulBinPointer)++;	 ulNrExcess--;	  /* reduce excess */
-		}
-	    }
-	    pulHisto++;		  /* restart redistributing on other bin location */
-	}
-    }
-}
-void MakeHistogram (kz_pixel_t* pImage, unsigned int uiXRes,
-		unsigned int uiSizeX, unsigned int uiSizeY,
-		unsigned long* pulHistogram,
-		unsigned int uiNrGreylevels, kz_pixel_t* pLookupTable)
-/* This function classifies the greylevels present in the array image into
- * a greylevel histogram. The pLookupTable specifies the relationship
- * between the greyvalue of the pixel (typically between 0 and 4095) and
- * the corresponding bin in the histogram (usually containing only 128 bins).
- */
-{
-    kz_pixel_t* pImagePointer;
-    unsigned int i;
-
-    for (i = 0; i < uiNrGreylevels; i++) pulHistogram[i] = 0L; /* clear histogram */
-
-    for (i = 0; i < uiSizeY; i++) {
-	pImagePointer = &pImage[uiSizeX];
-	while (pImage < pImagePointer) pulHistogram[pLookupTable[*pImage++]]++;
-	pImagePointer += uiXRes;
-	pImage = &pImagePointer[-uiSizeX];
-    }
-}
-
-void MapHistogram (unsigned long* pulHistogram, kz_pixel_t Min, kz_pixel_t Max,
-	       unsigned int uiNrGreylevels, unsigned long ulNrOfPixels)
-/* This function calculates the equalized lookup table (mapping) by
- * cumulating the input histogram. Note: lookup table is rescaled in range [Min..Max].
- */
-{
-    unsigned int i;  unsigned long ulSum = 0;
-    const float fScale = ((float)(Max - Min)) / ulNrOfPixels;
-    const unsigned long ulMin = (unsigned long) Min;
-
-    for (i = 0; i < uiNrGreylevels; i++) {
-	ulSum += pulHistogram[i]; pulHistogram[i]=(unsigned long)(ulMin+ulSum*fScale);
-	if (pulHistogram[i] > Max) pulHistogram[i] = Max;
-    }
-}
-
-void MakeLut (kz_pixel_t * pLUT, kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrBins)
-/* To speed up histogram clipping, the input image [Min,Max] is scaled down to
- * [0,uiNrBins-1]. This function calculates the LUT.
- */
-{
-    int i;
-    const kz_pixel_t BinSize = (kz_pixel_t) (1 + (Max - Min) / uiNrBins);
-
-    for (i = Min; i <= Max; i++)  pLUT[i] = (i - Min) / BinSize;
-}
-
-void Interpolate (kz_pixel_t * pImage, int uiXRes, unsigned long * pulMapLU,
-     unsigned long * pulMapRU, unsigned long * pulMapLB,  unsigned long * pulMapRB,
-     unsigned int uiXSize, unsigned int uiYSize, kz_pixel_t * pLUT)
-/* pImage      - pointer to input/output image
- * uiXRes      - resolution of image in x-direction
- * pulMap*     - mappings of greylevels from histograms
- * uiXSize     - uiXSize of image submatrix
- * uiYSize     - uiYSize of image submatrix
- * pLUT	       - lookup table containing mapping greyvalues to bins
- * This function calculates the new greylevel assignments of pixels within a submatrix
- * of the image with size uiXSize and uiYSize. This is done by a bilinear interpolation
- * between four different mappings in order to eliminate boundary artifacts.
- * It uses a division; since division is often an expensive operation, I added code to
- * perform a logical shift instead when feasible.
- */
-{
-    const unsigned int uiIncr = uiXRes-uiXSize; /* Pointer increment after processing row */
-    kz_pixel_t GreyValue; unsigned int uiNum = uiXSize*uiYSize; /* Normalization factor */
-
-    unsigned int uiXCoef, uiYCoef, uiXInvCoef, uiYInvCoef, uiShift = 0;
-
-    if (uiNum & (uiNum - 1))   /* If uiNum is not a power of two, use division */
-    for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
-	 uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {
-	for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
-	     uiXCoef++, uiXInvCoef--) {
-	    GreyValue = pLUT[*pImage];		   /* get histogram bin value */
-	    *pImage++ = (kz_pixel_t ) ((uiYInvCoef * (uiXInvCoef*pulMapLU[GreyValue]
-				      + uiXCoef * pulMapRU[GreyValue])
-				+ uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
-				      + uiXCoef * pulMapRB[GreyValue])) / uiNum);
-	}
-    }
-    else {			   /* avoid the division and use a right shift instead */
-	while (uiNum >>= 1) uiShift++;		   /* Calculate 2log of uiNum */
-	for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
-	     uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {
-	     for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
-	       uiXCoef++, uiXInvCoef--) {
-	       GreyValue = pLUT[*pImage];	  /* get histogram bin value */
-	       *pImage++ = (kz_pixel_t)((uiYInvCoef* (uiXInvCoef * pulMapLU[GreyValue]
-				      + uiXCoef * pulMapRU[GreyValue])
-				+ uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
-				      + uiXCoef * pulMapRB[GreyValue])) >> uiShift);
-	    }
-	}
-    }
-}
-
diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index 32231b00..8899f319 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -5,10 +5,9 @@ from skimage.util.dtype import dtype_range
 import skimage.color as color
 from skimage.util.dtype import convert
 
-from _adapthist import _adapthist
 
 __all__ = ['histogram', 'cumulative_distribution', 'equalize',
-           'rescale_intensity', 'adapthist']
+           'rescale_intensity']
 
 
 def histogram(image, nbins=256):
@@ -185,91 +184,4 @@ def rescale_intensity(image, in_range=None, out_range=None):
     image = np.clip(image, imin, imax)
 
     image = (image - imin) / float(imax - imin)
-    return dtype(image * (omax - omin) + omin)
-
-
-def adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, nbins=256,
-               out_range='full'):
-    '''Contrast Limited Adaptive Histogram Equalization
-
-    Parameters
-    ----------
-    image : array-like
-        original image
-    ntiles_x : int, optional
-        Tile regions in the X direction (2, 16)
-    ntiles_y : int, optional
-        Tile regions in the Y direction (2, 16)
-    clip_limit : float: optional
-        Normalized cliplimit (higher values give more contrast)
-    nbins : int, optional
-        Greybins for histogram ("dynamic range")
-    out_range : str, optional
-        Range of the output image data.
-           - 'original' - Use original image limits
-           - 'full' - Use full range of image data type
-
-    Returns
-    -------
-    out : np.ndarray
-        equalized image - may be a different shape than the original
-
-    Notes
-    -----
-    * The underlying 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 grayscale images, CLAHE is performed on one channel, and a grayscale is returned
-    * 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/
-    .. [2] https://en.wikipedia.org/wiki/CLAHE#CLAHE
-    '''
-    in_type = image.dtype.type
-    if out_range == 'full':
-        out_range = None
-    else:
-        out_range = (image.min(), image.max())
-    # must be converted to 12 bit for CLAHE
-    int_image = skimage.img_as_uint(image)
-    max_val = 2 ** 12 - 1
-    int_image = rescale_intensity(int_image, out_range=(0, max_val))
-    # handle color images - CLAHE accepts scalar images only
-    args = [int_image.copy(), 0, max_val, ntiles_x, ntiles_y, nbins,
-            clip_limit]
-    if image.ndim == 3:
-        # check for grayscale
-        if (np.allclose(image[:, :, 0], image[:, :, 1]) and
-            np.allclose(image[:, :, 1], image[:, :, 2])):
-            args[0] = int_image[:, :, 0]
-            out = _adapthist(*args)
-            image = int_image[:, :, :3]
-            for channel in range(3):
-                image[:out.shape[0], :out.shape[1], channel] = out
-        # for color images, convert to LAB space for processing
-        else:
-            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, max_val))
-            new_l = _adapthist(*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:
-        out = _adapthist(*args)
-        image = int_image
-        image[:out.shape[0], :out.shape[1]] = out
-    # restore to desired output type and output limits
-    image = rescale_intensity(image)
-    image = convert(image, in_type)
-    image = rescale_intensity(image, out_range=out_range)
-    return image
+    return dtype(image * (omax - omin) + omin)
\ No newline at end of file
diff --git a/skimage/exposure/setup.py b/skimage/exposure/setup.py
deleted file mode 100644
index 43bab23f..00000000
--- a/skimage/exposure/setup.py
+++ /dev/null
@@ -1,31 +0,0 @@
-#!/usr/bin/env python
-
-import os
-from skimage._build import cython
-
-base_path = os.path.abspath(os.path.dirname(__file__))
-
-
-def configuration(parent_package='', top_path=None):
-    from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
-
-    config = Configuration('exposure', parent_package, top_path)
-    config.add_data_dir('tests')
-
-    cython(['_adapthist.pyx'], working_path=base_path)
-
-    config.add_extension('_adapthist', sources=['_adapthist.c'],
-                         include_dirs=[get_numpy_include_dirs()])
-
-    return config
-
-if __name__ == '__main__':
-    from numpy.distutils.core import setup
-    setup(maintainer='scikits-image developers',
-          author='scikits-image developers',
-          maintainer_email='scikits-image@googlegroups.com',
-          description='Exposure',
-          url='https://github.com/scikits-image/scikits-image',
-          license='SciPy License (BSD Style)',
-          **(configuration(top_path='').todict())
-          )
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 2e64e741..3cf42a7f 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -1,6 +1,7 @@
 import numpy as np
 from numpy.testing import assert_array_almost_equal as assert_close
-
+import sys
+sys.path.insert(0, '../..')
 import skimage
 from skimage import data
 from skimage import exposure
@@ -82,30 +83,27 @@ def test_adapthist_scalar():
     img = skimage.img_as_ubyte(data.moon())
     adapted = exposure.adapthist(img, clip_limit=0.02)
     assert adapted.min() == 0
-    assert adapted.max() == 255
+    assert adapted.max() == (1<<16) - 1
     assert img.shape == adapted.shape
-    full_scale = skimage.exposure.rescale_intensity(img)
+    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), 22.19920073)
+    assert_almost_equal(peak_snr(full_scale, adapted), 28.6262034)
     assert_almost_equal(norm_brightness_err(full_scale, adapted),
-                        0.04161278)
+                        0.0410010)
     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.adapthist(img, nx=10, ny=9, clip_limit=0.01,
-                        nbins=128, out_range='original')
+    adapted = exposure.adapthist(img, 10, 9, clip_limit=0.01,
+                        nbins=128)
     assert_almost_equal = np.testing.assert_almost_equal
-    assert_almost_equal(adapted.min(), img.min())
-    assert_almost_equal(adapted.max(), img.max())
     assert img.shape == adapted.shape
-    assert_almost_equal(peak_snr(img, adapted), 131.4962063)
-    assert_almost_equal(norm_brightness_err(img, adapted), 0.0208805)
+    assert_almost_equal(peak_snr(img, adapted), 106.3020173)
+    assert_almost_equal(norm_brightness_err(img, adapted), 0.0218686)
     return data, adapted
 
 
@@ -116,12 +114,12 @@ def test_adapthist_color():
     adapted = exposure.adapthist(img, clip_limit=0.01)
     assert_almost_equal = np.testing.assert_almost_equal
     assert adapted.min() == 0
-    assert adapted.max() == 65535
+    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.29546231)
+    assert_almost_equal(peak_snr(full_scale, adapted), 64.7167515)
     assert_almost_equal(norm_brightness_err(full_scale, adapted),
-                        0.181473754)
+                        0.17922429)
     return data, adapted
 
 

From 79fb543256ac236ca59d93134935b41d52f5774e Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sun, 21 Oct 2012 13:49:35 -0500
Subject: [PATCH 09/31] Removed sys.path hack in test, fixed spelling of import
 in __init__.py

---
 skimage/exposure/__init__.py            | 2 +-
 skimage/exposure/tests/test_exposure.py | 2 --
 2 files changed, 1 insertion(+), 3 deletions(-)

diff --git a/skimage/exposure/__init__.py b/skimage/exposure/__init__.py
index f8a4f046..d495686a 100644
--- a/skimage/exposure/__init__.py
+++ b/skimage/exposure/__init__.py
@@ -1,3 +1,3 @@
 from .exposure import histogram, equalize, cumulative_distribution
 from .exposure import rescale_intensity
-from ._adaphist import adapthist
+from ._adapthist import adapthist
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 3cf42a7f..04c0cf0e 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -1,7 +1,5 @@
 import numpy as np
 from numpy.testing import assert_array_almost_equal as assert_close
-import sys
-sys.path.insert(0, '../..')
 import skimage
 from skimage import data
 from skimage import exposure

From 2e1729a9fbbc21fc0b04df8e68efbab9cfd6dada Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sat, 1 Dec 2012 23:46:52 -0600
Subject: [PATCH 10/31] Cleaned up adapthist func to use LAB for all
 RGB->required rest result change

---
 skimage/exposure/_adapthist.py          | 43 +++++++++----------------
 skimage/exposure/tests/test_exposure.py |  4 +--
 2 files changed, 18 insertions(+), 29 deletions(-)

diff --git a/skimage/exposure/_adapthist.py b/skimage/exposure/_adapthist.py
index f35a5da5..091514b6 100644
--- a/skimage/exposure/_adapthist.py
+++ b/skimage/exposure/_adapthist.py
@@ -71,37 +71,26 @@ def adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, nbins=256):
     .. [1] http://tog.acm.org/resources/GraphicsGems/
     .. [2] https://en.wikipedia.org/wiki/CLAHE#CLAHE
     '''
-    # convert to uint image
-    int_image = skimage.img_as_uint(image)
-    int_image = rescale_intensity(int_image, out_range=(0, NR_OF_GREY - 1))
     # handle color images - CLAHE accepts scalar images only
-    args = [int_image.copy(), ntiles_x, ntiles_y, clip_limit * nbins, nbins]
-    if image.ndim == 3:
-        # check for grayscale
-        if (np.allclose(image[:, :, 0], image[:, :, 1]) and
-            np.allclose(image[:, :, 1], image[:, :, 2])):
-            args[0] = int_image[:, :, 0]
-            out = _clahe(*args)
-            image = int_image[:, :, :3]
-            for channel in range(3):
-                image[:out.shape[0], :out.shape[1], channel] = out
-        # for color images, convert to LAB space for processing
-        else:
-            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))
+    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 = int_image
         image[:out.shape[0], :out.shape[1]] = out
-    return rescale_intensity(image)
+        image = rescale_intensity(image)
+    return image
 
 
 def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 04c0cf0e..dc48a3d8 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -100,8 +100,8 @@ def test_adapthist_grayscale():
                         nbins=128)
     assert_almost_equal = np.testing.assert_almost_equal
     assert img.shape == adapted.shape
-    assert_almost_equal(peak_snr(img, adapted), 106.3020173)
-    assert_almost_equal(norm_brightness_err(img, adapted), 0.0218686)
+    assert_almost_equal(peak_snr(img, adapted), 77.5841845)
+    assert_almost_equal(norm_brightness_err(img, adapted), 0.0376632)
     return data, adapted
 
 

From c6d0acfdabd5f14d2da353e76571638f3cc1fb0e Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sun, 5 Aug 2012 23:07:24 -0500
Subject: [PATCH 11/31] Initial commit of adapthist modifications.

---
 doc/examples/plot_equalize.py           |   8 +-
 skimage/exposure/_adapthist.pyx         |  40 +++
 skimage/exposure/clahe.h                | 329 ++++++++++++++++++++++++
 skimage/exposure/exposure.py            |  58 ++++-
 skimage/exposure/setup.py               |  31 +++
 skimage/exposure/tests/test_exposure.py |  71 +++++
 6 files changed, 535 insertions(+), 2 deletions(-)
 create mode 100644 skimage/exposure/_adapthist.pyx
 create mode 100644 skimage/exposure/clahe.h
 create mode 100644 skimage/exposure/setup.py

diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py
index 1e1291d8..2cb9cbea 100644
--- a/doc/examples/plot_equalize.py
+++ b/doc/examples/plot_equalize.py
@@ -63,9 +63,12 @@ img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))
 # Equalization
 img_eq = exposure.equalize(img)
 
+# Adaptive Equalization
+img_adapteq = exposure.adapthist(img)
+
 
 # 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')
@@ -78,6 +81,9 @@ 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_eq, axes[:, 3])
+ax_img.set_title('Adaptive equalization')
+ax_cdf.set_ylabel('Fraction of total intensity')
 
 # prevent overlap of y-axis labels
 plt.subplots_adjust(wspace=0.4)
diff --git a/skimage/exposure/_adapthist.pyx b/skimage/exposure/_adapthist.pyx
new file mode 100644
index 00000000..16e84336
--- /dev/null
+++ b/skimage/exposure/_adapthist.pyx
@@ -0,0 +1,40 @@
+import numpy as np
+cimport numpy as np
+
+cdef extern from "clahe.h":
+    ctypedef unsigned int kz_pixel_t
+    int CLAHE(kz_pixel_t * pImage, unsigned int uiXRes, unsigned int uiYRes,
+        kz_pixel_t Min,
+        kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,
+        unsigned int uiNrBins, float fCliplimit)
+
+
+def _adapthist(np.ndarray[kz_pixel_t, ndim=2] image, min, max, nx, ny, nbins,
+                 clip_limit=0):
+    '''
+    image -  to the input/output image
+    min - Minimum greyvalue of input image (also becomes min of output image)
+    max - Maximum greyvalue of input image (also becomes max of output image)
+    nx - Number of tile regions in the X direction (min 2, max uiMAX_REG_X)
+    ny - Number of tile regions in the Y direction (min 2, max uiMAX_REG_Y)
+    nbins - Number of greybins for histogram ("dynamic range")
+    clip_limit - Normalized cliplimit (higher values give more contrast)
+
+    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.
+    '''
+    # we need the image to be divisible by nx and ny
+    uiYRes = image.shape[0] - image.shape[0] % ny
+    uiXRes = image.shape[1] - image.shape[1] % nx
+
+    cdef np.ndarray[kz_pixel_t, ndim = 1] flattened
+    flattened = image[:uiYRes, :uiXRes].ravel()
+
+    ret = CLAHE(< kz_pixel_t * > flattened.data, uiXRes, uiYRes, min,
+                max, nx, ny, nbins, clip_limit * nbins)
+    if ret < 0:
+        print 'clahe error: ', ret
+    return flattened.reshape((uiYRes, uiXRes))
diff --git a/skimage/exposure/clahe.h b/skimage/exposure/clahe.h
new file mode 100644
index 00000000..8b115aa1
--- /dev/null
+++ b/skimage/exposure/clahe.h
@@ -0,0 +1,329 @@
+/*
+ANSI C code from the article
+ * "Contrast Limited Adaptive Histogram Equalization"
+ * by Karel Zuiderveld, karel@cv.ruu.nl
+ * in "Graphics Gems IV", Academic Press, 1994
+=============
+
+ACM Software License Agreement
+
+All software, both binary and source published by the Association for Computing Machinery (hereafter, Software) is copyrighted by the Association (hereafter, ACM) and ownership of all right, title and interest in and to the Software remains with ACM. By using or copying the Software, User agrees to abide by the terms of this Agreement.
+Noncommercial Use
+
+The ACM grants to you (hereafter, User) a royalty-free, nonexclusive right to execute, copy, modify and distribute both the binary and source code solely for academic, research and other similar noncommercial uses, subject to the following conditions:
+
+    User acknowledges that the Software is still in the development stage and that it is being supplied "as is," without any support services from ACM. Neither ACM nor the author makes any representations or warranties, express or implied, including, without limitation, any representations or warranties of the merchantability or fitness for any particular purpose, or that the application of the software, will not infringe on any patents or other proprietary rights of others.
+    ACM shall not be held liable for direct, indirect, incidental or consequential damages arising from any claim by User or any third party with respect to uses allowed under this Agreement, or from any use of the Software.
+    User agrees to fully indemnify and hold harmless ACM and/or the author(s) of the original work from and against any and all claims, demands, suits, losses, damages, costs and expenses arising out of the User's use of the Software, including, without limitation, arising out of the User's modification of the Software.
+    User may modify the Software and distribute that modified work to third parties provided that: (a) if posted separately, it clearly acknowledges that it contains material copyrighted by ACM (b) no charge is associated with such copies, (c) User agrees to notify ACM and the Author(s) of the distribution, and (d) User clearly notifies secondary users that such modified work is not the original Software.
+    User agrees that ACM, the authors of the original work and others may enjoy a royalty-free, non-exclusive license to use, copy, modify and redistribute these modifications to the Software made by the User and distributed to third parties as a derivative work under this agreement.
+    This agreement will terminate immediately upon User's breach of, or non-compliance with, any of its terms. User may be held liable for any copyright infringement or the infringement of any other proprietary rights in the Software that is caused or facilitated by the User's failure to abide by the terms of this agreement.
+    This agreement will be construed and enforced in accordance with the law of the state of New York applicable to contracts performed entirely within the State. The parties irrevocably consent to the exclusive jurisdiction of the state or federal courts located in the City of New York for all disputes concerning this agreement. 
+
+Commercial Use
+
+Any User wishing to make a commercial use of the Software must contact ACM at permissions@acm.org to arrange an appropriate license. Commercial use includes (1) integrating or incorporating all or part of the source code into a product for sale or license by, or on behalf of, User to third parties, or (2) distribution of the binary or source code to third parties for use with a commercial product sold or licensed by, or on behalf of, User.
+
+Revised 6/98
+*/
+
+#ifdef BYTE_IMAGE
+typedef unsigned char kz_pixel_t;	 /* for 8 bit-per-pixel images */
+#define uiNR_OF_GREY (256)
+#else
+typedef unsigned short kz_pixel_t;	 /* for 12 bit-per-pixel images (default) */
+# define uiNR_OF_GREY (4096)
+#endif
+
+/*
+ * ANSI C code from the article
+ * "Contrast Limited Adaptive Histogram Equalization"
+ * by Karel Zuiderveld, karel@cv.ruu.nl
+ * in "Graphics Gems IV", Academic Press, 1994
+ *
+ *
+ *  These functions implement Contrast Limited Adaptive Histogram Equalization.
+ *  The main routine (CLAHE) expects an input image that is stored contiguously in
+ *  memory;  the CLAHE output image overwrites the original input image and has the
+ *  same minimum and maximum values (which must be provided by the user).
+ *  This implementation assumes that the X- and Y image resolutions are an integer
+ *  multiple of the X- and Y sizes of the contextual regions. A check on various other
+ *  error conditions is performed.
+ *
+ *  #define the symbol BYTE_IMAGE to make this implementation suitable for
+ *  8-bit images. The maximum number of contextual regions can be redefined
+ *  by changing uiMAX_REG_X and/or uiMAX_REG_Y; the use of more than 256
+ *  contextual regions is not recommended.
+ *
+ *  The code is ANSI-C and is also C++ compliant.
+ *
+ *  Author: Karel Zuiderveld, Computer Vision Research Group,
+ *	     Utrecht, The Netherlands (karel@cv.ruu.nl)
+ */
+
+/*********************** Local prototypes ************************/
+static void ClipHistogram (unsigned long*, unsigned int, unsigned long);
+static void MakeHistogram (kz_pixel_t*, unsigned int, unsigned int, unsigned int,
+		unsigned long*, unsigned int, kz_pixel_t*);
+static void MapHistogram (unsigned long*, kz_pixel_t, kz_pixel_t,
+	       unsigned int, unsigned long);
+static void MakeLut (kz_pixel_t*, kz_pixel_t, kz_pixel_t, unsigned int);
+static void Interpolate (kz_pixel_t*, int, unsigned long*, unsigned long*,
+	unsigned long*, unsigned long*, unsigned int, unsigned int, kz_pixel_t*);
+
+/**************	 Start of actual code **************/
+#include <stdlib.h>			 /* To get prototypes of malloc() and free() */
+
+
+const unsigned int uiMAX_REG_X = 16;	  /* max. # contextual regions in x-direction */
+const unsigned int uiMAX_REG_Y = 16;	  /* max. # contextual regions in y-direction */
+
+
+
+/************************** main function CLAHE ******************/
+int CLAHE (kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes,
+	 kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,
+	      unsigned int uiNrBins, float fCliplimit)
+/*   pImage - Pointer to the input/output image
+ *   uiXRes - Image resolution in the X direction
+ *   uiYRes - Image resolution in the Y direction
+ *   Min - Minimum greyvalue of input image (also becomes minimum of output image)
+ *   Max - Maximum greyvalue of input image (also becomes maximum of output image)
+ *   uiNrX - Number of contextial regions in the X direction (min 2, max uiMAX_REG_X)
+ *   uiNrY - Number of contextial regions in the Y direction (min 2, max uiMAX_REG_Y)
+ *   uiNrBins - Number of greybins for histogram ("dynamic range")
+ *   float fCliplimit - Normalized cliplimit (higher values give more contrast)
+ * The number of "effective" greylevels in the output image is set by uiNrBins; 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.
+ */
+{
+    unsigned int uiX, uiY;		  /* counters */
+    unsigned int uiXSize, uiYSize, uiSubX, uiSubY; /* size of context. reg. and subimages */
+    unsigned int uiXL, uiXR, uiYU, uiYB;  /* auxiliary variables interpolation routine */
+    unsigned long ulClipLimit, ulNrPixels;/* clip limit and region pixel count */
+    kz_pixel_t* pImPointer;		   /* pointer to image */
+    kz_pixel_t aLUT[uiNR_OF_GREY];	    /* lookup table used for scaling of input image */
+    unsigned long* pulHist, *pulMapArray; /* pointer to histogram and mappings*/
+    unsigned long* pulLU, *pulLB, *pulRU, *pulRB; /* auxiliary pointers interpolation */
+    
+    if (uiNrX > uiMAX_REG_X) return -1;	   /* # of regions x-direction too large */
+    if (uiNrY > uiMAX_REG_Y) return -2;	   /* # of regions y-direction too large */
+    if (uiXRes % uiNrX) return -3;	  /* x-resolution no multiple of uiNrX */
+    if (uiYRes % uiNrY) return -4;	  /* y-resolution no multiple of uiNrY */
+    if (Max >= uiNR_OF_GREY) return -5;	   /* maximum too large */
+    if (Min >= Max) return -6;		  /* minimum equal or larger than maximum */
+    if (uiNrX < 2 || uiNrY < 2) return -7;/* at least 4 contextual regions required */
+    if (fCliplimit == 1.0) return 0;	  /* is OK, immediately returns original image. */
+    if (uiNrBins == 0) uiNrBins = 128;	  /* default value when not specified */
+
+    pulMapArray=(unsigned long *)malloc(sizeof(unsigned long)*uiNrX*uiNrY*uiNrBins);
+    if (pulMapArray == 0) return -8;	  /* Not enough memory! (try reducing uiNrBins) */
+
+    uiXSize = uiXRes/uiNrX; uiYSize = uiYRes/uiNrY;  /* Actual size of contextual regions */
+    ulNrPixels = (unsigned long)uiXSize * (unsigned long)uiYSize;
+
+    if(fCliplimit > 0.0) {		  /* Calculate actual cliplimit	 */
+       ulClipLimit = (unsigned long) (fCliplimit * (uiXSize * uiYSize) / uiNrBins);
+       ulClipLimit = (ulClipLimit < 1UL) ? 1UL : ulClipLimit;
+    }
+    else ulClipLimit = 1UL<<14;		  /* Large value, do not clip (AHE) */
+    MakeLut(aLUT, Min, Max, uiNrBins);	  /* Make lookup table for mapping of greyvalues */
+    /* Calculate greylevel mappings for each contextual region */
+    for (uiY = 0, pImPointer = pImage; uiY < uiNrY; uiY++) {
+	for (uiX = 0; uiX < uiNrX; uiX++, pImPointer += uiXSize) {
+	    pulHist = &pulMapArray[uiNrBins * (uiY * uiNrX + uiX)];
+	    MakeHistogram(pImPointer,uiXRes,uiXSize,uiYSize,pulHist,uiNrBins,aLUT);
+	    ClipHistogram(pulHist, uiNrBins, ulClipLimit);
+	    MapHistogram(pulHist, Min, Max, uiNrBins, ulNrPixels);
+	}
+	pImPointer += (uiYSize - 1) * uiXRes;		  /* skip lines, set pointer */
+    }
+    /* Interpolate greylevel mappings to get CLAHE image */
+    for (pImPointer = pImage, uiY = 0; uiY <= uiNrY; uiY++) {
+	if (uiY == 0) {					  /* special case: top row */
+	    uiSubY = uiYSize >> 1;  uiYU = 0; uiYB = 0;
+	}
+	else {
+	    if (uiY == uiNrY) {				  /* special case: bottom row */
+		uiSubY = uiYSize >> 1;	uiYU = uiNrY-1;	 uiYB = uiYU;
+	    }
+	    else {					  /* default values */
+		uiSubY = uiYSize; uiYU = uiY - 1; uiYB = uiYU + 1;
+	    }
+	}
+	for (uiX = 0; uiX <= uiNrX; uiX++) {
+	    if (uiX == 0) {				  /* special case: left column */
+		uiSubX = uiXSize >> 1; uiXL = 0; uiXR = 0;
+	    }
+	    else {
+		if (uiX == uiNrX) {			  /* special case: right column */
+		    uiSubX = uiXSize >> 1;  uiXL = uiNrX - 1; uiXR = uiXL;
+		}
+		else {					  /* default values */
+		    uiSubX = uiXSize; uiXL = uiX - 1; uiXR = uiXL + 1;
+		}
+	    }
+
+	    pulLU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXL)];
+	    pulRU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXR)];
+	    pulLB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXL)];
+	    pulRB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXR)];
+	    Interpolate(pImPointer,uiXRes,pulLU,pulRU,pulLB,pulRB,uiSubX,uiSubY,aLUT);
+	    pImPointer += uiSubX;			  /* set pointer on next matrix */
+	}
+	pImPointer += (uiSubY - 1) * uiXRes;
+    }
+    free(pulMapArray);					  /* free space for histograms */
+    return 0;						  /* return status OK */
+}
+void ClipHistogram (unsigned long* pulHistogram, unsigned int
+		     uiNrGreylevels, unsigned long ulClipLimit)
+/* This function performs 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).
+ */
+{
+    unsigned long* pulBinPointer, *pulEndPointer, *pulHisto;
+    unsigned long ulNrExcess, ulUpper, ulBinIncr, ulStepSize, i;
+    long lBinExcess;
+
+    ulNrExcess = 0;  pulBinPointer = pulHistogram;
+    for (i = 0; i < uiNrGreylevels; i++) { /* calculate total number of excess pixels */
+	lBinExcess = (long) pulBinPointer[i] - (long) ulClipLimit;
+	if (lBinExcess > 0) ulNrExcess += lBinExcess;	  /* excess in current bin */
+    };
+
+    /* Second part: clip histogram and redistribute excess pixels in each bin */
+    ulBinIncr = ulNrExcess / uiNrGreylevels;		  /* average binincrement */
+    ulUpper =  ulClipLimit - ulBinIncr;	 /* Bins larger than ulUpper set to cliplimit */
+
+    for (i = 0; i < uiNrGreylevels; i++) {
+      if (pulHistogram[i] > ulClipLimit) pulHistogram[i] = ulClipLimit; /* clip bin */
+      else {
+	  if (pulHistogram[i] > ulUpper) {		/* high bin count */
+	      ulNrExcess -= pulHistogram[i] - ulUpper; pulHistogram[i]=ulClipLimit;
+	  }
+	  else {					/* low bin count */
+	      ulNrExcess -= ulBinIncr; pulHistogram[i] += ulBinIncr;
+	  }
+       }
+    }
+
+    while (ulNrExcess) {   /* Redistribute remaining excess  */
+	pulEndPointer = &pulHistogram[uiNrGreylevels]; pulHisto = pulHistogram;
+
+	while (ulNrExcess && pulHisto < pulEndPointer) {
+	    ulStepSize = uiNrGreylevels / ulNrExcess;
+	    if (ulStepSize < 1) ulStepSize = 1;		  /* stepsize at least 1 */
+	    for (pulBinPointer=pulHisto; pulBinPointer < pulEndPointer && ulNrExcess;
+		 pulBinPointer += ulStepSize) {
+		if (*pulBinPointer < ulClipLimit) {
+		    (*pulBinPointer)++;	 ulNrExcess--;	  /* reduce excess */
+		}
+	    }
+	    pulHisto++;		  /* restart redistributing on other bin location */
+	}
+    }
+}
+void MakeHistogram (kz_pixel_t* pImage, unsigned int uiXRes,
+		unsigned int uiSizeX, unsigned int uiSizeY,
+		unsigned long* pulHistogram,
+		unsigned int uiNrGreylevels, kz_pixel_t* pLookupTable)
+/* This function classifies the greylevels present in the array image into
+ * a greylevel histogram. The pLookupTable specifies the relationship
+ * between the greyvalue of the pixel (typically between 0 and 4095) and
+ * the corresponding bin in the histogram (usually containing only 128 bins).
+ */
+{
+    kz_pixel_t* pImagePointer;
+    unsigned int i;
+
+    for (i = 0; i < uiNrGreylevels; i++) pulHistogram[i] = 0L; /* clear histogram */
+
+    for (i = 0; i < uiSizeY; i++) {
+	pImagePointer = &pImage[uiSizeX];
+	while (pImage < pImagePointer) pulHistogram[pLookupTable[*pImage++]]++;
+	pImagePointer += uiXRes;
+	pImage = &pImagePointer[-uiSizeX];
+    }
+}
+
+void MapHistogram (unsigned long* pulHistogram, kz_pixel_t Min, kz_pixel_t Max,
+	       unsigned int uiNrGreylevels, unsigned long ulNrOfPixels)
+/* This function calculates the equalized lookup table (mapping) by
+ * cumulating the input histogram. Note: lookup table is rescaled in range [Min..Max].
+ */
+{
+    unsigned int i;  unsigned long ulSum = 0;
+    const float fScale = ((float)(Max - Min)) / ulNrOfPixels;
+    const unsigned long ulMin = (unsigned long) Min;
+
+    for (i = 0; i < uiNrGreylevels; i++) {
+	ulSum += pulHistogram[i]; pulHistogram[i]=(unsigned long)(ulMin+ulSum*fScale);
+	if (pulHistogram[i] > Max) pulHistogram[i] = Max;
+    }
+}
+
+void MakeLut (kz_pixel_t * pLUT, kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrBins)
+/* To speed up histogram clipping, the input image [Min,Max] is scaled down to
+ * [0,uiNrBins-1]. This function calculates the LUT.
+ */
+{
+    int i;
+    const kz_pixel_t BinSize = (kz_pixel_t) (1 + (Max - Min) / uiNrBins);
+
+    for (i = Min; i <= Max; i++)  pLUT[i] = (i - Min) / BinSize;
+}
+
+void Interpolate (kz_pixel_t * pImage, int uiXRes, unsigned long * pulMapLU,
+     unsigned long * pulMapRU, unsigned long * pulMapLB,  unsigned long * pulMapRB,
+     unsigned int uiXSize, unsigned int uiYSize, kz_pixel_t * pLUT)
+/* pImage      - pointer to input/output image
+ * uiXRes      - resolution of image in x-direction
+ * pulMap*     - mappings of greylevels from histograms
+ * uiXSize     - uiXSize of image submatrix
+ * uiYSize     - uiYSize of image submatrix
+ * pLUT	       - lookup table containing mapping greyvalues to bins
+ * This function calculates the new greylevel assignments of pixels within a submatrix
+ * of the image with size uiXSize and uiYSize. This is done by a bilinear interpolation
+ * between four different mappings in order to eliminate boundary artifacts.
+ * It uses a division; since division is often an expensive operation, I added code to
+ * perform a logical shift instead when feasible.
+ */
+{
+    const unsigned int uiIncr = uiXRes-uiXSize; /* Pointer increment after processing row */
+    kz_pixel_t GreyValue; unsigned int uiNum = uiXSize*uiYSize; /* Normalization factor */
+
+    unsigned int uiXCoef, uiYCoef, uiXInvCoef, uiYInvCoef, uiShift = 0;
+
+    if (uiNum & (uiNum - 1))   /* If uiNum is not a power of two, use division */
+    for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
+	 uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {
+	for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
+	     uiXCoef++, uiXInvCoef--) {
+	    GreyValue = pLUT[*pImage];		   /* get histogram bin value */
+	    *pImage++ = (kz_pixel_t ) ((uiYInvCoef * (uiXInvCoef*pulMapLU[GreyValue]
+				      + uiXCoef * pulMapRU[GreyValue])
+				+ uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
+				      + uiXCoef * pulMapRB[GreyValue])) / uiNum);
+	}
+    }
+    else {			   /* avoid the division and use a right shift instead */
+	while (uiNum >>= 1) uiShift++;		   /* Calculate 2log of uiNum */
+	for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
+	     uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {
+	     for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
+	       uiXCoef++, uiXInvCoef--) {
+	       GreyValue = pLUT[*pImage];	  /* get histogram bin value */
+	       *pImage++ = (kz_pixel_t)((uiYInvCoef* (uiXInvCoef * pulMapLU[GreyValue]
+				      + uiXCoef * pulMapRU[GreyValue])
+				+ uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
+				      + uiXCoef * pulMapRB[GreyValue])) >> uiShift);
+	    }
+	}
+    }
+}
+
diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index bffe660a..1c4c69da 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -2,10 +2,13 @@ import numpy as np
 
 from skimage import img_as_float
 from skimage.util.dtype import dtype_range
+from skimage.color import rgb2gray
+from skimage.util.dtype import convert
 
+from _adapthist import _adapthist
 
 __all__ = ['histogram', 'cumulative_distribution', 'equalize',
-           'rescale_intensity']
+           'rescale_intensity', 'adapthist']
 
 
 def histogram(image, nbins=256):
@@ -189,3 +192,56 @@ def rescale_intensity(image, in_range=None, out_range=None):
 
     image = (image - imin) / float(imax - imin)
     return dtype(image * (omax - omin) + omin)
+
+
+def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
+    '''Contrast Limited Adaptive Histogram Equalization
+
+    Parameters
+    ----------
+    image : array-like
+        original image
+    nx : int, optional
+        Tile regions in the X direction (2, 16)
+    ny : int, optional
+        Tile regions in the Y direction (2, 16)
+    clip_limit : float: optional
+        Normalized cliplimit (higher values give more contrast)
+    nbins : int, optional
+        Greybins for histogram ("dynamic range")
+    out_range : str, optional
+        Range of the output image data.
+           - 'original' - Use original image limits
+           - 'full' - Use full range of image data type
+
+    Returns
+    -------
+    out - np.ndarray :
+        equalized image - may be a different shape than the original
+    '''
+    in_type = image.dtype.type
+    if out_range == 'full':
+        out_range = None
+    else:
+        out_range = (image.min(), image.max())
+    # must be converted to 12 bit for CLAHE
+    image = skimage.img_as_uint(image)
+    MAX_VAL = 2 ** 12 - 1
+    image = rescale_intensity(image, out_range=(0, MAX_VAL))
+    # handle color images - CLAHE accepts scalar images only
+    args = [image.copy(), 0, MAX_VAL, nx, ny, nbins, clip_limit]
+    if image.ndim == 3:
+        image = image[:, :, :3]
+        for channel in range(3):
+            args[0] = image[:, :, channel]
+            out = _adapthist(*args)
+            image[:out.shape[0], :out.shape[1], channel] = out
+    else:
+        out = _adapthist(*args)
+        image[:out.shape[0], :out.shape[1]] = out
+    # restore to desired output type and output limits
+    image = rescale_intensity(image)
+    if in_type != np.uint16:
+        image = convert(image, in_type)
+    image = rescale_intensity(image, out_range=out_range)
+    return image
diff --git a/skimage/exposure/setup.py b/skimage/exposure/setup.py
new file mode 100644
index 00000000..43bab23f
--- /dev/null
+++ b/skimage/exposure/setup.py
@@ -0,0 +1,31 @@
+#!/usr/bin/env python
+
+import os
+from skimage._build import cython
+
+base_path = os.path.abspath(os.path.dirname(__file__))
+
+
+def configuration(parent_package='', top_path=None):
+    from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
+
+    config = Configuration('exposure', parent_package, top_path)
+    config.add_data_dir('tests')
+
+    cython(['_adapthist.pyx'], working_path=base_path)
+
+    config.add_extension('_adapthist', sources=['_adapthist.c'],
+                         include_dirs=[get_numpy_include_dirs()])
+
+    return config
+
+if __name__ == '__main__':
+    from numpy.distutils.core import setup
+    setup(maintainer='scikits-image developers',
+          author='scikits-image developers',
+          maintainer_email='scikits-image@googlegroups.com',
+          description='Exposure',
+          url='https://github.com/scikits-image/scikits-image',
+          license='SciPy License (BSD Style)',
+          **(configuration(top_path='').todict())
+          )
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index b6ed0b12..fbad207b 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -4,6 +4,7 @@ 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
 
 
 # Test histogram equalization
@@ -71,6 +72,76 @@ def test_rescale_out_range():
     assert_close(out, [0, 63, 127])
 
 
+# Test rescale intensity
+# ======================
+
+def test_adapthist_ubyte():
+    '''Test a scalar uint8 image
+    '''
+    img = skimage.img_as_ubyte(data.moon())
+    adapted = exposure.adapthist(img, clip_limit=0.02)
+    assert adapted.min() == 0
+    assert adapted.max() == 255
+    assert img.shape == adapted.shape
+    assert peak_snr(img, adapted) > 22
+    assert norm_brightness_err(img, adapted) < 0.05
+    return img, adapted
+
+
+def test_adapthist_float():
+    '''Test an RGB float image
+    '''
+    img = skimage.img_as_float(data.lena())
+    adapted = exposure.adapthist(img, nx=10, ny=9, clip_limit=0.01,
+                        nbins=128, out_range='original')
+    assert_almost_equal(adapted.min() , img.min())
+    assert_almost_equal(adapted.min(), img.min())
+    assert img.shape == adapted.shape
+    assert peak_snr(img, adapted) > 136
+    assert norm_brightness_err(img, adapted) < 0.02
+    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())
+    mse = 1. / img1.size * np.square(img1 - img2).sum()
+    _, max_ = dtype_range[img1.dtype.type]
+    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
+        Normalize 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()

From 9eeb99089821a5be54bfe75e726414cdfd4a23b6 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sun, 5 Aug 2012 23:27:24 -0500
Subject: [PATCH 12/31] Tests pass and example works on build.

---
 doc/examples/plot_equalize.py           | 4 ++--
 skimage/exposure/__init__.py            | 2 +-
 skimage/exposure/tests/test_exposure.py | 2 ++
 3 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py
index 2cb9cbea..8b06890a 100644
--- a/doc/examples/plot_equalize.py
+++ b/doc/examples/plot_equalize.py
@@ -64,7 +64,7 @@ img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))
 img_eq = exposure.equalize(img)
 
 # Adaptive Equalization
-img_adapteq = exposure.adapthist(img)
+img_adapteq = exposure.adapthist(img, clip_limit=0.03)
 
 
 # Display results
@@ -81,7 +81,7 @@ 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_eq, axes[:, 3])
+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')
 
diff --git a/skimage/exposure/__init__.py b/skimage/exposure/__init__.py
index 7e19d317..a5e6c25a 100644
--- a/skimage/exposure/__init__.py
+++ b/skimage/exposure/__init__.py
@@ -1,2 +1,2 @@
-from .exposure import histogram, equalize, cumulative_distribution
+from .exposure import histogram, equalize, cumulative_distribution, adapthist
 from .exposure import rescale_intensity
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index fbad207b..963758f3 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -5,6 +5,7 @@ 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
@@ -94,6 +95,7 @@ def test_adapthist_float():
     img = skimage.img_as_float(data.lena())
     adapted = exposure.adapthist(img, nx=10, ny=9, clip_limit=0.01,
                         nbins=128, out_range='original')
+    assert_almost_equal = np.testing.assert_almost_equal
     assert_almost_equal(adapted.min() , img.min())
     assert_almost_equal(adapted.min(), img.min())
     assert img.shape == adapted.shape

From ea850757c8232e575de5c00031856f6dfedc04b7 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sun, 5 Aug 2012 23:36:39 -0500
Subject: [PATCH 13/31] Added Notes to the docstring.

---
 skimage/exposure/exposure.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index 1c4c69da..c37afcc3 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -218,6 +218,14 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
     -------
     out - np.ndarray :
         equalized image - may be a different shape than the original
+        
+    Notes
+    -----
+    * The underlying 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 RGB or RGBA images, the algorithm is run on each channel.
+    * For RGBA images, the original alpha channel is removed.
     '''
     in_type = image.dtype.type
     if out_range == 'full':

From 715271bde9edab2747ebfac6ab8e5c500ced67a2 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sun, 5 Aug 2012 23:42:44 -0500
Subject: [PATCH 14/31] Pep8 fixes.

---
 doc/examples/plot_equalize.py           | 2 +-
 skimage/exposure/exposure.py            | 2 +-
 skimage/exposure/tests/test_exposure.py | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py
index 8b06890a..e09b5a07 100644
--- a/doc/examples/plot_equalize.py
+++ b/doc/examples/plot_equalize.py
@@ -26,6 +26,7 @@ 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.
 
@@ -88,4 +89,3 @@ ax_cdf.set_ylabel('Fraction of total intensity')
 # prevent overlap of y-axis labels
 plt.subplots_adjust(wspace=0.4)
 plt.show()
-
diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index c37afcc3..0a2a253d 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -218,7 +218,7 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
     -------
     out - np.ndarray :
         equalized image - may be a different shape than the original
-        
+
     Notes
     -----
     * The underlying algorithm relies on an image whose rows and columns are even multiples of
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 963758f3..9d70969a 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -96,7 +96,7 @@ def test_adapthist_float():
     adapted = exposure.adapthist(img, nx=10, ny=9, clip_limit=0.01,
                         nbins=128, out_range='original')
     assert_almost_equal = np.testing.assert_almost_equal
-    assert_almost_equal(adapted.min() , img.min())
+    assert_almost_equal(adapted.min(), img.min())
     assert_almost_equal(adapted.min(), img.min())
     assert img.shape == adapted.shape
     assert peak_snr(img, adapted) > 136

From 20f4cee5da839edf68bded25d198f809e81cab05 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Mon, 6 Aug 2012 19:06:57 -0500
Subject: [PATCH 15/31] Replaced generic ACM license with the Graphics Gems IV
 license.

---
 skimage/exposure/clahe.h | 31 ++++++++++++-------------------
 1 file changed, 12 insertions(+), 19 deletions(-)

diff --git a/skimage/exposure/clahe.h b/skimage/exposure/clahe.h
index 8b115aa1..a2721b98 100644
--- a/skimage/exposure/clahe.h
+++ b/skimage/exposure/clahe.h
@@ -5,26 +5,19 @@ ANSI C code from the article
  * in "Graphics Gems IV", Academic Press, 1994
 =============
 
-ACM Software License Agreement
+http://tog.acm.org/resources/GraphicsGems/
 
-All software, both binary and source published by the Association for Computing Machinery (hereafter, Software) is copyrighted by the Association (hereafter, ACM) and ownership of all right, title and interest in and to the Software remains with ACM. By using or copying the Software, User agrees to abide by the terms of this Agreement.
-Noncommercial Use
+EULA: 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.
 
-The ACM grants to you (hereafter, User) a royalty-free, nonexclusive right to execute, copy, modify and distribute both the binary and source code solely for academic, research and other similar noncommercial uses, subject to the following conditions:
-
-    User acknowledges that the Software is still in the development stage and that it is being supplied "as is," without any support services from ACM. Neither ACM nor the author makes any representations or warranties, express or implied, including, without limitation, any representations or warranties of the merchantability or fitness for any particular purpose, or that the application of the software, will not infringe on any patents or other proprietary rights of others.
-    ACM shall not be held liable for direct, indirect, incidental or consequential damages arising from any claim by User or any third party with respect to uses allowed under this Agreement, or from any use of the Software.
-    User agrees to fully indemnify and hold harmless ACM and/or the author(s) of the original work from and against any and all claims, demands, suits, losses, damages, costs and expenses arising out of the User's use of the Software, including, without limitation, arising out of the User's modification of the Software.
-    User may modify the Software and distribute that modified work to third parties provided that: (a) if posted separately, it clearly acknowledges that it contains material copyrighted by ACM (b) no charge is associated with such copies, (c) User agrees to notify ACM and the Author(s) of the distribution, and (d) User clearly notifies secondary users that such modified work is not the original Software.
-    User agrees that ACM, the authors of the original work and others may enjoy a royalty-free, non-exclusive license to use, copy, modify and redistribute these modifications to the Software made by the User and distributed to third parties as a derivative work under this agreement.
-    This agreement will terminate immediately upon User's breach of, or non-compliance with, any of its terms. User may be held liable for any copyright infringement or the infringement of any other proprietary rights in the Software that is caused or facilitated by the User's failure to abide by the terms of this agreement.
-    This agreement will be construed and enforced in accordance with the law of the state of New York applicable to contracts performed entirely within the State. The parties irrevocably consent to the exclusive jurisdiction of the state or federal courts located in the City of New York for all disputes concerning this agreement. 
-
-Commercial Use
-
-Any User wishing to make a commercial use of the Software must contact ACM at permissions@acm.org to arrange an appropriate license. Commercial use includes (1) integrating or incorporating all or part of the source code into a product for sale or license by, or on behalf of, User to third parties, or (2) distribution of the binary or source code to third parties for use with a commercial product sold or licensed by, or on behalf of, User.
-
-Revised 6/98
 */
 
 #ifdef BYTE_IMAGE
@@ -107,7 +100,7 @@ int CLAHE (kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes,
     kz_pixel_t aLUT[uiNR_OF_GREY];	    /* lookup table used for scaling of input image */
     unsigned long* pulHist, *pulMapArray; /* pointer to histogram and mappings*/
     unsigned long* pulLU, *pulLB, *pulRU, *pulRB; /* auxiliary pointers interpolation */
-    
+
     if (uiNrX > uiMAX_REG_X) return -1;	   /* # of regions x-direction too large */
     if (uiNrY > uiMAX_REG_Y) return -2;	   /* # of regions y-direction too large */
     if (uiXRes % uiNrX) return -3;	  /* x-resolution no multiple of uiNrX */

From c7c4db6e3c2556991695f81c679113e111139254 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sat, 11 Aug 2012 15:48:18 -0500
Subject: [PATCH 16/31] Improved support for color images using lab transform.

---
 skimage/exposure/exposure.py            | 60 ++++++++++++++++++++-----
 skimage/exposure/tests/test_exposure.py | 38 ++++++++++++----
 2 files changed, 77 insertions(+), 21 deletions(-)

diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index 0a2a253d..d34068f6 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -2,7 +2,7 @@ import numpy as np
 
 from skimage import img_as_float
 from skimage.util.dtype import dtype_range
-from skimage.color import rgb2gray
+import skimage.color as color
 from skimage.util.dtype import convert
 
 from _adapthist import _adapthist
@@ -216,7 +216,7 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
 
     Returns
     -------
-    out - np.ndarray :
+    out : np.ndarray
         equalized image - may be a different shape than the original
 
     Notes
@@ -224,8 +224,17 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
     * The underlying 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 RGB or RGBA images, the algorithm is run on each channel.
+    * For grayscale images, CLAHE is performed on one channel, and a grayscale is returned
+    * 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/
+    .. [2] https://en.wikipedia.org/wiki/CLAHE#CLAHE
     '''
     in_type = image.dtype.type
     if out_range == 'full':
@@ -233,23 +242,50 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
     else:
         out_range = (image.min(), image.max())
     # must be converted to 12 bit for CLAHE
-    image = skimage.img_as_uint(image)
+    int_image = skimage.img_as_uint(image)
     MAX_VAL = 2 ** 12 - 1
-    image = rescale_intensity(image, out_range=(0, MAX_VAL))
+    int_image = rescale_intensity(int_image, out_range=(0, MAX_VAL))
     # handle color images - CLAHE accepts scalar images only
-    args = [image.copy(), 0, MAX_VAL, nx, ny, nbins, clip_limit]
+    args = [int_image.copy(), 0, MAX_VAL, nx, ny, nbins, clip_limit]
     if image.ndim == 3:
-        image = image[:, :, :3]
-        for channel in range(3):
-            args[0] = image[:, :, channel]
+        # check for grayscale
+        if (np.allclose(image[:, :, 0], image[:, :, 1]) and
+            np.allclose(image[:, :, 2], image[:, :, 3])):
+            args[0] = image[:, :, 0]
             out = _adapthist(*args)
-            image[:out.shape[0], :out.shape[1], channel] = out
+            image = int_image[:, :, :3]
+            for channel in range(3):
+                image[:out.shape[0], :out.shape[1], channel] = out
+        # for color images, convert to LAB space for processing
+        else:
+            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, MAX_VAL))
+            new_L = _adapthist(*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:
         out = _adapthist(*args)
+        image = int_image
         image[:out.shape[0], :out.shape[1]] = out
     # restore to desired output type and output limits
     image = rescale_intensity(image)
-    if in_type != np.uint16:
-        image = convert(image, in_type)
+    image = convert(image, in_type)
     image = rescale_intensity(image, out_range=out_range)
     return image
+
+if __name__ == '__main__':
+    from skimage import data
+    import matplotlib.pyplot as plt
+    img = skimage.img_as_uint(data.lena())
+    adapted = adapthist(img, nx=10, ny=9, clip_limit=0.01,
+                            nbins=128, out_range='original')
+    plt.imshow(img)
+    plt.figure(); plt.imshow(skimage.img_as_ubyte(adapted))
+    plt.figure(); plt.imshow(color.lab2rgb(color.rgb2lab(img)))
+    plt.show()
+    print 'Done'
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 9d70969a..8ce7c8c6 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -76,7 +76,7 @@ def test_rescale_out_range():
 # Test rescale intensity
 # ======================
 
-def test_adapthist_ubyte():
+def test_adapthist_scalar():
     '''Test a scalar uint8 image
     '''
     img = skimage.img_as_ubyte(data.moon())
@@ -84,23 +84,43 @@ def test_adapthist_ubyte():
     assert adapted.min() == 0
     assert adapted.max() == 255
     assert img.shape == adapted.shape
-    assert peak_snr(img, adapted) > 22
-    assert norm_brightness_err(img, adapted) < 0.05
+    full_scale = skimage.exposure.rescale_intensity(img)
+    assert_almost_equal = np.testing.assert_almost_equal
+    assert_almost_equal(peak_snr(full_scale, adapted), 22.19920073)
+    assert_almost_equal(norm_brightness_err(full_scale, adapted),
+                        0.04161278)
     return img, adapted
 
 
-def test_adapthist_float():
-    '''Test an RGB float image
+def test_adapthist_grayscale():
+    '''Test a grayscale float image
     '''
     img = skimage.img_as_float(data.lena())
+    img = rgb2gray(img)
     adapted = exposure.adapthist(img, nx=10, ny=9, clip_limit=0.01,
                         nbins=128, out_range='original')
     assert_almost_equal = np.testing.assert_almost_equal
     assert_almost_equal(adapted.min(), img.min())
-    assert_almost_equal(adapted.min(), img.min())
+    assert_almost_equal(adapted.max(), img.max())
     assert img.shape == adapted.shape
-    assert peak_snr(img, adapted) > 136
-    assert norm_brightness_err(img, adapted) < 0.02
+    assert_almost_equal(peak_snr(img, adapted), 131.4962063)
+    assert_almost_equal(norm_brightness_err(img, adapted), 0.0208805)
+    return data, adapted
+
+
+def test_adapthist_color():
+    '''Test a color uint16 image
+    '''
+    img = skimage.img_as_uint(data.lena())
+    adapted = exposure.adapthist(img, clip_limit=0.01)
+    assert_almost_equal = np.testing.assert_almost_equal
+    assert adapted.min() == 0
+    assert adapted.max() == 65535
+    assert img.shape == adapted.shape
+    full_scale = skimage.exposure.rescale_intensity(img)
+    assert_almost_equal(peak_snr(full_scale, adapted), 64.29546231)
+    assert_almost_equal(norm_brightness_err(full_scale, adapted),
+                        0.181473754)
     return data, adapted
 
 
@@ -135,7 +155,7 @@ def norm_brightness_err(img1, img2):
     Returns
     -------
     norm_brightness_error : float
-        Normalize absolute mean brightness error
+        Normalized absolute mean brightness error
     '''
     if img1.ndim == 3:
         img1, img2 = rgb2gray(img1), rgb2gray(img2)

From 269929333b3c72e9b5a7da34fb2e8a99607f1aab Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sat, 11 Aug 2012 21:28:08 -0500
Subject: [PATCH 17/31] Passed all tests with full coverage, PEP8 and Pylint
 performed.

---
 skimage/exposure/exposure.py            | 44 ++++++++++---------------
 skimage/exposure/tests/test_exposure.py |  1 +
 2 files changed, 18 insertions(+), 27 deletions(-)

diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index d34068f6..a85a209d 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -194,16 +194,17 @@ def rescale_intensity(image, in_range=None, out_range=None):
     return dtype(image * (omax - omin) + omin)
 
 
-def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
+def adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, nbins=256,
+               out_range='full'):
     '''Contrast Limited Adaptive Histogram Equalization
 
     Parameters
     ----------
     image : array-like
         original image
-    nx : int, optional
+    ntiles_x : int, optional
         Tile regions in the X direction (2, 16)
-    ny : int, optional
+    ntiles_y : int, optional
         Tile regions in the Y direction (2, 16)
     clip_limit : float: optional
         Normalized cliplimit (higher values give more contrast)
@@ -243,15 +244,16 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
         out_range = (image.min(), image.max())
     # must be converted to 12 bit for CLAHE
     int_image = skimage.img_as_uint(image)
-    MAX_VAL = 2 ** 12 - 1
-    int_image = rescale_intensity(int_image, out_range=(0, MAX_VAL))
+    max_val = 2 ** 12 - 1
+    int_image = rescale_intensity(int_image, out_range=(0, max_val))
     # handle color images - CLAHE accepts scalar images only
-    args = [int_image.copy(), 0, MAX_VAL, nx, ny, nbins, clip_limit]
+    args = [int_image.copy(), 0, max_val, ntiles_x, ntiles_y, nbins,
+            clip_limit]
     if image.ndim == 3:
         # check for grayscale
         if (np.allclose(image[:, :, 0], image[:, :, 1]) and
-            np.allclose(image[:, :, 2], image[:, :, 3])):
-            args[0] = image[:, :, 0]
+            np.allclose(image[:, :, 1], image[:, :, 2])):
+            args[0] = int_image[:, :, 0]
             out = _adapthist(*args)
             image = int_image[:, :, :3]
             for channel in range(3):
@@ -259,13 +261,13 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
         # for color images, convert to LAB space for processing
         else:
             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, MAX_VAL))
-            new_L = _adapthist(*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
+            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, max_val))
+            new_l = _adapthist(*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:
@@ -277,15 +279,3 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
     image = convert(image, in_type)
     image = rescale_intensity(image, out_range=out_range)
     return image
-
-if __name__ == '__main__':
-    from skimage import data
-    import matplotlib.pyplot as plt
-    img = skimage.img_as_uint(data.lena())
-    adapted = adapthist(img, nx=10, ny=9, clip_limit=0.01,
-                            nbins=128, out_range='original')
-    plt.imshow(img)
-    plt.figure(); plt.imshow(skimage.img_as_ubyte(adapted))
-    plt.figure(); plt.imshow(color.lab2rgb(color.rgb2lab(img)))
-    plt.show()
-    print 'Done'
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 8ce7c8c6..2e64e741 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -97,6 +97,7 @@ def test_adapthist_grayscale():
     '''
     img = skimage.img_as_float(data.lena())
     img = rgb2gray(img)
+    img = np.dstack((img, img, img))
     adapted = exposure.adapthist(img, nx=10, ny=9, clip_limit=0.01,
                         nbins=128, out_range='original')
     assert_almost_equal = np.testing.assert_almost_equal

From 94af01cdbde63b6d4a55b4235732c084fd696241 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sun, 21 Oct 2012 13:25:18 -0500
Subject: [PATCH 18/31] Translated clahe C code to python, added to list of
 contributors.

---
 CONTRIBUTORS.txt                        |   7 +
 skimage/exposure/__init__.py            |   3 +-
 skimage/exposure/_adapthist.py          | 337 ++++++++++++++++++++++++
 skimage/exposure/_adapthist.pyx         |  40 ---
 skimage/exposure/clahe.h                | 322 ----------------------
 skimage/exposure/exposure.py            |  90 +------
 skimage/exposure/setup.py               |  31 ---
 skimage/exposure/tests/test_exposure.py |  28 +-
 8 files changed, 360 insertions(+), 498 deletions(-)
 create mode 100644 skimage/exposure/_adapthist.py
 delete mode 100644 skimage/exposure/_adapthist.pyx
 delete mode 100644 skimage/exposure/clahe.h
 delete mode 100644 skimage/exposure/setup.py

diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index 012fd53d..5f796b3f 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -123,3 +123,10 @@
 
 - Luis Pedro Coelho
   imread plugin
+=======
+  Polygon, circle and ellipse drawing functions
+  Adaptive thresholding
+  Implementation of Matlab's `regionprops`
+
+- Steven Silvester, Karel Zuiderveld
+  Adaptive Histogram Equalization
diff --git a/skimage/exposure/__init__.py b/skimage/exposure/__init__.py
index a5e6c25a..f8a4f046 100644
--- a/skimage/exposure/__init__.py
+++ b/skimage/exposure/__init__.py
@@ -1,2 +1,3 @@
-from .exposure import histogram, equalize, cumulative_distribution, adapthist
+from .exposure import histogram, equalize, cumulative_distribution
 from .exposure import rescale_intensity
+from ._adaphist import adapthist
diff --git a/skimage/exposure/_adapthist.py b/skimage/exposure/_adapthist.py
new file mode 100644
index 00000000..f35a5da5
--- /dev/null
+++ b/skimage/exposure/_adapthist.py
@@ -0,0 +1,337 @@
+'''
+Adapted code from the article
+ * "Contrast Limited Adaptive Histogram Equalization"
+ * by Karel Zuiderveld, karel@cv.ruu.nl
+ * in "Graphics Gems IV", Academic Press, 1994
+=============
+
+http://tog.acm.org/resources/GraphicsGems/
+
+EULA: 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.
+ *
+ *  Author: Karel Zuiderveld, Computer Vision Research Group,
+ *	     Utrecht, The Netherlands (karel@cv.ruu.nl)
+'''
+import numpy as np
+import skimage
+from skimage import color
+from skimage.exposure import rescale_intensity
+
+
+MAX_REG_X = 16  # max. # contextual regions in x-direction */
+MAX_REG_Y = 16  # max. # contextual regions in y-direction */
+NR_OF_GREY = 1 << 14  # number of grayscale levels to use in CLAHE algorithm
+
+
+def adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, nbins=256):
+    '''Contrast Limited Adaptive Histogram Equalization
+
+    Parameters
+    ----------
+    image : array-like
+        original image
+    ntiles_x : int, optional
+        Tile regions in the X direction (2, 16)
+    ntiles_y : int, optional
+        Tile regions in the Y direction (2, 16)
+    clip_limit : float: optional
+        Normalized cliplimit (higher values give more contrast)
+    nbins : int, optional
+        Greybins for histogram ("dynamic range")
+
+    Returns
+    -------
+    out : np.ndarray
+        equalized image - grayscale images are uint16, color images are float
+
+    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 grayscale images, CLAHE is performed on one channel,
+      and a grayscale is returned
+    * 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/
+    .. [2] https://en.wikipedia.org/wiki/CLAHE#CLAHE
+    '''
+    # convert to uint image
+    int_image = skimage.img_as_uint(image)
+    int_image = rescale_intensity(int_image, out_range=(0, NR_OF_GREY - 1))
+    # handle color images - CLAHE accepts scalar images only
+    args = [int_image.copy(), ntiles_x, ntiles_y, clip_limit * nbins, nbins]
+    if image.ndim == 3:
+        # check for grayscale
+        if (np.allclose(image[:, :, 0], image[:, :, 1]) and
+            np.allclose(image[:, :, 1], image[:, :, 2])):
+            args[0] = int_image[:, :, 0]
+            out = _clahe(*args)
+            image = int_image[:, :, :3]
+            for channel in range(3):
+                image[:out.shape[0], :out.shape[1], channel] = out
+        # for color images, convert to LAB space for processing
+        else:
+            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:
+        out = _clahe(*args)
+        image = int_image
+        image[:out.shape[0], :out.shape[1]] = out
+    return rescale_intensity(image)
+
+
+def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
+    '''Contrast Limited Adaptive Histogram Equalization
+
+    Parameters
+    ----------
+    image : array-like
+        original image
+    ntiles_x : int, optional
+        Tile regions in the X direction (2, 16)
+    ntiles_y : int, optional
+        Tile regions in the Y direction (2, 16)
+    clip_limit : float: optional
+        Normalized cliplimit (higher values give more contrast)
+    nbins : int, optional
+        Greybins for histogram ("dynamic range")
+
+    Returns
+    -------
+    out : np.ndarray
+
+    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_x * ntiles_y, 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
+    # Calculate greylevel mappings for each contextual region
+    ystart = 0
+    for y in range(ntiles_y):
+        xstart = 0
+        for x in range(ntiles_x):
+            sub_img = image[ystart: ystart + y_size,
+                            xstart: xstart + x_size]
+            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 * ntiles_x + x] = hist
+            xstart += x_size
+        ystart += y_size
+    # 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 * ntiles_x + xL]
+            mapRU = map_array[yU * ntiles_x + xR]
+            mapLB = map_array[yB * ntiles_x + xL]
+            mapRB = map_array[yB * ntiles_x + xR]
+            interpolate(image, xstart, xstep, ystart, ystep,
+                        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 : np.ndarray
+        histogram array
+    clip_limit : int
+        maximum allowed bin count
+
+    Returns
+    -------
+    hist : np.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):
+    '''Calculates the equalized lookup table (mapping)
+
+    It does so by cumulating the input histogram.
+
+    hist : np.ndarray
+        clipped histogram
+    min_val : int
+        min value for mapping
+    max_val : int
+        max value for mapping
+    n_pixels : int
+        number of pixels in the region
+
+    Returns
+    -------
+    out : np.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, xstart, xstep, ystart, ystep,
+                mapLU, mapRU, mapLB, mapRB, aLUT):
+    '''Find the new grayscale level for a region using bilinear interpolation
+
+    Parameters
+    ----------
+    image : np.ndarray
+        full image
+    xstart, xstop : int
+       indices of xslice
+    ystart, ystop : int
+        indices of yslice
+    map* : np.ndarray
+        mappings of greylevels from histograms
+    aLUT : np.ndarray
+        maps grayscale levels in image to histogram levels
+
+    Returns
+    -------
+    out : np.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 = xstep * ystep  # Normalization factor
+
+    # interpolation weight matrices
+    x_coef, y_coef = np.meshgrid(np.arange(xstep),
+                                 np.arange(ystep))
+    x_inv_coef, y_inv_coef = x_coef[:, ::-1] + 1, y_coef[::-1] + 1
+
+    im_slice = image[ystart: ystart + ystep, xstart: xstart + xstep]
+    im_slice = aLUT[im_slice]
+    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)
+    image[ystart: ystart + ystep, xstart: xstart + xstep] = new
diff --git a/skimage/exposure/_adapthist.pyx b/skimage/exposure/_adapthist.pyx
deleted file mode 100644
index 16e84336..00000000
--- a/skimage/exposure/_adapthist.pyx
+++ /dev/null
@@ -1,40 +0,0 @@
-import numpy as np
-cimport numpy as np
-
-cdef extern from "clahe.h":
-    ctypedef unsigned int kz_pixel_t
-    int CLAHE(kz_pixel_t * pImage, unsigned int uiXRes, unsigned int uiYRes,
-        kz_pixel_t Min,
-        kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,
-        unsigned int uiNrBins, float fCliplimit)
-
-
-def _adapthist(np.ndarray[kz_pixel_t, ndim=2] image, min, max, nx, ny, nbins,
-                 clip_limit=0):
-    '''
-    image -  to the input/output image
-    min - Minimum greyvalue of input image (also becomes min of output image)
-    max - Maximum greyvalue of input image (also becomes max of output image)
-    nx - Number of tile regions in the X direction (min 2, max uiMAX_REG_X)
-    ny - Number of tile regions in the Y direction (min 2, max uiMAX_REG_Y)
-    nbins - Number of greybins for histogram ("dynamic range")
-    clip_limit - Normalized cliplimit (higher values give more contrast)
-
-    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.
-    '''
-    # we need the image to be divisible by nx and ny
-    uiYRes = image.shape[0] - image.shape[0] % ny
-    uiXRes = image.shape[1] - image.shape[1] % nx
-
-    cdef np.ndarray[kz_pixel_t, ndim = 1] flattened
-    flattened = image[:uiYRes, :uiXRes].ravel()
-
-    ret = CLAHE(< kz_pixel_t * > flattened.data, uiXRes, uiYRes, min,
-                max, nx, ny, nbins, clip_limit * nbins)
-    if ret < 0:
-        print 'clahe error: ', ret
-    return flattened.reshape((uiYRes, uiXRes))
diff --git a/skimage/exposure/clahe.h b/skimage/exposure/clahe.h
deleted file mode 100644
index a2721b98..00000000
--- a/skimage/exposure/clahe.h
+++ /dev/null
@@ -1,322 +0,0 @@
-/*
-ANSI C code from the article
- * "Contrast Limited Adaptive Histogram Equalization"
- * by Karel Zuiderveld, karel@cv.ruu.nl
- * in "Graphics Gems IV", Academic Press, 1994
-=============
-
-http://tog.acm.org/resources/GraphicsGems/
-
-EULA: 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.
-
-*/
-
-#ifdef BYTE_IMAGE
-typedef unsigned char kz_pixel_t;	 /* for 8 bit-per-pixel images */
-#define uiNR_OF_GREY (256)
-#else
-typedef unsigned short kz_pixel_t;	 /* for 12 bit-per-pixel images (default) */
-# define uiNR_OF_GREY (4096)
-#endif
-
-/*
- * ANSI C code from the article
- * "Contrast Limited Adaptive Histogram Equalization"
- * by Karel Zuiderveld, karel@cv.ruu.nl
- * in "Graphics Gems IV", Academic Press, 1994
- *
- *
- *  These functions implement Contrast Limited Adaptive Histogram Equalization.
- *  The main routine (CLAHE) expects an input image that is stored contiguously in
- *  memory;  the CLAHE output image overwrites the original input image and has the
- *  same minimum and maximum values (which must be provided by the user).
- *  This implementation assumes that the X- and Y image resolutions are an integer
- *  multiple of the X- and Y sizes of the contextual regions. A check on various other
- *  error conditions is performed.
- *
- *  #define the symbol BYTE_IMAGE to make this implementation suitable for
- *  8-bit images. The maximum number of contextual regions can be redefined
- *  by changing uiMAX_REG_X and/or uiMAX_REG_Y; the use of more than 256
- *  contextual regions is not recommended.
- *
- *  The code is ANSI-C and is also C++ compliant.
- *
- *  Author: Karel Zuiderveld, Computer Vision Research Group,
- *	     Utrecht, The Netherlands (karel@cv.ruu.nl)
- */
-
-/*********************** Local prototypes ************************/
-static void ClipHistogram (unsigned long*, unsigned int, unsigned long);
-static void MakeHistogram (kz_pixel_t*, unsigned int, unsigned int, unsigned int,
-		unsigned long*, unsigned int, kz_pixel_t*);
-static void MapHistogram (unsigned long*, kz_pixel_t, kz_pixel_t,
-	       unsigned int, unsigned long);
-static void MakeLut (kz_pixel_t*, kz_pixel_t, kz_pixel_t, unsigned int);
-static void Interpolate (kz_pixel_t*, int, unsigned long*, unsigned long*,
-	unsigned long*, unsigned long*, unsigned int, unsigned int, kz_pixel_t*);
-
-/**************	 Start of actual code **************/
-#include <stdlib.h>			 /* To get prototypes of malloc() and free() */
-
-
-const unsigned int uiMAX_REG_X = 16;	  /* max. # contextual regions in x-direction */
-const unsigned int uiMAX_REG_Y = 16;	  /* max. # contextual regions in y-direction */
-
-
-
-/************************** main function CLAHE ******************/
-int CLAHE (kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes,
-	 kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,
-	      unsigned int uiNrBins, float fCliplimit)
-/*   pImage - Pointer to the input/output image
- *   uiXRes - Image resolution in the X direction
- *   uiYRes - Image resolution in the Y direction
- *   Min - Minimum greyvalue of input image (also becomes minimum of output image)
- *   Max - Maximum greyvalue of input image (also becomes maximum of output image)
- *   uiNrX - Number of contextial regions in the X direction (min 2, max uiMAX_REG_X)
- *   uiNrY - Number of contextial regions in the Y direction (min 2, max uiMAX_REG_Y)
- *   uiNrBins - Number of greybins for histogram ("dynamic range")
- *   float fCliplimit - Normalized cliplimit (higher values give more contrast)
- * The number of "effective" greylevels in the output image is set by uiNrBins; 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.
- */
-{
-    unsigned int uiX, uiY;		  /* counters */
-    unsigned int uiXSize, uiYSize, uiSubX, uiSubY; /* size of context. reg. and subimages */
-    unsigned int uiXL, uiXR, uiYU, uiYB;  /* auxiliary variables interpolation routine */
-    unsigned long ulClipLimit, ulNrPixels;/* clip limit and region pixel count */
-    kz_pixel_t* pImPointer;		   /* pointer to image */
-    kz_pixel_t aLUT[uiNR_OF_GREY];	    /* lookup table used for scaling of input image */
-    unsigned long* pulHist, *pulMapArray; /* pointer to histogram and mappings*/
-    unsigned long* pulLU, *pulLB, *pulRU, *pulRB; /* auxiliary pointers interpolation */
-
-    if (uiNrX > uiMAX_REG_X) return -1;	   /* # of regions x-direction too large */
-    if (uiNrY > uiMAX_REG_Y) return -2;	   /* # of regions y-direction too large */
-    if (uiXRes % uiNrX) return -3;	  /* x-resolution no multiple of uiNrX */
-    if (uiYRes % uiNrY) return -4;	  /* y-resolution no multiple of uiNrY */
-    if (Max >= uiNR_OF_GREY) return -5;	   /* maximum too large */
-    if (Min >= Max) return -6;		  /* minimum equal or larger than maximum */
-    if (uiNrX < 2 || uiNrY < 2) return -7;/* at least 4 contextual regions required */
-    if (fCliplimit == 1.0) return 0;	  /* is OK, immediately returns original image. */
-    if (uiNrBins == 0) uiNrBins = 128;	  /* default value when not specified */
-
-    pulMapArray=(unsigned long *)malloc(sizeof(unsigned long)*uiNrX*uiNrY*uiNrBins);
-    if (pulMapArray == 0) return -8;	  /* Not enough memory! (try reducing uiNrBins) */
-
-    uiXSize = uiXRes/uiNrX; uiYSize = uiYRes/uiNrY;  /* Actual size of contextual regions */
-    ulNrPixels = (unsigned long)uiXSize * (unsigned long)uiYSize;
-
-    if(fCliplimit > 0.0) {		  /* Calculate actual cliplimit	 */
-       ulClipLimit = (unsigned long) (fCliplimit * (uiXSize * uiYSize) / uiNrBins);
-       ulClipLimit = (ulClipLimit < 1UL) ? 1UL : ulClipLimit;
-    }
-    else ulClipLimit = 1UL<<14;		  /* Large value, do not clip (AHE) */
-    MakeLut(aLUT, Min, Max, uiNrBins);	  /* Make lookup table for mapping of greyvalues */
-    /* Calculate greylevel mappings for each contextual region */
-    for (uiY = 0, pImPointer = pImage; uiY < uiNrY; uiY++) {
-	for (uiX = 0; uiX < uiNrX; uiX++, pImPointer += uiXSize) {
-	    pulHist = &pulMapArray[uiNrBins * (uiY * uiNrX + uiX)];
-	    MakeHistogram(pImPointer,uiXRes,uiXSize,uiYSize,pulHist,uiNrBins,aLUT);
-	    ClipHistogram(pulHist, uiNrBins, ulClipLimit);
-	    MapHistogram(pulHist, Min, Max, uiNrBins, ulNrPixels);
-	}
-	pImPointer += (uiYSize - 1) * uiXRes;		  /* skip lines, set pointer */
-    }
-    /* Interpolate greylevel mappings to get CLAHE image */
-    for (pImPointer = pImage, uiY = 0; uiY <= uiNrY; uiY++) {
-	if (uiY == 0) {					  /* special case: top row */
-	    uiSubY = uiYSize >> 1;  uiYU = 0; uiYB = 0;
-	}
-	else {
-	    if (uiY == uiNrY) {				  /* special case: bottom row */
-		uiSubY = uiYSize >> 1;	uiYU = uiNrY-1;	 uiYB = uiYU;
-	    }
-	    else {					  /* default values */
-		uiSubY = uiYSize; uiYU = uiY - 1; uiYB = uiYU + 1;
-	    }
-	}
-	for (uiX = 0; uiX <= uiNrX; uiX++) {
-	    if (uiX == 0) {				  /* special case: left column */
-		uiSubX = uiXSize >> 1; uiXL = 0; uiXR = 0;
-	    }
-	    else {
-		if (uiX == uiNrX) {			  /* special case: right column */
-		    uiSubX = uiXSize >> 1;  uiXL = uiNrX - 1; uiXR = uiXL;
-		}
-		else {					  /* default values */
-		    uiSubX = uiXSize; uiXL = uiX - 1; uiXR = uiXL + 1;
-		}
-	    }
-
-	    pulLU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXL)];
-	    pulRU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXR)];
-	    pulLB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXL)];
-	    pulRB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXR)];
-	    Interpolate(pImPointer,uiXRes,pulLU,pulRU,pulLB,pulRB,uiSubX,uiSubY,aLUT);
-	    pImPointer += uiSubX;			  /* set pointer on next matrix */
-	}
-	pImPointer += (uiSubY - 1) * uiXRes;
-    }
-    free(pulMapArray);					  /* free space for histograms */
-    return 0;						  /* return status OK */
-}
-void ClipHistogram (unsigned long* pulHistogram, unsigned int
-		     uiNrGreylevels, unsigned long ulClipLimit)
-/* This function performs 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).
- */
-{
-    unsigned long* pulBinPointer, *pulEndPointer, *pulHisto;
-    unsigned long ulNrExcess, ulUpper, ulBinIncr, ulStepSize, i;
-    long lBinExcess;
-
-    ulNrExcess = 0;  pulBinPointer = pulHistogram;
-    for (i = 0; i < uiNrGreylevels; i++) { /* calculate total number of excess pixels */
-	lBinExcess = (long) pulBinPointer[i] - (long) ulClipLimit;
-	if (lBinExcess > 0) ulNrExcess += lBinExcess;	  /* excess in current bin */
-    };
-
-    /* Second part: clip histogram and redistribute excess pixels in each bin */
-    ulBinIncr = ulNrExcess / uiNrGreylevels;		  /* average binincrement */
-    ulUpper =  ulClipLimit - ulBinIncr;	 /* Bins larger than ulUpper set to cliplimit */
-
-    for (i = 0; i < uiNrGreylevels; i++) {
-      if (pulHistogram[i] > ulClipLimit) pulHistogram[i] = ulClipLimit; /* clip bin */
-      else {
-	  if (pulHistogram[i] > ulUpper) {		/* high bin count */
-	      ulNrExcess -= pulHistogram[i] - ulUpper; pulHistogram[i]=ulClipLimit;
-	  }
-	  else {					/* low bin count */
-	      ulNrExcess -= ulBinIncr; pulHistogram[i] += ulBinIncr;
-	  }
-       }
-    }
-
-    while (ulNrExcess) {   /* Redistribute remaining excess  */
-	pulEndPointer = &pulHistogram[uiNrGreylevels]; pulHisto = pulHistogram;
-
-	while (ulNrExcess && pulHisto < pulEndPointer) {
-	    ulStepSize = uiNrGreylevels / ulNrExcess;
-	    if (ulStepSize < 1) ulStepSize = 1;		  /* stepsize at least 1 */
-	    for (pulBinPointer=pulHisto; pulBinPointer < pulEndPointer && ulNrExcess;
-		 pulBinPointer += ulStepSize) {
-		if (*pulBinPointer < ulClipLimit) {
-		    (*pulBinPointer)++;	 ulNrExcess--;	  /* reduce excess */
-		}
-	    }
-	    pulHisto++;		  /* restart redistributing on other bin location */
-	}
-    }
-}
-void MakeHistogram (kz_pixel_t* pImage, unsigned int uiXRes,
-		unsigned int uiSizeX, unsigned int uiSizeY,
-		unsigned long* pulHistogram,
-		unsigned int uiNrGreylevels, kz_pixel_t* pLookupTable)
-/* This function classifies the greylevels present in the array image into
- * a greylevel histogram. The pLookupTable specifies the relationship
- * between the greyvalue of the pixel (typically between 0 and 4095) and
- * the corresponding bin in the histogram (usually containing only 128 bins).
- */
-{
-    kz_pixel_t* pImagePointer;
-    unsigned int i;
-
-    for (i = 0; i < uiNrGreylevels; i++) pulHistogram[i] = 0L; /* clear histogram */
-
-    for (i = 0; i < uiSizeY; i++) {
-	pImagePointer = &pImage[uiSizeX];
-	while (pImage < pImagePointer) pulHistogram[pLookupTable[*pImage++]]++;
-	pImagePointer += uiXRes;
-	pImage = &pImagePointer[-uiSizeX];
-    }
-}
-
-void MapHistogram (unsigned long* pulHistogram, kz_pixel_t Min, kz_pixel_t Max,
-	       unsigned int uiNrGreylevels, unsigned long ulNrOfPixels)
-/* This function calculates the equalized lookup table (mapping) by
- * cumulating the input histogram. Note: lookup table is rescaled in range [Min..Max].
- */
-{
-    unsigned int i;  unsigned long ulSum = 0;
-    const float fScale = ((float)(Max - Min)) / ulNrOfPixels;
-    const unsigned long ulMin = (unsigned long) Min;
-
-    for (i = 0; i < uiNrGreylevels; i++) {
-	ulSum += pulHistogram[i]; pulHistogram[i]=(unsigned long)(ulMin+ulSum*fScale);
-	if (pulHistogram[i] > Max) pulHistogram[i] = Max;
-    }
-}
-
-void MakeLut (kz_pixel_t * pLUT, kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrBins)
-/* To speed up histogram clipping, the input image [Min,Max] is scaled down to
- * [0,uiNrBins-1]. This function calculates the LUT.
- */
-{
-    int i;
-    const kz_pixel_t BinSize = (kz_pixel_t) (1 + (Max - Min) / uiNrBins);
-
-    for (i = Min; i <= Max; i++)  pLUT[i] = (i - Min) / BinSize;
-}
-
-void Interpolate (kz_pixel_t * pImage, int uiXRes, unsigned long * pulMapLU,
-     unsigned long * pulMapRU, unsigned long * pulMapLB,  unsigned long * pulMapRB,
-     unsigned int uiXSize, unsigned int uiYSize, kz_pixel_t * pLUT)
-/* pImage      - pointer to input/output image
- * uiXRes      - resolution of image in x-direction
- * pulMap*     - mappings of greylevels from histograms
- * uiXSize     - uiXSize of image submatrix
- * uiYSize     - uiYSize of image submatrix
- * pLUT	       - lookup table containing mapping greyvalues to bins
- * This function calculates the new greylevel assignments of pixels within a submatrix
- * of the image with size uiXSize and uiYSize. This is done by a bilinear interpolation
- * between four different mappings in order to eliminate boundary artifacts.
- * It uses a division; since division is often an expensive operation, I added code to
- * perform a logical shift instead when feasible.
- */
-{
-    const unsigned int uiIncr = uiXRes-uiXSize; /* Pointer increment after processing row */
-    kz_pixel_t GreyValue; unsigned int uiNum = uiXSize*uiYSize; /* Normalization factor */
-
-    unsigned int uiXCoef, uiYCoef, uiXInvCoef, uiYInvCoef, uiShift = 0;
-
-    if (uiNum & (uiNum - 1))   /* If uiNum is not a power of two, use division */
-    for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
-	 uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {
-	for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
-	     uiXCoef++, uiXInvCoef--) {
-	    GreyValue = pLUT[*pImage];		   /* get histogram bin value */
-	    *pImage++ = (kz_pixel_t ) ((uiYInvCoef * (uiXInvCoef*pulMapLU[GreyValue]
-				      + uiXCoef * pulMapRU[GreyValue])
-				+ uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
-				      + uiXCoef * pulMapRB[GreyValue])) / uiNum);
-	}
-    }
-    else {			   /* avoid the division and use a right shift instead */
-	while (uiNum >>= 1) uiShift++;		   /* Calculate 2log of uiNum */
-	for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;
-	     uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {
-	     for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;
-	       uiXCoef++, uiXInvCoef--) {
-	       GreyValue = pLUT[*pImage];	  /* get histogram bin value */
-	       *pImage++ = (kz_pixel_t)((uiYInvCoef* (uiXInvCoef * pulMapLU[GreyValue]
-				      + uiXCoef * pulMapRU[GreyValue])
-				+ uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]
-				      + uiXCoef * pulMapRB[GreyValue])) >> uiShift);
-	    }
-	}
-    }
-}
-
diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index a85a209d..6ef7ff66 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -5,10 +5,9 @@ from skimage.util.dtype import dtype_range
 import skimage.color as color
 from skimage.util.dtype import convert
 
-from _adapthist import _adapthist
 
 __all__ = ['histogram', 'cumulative_distribution', 'equalize',
-           'rescale_intensity', 'adapthist']
+           'rescale_intensity']
 
 
 def histogram(image, nbins=256):
@@ -192,90 +191,3 @@ def rescale_intensity(image, in_range=None, out_range=None):
 
     image = (image - imin) / float(imax - imin)
     return dtype(image * (omax - omin) + omin)
-
-
-def adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, nbins=256,
-               out_range='full'):
-    '''Contrast Limited Adaptive Histogram Equalization
-
-    Parameters
-    ----------
-    image : array-like
-        original image
-    ntiles_x : int, optional
-        Tile regions in the X direction (2, 16)
-    ntiles_y : int, optional
-        Tile regions in the Y direction (2, 16)
-    clip_limit : float: optional
-        Normalized cliplimit (higher values give more contrast)
-    nbins : int, optional
-        Greybins for histogram ("dynamic range")
-    out_range : str, optional
-        Range of the output image data.
-           - 'original' - Use original image limits
-           - 'full' - Use full range of image data type
-
-    Returns
-    -------
-    out : np.ndarray
-        equalized image - may be a different shape than the original
-
-    Notes
-    -----
-    * The underlying 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 grayscale images, CLAHE is performed on one channel, and a grayscale is returned
-    * 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/
-    .. [2] https://en.wikipedia.org/wiki/CLAHE#CLAHE
-    '''
-    in_type = image.dtype.type
-    if out_range == 'full':
-        out_range = None
-    else:
-        out_range = (image.min(), image.max())
-    # must be converted to 12 bit for CLAHE
-    int_image = skimage.img_as_uint(image)
-    max_val = 2 ** 12 - 1
-    int_image = rescale_intensity(int_image, out_range=(0, max_val))
-    # handle color images - CLAHE accepts scalar images only
-    args = [int_image.copy(), 0, max_val, ntiles_x, ntiles_y, nbins,
-            clip_limit]
-    if image.ndim == 3:
-        # check for grayscale
-        if (np.allclose(image[:, :, 0], image[:, :, 1]) and
-            np.allclose(image[:, :, 1], image[:, :, 2])):
-            args[0] = int_image[:, :, 0]
-            out = _adapthist(*args)
-            image = int_image[:, :, :3]
-            for channel in range(3):
-                image[:out.shape[0], :out.shape[1], channel] = out
-        # for color images, convert to LAB space for processing
-        else:
-            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, max_val))
-            new_l = _adapthist(*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:
-        out = _adapthist(*args)
-        image = int_image
-        image[:out.shape[0], :out.shape[1]] = out
-    # restore to desired output type and output limits
-    image = rescale_intensity(image)
-    image = convert(image, in_type)
-    image = rescale_intensity(image, out_range=out_range)
-    return image
diff --git a/skimage/exposure/setup.py b/skimage/exposure/setup.py
deleted file mode 100644
index 43bab23f..00000000
--- a/skimage/exposure/setup.py
+++ /dev/null
@@ -1,31 +0,0 @@
-#!/usr/bin/env python
-
-import os
-from skimage._build import cython
-
-base_path = os.path.abspath(os.path.dirname(__file__))
-
-
-def configuration(parent_package='', top_path=None):
-    from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
-
-    config = Configuration('exposure', parent_package, top_path)
-    config.add_data_dir('tests')
-
-    cython(['_adapthist.pyx'], working_path=base_path)
-
-    config.add_extension('_adapthist', sources=['_adapthist.c'],
-                         include_dirs=[get_numpy_include_dirs()])
-
-    return config
-
-if __name__ == '__main__':
-    from numpy.distutils.core import setup
-    setup(maintainer='scikits-image developers',
-          author='scikits-image developers',
-          maintainer_email='scikits-image@googlegroups.com',
-          description='Exposure',
-          url='https://github.com/scikits-image/scikits-image',
-          license='SciPy License (BSD Style)',
-          **(configuration(top_path='').todict())
-          )
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 2e64e741..3cf42a7f 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -1,6 +1,7 @@
 import numpy as np
 from numpy.testing import assert_array_almost_equal as assert_close
-
+import sys
+sys.path.insert(0, '../..')
 import skimage
 from skimage import data
 from skimage import exposure
@@ -82,30 +83,27 @@ def test_adapthist_scalar():
     img = skimage.img_as_ubyte(data.moon())
     adapted = exposure.adapthist(img, clip_limit=0.02)
     assert adapted.min() == 0
-    assert adapted.max() == 255
+    assert adapted.max() == (1<<16) - 1
     assert img.shape == adapted.shape
-    full_scale = skimage.exposure.rescale_intensity(img)
+    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), 22.19920073)
+    assert_almost_equal(peak_snr(full_scale, adapted), 28.6262034)
     assert_almost_equal(norm_brightness_err(full_scale, adapted),
-                        0.04161278)
+                        0.0410010)
     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.adapthist(img, nx=10, ny=9, clip_limit=0.01,
-                        nbins=128, out_range='original')
+    adapted = exposure.adapthist(img, 10, 9, clip_limit=0.01,
+                        nbins=128)
     assert_almost_equal = np.testing.assert_almost_equal
-    assert_almost_equal(adapted.min(), img.min())
-    assert_almost_equal(adapted.max(), img.max())
     assert img.shape == adapted.shape
-    assert_almost_equal(peak_snr(img, adapted), 131.4962063)
-    assert_almost_equal(norm_brightness_err(img, adapted), 0.0208805)
+    assert_almost_equal(peak_snr(img, adapted), 106.3020173)
+    assert_almost_equal(norm_brightness_err(img, adapted), 0.0218686)
     return data, adapted
 
 
@@ -116,12 +114,12 @@ def test_adapthist_color():
     adapted = exposure.adapthist(img, clip_limit=0.01)
     assert_almost_equal = np.testing.assert_almost_equal
     assert adapted.min() == 0
-    assert adapted.max() == 65535
+    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.29546231)
+    assert_almost_equal(peak_snr(full_scale, adapted), 64.7167515)
     assert_almost_equal(norm_brightness_err(full_scale, adapted),
-                        0.181473754)
+                        0.17922429)
     return data, adapted
 
 

From 7f86b0c40068e22bdb6ef65d1c9b03686fc4a47e Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sun, 21 Oct 2012 13:49:35 -0500
Subject: [PATCH 19/31] Removed sys.path hack in test, fixed spelling of import
 in __init__.py

---
 skimage/exposure/__init__.py            | 2 +-
 skimage/exposure/tests/test_exposure.py | 2 --
 2 files changed, 1 insertion(+), 3 deletions(-)

diff --git a/skimage/exposure/__init__.py b/skimage/exposure/__init__.py
index f8a4f046..d495686a 100644
--- a/skimage/exposure/__init__.py
+++ b/skimage/exposure/__init__.py
@@ -1,3 +1,3 @@
 from .exposure import histogram, equalize, cumulative_distribution
 from .exposure import rescale_intensity
-from ._adaphist import adapthist
+from ._adapthist import adapthist
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 3cf42a7f..04c0cf0e 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -1,7 +1,5 @@
 import numpy as np
 from numpy.testing import assert_array_almost_equal as assert_close
-import sys
-sys.path.insert(0, '../..')
 import skimage
 from skimage import data
 from skimage import exposure

From b8844beb793f2341677ba4411307ea1624fb6c80 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sat, 1 Dec 2012 23:46:52 -0600
Subject: [PATCH 20/31] Cleaned up adapthist func to use LAB for all
 RGB->required rest result change

---
 skimage/exposure/_adapthist.py          | 43 +++++++++----------------
 skimage/exposure/tests/test_exposure.py |  4 +--
 2 files changed, 18 insertions(+), 29 deletions(-)

diff --git a/skimage/exposure/_adapthist.py b/skimage/exposure/_adapthist.py
index f35a5da5..091514b6 100644
--- a/skimage/exposure/_adapthist.py
+++ b/skimage/exposure/_adapthist.py
@@ -71,37 +71,26 @@ def adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, nbins=256):
     .. [1] http://tog.acm.org/resources/GraphicsGems/
     .. [2] https://en.wikipedia.org/wiki/CLAHE#CLAHE
     '''
-    # convert to uint image
-    int_image = skimage.img_as_uint(image)
-    int_image = rescale_intensity(int_image, out_range=(0, NR_OF_GREY - 1))
     # handle color images - CLAHE accepts scalar images only
-    args = [int_image.copy(), ntiles_x, ntiles_y, clip_limit * nbins, nbins]
-    if image.ndim == 3:
-        # check for grayscale
-        if (np.allclose(image[:, :, 0], image[:, :, 1]) and
-            np.allclose(image[:, :, 1], image[:, :, 2])):
-            args[0] = int_image[:, :, 0]
-            out = _clahe(*args)
-            image = int_image[:, :, :3]
-            for channel in range(3):
-                image[:out.shape[0], :out.shape[1], channel] = out
-        # for color images, convert to LAB space for processing
-        else:
-            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))
+    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 = int_image
         image[:out.shape[0], :out.shape[1]] = out
-    return rescale_intensity(image)
+        image = rescale_intensity(image)
+    return image
 
 
 def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 04c0cf0e..dc48a3d8 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -100,8 +100,8 @@ def test_adapthist_grayscale():
                         nbins=128)
     assert_almost_equal = np.testing.assert_almost_equal
     assert img.shape == adapted.shape
-    assert_almost_equal(peak_snr(img, adapted), 106.3020173)
-    assert_almost_equal(norm_brightness_err(img, adapted), 0.0218686)
+    assert_almost_equal(peak_snr(img, adapted), 77.5841845)
+    assert_almost_equal(norm_brightness_err(img, adapted), 0.0376632)
     return data, adapted
 
 

From 08ef1e9b7541fb1465f73539400442ebaf929deb Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sat, 8 Dec 2012 14:29:36 -0600
Subject: [PATCH 21/31] Rebased to master, incorporated view_as_blocks, cleaned
 up interpolate method and example

---
 CONTRIBUTORS.txt                        |  4 --
 doc/examples/plot_equalize.py           |  7 ++--
 skimage/exposure/_adapthist.py          | 49 +++++++++++--------------
 skimage/exposure/tests/test_exposure.py |  5 ++-
 4 files changed, 30 insertions(+), 35 deletions(-)

diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index 5f796b3f..fd3ff235 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -123,10 +123,6 @@
 
 - Luis Pedro Coelho
   imread plugin
-=======
-  Polygon, circle and ellipse drawing functions
-  Adaptive thresholding
-  Implementation of Matlab's `regionprops`
 
 - Steven Silvester, Karel Zuiderveld
   Adaptive Histogram Equalization
diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py
index e09b5a07..14e4b9f2 100644
--- a/doc/examples/plot_equalize.py
+++ b/doc/examples/plot_equalize.py
@@ -18,7 +18,7 @@ that fall within the 2nd and 98th percentiles [2]_.
 
 """
 
-from skimage import data
+from skimage import data, img_as_ubyte
 from skimage.util.dtype import dtype_range
 from skimage import exposure
 
@@ -39,7 +39,7 @@ 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='stepfilled')
     ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
     ax_hist.set_xlabel('Pixel intensity')
 
@@ -63,10 +63,11 @@ img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))
 
 # Equalization
 img_eq = exposure.equalize(img)
+img_eq = img_as_ubyte(img_eq)
 
 # Adaptive Equalization
 img_adapteq = exposure.adapthist(img, clip_limit=0.03)
-
+img_adapteq = img_as_ubyte(img_adapteq)
 
 # Display results
 f, axes = plt.subplots(2, 4, figsize=(8, 4))
diff --git a/skimage/exposure/_adapthist.py b/skimage/exposure/_adapthist.py
index 091514b6..ff144700 100644
--- a/skimage/exposure/_adapthist.py
+++ b/skimage/exposure/_adapthist.py
@@ -25,6 +25,7 @@ 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 */
@@ -58,8 +59,6 @@ def adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, nbins=256):
     * 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 grayscale images, CLAHE is performed on one channel,
-      and a grayscale is returned
     * 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
@@ -127,7 +126,7 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
     if clip_limit == 1.0:
         return image  # is OK, immediately returns original image.
 
-    map_array = np.zeros((ntiles_x * ntiles_y, nbins), dtype=int)
+    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
@@ -146,21 +145,17 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
     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
-    ystart = 0
     for y in range(ntiles_y):
-        xstart = 0
         for x in range(ntiles_x):
-            sub_img = image[ystart: ystart + y_size,
-                            xstart: xstart + x_size]
+            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 * ntiles_x + x] = hist
-            xstart += x_size
-        ystart += y_size
+            map_array[y, x] = hist
     # Interpolate greylevel mappings to get CLAHE image
     ystart = 0
     for y in range(ntiles_y + 1):
@@ -190,11 +185,13 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
                 xstep = x_size
                 xL = x - 1
                 xR = xL + 1
-            mapLU = map_array[yU * ntiles_x + xL]
-            mapRU = map_array[yU * ntiles_x + xR]
-            mapLB = map_array[yB * ntiles_x + xL]
-            mapRB = map_array[yB * ntiles_x + xR]
-            interpolate(image, xstart, xstep, ystart, ystep,
+            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
@@ -280,7 +277,7 @@ def map_histogram(hist, min_val, max_val, n_pixels):
     return out.astype(int)
 
 
-def interpolate(image, xstart, xstep, ystart, ystep,
+def interpolate(image, xslice, yslice,
                 mapLU, mapRU, mapLB, mapRB, aLUT):
     '''Find the new grayscale level for a region using bilinear interpolation
 
@@ -288,10 +285,8 @@ def interpolate(image, xstart, xstep, ystart, ystep,
     ----------
     image : np.ndarray
         full image
-    xstart, xstop : int
-       indices of xslice
-    ystart, ystop : int
-        indices of yslice
+    xslice, yslice : array-like
+       indices of the region
     map* : np.ndarray
         mappings of greylevels from histograms
     aLUT : np.ndarray
@@ -309,18 +304,18 @@ def interpolate(image, xstart, xstep, ystart, ystep,
     This is done by a bilinear interpolation between four different
     mappings in order to eliminate boundary artifacts.
     '''
-    norm = xstep * ystep  # Normalization factor
-
+    norm = xslice.size * yslice.size  # Normalization factor
     # interpolation weight matrices
-    x_coef, y_coef = np.meshgrid(np.arange(xstep),
-                                 np.arange(ystep))
+    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
 
-    im_slice = image[ystart: ystart + ystep, xstart: xstart + xstep]
-    im_slice = aLUT[im_slice]
+    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)
-    image[ystart: ystart + ystep, xstart: xstart + xstep] = new
+    view[:, :] = new
+    return image
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index dc48a3d8..91b15b09 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -1,11 +1,13 @@
 import numpy as np
 from numpy.testing import assert_array_almost_equal as assert_close
 import skimage
+from skimage import io
 from skimage import data
 from skimage import exposure
 from skimage.color import rgb2gray
 from skimage.util.dtype import dtype_range
 
+io.use_plugin('qt')
 
 # Test histogram equalization
 # ===========================
@@ -81,7 +83,7 @@ def test_adapthist_scalar():
     img = skimage.img_as_ubyte(data.moon())
     adapted = exposure.adapthist(img, clip_limit=0.02)
     assert adapted.min() == 0
-    assert adapted.max() == (1<<16) - 1
+    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
@@ -90,6 +92,7 @@ def test_adapthist_scalar():
                         0.0410010)
     return img, adapted
 
+
 def test_adapthist_grayscale():
     '''Test a grayscale float image
     '''

From e91129de537088e616178ae12084ea380b50bac3 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sat, 8 Dec 2012 14:53:34 -0600
Subject: [PATCH 22/31] Renamed functions to equalize_hist/adapthist, with a
 deprecation warning for equalize

---
 skimage/exposure/__init__.py   | 6 +++---
 skimage/exposure/_adapthist.py | 3 ++-
 skimage/exposure/exposure.py   | 6 ++++++
 3 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/skimage/exposure/__init__.py b/skimage/exposure/__init__.py
index d495686a..ae75c982 100644
--- a/skimage/exposure/__init__.py
+++ b/skimage/exposure/__init__.py
@@ -1,3 +1,3 @@
-from .exposure import histogram, equalize, cumulative_distribution
-from .exposure import rescale_intensity
-from ._adapthist import adapthist
+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
index 57541de1..30d92d48 100644
--- a/skimage/exposure/_adapthist.py
+++ b/skimage/exposure/_adapthist.py
@@ -33,7 +33,8 @@ MAX_REG_Y = 16  # max. # contextual regions in y-direction */
 NR_OF_GREY = 1 << 14  # number of grayscale levels to use in CLAHE algorithm
 
 
-def adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, nbins=256):
+def equalize_adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01,
+                       nbins=256):
     '''Contrast Limited Adaptive Histogram Equalization
 
     Parameters
diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py
index 6ef7ff66..465b9182 100644
--- a/skimage/exposure/exposure.py
+++ b/skimage/exposure/exposure.py
@@ -4,6 +4,7 @@ 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',
@@ -78,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

From 71f505f71c3f6bbf2b9e8dd486586e5a17b9803a Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Sat, 8 Dec 2012 15:03:11 -0600
Subject: [PATCH 23/31] Used the new function names in example and test.

---
 doc/examples/plot_equalize.py           |  5 ++---
 skimage/exposure/tests/test_exposure.py | 14 ++++++--------
 2 files changed, 8 insertions(+), 11 deletions(-)

diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py
index 14e4b9f2..13ffa320 100644
--- a/doc/examples/plot_equalize.py
+++ b/doc/examples/plot_equalize.py
@@ -23,7 +23,6 @@ from skimage.util.dtype import dtype_range
 from skimage import exposure
 
 import matplotlib.pyplot as plt
-
 import numpy as np
 
 
@@ -62,11 +61,11 @@ 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)
 img_eq = img_as_ubyte(img_eq)
 
 # Adaptive Equalization
-img_adapteq = exposure.adapthist(img, clip_limit=0.03)
+img_adapteq = exposure.equalize_adapthist(img, clip_limit=0.03)
 img_adapteq = img_as_ubyte(img_adapteq)
 
 # Display results
diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 17206d1b..17633cf9 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -1,13 +1,11 @@
 import numpy as np
 from numpy.testing import assert_array_almost_equal as assert_close
 import skimage
-from skimage import io
 from skimage import data
 from skimage import exposure
 from skimage.color import rgb2gray
 from skimage.util.dtype import dtype_range
 
-io.use_plugin('qt')
 
 # Test histogram equalization
 # ===========================
@@ -18,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)
@@ -26,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)
@@ -81,8 +79,8 @@ def test_adapthist_scalar():
     '''Test a scalar uint8 image
     '''
     img = skimage.img_as_ubyte(data.moon())
-    adapted = exposure.adapthist(img, clip_limit=0.02)
-    assert adapted.min() == 0=
+    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))
@@ -99,7 +97,7 @@ def test_adapthist_grayscale():
     img = skimage.img_as_float(data.lena())
     img = rgb2gray(img)
     img = np.dstack((img, img, img))
-    adapted = exposure.adapthist(img, 10, 9, clip_limit=0.01,
+    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
@@ -112,7 +110,7 @@ def test_adapthist_color():
     '''Test a color uint16 image
     '''
     img = skimage.img_as_uint(data.lena())
-    adapted = exposure.adapthist(img, clip_limit=0.01)
+    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

From ac09463cf8b7571e65f51fc892aaa3cc0d9fde70 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Mon, 10 Dec 2012 20:19:23 -0600
Subject: [PATCH 24/31] Fixed comment header in test file.

---
 skimage/exposure/tests/test_exposure.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 17633cf9..1517eb11 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -72,8 +72,8 @@ def test_rescale_out_range():
     assert_close(out, [0, 63, 127])
 
 
-# Test rescale intensity
-# ======================
+# Test adaptive histogram equalization
+# ====================================
 
 def test_adapthist_scalar():
     '''Test a scalar uint8 image

From bb4f2b96dfe25e07558d6d1836c10e0563402416 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Mon, 10 Dec 2012 20:37:01 -0600
Subject: [PATCH 25/31] Modified test to pass with matplotlib plugin.

---
 skimage/exposure/tests/test_exposure.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 1517eb11..f0bdde37 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -5,7 +5,9 @@ from skimage import data
 from skimage import exposure
 from skimage.color import rgb2gray
 from skimage.util.dtype import dtype_range
+from skimage.io import use_plugin
 
+use_plugin('matplotlib')
 
 # Test histogram equalization
 # ===========================
@@ -101,7 +103,7 @@ def test_adapthist_grayscale():
                         nbins=128)
     assert_almost_equal = np.testing.assert_almost_equal
     assert img.shape == adapted.shape
-    assert_almost_equal(peak_snr(img, adapted), 77.5841845)
+    assert_almost_equal(peak_snr(img, adapted), 77.58419, 5)
     assert_almost_equal(norm_brightness_err(img, adapted), 0.0376632)
     return data, adapted
 

From 5444867fd09e19d39ffd200a89ccc43f82994a30 Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Mon, 10 Dec 2012 20:38:42 -0600
Subject: [PATCH 26/31] Removed explicit io_plugin from test

---
 skimage/exposure/tests/test_exposure.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index f0bdde37..30f51fc5 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -5,9 +5,7 @@ from skimage import data
 from skimage import exposure
 from skimage.color import rgb2gray
 from skimage.util.dtype import dtype_range
-from skimage.io import use_plugin
 
-use_plugin('matplotlib')
 
 # Test histogram equalization
 # ===========================

From 0e34c91e37f7ffbe401012a7f8a161bc95ce7fea Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Tue, 11 Dec 2012 18:08:27 -0600
Subject: [PATCH 27/31] Changed all test to use 3 significant digits.

---
 skimage/exposure/tests/test_exposure.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 30f51fc5..26607f61 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -85,9 +85,9 @@ def test_adapthist_scalar():
     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), 28.6262034)
+    assert_almost_equal(peak_snr(full_scale, adapted), 28.626, 3)
     assert_almost_equal(norm_brightness_err(full_scale, adapted),
-                        0.0410010)
+                        0.041, 3)
     return img, adapted
 
 
@@ -101,8 +101,8 @@ def test_adapthist_grayscale():
                         nbins=128)
     assert_almost_equal = np.testing.assert_almost_equal
     assert img.shape == adapted.shape
-    assert_almost_equal(peak_snr(img, adapted), 77.58419, 5)
-    assert_almost_equal(norm_brightness_err(img, adapted), 0.0376632)
+    assert_almost_equal(peak_snr(img, adapted), 77.584, 5)
+    assert_almost_equal(norm_brightness_err(img, adapted), 0.038, 3)
     return data, adapted
 
 
@@ -116,9 +116,9 @@ def test_adapthist_color():
     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.7167515)
+    assert_almost_equal(peak_snr(full_scale, adapted), 64.717, 3)
     assert_almost_equal(norm_brightness_err(full_scale, adapted),
-                        0.17922429)
+                        0.179, 3)
     return data, adapted
 
 

From 9adae8b0c231fe79a6599946b2ebc4f2aebc4acf Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Tue, 11 Dec 2012 18:57:24 -0600
Subject: [PATCH 28/31] Fixed last remaining sig figure problem.

---
 skimage/exposure/tests/test_exposure.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index 26607f61..b41f53fe 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -101,7 +101,7 @@ def test_adapthist_grayscale():
                         nbins=128)
     assert_almost_equal = np.testing.assert_almost_equal
     assert img.shape == adapted.shape
-    assert_almost_equal(peak_snr(img, adapted), 77.584, 5)
+    assert_almost_equal(peak_snr(img, adapted), 77.584, 3)
     assert_almost_equal(norm_brightness_err(img, adapted), 0.038, 3)
     return data, adapted
 

From 530fbbb2fc4cb739c8af075570246f2db3e0398a Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Tue, 27 Nov 2012 16:26:47 -0800
Subject: [PATCH 29/31] DOC: Minor clean-ups of equalization demo.

---
 doc/examples/plot_equalize.py | 19 +++++++++++--------
 1 file changed, 11 insertions(+), 8 deletions(-)

diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py
index 13ffa320..f57cea9c 100644
--- a/doc/examples/plot_equalize.py
+++ b/doc/examples/plot_equalize.py
@@ -18,7 +18,7 @@ that fall within the 2nd and 98th percentiles [2]_.
 
 """
 
-from skimage import data, img_as_ubyte
+from skimage import data, img_as_float
 from skimage.util.dtype import dtype_range
 from skimage import exposure
 
@@ -30,6 +30,7 @@ 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, histtype='stepfilled')
+    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
 
@@ -62,29 +63,31 @@ img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))
 
 # Equalization
 img_eq = exposure.equalize_hist(img)
-img_eq = img_as_ubyte(img_eq)
 
 # Adaptive Equalization
 img_adapteq = exposure.equalize_adapthist(img, clip_limit=0.03)
-img_adapteq = img_as_ubyte(img_adapteq)
 
 # Display results
 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)

From 7ebbdfd75fb0cd38153698815746974c12588742 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Tue, 11 Dec 2012 18:14:45 -0800
Subject: [PATCH 30/31] Minor code cleanup.

---
 skimage/exposure/_adapthist.py | 147 ++++++++++++++++-----------------
 1 file changed, 72 insertions(+), 75 deletions(-)

diff --git a/skimage/exposure/_adapthist.py b/skimage/exposure/_adapthist.py
index 30d92d48..667d2cb2 100644
--- a/skimage/exposure/_adapthist.py
+++ b/skimage/exposure/_adapthist.py
@@ -1,26 +1,18 @@
-'''
-Adapted code from the article
- * "Contrast Limited Adaptive Histogram Equalization"
- * by Karel Zuiderveld, karel@cv.ruu.nl
- * in "Graphics Gems IV", Academic Press, 1994
-=============
+"""
+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/
+http://tog.acm.org/resources/GraphicsGems/gems.html#gemsvi
 
-EULA: 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
+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.
- *
- *  Author: Karel Zuiderveld, Computer Vision Research Group,
- *	     Utrecht, The Netherlands (karel@cv.ruu.nl)
-'''
+"""
 import numpy as np
 import skimage
 from skimage import color
@@ -30,41 +22,37 @@ 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 = 1 << 14  # number of grayscale levels to use in CLAHE algorithm
+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
+    """Contrast Limited Adaptive Histogram Equalization.
 
     Parameters
     ----------
     image : array-like
-        original image
+        Input image.
     ntiles_x : int, optional
-        Tile regions in the X direction (2, 16)
+        Number of tile regions in the X direction.  Ranges between 2 and 16.
     ntiles_y : int, optional
-        Tile regions in the Y direction (2, 16)
+        Number of tile regions in the Y direction.  Ranges between 2 and 16.
     clip_limit : float: optional
-        Normalized cliplimit (higher values give more contrast)
+        Clipping limit, normalized between 0 and 1 (higher values give more
+        contrast).
     nbins : int, optional
-        Greybins for histogram ("dynamic range")
+        Number of gray bins for histogram ("dynamic range").
 
     Returns
     -------
-    out : np.ndarray
-        equalized image - grayscale images are uint16, color images are float
+    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.
-<<<<<<< HEAD
-=======
-    * For grayscale images, CLAHE is performed on one channel,
-      and a grayscale is returned
->>>>>>> 2e1729a9fbbc21fc0b04df8e68efbab9cfd6dada
     * 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
@@ -73,10 +61,9 @@ def equalize_adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01,
 
     References
     ----------
-    .. [1] http://tog.acm.org/resources/GraphicsGems/
+    .. [1] http://tog.acm.org/resources/GraphicsGems/gems.html#gemsvi
     .. [2] https://en.wikipedia.org/wiki/CLAHE#CLAHE
-    '''
-    # handle color images - CLAHE accepts scalar images only
+    """
     args = [None, ntiles_x, ntiles_y, clip_limit * nbins, nbins]
     if image.ndim > 2:
         lab_img = color.rgb2lab(skimage.img_as_float(image))
@@ -99,31 +86,32 @@ def equalize_adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01,
 
 
 def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
-    '''Contrast Limited Adaptive Histogram Equalization
+    """Contrast Limited Adaptive Histogram Equalization.
 
     Parameters
     ----------
     image : array-like
-        original image
+        Input image.
     ntiles_x : int, optional
-        Tile regions in the X direction (2, 16)
+        Number of tile regions in the X direction.  Ranges between 2 and 16.
     ntiles_y : int, optional
-        Tile regions in the Y direction (2, 16)
-    clip_limit : float: optional
-        Normalized cliplimit (higher values give more contrast)
+        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
-        Greybins for histogram ("dynamic range")
+        Number of gray bins for histogram ("dynamic range").
 
     Returns
     -------
-    out : np.ndarray
+    out : ndarray
+        Equalized image.
 
-    The number of "effective" greylevels in the output image is set by nbins;
+    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)
@@ -148,10 +136,12 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
             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):
@@ -162,6 +152,7 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
             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):
@@ -178,6 +169,7 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
             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
@@ -191,21 +183,26 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
                 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.
+    """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
@@ -213,16 +210,16 @@ def clip_histogram(hist, clip_limit):
 
     Parameters
     ----------
-    hist : np.ndarray
-        histogram array
+    hist : ndarray
+        Histogram array.
     clip_limit : int
-        maximum allowed bin count
+        Maximum allowed bin count.
 
     Returns
     -------
-    hist : np.ndarray
-        clipped histogram
-    '''
+    hist : ndarray
+        Clipped histogram.
+    """
     # calculate total number of excess pixels
     excess_mask = hist > clip_limit
     excess = hist[excess_mask]
@@ -253,28 +250,29 @@ def clip_histogram(hist, clip_limit):
             hist[under] += 1
             n_excess -= hist[under].size
             index += 1
+
     return hist
 
 
 def map_histogram(hist, min_val, max_val, n_pixels):
-    '''Calculates the equalized lookup table (mapping)
+    """Calculate the equalized lookup table (mapping).
 
     It does so by cumulating the input histogram.
 
-    hist : np.ndarray
-        clipped histogram
+    hist : ndarray
+        Clipped histogram.
     min_val : int
-        min value for mapping
+        Minimum value for mapping.
     max_val : int
-        max value for mapping
+        Maximum value for mapping.
     n_pixels : int
-        number of pixels in the region
+        Number of pixels in the region.
 
     Returns
     -------
-    out : np.ndarray
-       mapped intensity LUT
-    '''
+    out : ndarray
+       Mapped intensity LUT.
+    """
     out = np.cumsum(hist).astype(float)
     scale = ((float)(max_val - min_val)) / n_pixels
     out *= scale
@@ -285,23 +283,23 @@ def map_histogram(hist, min_val, max_val, n_pixels):
 
 def interpolate(image, xslice, yslice,
                 mapLU, mapRU, mapLB, mapRB, aLUT):
-    '''Find the new grayscale level for a region using bilinear interpolation
+    """Find the new grayscale level for a region using bilinear interpolation.
 
     Parameters
     ----------
-    image : np.ndarray
-        full image
+    image : ndarray
+        Full image.
     xslice, yslice : array-like
-       indices of the region
-    map* : np.ndarray
-        mappings of greylevels from histograms
-    aLUT : np.ndarray
-        maps grayscale levels in image to histogram levels
+       Indices of the region.
+    map* : ndarray
+        Mappings of greylevels from histograms.
+    aLUT : ndarray
+        Maps grayscale levels in image to histogram levels.
 
     Returns
     -------
-    out : np.ndarray
-        original image with the subregion replaced
+    out : ndarray
+        Original image with the subregion replaced.
 
     Note
     ----
@@ -309,7 +307,7 @@ def interpolate(image, xslice, yslice,
     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),
@@ -325,4 +323,3 @@ def interpolate(image, xslice, yslice,
            / norm)
     view[:, :] = new
     return image
-

From 73f5e9dccc375f930642a9166abb6d46a3d366fa Mon Sep 17 00:00:00 2001
From: Steven Silvester <steven.silvester@ieee.org>
Date: Tue, 11 Dec 2012 23:14:30 -0600
Subject: [PATCH 31/31] Fixed error in peak_snr, now uses float arithmetic
 throughout.

---
 skimage/exposure/tests/test_exposure.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py
index b41f53fe..258b4224 100644
--- a/skimage/exposure/tests/test_exposure.py
+++ b/skimage/exposure/tests/test_exposure.py
@@ -85,7 +85,7 @@ def test_adapthist_scalar():
     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), 28.626, 3)
+    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
@@ -137,8 +137,11 @@ def peak_snr(img1, img2):
     '''
     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)