diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py
index 8e5609bc..2f0bb13e 100644
--- a/doc/examples/plot_equalize.py
+++ b/doc/examples/plot_equalize.py
@@ -3,70 +3,81 @@
 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.
 
-For comparison, this example also shows an image after its intensity levels are
-uniformly stretched.
+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 matplotlib import ticker
 
 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.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)
+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)
+
+    # Display cumulative distribution
+    img_cdf, bins = exposure.cumulative_distribution(img, bins)
+    ax_cdf.plot(bins, img_cdf, 'r')
+
+    return ax_img, ax_hist, ax_cdf
 
 
-img_orig = data.camera()
-# squeeze image intensities to lower image contrast
-img = img_orig / 5 + 100
-img_rescale = exposure.rescale_intensity(img)
+# Load an example image
+img = data.moon()
+
+# Contrast stretching
+p5 = np.percentile(img, 2)
+p95 = np.percentile(img, 98)
+img_rescale = exposure.rescale_intensity(img, in_range=(p5, p95))
+
+# Equalization
 img_eq = exposure.equalize(img)
 
-plt.subplot(2, 3, 1)
-plt.imshow(img, cmap=plt.cm.gray, vmin=0, vmax=255)
-plt.title('Low contrast image')
-plt.axis('off')
-plt.subplot(2, 3, 4)
-plt.ylabel('Number of pixels')
-plot_hist(img)
 
-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, 3, 5)
-plot_hist(img_rescale)
+# Display results
+f, axes = plt.subplots(2, 3, figsize=(11, 5))
+
+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_cdf = plot_img_and_hist(img_eq, axes[:, 2])
+ax_img.set_title('Histogram equalization')
+ax_cdf.set_ylabel('Fraction of total intensity')
 
-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)
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