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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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/10] 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