wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-23 13:10:18 +08:00
Code Issues Packages Projects Releases Wiki Activity

Translated clahe C code to python, added to list of contributors.

Browse Source
This commit is contained in:
Steven Silvester
2012-12-07 17:25:04 -06:00
parent 269929333b
commit 94af01cdbd
8 changed files with 360 additions and 498 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CONTRIBUTORS.txt
+7
View File
@@ -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
skimage/exposure/__init__.py
+2 -1
View File
@@ -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
skimage/exposure/_adapthist.py
+337
View File
@@ -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
skimage/exposure/_adapthist.pyx
-40
View File
@@ -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))
skimage/exposure/clahe.h
-322
View File
@@ -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);
}
}
}
}
skimage/exposure/exposure.py
+1 -89
View File
@@ -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
skimage/exposure/setup.py
-31
View File
@@ -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())
)
skimage/exposure/tests/test_exposure.py
+13 -15
View File
@@ -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