From 08ef1e9b7541fb1465f73539400442ebaf929deb Mon Sep 17 00:00:00 2001 From: Steven Silvester Date: Sat, 8 Dec 2012 14:29:36 -0600 Subject: [PATCH] Rebased to master, incorporated view_as_blocks, cleaned up interpolate method and example --- CONTRIBUTORS.txt | 4 -- doc/examples/plot_equalize.py | 7 ++-- skimage/exposure/_adapthist.py | 49 +++++++++++-------------- skimage/exposure/tests/test_exposure.py | 5 ++- 4 files changed, 30 insertions(+), 35 deletions(-) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 5f796b3f..fd3ff235 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -123,10 +123,6 @@ - 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/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py index e09b5a07..14e4b9f2 100644 --- a/doc/examples/plot_equalize.py +++ b/doc/examples/plot_equalize.py @@ -18,7 +18,7 @@ that fall within the 2nd and 98th percentiles [2]_. """ -from skimage import data +from skimage import data, img_as_ubyte from skimage.util.dtype import dtype_range from skimage import exposure @@ -39,7 +39,7 @@ def plot_img_and_hist(img, axes, bins=256): ax_img.set_axis_off() # Display histogram - ax_hist.hist(img.ravel(), bins=bins) + ax_hist.hist(img.ravel(), bins=bins, histtype='stepfilled') ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0)) ax_hist.set_xlabel('Pixel intensity') @@ -63,10 +63,11 @@ img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98)) # Equalization img_eq = exposure.equalize(img) +img_eq = img_as_ubyte(img_eq) # Adaptive Equalization img_adapteq = exposure.adapthist(img, clip_limit=0.03) - +img_adapteq = img_as_ubyte(img_adapteq) # Display results f, axes = plt.subplots(2, 4, figsize=(8, 4)) diff --git a/skimage/exposure/_adapthist.py b/skimage/exposure/_adapthist.py index 091514b6..ff144700 100644 --- a/skimage/exposure/_adapthist.py +++ b/skimage/exposure/_adapthist.py @@ -25,6 +25,7 @@ import numpy as np import skimage from skimage import color from skimage.exposure import rescale_intensity +from skimage.util import view_as_blocks MAX_REG_X = 16 # max. # contextual regions in x-direction */ @@ -58,8 +59,6 @@ def adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, nbins=256): * 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 @@ -127,7 +126,7 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128): if clip_limit == 1.0: return image # is OK, immediately returns original image. - map_array = np.zeros((ntiles_x * ntiles_y, nbins), dtype=int) + 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 @@ -146,21 +145,17 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128): 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 - 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] + sub_img = img_blocks[y, x] 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 + map_array[y, x] = hist # Interpolate greylevel mappings to get CLAHE image ystart = 0 for y in range(ntiles_y + 1): @@ -190,11 +185,13 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128): 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 = map_array[yU, xL] + mapRU = map_array[yU, xR] + mapLB = map_array[yB, xL] + mapRB = map_array[yB, xR] + xslice = np.arange(xstart, xstart + xstep) + yslice = np.arange(ystart, ystart + ystep) + interpolate(image, xslice, yslice, mapLU, mapRU, mapLB, mapRB, aLUT) xstart += xstep # set pointer on next matrix */ ystart += ystep @@ -280,7 +277,7 @@ def map_histogram(hist, min_val, max_val, n_pixels): return out.astype(int) -def interpolate(image, xstart, xstep, ystart, ystep, +def interpolate(image, xslice, yslice, mapLU, mapRU, mapLB, mapRB, aLUT): '''Find the new grayscale level for a region using bilinear interpolation @@ -288,10 +285,8 @@ def interpolate(image, xstart, xstep, ystart, ystep, ---------- image : np.ndarray full image - xstart, xstop : int - indices of xslice - ystart, ystop : int - indices of yslice + xslice, yslice : array-like + indices of the region map* : np.ndarray mappings of greylevels from histograms aLUT : np.ndarray @@ -309,18 +304,18 @@ def interpolate(image, xstart, xstep, ystart, ystep, This is done by a bilinear interpolation between four different mappings in order to eliminate boundary artifacts. ''' - norm = xstep * ystep # Normalization factor - + norm = xslice.size * yslice.size # Normalization factor # interpolation weight matrices - x_coef, y_coef = np.meshgrid(np.arange(xstep), - np.arange(ystep)) + x_coef, y_coef = np.meshgrid(np.arange(xslice.size), + np.arange(yslice.size)) 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] + view = image[yslice[0]: yslice[-1] + 1, xslice[0]: xslice[-1] + 1] + im_slice = aLUT[view] 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 + view[:, :] = new + return image diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py index dc48a3d8..91b15b09 100644 --- a/skimage/exposure/tests/test_exposure.py +++ b/skimage/exposure/tests/test_exposure.py @@ -1,11 +1,13 @@ import numpy as np from numpy.testing import assert_array_almost_equal as assert_close import skimage +from skimage import io from skimage import data from skimage import exposure from skimage.color import rgb2gray from skimage.util.dtype import dtype_range +io.use_plugin('qt') # Test histogram equalization # =========================== @@ -81,7 +83,7 @@ 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() == (1<<16) - 1 + assert adapted.max() == (1 << 16) - 1 assert img.shape == adapted.shape full_scale = skimage.exposure.rescale_intensity(skimage.img_as_uint(img)) assert_almost_equal = np.testing.assert_almost_equal @@ -90,6 +92,7 @@ def test_adapthist_scalar(): 0.0410010) return img, adapted + def test_adapthist_grayscale(): '''Test a grayscale float image '''