wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-06-30 03:13:20 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #893 from blink1073/fix-adapthist-757

Browse Source 
Add option for Adaptive Histogram border pixel handling.
This commit is contained in:
Tony S Yu
2014-07-15 22:50:50 -05:00
parent 74988b88d5 70a3760bca
commit 7cf80e6ae5
2 changed files with 101 additions and 60 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/exposure/_adapthist.py
+68 -48
View File
@@ -50,37 +50,38 @@ def equalize_adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01,
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 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 to HSV color space
- The CLAHE algorithm is run on the V (Value) channel
- The image is converted back to RGB space and returned
* For RGBA images, the original alpha channel is removed.
* The CLAHE algorithm relies on image blocks of equal size. This may
result in extra border pixels that would not be handled. In that case,
we pad the image with a repeat of the border pixels, apply the
algorithm, and then trim the image to original size.
References
----------
.. [1] http://tog.acm.org/resources/GraphicsGems/gems.html#gemsvi
.. [2] https://en.wikipedia.org/wiki/CLAHE#CLAHE
"""
args = [None, ntiles_x, ntiles_y, clip_limit * nbins, nbins]
if image.ndim > 2:
lab_img = color.rgb2lab(skimage.img_as_float(image))
l_chan = lab_img[:, :, 0]
l_chan /= np.max(np.abs(l_chan))
l_chan = skimage.img_as_uint(l_chan)
args[0] = rescale_intensity(l_chan, out_range=(0, NR_OF_GREY - 1))
new_l = _clahe(*args).astype(float)
new_l = rescale_intensity(new_l, out_range=(0, 100))
lab_img[:new_l.shape[0], :new_l.shape[1], 0] = new_l
image = color.lab2rgb(lab_img)
image = rescale_intensity(image, out_range=(0, 1))
ndim = image.ndim
if ndim == 3:
if image.shape[2] == 4:
image = image[:, :, :3]
image = skimage.img_as_float(image)
image = rescale_intensity(image)
hsv_img = color.rgb2hsv(image)
image = hsv_img[:, :, 2].copy()
image = skimage.img_as_uint(image)
image = rescale_intensity(image, out_range=(0, NR_OF_GREY - 1))
out = _clahe(image, ntiles_x, ntiles_y, clip_limit * nbins, nbins)
image[:out.shape[0], :out.shape[1]] = out
image = skimage.img_as_float(image)
if ndim == 3:
hsv_img[:, :, 2] = rescale_intensity(image)
image = color.hsv2rgb(hsv_img)
else:
image = skimage.img_as_uint(image)
args[0] = rescale_intensity(image, out_range=(0, NR_OF_GREY - 1))
out = _clahe(*args)
image[:out.shape[0], :out.shape[1]] = out
image = rescale_intensity(image)
return image
@@ -114,39 +115,55 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
"""
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)
ntiles_y = max(ntiles_y, 2)
ntiles_x = max(ntiles_x, 2)
if clip_limit == 1.0:
return image # is OK, immediately returns original image.
h_inner = image.shape[0] - image.shape[0] % ntiles_y
w_inner = image.shape[1] - image.shape[1] % ntiles_x
# make the tile size divisible by 2
while h_inner % (2 * ntiles_y):
h_inner -= 1
while w_inner % (2 * ntiles_x):
w_inner -= 1
orig_shape = image.shape
width = w_inner // ntiles_x # Actual size of contextual regions
height = h_inner // ntiles_y
if h_inner != image.shape[0]:
ntiles_y += 1
if w_inner != image.shape[1]:
ntiles_x += 1
if h_inner != image.shape[1] or w_inner != image.shape[0]:
h_pad = height * ntiles_y - image.shape[0]
w_pad = width * ntiles_x - image.shape[1]
image = np.pad(image, ((0, h_pad), (0, w_pad)), 'reflect')
h_inner, w_inner = image.shape
bin_size = 1 + NR_OF_GREY / nbins
lut = np.arange(NR_OF_GREY)
lut //= bin_size
img_blocks = view_as_blocks(image, (height, width))
map_array = np.zeros((ntiles_y, ntiles_x, nbins), dtype=int)
y_res = image.shape[0] - image.shape[0] % ntiles_y
x_res = image.shape[1] - image.shape[1] % ntiles_x
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
n_pixels = width * height
if clip_limit > 0.0: # Calculate actual cliplimit
clip_limit = int(clip_limit * (x_size * y_size) / nbins)
clip_limit = int(clip_limit * (width * height) / 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
img_blocks = view_as_blocks(image, (y_size, x_size))
# Calculate greylevel mappings for each contextual region
for y in range(ntiles_y):
for x in range(ntiles_x):
sub_img = img_blocks[y, x]
hist = aLUT[sub_img.ravel()]
hist = lut[sub_img.ravel()]
hist = np.bincount(hist)
hist = np.append(hist, np.zeros(nbins - hist.size, dtype=int))
hist = clip_histogram(hist, clip_limit)
@@ -158,29 +175,29 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
for y in range(ntiles_y + 1):
xstart = 0
if y == 0: # special case: top row
ystep = y_size / 2.0
ystep = height / 2.0
yU = 0
yB = 0
elif y == ntiles_y: # special case: bottom row
ystep = y_size / 2.0
ystep = height / 2.0
yU = ntiles_y - 1
yB = yU
else: # default values
ystep = y_size
ystep = height
yU = y - 1
yB = yB + 1
for x in range(ntiles_x + 1):
if x == 0: # special case: left column
xstep = x_size / 2.0
xstep = width / 2.0
xL = 0
xR = 0
elif x == ntiles_x: # special case: right column
xstep = x_size / 2.0
xstep = width / 2.0
xL = ntiles_x - 1
xR = xL
else: # default values
xstep = x_size
xstep = width
xL = x - 1
xR = xL + 1
@@ -192,12 +209,15 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128):
xslice = np.arange(xstart, xstart + xstep)
yslice = np.arange(ystart, ystart + ystep)
interpolate(image, xslice, yslice,
mapLU, mapRU, mapLB, mapRB, aLUT)
mapLU, mapRU, mapLB, mapRB, lut)
xstart += xstep # set pointer on next matrix */
ystart += ystep
if image.shape != orig_shape:
image = image[:orig_shape[0], :orig_shape[1]]
return image
@@ -284,7 +304,7 @@ def map_histogram(hist, min_val, max_val, n_pixels):
def interpolate(image, xslice, yslice,
mapLU, mapRU, mapLB, mapRB, aLUT):
mapLU, mapRU, mapLB, mapRB, lut):
"""Find the new grayscale level for a region using bilinear interpolation.
Parameters
@@ -295,7 +315,7 @@ def interpolate(image, xslice, yslice,
Indices of the region.
map* : ndarray
Mappings of greylevels from histograms.
aLUT : ndarray
lut : ndarray
Maps grayscale levels in image to histogram levels.
Returns
@@ -317,7 +337,7 @@ def interpolate(image, xslice, yslice,
view = image[int(yslice[0]):int(yslice[-1] + 1),
int(xslice[0]):int(xslice[-1] + 1)]
im_slice = aLUT[view]
im_slice = lut[view]
new = ((y_inv_coef * (x_inv_coef * mapLU[im_slice]
+ x_coef * mapRU[im_slice])
+ y_coef * (x_inv_coef * mapLB[im_slice]
skimage/exposure/tests/test_exposure.py
+33 -12
View File
@@ -11,10 +11,13 @@ from skimage.exposure.exposure import intensity_range
from skimage.color import rgb2gray
from skimage.util.dtype import dtype_range
import matplotlib.pyplot as plt
# Test histogram equalization
# ===========================
np.random.seed(0)
# squeeze image intensities to lower image contrast
test_img = skimage.img_as_float(data.camera())
test_img = exposure.rescale_intensity(test_img / 5. + 100)
@@ -147,13 +150,13 @@ def test_adapthist_scalar():
'''
img = skimage.img_as_ubyte(data.moon())
adapted = exposure.equalize_adapthist(img, clip_limit=0.02)
assert adapted.min() == 0
assert adapted.max() == (1 << 16) - 1
assert adapted.min() == 0.0
assert adapted.max() == 1.0
assert img.shape == adapted.shape
full_scale = skimage.exposure.rescale_intensity(skimage.img_as_uint(img))
full_scale = skimage.exposure.rescale_intensity(skimage.img_as_float(img))
assert_almost_equal = np.testing.assert_almost_equal
assert_almost_equal(peak_snr(full_scale, adapted), 101.231, 3)
assert_almost_equal(peak_snr(full_scale, adapted), 101.2295, 3)
assert_almost_equal(norm_brightness_err(full_scale, adapted),
0.041, 3)
return img, adapted
@@ -169,8 +172,8 @@ def test_adapthist_grayscale():
nbins=128)
assert_almost_equal = np.testing.assert_almost_equal
assert img.shape == adapted.shape
assert_almost_equal(peak_snr(img, adapted), 97.531, 3)
assert_almost_equal(norm_brightness_err(img, adapted), 0.0313, 3)
assert_almost_equal(peak_snr(img, adapted), 104.307, 3)
assert_almost_equal(norm_brightness_err(img, adapted), 0.0265, 3)
return data, adapted
@@ -183,17 +186,34 @@ def test_adapthist_color():
hist, bin_centers = exposure.histogram(img)
assert len(w) > 0
adapted = exposure.equalize_adapthist(img, clip_limit=0.01)
assert_almost_equal = np.testing.assert_almost_equal
assert adapted.min() == 0
assert adapted.max() == 1.0
assert img.shape == adapted.shape
full_scale = skimage.exposure.rescale_intensity(img)
assert_almost_equal(peak_snr(full_scale, adapted), 102.940, 3)
assert_almost_equal(peak_snr(full_scale, adapted), 105.50517, 3)
assert_almost_equal(norm_brightness_err(full_scale, adapted),
0.0110, 3)
0.0544, 3)
return data, adapted
def test_adapthist_alpha():
'''Test an RGBA color image
'''
img = skimage.img_as_float(data.lena())
alpha = np.ones((img.shape[0], img.shape[1]), dtype=float)
img = np.dstack((img, alpha))
adapted = exposure.equalize_adapthist(img)
assert adapted.shape != img.shape
img = img[:, :, :3]
full_scale = skimage.exposure.rescale_intensity(img)
assert img.shape == adapted.shape
assert_almost_equal = np.testing.assert_almost_equal
assert_almost_equal(peak_snr(full_scale, adapted), 105.50198, 3)
assert_almost_equal(norm_brightness_err(full_scale, adapted), 0.0544, 3)
def peak_snr(img1, img2):
'''Peak signal to noise ratio of two images
@@ -236,10 +256,6 @@ def norm_brightness_err(img1, img2):
return nbe
if __name__ == '__main__':
from numpy import testing
testing.run_module_suite()
# Test Gamma Correction
# =====================
@@ -409,3 +425,8 @@ def test_adjust_inv_sigmoid_cutoff_half():
def test_negative():
image = np.arange(-10, 245, 4).reshape(8, 8).astype(np.double)
assert_raises(ValueError, exposure.adjust_gamma, image)
if __name__ == '__main__':
from numpy import testing
testing.run_module_suite()