diff --git a/TODO.txt b/TODO.txt index 98e8653b..a96b3731 100644 --- a/TODO.txt +++ b/TODO.txt @@ -1,5 +1,10 @@ Remember to list any API changes below in `doc/source/api_changes.txt`. +Version 0.14 +------------ +* Remove deprecated ``ntiles_*` kwargs in ``equalize_adapthist``. + + Version 0.13 ------------ * Remove deprecated `None` defaults for `skimage.exposure.rescale_intensity` @@ -12,6 +17,8 @@ Version 0.13 `hprewitt`, `vprewitt`, `roberts_positive_diagonal`, `roberts_negative_diagonal` in `skimage/filters/edges.py` + + Version 0.12 ------------ * Change `label` to mark background as 0, not -1, which is consistent with diff --git a/doc/source/api_changes.txt b/doc/source/api_changes.txt index a8edfc3d..e87445b3 100644 --- a/doc/source/api_changes.txt +++ b/doc/source/api_changes.txt @@ -1,5 +1,7 @@ Version 0.12 ------------ +- ``equalize_adapthist`` now takes a ``kernel_size`` keyword argument, replacing + the ``ntiles_*`` arguments. - The functions ``blob_dog``, ``blob_log`` and ``blob_doh`` now return float arrays instead of integer arrays. diff --git a/skimage/exposure/_adapthist.py b/skimage/exposure/_adapthist.py index b71916f5..343f71b3 100644 --- a/skimage/exposure/_adapthist.py +++ b/skimage/exposure/_adapthist.py @@ -13,21 +13,20 @@ 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. """ +from __future__ import division +import numbers import numpy as np from .. import img_as_float, img_as_uint from ..color.adapt_rgb import adapt_rgb, hsv_value from ..exposure import rescale_intensity -from ..util import view_as_blocks, pad +from .._shared.utils import skimage_deprecation, warnings - -MAX_REG_X = 16 # max. # contextual regions in x-direction */ -MAX_REG_Y = 16 # max. # contextual regions in y-direction */ NR_OF_GREY = 2 ** 14 # number of grayscale levels to use in CLAHE algorithm @adapt_rgb(hsv_value) def equalize_adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, - nbins=256): + nbins=256, kernel_size=None): """Contrast Limited Adaptive Histogram Equalization (CLAHE). An algorithm for local contrast enhancement, that uses histograms computed @@ -38,10 +37,14 @@ def equalize_adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, ---------- image : array-like Input image. - ntiles_x : int, optional - Number of tile regions in the X direction. Ranges between 1 and 16. - ntiles_y : int, optional - Number of tile regions in the Y direction. Ranges between 1 and 16. + kernel_size: integer or 2-tuple + Defines the shape of contextual regions used in the algorithm. + If an integer is given, the shape will be a square of + sidelength given by this value. + ntiles_x : int, optional (deprecated in favor of ``kernel_size``) + Number of tile regions in the X direction (horizontal). + ntiles_y : int, optional (deprecated in favor of ``kernel_size``) + Number of tile regions in the Y direction (vertical). clip_limit : float: optional Clipping limit, normalized between 0 and 1 (higher values give more contrast). @@ -64,10 +67,6 @@ def equalize_adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, - 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 ---------- @@ -76,23 +75,35 @@ def equalize_adapthist(image, ntiles_x=8, ntiles_y=8, clip_limit=0.01, """ image = 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 + + if kernel_size is None: + warnings.warn('`ntiles_*` have been deprecated in favor of ' + '`kernel_size`. The `ntiles_*` keyword arguments ' + 'will be removed in v0.14', skimage_deprecation) + ntiles_x = ntiles_x or 8 + ntiles_y = ntiles_y or 8 + kernel_size = (np.round(image.shape[0] / ntiles_y), + np.round(image.shape[1] / ntiles_x)) + + if isinstance(kernel_size, numbers.Number): + kernel_size = (kernel_size, kernel_size) + + kernel_size = [int(k) for k in kernel_size] + + image = _clahe(image, kernel_size, clip_limit * nbins, nbins) image = img_as_float(image) return rescale_intensity(image) -def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128): +def _clahe(image, kernel_size, clip_limit, nbins=128): """Contrast Limited Adaptive Histogram Equalization. Parameters ---------- image : array-like Input image. - ntiles_x : int, optional - Number of tile regions in the X direction. Ranges between 2 and 16. - ntiles_y : int, optional - Number of tile regions in the Y direction. Ranges between 2 and 16. + kernel_size: 2-tuple + Defines the shape of contextual regions used in the algorithm. clip_limit : float, optional Normalized clipping limit (higher values give more contrast). nbins : int, optional @@ -109,106 +120,87 @@ def _clahe(image, ntiles_x, ntiles_y, clip_limit, nbins=128): 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) 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 + nr = int(np.ceil(image.shape[0] / kernel_size[0])) + nc = int(np.ceil(image.shape[1] / kernel_size[1])) - # make the tile size divisible by 2 - h_inner -= h_inner % (2 * ntiles_y) - w_inner -= w_inner % (2 * ntiles_x) - - 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 = pad(image, ((0, h_pad), (0, w_pad)), mode='reflect') - h_inner, w_inner = image.shape + row_step = int(np.floor(image.shape[0] / nr)) + col_step = int(np.floor(image.shape[1] / nc)) 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) - n_pixels = width * height - - if clip_limit > 0.0: # Calculate actual cliplimit - 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) + map_array = np.zeros((nr, nc, nbins), dtype=int) # 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] + for r in range(nr): + for c in range(nc): + sub_img = image[r * row_step: (r + 1) * row_step, + c * col_step: (c + 1) * col_step] + + if clip_limit > 0.0: # Calculate actual cliplimit + clim = int(clip_limit * sub_img.size / nbins) + if clim < 1: + clim = 1 + else: + clim = NR_OF_GREY # Large value, do not clip (AHE) + 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) - hist = map_histogram(hist, 0, NR_OF_GREY - 1, n_pixels) - map_array[y, x] = hist + hist = clip_histogram(hist, clim) + hist = map_histogram(hist, 0, NR_OF_GREY - 1, sub_img.size) + map_array[r, c] = hist # 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 = height / 2.0 - yU = 0 - yB = 0 - elif y == ntiles_y: # special case: bottom row - ystep = height / 2.0 - yU = ntiles_y - 1 - yB = yU + rstart = 0 + for r in range(nr + 1): + cstart = 0 + if r == 0: # special case: top row + r_offset = row_step / 2.0 + rU = 0 + rB = 0 + elif r == nr: # special case: bottom row + r_offset = row_step / 2.0 + rU = nr - 1 + rB = rU else: # default values - ystep = height - yU = y - 1 - yB = yB + 1 + r_offset = row_step + rU = r - 1 + rB = rB + 1 - for x in range(ntiles_x + 1): - if x == 0: # special case: left column - xstep = width / 2.0 - xL = 0 - xR = 0 - elif x == ntiles_x: # special case: right column - xstep = width / 2.0 - xL = ntiles_x - 1 - xR = xL + for c in range(nc + 1): + if c == 0: # special case: left column + c_offset = col_step / 2.0 + cL = 0 + cR = 0 + elif c == nc: # special case: right column + c_offset = col_step / 2.0 + cL = nc - 1 + cR = cL else: # default values - xstep = width - xL = x - 1 - xR = xL + 1 + c_offset = col_step + cL = c - 1 + cR = cL + 1 - mapLU = map_array[yU, xL] - mapRU = map_array[yU, xR] - mapLB = map_array[yB, xL] - mapRB = map_array[yB, xR] + mapLU = map_array[rU, cL] + mapRU = map_array[rU, cR] + mapLB = map_array[rB, cL] + mapRB = map_array[rB, cR] - xslice = np.arange(xstart, xstart + xstep) - yslice = np.arange(ystart, ystart + ystep) - interpolate(image, xslice, yslice, + cslice = np.arange(cstart, cstart + c_offset) + rslice = np.arange(rstart, rstart + r_offset) + + interpolate(image, cslice, rslice, mapLU, mapRU, mapLB, mapRB, lut) - xstart += xstep # set pointer on next matrix */ + cstart += c_offset # set pointer on next matrix */ - ystart += ystep - - if image.shape != orig_shape: - image = image[:orig_shape[0], :orig_shape[1]] + rstart += r_offset return image diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py index 265ce246..8fb6d621 100644 --- a/skimage/exposure/tests/test_exposure.py +++ b/skimage/exposure/tests/test_exposure.py @@ -192,7 +192,7 @@ def test_adapthist_scalar(): """Test a scalar uint8 image """ img = skimage.img_as_ubyte(data.moon()) - adapted = exposure.equalize_adapthist(img, clip_limit=0.02) + adapted = exposure.equalize_adapthist(img, kernel_size=64, clip_limit=0.02) assert adapted.min() == 0.0 assert adapted.max() == 1.0 assert img.shape == adapted.shape @@ -211,13 +211,14 @@ def test_adapthist_grayscale(): img = skimage.img_as_float(data.astronaut()) img = rgb2gray(img) img = np.dstack((img, img, img)) - with expected_warnings(['precision loss|non-contiguous input']): - adapted = exposure.equalize_adapthist(img, 10, 9, clip_limit=0.01, - nbins=128) - assert_almost_equal = np.testing.assert_almost_equal + with expected_warnings(['precision loss|non-contiguous input', + 'deprecated']): + adapted_old = exposure.equalize_adapthist(img, 10, 9, clip_limit=0.01, + nbins=128) + adapted = exposure.equalize_adapthist(img, kernel_size=(57, 51), clip_limit=0.01, nbins=128) assert img.shape == adapted.shape - assert_almost_equal(peak_snr(img, adapted), 97.6876, 3) - assert_almost_equal(norm_brightness_err(img, adapted), 0.0591, 3) + assert_almost_equal(peak_snr(img, adapted), 101.750, 3) + assert_almost_equal(norm_brightness_err(img, adapted), 0.0540, 3) return data, adapted @@ -229,7 +230,7 @@ def test_adapthist_color(): warnings.simplefilter('always') hist, bin_centers = exposure.histogram(img) assert len(w) > 0 - with expected_warnings(['precision loss']): + with expected_warnings(['precision loss', 'deprecated']): adapted = exposure.equalize_adapthist(img, clip_limit=0.01) assert_almost_equal = np.testing.assert_almost_equal @@ -248,7 +249,7 @@ def test_adapthist_alpha(): img = skimage.img_as_float(data.astronaut()) alpha = np.ones((img.shape[0], img.shape[1]), dtype=float) img = np.dstack((img, alpha)) - with expected_warnings(['precision loss']): + with expected_warnings(['precision loss', 'deprecated']): adapted = exposure.equalize_adapthist(img) assert adapted.shape != img.shape img = img[:, :, :3] diff --git a/skimage/util/shape.py b/skimage/util/shape.py index 48e9a5ef..b786ec1c 100644 --- a/skimage/util/shape.py +++ b/skimage/util/shape.py @@ -104,7 +104,7 @@ def view_as_blocks(arr_in, block_shape): return arr_out -def view_as_windows(arr_in, window_shape, step=1, optimal_step=False): +def view_as_windows(arr_in, window_shape, step=1): """Rolling window view of the input n-dimensional array. Windows are overlapping views of the input array, with adjacent windows @@ -122,9 +122,6 @@ def view_as_windows(arr_in, window_shape, step=1, optimal_step=False): step : integer or tuple of length arr_in.ndim Indicates step size at which extraction shall be performed. If integer is given, then the step is uniform in all dimensions. - optimal_step: bool, optional - When True, selects a ``step`` that will give full coverage of - ``arr_in`` with minimal overlap. Returns ------- @@ -229,15 +226,6 @@ def view_as_windows(arr_in, window_shape, step=1, optimal_step=False): if not (len(window_shape) == ndim): raise ValueError("`window_shape` is incompatible with `arr_in.shape`") - if optimal_step: - rem = np.array(arr_in.shape) - np.array(window_shape) - step = list(window_shape) - - for (ind, size) in enumerate(window_shape): - ns = int(np.ceil(arr_in.shape[ind] / size)) - while step[ind] * (ns - 1) > rem[ind]: - step[ind] -= 1 - if isinstance(step, numbers.Number): if step < 1: raise ValueError("`step` must be >= 1") diff --git a/skimage/util/tests/test_shape.py b/skimage/util/tests/test_shape.py index 5c59edb0..7a6a785d 100644 --- a/skimage/util/tests/test_shape.py +++ b/skimage/util/tests/test_shape.py @@ -176,21 +176,5 @@ def test_view_as_windows_step_tuple(): [22, 23]]]]) -def test_view_as_windows_optimal_step(): - A = np.arange(24).reshape((6, 4)) - B = view_as_windows(A, (3, 2), optimal_step=True) - assert B.shape == (2, 2, 3, 2) - assert B.size == A.size - - A = np.arange(512 * 512).reshape((512, 512)) - B = view_as_windows(A, (10, 10), optimal_step=True) - assert B.shape == (56, 56, 10, 10) - assert B.size >= A.size - - C = view_as_windows(A, (11, 9), optimal_step=True) - assert C.shape == (51, 63, 11, 9) - assert C.size >= A.size - - if __name__ == '__main__': np.testing.run_module_suite()