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