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