Merge pull request #107 from tonysyu/rescale-intensity

ENH: Add rescale_intensity function with test and example.

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()
 

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")

skimage/exposure/__init__.py

@@ -1 +1,2 @@
 from exposure import histogram, equalize, cumulative_distribution
+from exposure import rescale_intensity

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)
+

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()