diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py index 9cb30dd7..8e5609bc 100644 --- a/doc/examples/plot_equalize.py +++ b/doc/examples/plot_equalize.py @@ -8,10 +8,14 @@ histogram equalization. Histogram equalization enhances contrast by "spreading out the most frequent intensity values" in an image [1]_. The equalized image has a roughly linear cumulative distribution function, as shown in this example. +For comparison, this example also shows an image after its intensity levels are +uniformly stretched. + .. [1] http://en.wikipedia.org/wiki/Histogram_equalization """ import matplotlib.pyplot as plt +from matplotlib import ticker from skimage import data from skimage.util.dtype import dtype_range @@ -22,35 +26,49 @@ def plot_hist(img, bins=256): """Plot histogram and cumulative histogram for image""" img_cdf, bins = exposure.cumulative_distribution(img, bins) plt.hist(img.ravel(), bins=bins) - plt.ylabel('Number of pixels') plt.xlabel('Pixel intensiy') + # Shorten y-tick labels using scientific notation + y_formatter = ticker.ScalarFormatter(useOffset=True) + y_formatter.set_powerlimits((0, 0)) # force use of scientific notation + ax = plt.gca() + ax.yaxis.set_major_formatter(y_formatter) ax_cdf = plt.twinx() ax_cdf.plot(bins, img_cdf, 'r') xmin, xmax = dtype_range[img.dtype.type] plt.xlim(xmin, xmax) - ax_cdf.set_ylabel('Fraction of total intensity') img_orig = data.camera() # squeeze image intensities to lower image contrast img = img_orig / 5 + 100 +img_rescale = exposure.rescale_intensity(img) img_eq = exposure.equalize(img) -plt.subplot(2, 2, 1) +plt.subplot(2, 3, 1) plt.imshow(img, cmap=plt.cm.gray, vmin=0, vmax=255) -plt.title('Low contrast input image') +plt.title('Low contrast image') plt.axis('off') -plt.subplot(2, 2, 2) +plt.subplot(2, 3, 4) +plt.ylabel('Number of pixels') plot_hist(img) -plt.subplot(2, 2, 3) -plt.imshow(img_eq, cmap=plt.cm.gray, vmin=0, vmax=1) -plt.title('After\nhistogram equalization') +plt.subplot(2, 3, 2) +plt.imshow(img_rescale, cmap=plt.cm.gray, vmin=0, vmax=255) +plt.title('Rescale intensities') plt.axis('off') -plt.subplot(2, 2, 4) -plot_hist(img_eq) +plt.subplot(2, 3, 5) +plot_hist(img_rescale) -plt.subplots_adjust(left=0.05, hspace=0.25, wspace=0.3, top=0.95, bottom=0.1) +plt.subplot(2, 3, 3) +plt.imshow(img_eq, cmap=plt.cm.gray, vmin=0, vmax=1) +plt.title('Histogram equalization') +plt.axis('off') +plt.subplot(2, 3, 6) +plot_hist(img_eq) +plt.ylabel('Fraction of total intensity') + +# prevent overlap of y-axis labels +plt.subplots_adjust(wspace=0.4) plt.show() diff --git a/skimage/exposure/__init__.py b/skimage/exposure/__init__.py index ee5969ab..83933f62 100644 --- a/skimage/exposure/__init__.py +++ b/skimage/exposure/__init__.py @@ -1 +1,2 @@ from exposure import histogram, equalize, cumulative_distribution +from exposure import rescale_intensity diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py index 42a60dbf..2d3b82f5 100644 --- a/skimage/exposure/exposure.py +++ b/skimage/exposure/exposure.py @@ -1,9 +1,11 @@ import numpy as np import skimage +from skimage.util.dtype import dtype_range -__all__ = ['histogram', 'cumulative_distribution', 'equalize'] +__all__ = ['histogram', 'cumulative_distribution', 'equalize', + 'rescale_intensity'] def histogram(image, nbins=256): @@ -104,3 +106,79 @@ def equalize(image, nbins=256): out = np.interp(image.flat, bin_centers, cdf) return out.reshape(image.shape) + +def rescale_intensity(image, in_range=None, out_range=None): + """Return image after stretching or shrinking its intensity levels. + + The image intensities are uniformly rescaled such that the minimum and + maximum values given by `in_range` match those given by `out_range`. + + Parameters + ---------- + image : array + Image array. + in_range : 2-tuple (float, float) + Min and max *allowed* intensity values of input image. If None, the + *allowed* min/max values are set to the *actual* min/max values in the + input image. + out_range : 2-tuple (float, float) + Min and max intensity values of output image. If None, use the min/max + intensities of the image data type. See `skimage.util.dtype` for + details. + + Returns + ------- + out : array + Image array after rescaling its intensity. This image is the same dtype + as the input image. + + Examples + -------- + By default, intensities are stretched to the limits allowed by the dtype: + >>> image = np.array([51, 102, 153], dtype=np.uint8) + >>> rescale_intensity(image) + array([ 0, 127, 255], dtype=uint8) + + It's easy to accidentally convert an image dtype from uint8 to float: + >>> 1.0 * image + array([ 51., 102., 153.]) + + Use `rescale_intensity` to rescale to the proper range for float dtypes: + >>> image_float = 1.0 * image + >>> rescale_intensity(image_float) + array([ 0. , 0.5, 1. ]) + + To maintain the low contrast of the original, use the `in_range` parameter: + >>> rescale_intensity(image_float, in_range=(0, 255)) + array([ 0.2, 0.4, 0.6]) + + If the min/max value of `in_range` is more/less than the min/max image + intensity, then the intensity levels are clipped: + >>> rescale_intensity(image_float, in_range=(0, 102)) + array([ 0.5, 1. , 1. ]) + + If you have an image with signed integers but want to rescale the image to + just the positive range, use the `out_range` parameter: + >>> image = np.array([-10, 0, 10], dtype=np.int8) + >>> rescale_intensity(image, out_range=(0, 127)) + array([ 0, 63, 127], dtype=int8) + + """ + dtype = image.dtype.type + + if in_range is None: + imin = np.min(image) + imax = np.max(image) + else: + imin, imax = in_range + + if out_range is None: + omin, omax = dtype_range[dtype] + else: + omin, omax = out_range + + image = np.clip(image, imin, imax) + + image = (image - imin) / float(imax - imin) + return dtype(image * (omax - omin) + omin) + diff --git a/skimage/exposure/tests/test_exposure.py b/skimage/exposure/tests/test_exposure.py index 22339de7..6b961de1 100644 --- a/skimage/exposure/tests/test_exposure.py +++ b/skimage/exposure/tests/test_exposure.py @@ -1,10 +1,14 @@ import numpy as np +from numpy.testing import assert_array_almost_equal as assert_close import skimage from skimage import data from skimage import exposure +# Test histogram equalization +# =========================== + # squeeze image intensities to lower image contrast test_img = data.camera() / 5 + 100 @@ -31,6 +35,41 @@ def check_cdf_slope(cdf): assert 0.9 < slope < 1.1 +# Test rescale intensity +# ====================== + +def test_rescale_stretch(): + image = np.array([51, 102, 153], dtype=np.uint8) + out = exposure.rescale_intensity(image) + assert out.dtype == np.uint8 + assert_close(out, [0, 127, 255]) + + +def test_rescale_shrink(): + image = np.array([51., 102., 153.]) + out = exposure.rescale_intensity(image) + assert_close(out, [0, 0.5, 1]) + + +def test_rescale_in_range(): + image = np.array([51., 102., 153.]) + out = exposure.rescale_intensity(image, in_range=(0, 255)) + assert_close(out, [0.2, 0.4, 0.6]) + + +def test_rescale_in_range_clip(): + image = np.array([51., 102., 153.]) + out = exposure.rescale_intensity(image, in_range=(0, 102)) + assert_close(out, [0.5, 1, 1]) + + +def test_rescale_out_range(): + image = np.array([-10, 0, 10], dtype=np.int8) + out = exposure.rescale_intensity(image, out_range=(0, 127)) + assert out.dtype == np.int8 + assert_close(out, [0, 63, 127]) + + if __name__ == '__main__': from numpy import testing testing.run_module_suite()