diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py index 9cb30dd7..111cf891 100644 --- a/doc/examples/plot_equalize.py +++ b/doc/examples/plot_equalize.py @@ -3,54 +3,83 @@ Histogram Equalization ====================== -This examples takes an image with low contrast and enhances its contrast using -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. +This examples enhances an image with low contrast, using a method called +*histogram equalization*, which "spreads out the most frequent intensity +values" in an image [1]_. The equalized image has a roughly linear cumulative +distribution function. + +While histogram equalization has the advantage that it requires no parameters, +it sometimes yields unnatural looking images. An alternative method is +*contrast stretching*, where the image is rescaled to include all intensities +that fall within the 2nd and 98th percentiles [2]_. .. [1] http://en.wikipedia.org/wiki/Histogram_equalization +.. [2] http://homepages.inf.ed.ac.uk/rbf/HIPR2/stretch.htm """ -import matplotlib.pyplot as plt from skimage import data from skimage.util.dtype import dtype_range from skimage import exposure +import matplotlib.pyplot as plt -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') +import numpy as np + +def plot_img_and_hist(img, axes, bins=256): + """Plot an image along with its histogram and cumulative histogram. + + """ + ax_img, ax_hist = axes + ax_cdf = ax_hist.twinx() + + # Display image + ax_img.imshow(img, cmap=plt.cm.gray) + ax_img.set_axis_off() + + # Display histogram + ax_hist.hist(img.ravel(), bins=bins) + ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0)) + ax_hist.set_xlabel('Pixel intensity') - ax_cdf = plt.twinx() - ax_cdf.plot(bins, img_cdf, 'r') xmin, xmax = dtype_range[img.dtype.type] - plt.xlim(xmin, xmax) + ax_hist.set_xlim(xmin, xmax) - ax_cdf.set_ylabel('Fraction of total intensity') + # Display cumulative distribution + img_cdf, bins = exposure.cumulative_distribution(img, bins) + ax_cdf.plot(bins, img_cdf, 'r') -img_orig = data.camera() -# squeeze image intensities to lower image contrast -img = img_orig / 5 + 100 + return ax_img, ax_hist, ax_cdf + + +# Load an example image +img = data.moon() + +# Contrast stretching +p2 = np.percentile(img, 2) +p98 = np.percentile(img, 98) +img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98)) + +# Equalization img_eq = exposure.equalize(img) -plt.subplot(2, 2, 1) -plt.imshow(img, cmap=plt.cm.gray, vmin=0, vmax=255) -plt.title('Low contrast input image') -plt.axis('off') -plt.subplot(2, 2, 2) -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.axis('off') -plt.subplot(2, 2, 4) -plot_hist(img_eq) +# Display results +f, axes = plt.subplots(2, 3, figsize=(11, 5)) -plt.subplots_adjust(left=0.05, hspace=0.25, wspace=0.3, top=0.95, bottom=0.1) +ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0]) +ax_img.set_title('Low contrast image') +ax_hist.set_ylabel('Number of pixels') + +ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_rescale, axes[:, 1]) +ax_img.set_title('Contrast stretching') + +ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_eq, axes[:, 2]) +ax_img.set_title('Histogram equalization') +ax_cdf.set_ylabel('Fraction of total intensity') + + +# prevent overlap of y-axis labels +plt.subplots_adjust(wspace=0.4) plt.show() diff --git a/skimage/data/__init__.py b/skimage/data/__init__.py index 6d955f4f..c5159092 100644 --- a/skimage/data/__init__.py +++ b/skimage/data/__init__.py @@ -73,3 +73,12 @@ def coins(): """ return load("coins.png") + +def moon(): + """Surface of the moon. + + This low-contrast image of the surface of the moon is useful for + illustrating histogram equalization and contrast stretching. + + """ + return load("moon.png") diff --git a/skimage/data/moon.png b/skimage/data/moon.png new file mode 100644 index 00000000..6a6f077e Binary files /dev/null and b/skimage/data/moon.png differ 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()