wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-02 17:40:31 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1547 from blink1073/adapthist-window

Browse Source 
Update Adapthist to use more intelligent tiles
This commit is contained in:
Juan Nunez-Iglesias
2015-07-14 14:08:17 -04:00
parent 2690fef1db 80fe225fc8
commit 7124716d24
6 changed files with 107 additions and 133 deletions
Download Patch File Download Diff File Expand all files Collapse all files
TODO.txt
+7
View File
@@ -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
doc/source/api_changes.txt
+2
View File
@@ -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.
skimage/exposure/_adapthist.py
+87 -95
View File
@@ -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
skimage/exposure/tests/test_exposure.py
+10 -9
View File
@@ -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]
skimage/util/shape.py
+1 -13
View File
@@ -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")
skimage/util/tests/test_shape.py
-16
View File
@@ -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()