diff --git a/doc/examples/applications/plot_coins_segmentation.py b/doc/examples/applications/plot_coins_segmentation.py
index fa7aa1f4..49ff399f 100644
--- a/doc/examples/applications/plot_coins_segmentation.py
+++ b/doc/examples/applications/plot_coins_segmentation.py
@@ -16,13 +16,11 @@ from skimage import data
 coins = data.coins()
 hist = np.histogram(coins, bins=np.arange(0, 256))
 
-plt.figure(figsize=(8, 3))
-plt.subplot(121)
-plt.imshow(coins, cmap=plt.cm.gray, interpolation='nearest')
-plt.axis('off')
-plt.subplot(122)
-plt.plot(hist[1][:-1], hist[0], lw=2)
-plt.title('histogram of grey values')
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3))
+ax1.imshow(coins, cmap=plt.cm.gray, interpolation='nearest')
+ax1.axis('off')
+ax2.plot(hist[1][:-1], hist[0], lw=2)
+ax2.set_title('histogram of grey values')
 
 """
 .. image:: PLOT2RST.current_figure
@@ -37,17 +35,15 @@ background with the coins:
 
 """
 
-plt.figure(figsize=(6, 3))
-plt.subplot(121)
-plt.imshow(coins > 100, cmap=plt.cm.gray, interpolation='nearest')
-plt.title('coins > 100')
-plt.axis('off')
-plt.subplot(122)
-plt.imshow(coins > 150, cmap=plt.cm.gray, interpolation='nearest')
-plt.title('coins > 150')
-plt.axis('off')
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3))
+ax1.imshow(coins > 100, cmap=plt.cm.gray, interpolation='nearest')
+ax1.set_title('coins > 100')
+ax1.axis('off')
+ax2.imshow(coins > 150, cmap=plt.cm.gray, interpolation='nearest')
+ax2.set_title('coins > 150')
+ax2.axis('off')
 margins = dict(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
-plt.subplots_adjust(**margins)
+fig.subplots_adjust(**margins)
 
 """
 .. image:: PLOT2RST.current_figure
@@ -64,10 +60,10 @@ edge-detector.
 from skimage.filter import canny
 edges = canny(coins/255.)
 
-plt.figure(figsize=(4, 3))
-plt.imshow(edges, cmap=plt.cm.gray, interpolation='nearest')
-plt.axis('off')
-plt.title('Canny detector')
+fig, ax = plt.subplots(figsize=(4, 3))
+ax.imshow(edges, cmap=plt.cm.gray, interpolation='nearest')
+ax.axis('off')
+ax.set_title('Canny detector')
 
 """
 .. image:: PLOT2RST.current_figure
@@ -79,10 +75,10 @@ from scipy import ndimage
 
 fill_coins = ndimage.binary_fill_holes(edges)
 
-plt.figure(figsize=(4, 3))
-plt.imshow(fill_coins, cmap=plt.cm.gray, interpolation='nearest')
-plt.axis('off')
-plt.title('Filling the holes')
+fig, ax = plt.subplots(figsize=(4, 3))
+ax.imshow(fill_coins, cmap=plt.cm.gray, interpolation='nearest')
+ax.axis('off')
+ax.set_title('Filling the holes')
 
 """
 .. image:: PLOT2RST.current_figure
@@ -93,10 +89,10 @@ objects.
 from skimage import morphology
 coins_cleaned = morphology.remove_small_objects(fill_coins, 21)
 
-plt.figure(figsize=(4, 3))
-plt.imshow(coins_cleaned, cmap=plt.cm.gray, interpolation='nearest')
-plt.axis('off')
-plt.title('Removing small objects')
+fig, ax = plt.subplots(figsize=(4, 3))
+ax.imshow(coins_cleaned, cmap=plt.cm.gray, interpolation='nearest')
+ax.axis('off')
+ax.set_title('Removing small objects')
 
 """
 .. image:: PLOT2RST.current_figure
@@ -117,10 +113,10 @@ from skimage.filter import sobel
 
 elevation_map = sobel(coins)
 
-plt.figure(figsize=(4, 3))
-plt.imshow(elevation_map, cmap=plt.cm.jet, interpolation='nearest')
-plt.axis('off')
-plt.title('elevation_map')
+fig, ax = plt.subplots(figsize=(4, 3))
+ax.imshow(elevation_map, cmap=plt.cm.jet, interpolation='nearest')
+ax.axis('off')
+ax.set_title('elevation_map')
 
 """
 .. image:: PLOT2RST.current_figure
@@ -133,10 +129,10 @@ markers = np.zeros_like(coins)
 markers[coins < 30] = 1
 markers[coins > 150] = 2
 
-plt.figure(figsize=(4, 3))
-plt.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest')
-plt.axis('off')
-plt.title('markers')
+fig, ax = plt.subplots(figsize=(4, 3))
+ax.imshow(markers, cmap=plt.cm.spectral, interpolation='nearest')
+ax.axis('off')
+ax.set_title('markers')
 
 """
 .. image:: PLOT2RST.current_figure
@@ -147,10 +143,10 @@ starting from the markers determined above:
 """
 segmentation = morphology.watershed(elevation_map, markers)
 
-plt.figure(figsize=(4, 3))
-plt.imshow(segmentation, cmap=plt.cm.gray, interpolation='nearest')
-plt.axis('off')
-plt.title('segmentation')
+fig, ax = plt.subplots(figsize=(4, 3))
+ax.imshow(segmentation, cmap=plt.cm.gray, interpolation='nearest')
+ax.axis('off')
+ax.set_title('segmentation')
 
 """
 .. image:: PLOT2RST.current_figure
@@ -166,16 +162,14 @@ segmentation = ndimage.binary_fill_holes(segmentation - 1)
 labeled_coins, _ = ndimage.label(segmentation)
 image_label_overlay = label2rgb(labeled_coins, image=coins)
 
-plt.figure(figsize=(6, 3))
-plt.subplot(121)
-plt.imshow(coins, cmap=plt.cm.gray, interpolation='nearest')
-plt.contour(segmentation, [0.5], linewidths=1.2, colors='y')
-plt.axis('off')
-plt.subplot(122)
-plt.imshow(image_label_overlay, interpolation='nearest')
-plt.axis('off')
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3))
+ax1.imshow(coins, cmap=plt.cm.gray, interpolation='nearest')
+ax1.contour(segmentation, [0.5], linewidths=1.2, colors='y')
+ax1.axis('off')
+ax2.imshow(image_label_overlay, interpolation='nearest')
+ax2.axis('off')
 
-plt.subplots_adjust(**margins)
+fig.subplots_adjust(**margins)
 
 """
 .. image:: PLOT2RST.current_figure
diff --git a/doc/examples/applications/plot_morphology.py b/doc/examples/applications/plot_morphology.py
index 84cf2972..68ecc6d6 100644
--- a/doc/examples/applications/plot_morphology.py
+++ b/doc/examples/applications/plot_morphology.py
@@ -30,9 +30,9 @@ from skimage.data import data_dir
 from skimage.util import img_as_ubyte
 from skimage import io
 
-plt.gray()
 phantom = img_as_ubyte(io.imread(data_dir+'/phantom.png', as_grey=True))
-plt.imshow(phantom)
+fig, ax = plt.subplots()
+ax.imshow(phantom, cmap=plt.cm.gray)
 
 """
 .. image:: PLOT2RST.current_figure
@@ -43,10 +43,10 @@ Let's also define a convenience function for plotting comparisons:
 def plot_comparison(original, filtered, filter_name):
 
     fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4))
-    ax1.imshow(original)
+    ax1.imshow(original, cmap=plt.cm.gray)
     ax1.set_title('original')
     ax1.axis('off')
-    ax2.imshow(filtered)
+    ax2.imshow(filtered, cmap=plt.cm.gray)
     ax2.set_title(filter_name)
     ax2.axis('off')
 
diff --git a/doc/examples/applications/plot_rank_filters.py b/doc/examples/applications/plot_rank_filters.py
index 87dad0d0..5fa88aa0 100644
--- a/doc/examples/applications/plot_rank_filters.py
+++ b/doc/examples/applications/plot_rank_filters.py
@@ -44,13 +44,11 @@ from skimage import data
 noisy_image = img_as_ubyte(data.camera())
 hist = np.histogram(noisy_image, bins=np.arange(0, 256))
 
-plt.figure(figsize=(8, 3))
-plt.subplot(1, 2, 1)
-plt.imshow(noisy_image, interpolation='nearest')
-plt.axis('off')
-plt.subplot(1, 2, 2)
-plt.plot(hist[1][:-1], hist[0], lw=2)
-plt.title('Histogram of grey values')
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 3))
+ax1.imshow(noisy_image, interpolation='nearest', cmap=plt.cm.gray)
+ax1.axis('off')
+ax2.plot(hist[1][:-1], hist[0], lw=2)
+ax2.set_title('Histogram of grey values')
 
 """
 
@@ -77,27 +75,24 @@ noisy_image = img_as_ubyte(data.camera())
 noisy_image[noise > 0.99] = 255
 noisy_image[noise < 0.01] = 0
 
-fig = plt.figure(figsize=(10, 7))
+fig, ax = plt.subplots(2, 2, figsize=(10, 7))
+ax1, ax2, ax3, ax4 = ax.ravel()
 
-plt.subplot(2, 2, 1)
-plt.imshow(noisy_image, vmin=0, vmax=255)
-plt.title('Noisy image')
-plt.axis('off')
+ax1.imshow(noisy_image, vmin=0, vmax=255, cmap=plt.cm.gray)
+ax1.set_title('Noisy image')
+ax1.axis('off')
 
-plt.subplot(2, 2, 2)
-plt.imshow(median(noisy_image, disk(1)), vmin=0, vmax=255)
-plt.title('Median $r=1$')
-plt.axis('off')
+ax2.imshow(median(noisy_image, disk(1)), vmin=0, vmax=255, cmap=plt.cm.gray)
+ax2.set_title('Median $r=1$')
+ax2.axis('off')
 
-plt.subplot(2, 2, 3)
-plt.imshow(median(noisy_image, disk(5)), vmin=0, vmax=255)
-plt.title('Median $r=5$')
-plt.axis('off')
+ax3.imshow(median(noisy_image, disk(5)), vmin=0, vmax=255, cmap=plt.cm.gray)
+ax3.set_title('Median $r=5$')
+ax3.axis('off')
 
-plt.subplot(2, 2, 4)
-plt.imshow(median(noisy_image, disk(20)), vmin=0, vmax=255)
-plt.title('Median $r=20$')
-plt.axis('off')
+ax4.imshow(median(noisy_image, disk(20)), vmin=0, vmax=255, cmap=plt.cm.gray)
+ax4.set_title('Median $r=20$')
+ax4.axis('off')
 
 """
 
@@ -119,19 +114,17 @@ image.
 
 from skimage.filter.rank import mean
 
-fig = plt.figure(figsize=[10, 7])
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 7])
 
 loc_mean = mean(noisy_image, disk(10))
 
-plt.subplot(1, 2, 1)
-plt.imshow(noisy_image, vmin=0, vmax=255)
-plt.title('Original')
-plt.axis('off')
+ax1.imshow(noisy_image, vmin=0, vmax=255, cmap=plt.cm.gray)
+ax1.set_title('Original')
+ax1.axis('off')
 
-plt.subplot(1, 2, 2)
-plt.imshow(loc_mean, vmin=0, vmax=255)
-plt.title('Local mean $r=10$')
-plt.axis('off')
+ax2.imshow(loc_mean, vmin=0, vmax=255, cmap=plt.cm.gray)
+ax2.set_title('Local mean $r=10$')
+ax2.axis('off')
 
 """
 
@@ -155,25 +148,22 @@ noisy_image = img_as_ubyte(data.camera())
 
 bilat = bilateral_mean(noisy_image.astype(np.uint16), disk(20), s0=10, s1=10)
 
-fig = plt.figure(figsize=[10, 7])
+fig, ax = plt.subplots(2, 2, figsize=(10, 7))
+ax1, ax2, ax3, ax4 = ax.ravel()
 
-plt.subplot(2, 2, 1)
-plt.imshow(noisy_image, cmap=plt.cm.gray)
-plt.title('Original')
-plt.axis('off')
+ax1.imshow(noisy_image, cmap=plt.cm.gray)
+ax1.set_title('Original')
+ax1.axis('off')
 
-plt.subplot(2, 2, 3)
-plt.imshow(bilat, cmap=plt.cm.gray)
-plt.title('Bilateral mean')
-plt.axis('off')
+ax2.imshow(bilat, cmap=plt.cm.gray)
+ax2.set_title('Bilateral mean')
+ax2.axis('off')
 
-plt.subplot(2, 2, 2)
-plt.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)
-plt.axis('off')
+ax3.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)
+ax3.axis('off')
 
-plt.subplot(2, 2, 4)
-plt.imshow(bilat[200:350, 350:450], cmap=plt.cm.gray)
-plt.axis('off')
+ax4.imshow(bilat[200:350, 350:450], cmap=plt.cm.gray)
+ax4.axis('off')
 
 """
 
@@ -211,31 +201,26 @@ hist = np.histogram(noisy_image, bins=np.arange(0, 256))
 glob_hist = np.histogram(glob, bins=np.arange(0, 256))
 loc_hist = np.histogram(loc, bins=np.arange(0, 256))
 
-plt.figure(figsize=(10, 10))
+fig, ax = plt.subplots(3, 2, figsize=(10, 10))
+ax1, ax2, ax3, ax4, ax5, ax6 = ax.ravel()
 
-plt.subplot(321)
-plt.imshow(noisy_image, interpolation='nearest')
-plt.axis('off')
+ax1.imshow(noisy_image, interpolation='nearest', cmap=plt.cm.gray)
+ax1.axis('off')
 
-plt.subplot(322)
-plt.plot(hist[1][:-1], hist[0], lw=2)
-plt.title('Histogram of gray values')
+ax2.plot(hist[1][:-1], hist[0], lw=2)
+ax2.set_title('Histogram of gray values')
 
-plt.subplot(323)
-plt.imshow(glob, interpolation='nearest')
-plt.axis('off')
+ax3.imshow(glob, interpolation='nearest', cmap=plt.cm.gray)
+ax3.axis('off')
 
-plt.subplot(324)
-plt.plot(glob_hist[1][:-1], glob_hist[0], lw=2)
-plt.title('Histogram of gray values')
+ax4.plot(glob_hist[1][:-1], glob_hist[0], lw=2)
+ax4.set_title('Histogram of gray values')
 
-plt.subplot(325)
-plt.imshow(loc, interpolation='nearest')
-plt.axis('off')
+ax5.imshow(loc, interpolation='nearest', cmap=plt.cm.gray)
+ax5.axis('off')
 
-plt.subplot(326)
-plt.plot(loc_hist[1][:-1], loc_hist[0], lw=2)
-plt.title('Histogram of gray values')
+ax6.plot(loc_hist[1][:-1], loc_hist[0], lw=2)
+ax6.set_title('Histogram of gray values')
 
 """
 
@@ -256,17 +241,15 @@ noisy_image = img_as_ubyte(data.camera())
 
 auto = autolevel(noisy_image.astype(np.uint16), disk(20))
 
-fig = plt.figure(figsize=[10, 7])
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 7])
 
-plt.subplot(1, 2, 1)
-plt.imshow(noisy_image, cmap=plt.cm.gray)
-plt.title('Original')
-plt.axis('off')
+ax1.imshow(noisy_image, cmap=plt.cm.gray)
+ax1.set_title('Original')
+ax1.axis('off')
 
-plt.subplot(1, 2, 2)
-plt.imshow(auto, cmap=plt.cm.gray)
-plt.title('Local autolevel')
-plt.axis('off')
+ax2.imshow(auto, cmap=plt.cm.gray)
+ax2.set_title('Local autolevel')
+ax2.axis('off')
 
 """
 
@@ -297,7 +280,7 @@ fig, axes = plt.subplots(nrows=3, figsize=(7, 8))
 ax0, ax1, ax2 = axes
 plt.gray()
 
-ax0.imshow(np.hstack((image, loc_autolevel)))
+ax0.imshow(np.hstack((image, loc_autolevel)), cmap=plt.cm.gray)
 ax0.set_title('Original / auto-level')
 
 ax1.imshow(
@@ -326,24 +309,22 @@ noisy_image = img_as_ubyte(data.camera())
 
 enh = enhance_contrast(noisy_image, disk(5))
 
-fig = plt.figure(figsize=[10, 7])
-plt.subplot(2, 2, 1)
-plt.imshow(noisy_image, cmap=plt.cm.gray)
-plt.title('Original')
-plt.axis('off')
+fig, ax = plt.subplots(2, 2, figsize=[10, 7])
+ax1, ax2, ax3, ax4 = ax.ravel()
 
-plt.subplot(2, 2, 3)
-plt.imshow(enh, cmap=plt.cm.gray)
-plt.title('Local morphological contrast enhancement')
-plt.axis('off')
+ax1.imshow(noisy_image, cmap=plt.cm.gray)
+ax1.set_title('Original')
+ax1.axis('off')
 
-plt.subplot(2, 2, 2)
-plt.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)
-plt.axis('off')
+ax2.imshow(enh, cmap=plt.cm.gray)
+ax2.set_title('Local morphological contrast enhancement')
+ax2.axis('off')
 
-plt.subplot(2, 2, 4)
-plt.imshow(enh[200:350, 350:450], cmap=plt.cm.gray)
-plt.axis('off')
+ax3.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)
+ax3.axis('off')
+
+ax4.imshow(enh[200:350, 350:450], cmap=plt.cm.gray)
+ax4.axis('off')
 
 """
 
@@ -360,24 +341,22 @@ noisy_image = img_as_ubyte(data.camera())
 
 penh = enhance_contrast_percentile(noisy_image, disk(5), p0=.1, p1=.9)
 
-fig = plt.figure(figsize=[10, 7])
-plt.subplot(2, 2, 1)
-plt.imshow(noisy_image, cmap=plt.cm.gray)
-plt.title('Original')
-plt.axis('off')
+fig, ax = plt.subplots(2, 2, figsize=[10, 7])
+ax1, ax2, ax3, ax4 = ax.ravel()
 
-plt.subplot(2, 2, 3)
-plt.imshow(penh, cmap=plt.cm.gray)
-plt.title('Local percentile morphological\n contrast enhancement')
-plt.axis('off')
+ax1.imshow(noisy_image, cmap=plt.cm.gray)
+ax1.set_title('Original')
+ax1.axis('off')
 
-plt.subplot(2, 2, 2)
-plt.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)
-plt.axis('off')
+ax2.imshow(penh, cmap=plt.cm.gray)
+ax2.set_title('Local percentile morphological\n contrast enhancement')
+ax2.axis('off')
 
-plt.subplot(2, 2, 4)
-plt.imshow(penh[200:350, 350:450], cmap=plt.cm.gray)
-plt.axis('off')
+ax3.imshow(noisy_image[200:350, 350:450], cmap=plt.cm.gray)
+ax3.axis('off')
+
+ax4.imshow(penh[200:350, 350:450], cmap=plt.cm.gray)
+ax4.axis('off')
 
 """
 
@@ -419,29 +398,24 @@ loc_otsu = p8 >= t_loc_otsu
 t_glob_otsu = threshold_otsu(p8)
 glob_otsu = p8 >= t_glob_otsu
 
-plt.figure()
+fig, ax = plt.subplots(2, 2)
+ax1, ax2, ax3, ax4 = ax.ravel()
 
-plt.subplot(2, 2, 1)
-plt.imshow(p8, cmap=plt.cm.gray)
-plt.title('Original')
-plt.colorbar()
-plt.axis('off')
+fig.colorbar(ax1.imshow(p8, cmap=plt.cm.gray), ax=ax1)
+ax1.set_title('Original')
+ax1.axis('off')
 
-plt.subplot(2, 2, 2)
-plt.imshow(t_loc_otsu, cmap=plt.cm.gray)
-plt.title('Local Otsu ($r=%d$)' % radius)
-plt.colorbar()
-plt.axis('off')
+fig.colorbar(ax2.imshow(t_loc_otsu, cmap=plt.cm.gray), ax=ax2)
+ax2.set_title('Local Otsu ($r=%d$)' % radius)
+ax2.axis('off')
 
-plt.subplot(2, 2, 3)
-plt.imshow(p8 >= t_loc_otsu, cmap=plt.cm.gray)
-plt.title('Original >= local Otsu' % t_glob_otsu)
-plt.axis('off')
+ax3.imshow(p8 >= t_loc_otsu, cmap=plt.cm.gray)
+ax3.set_title('Original >= local Otsu' % t_glob_otsu)
+ax3.axis('off')
 
-plt.subplot(2, 2, 4)
-plt.imshow(glob_otsu, cmap=plt.cm.gray)
-plt.title('Global Otsu ($t=%d$)' % t_glob_otsu)
-plt.axis('off')
+ax4.imshow(glob_otsu, cmap=plt.cm.gray)
+ax4.set_title('Global Otsu ($t=%d$)' % t_glob_otsu)
+ax4.axis('off')
 
 """
 
@@ -460,17 +434,15 @@ m = (np.tile(x, (n, 1)) * np.linspace(0.1, 1, n) * 128 + 128).astype(np.uint8)
 radius = 10
 t = rank.otsu(m, disk(radius))
 
-plt.figure()
+fig, (ax1, ax2) = plt.subplots(1, 2)
 
-plt.subplot(1, 2, 1)
-plt.imshow(m)
-plt.title('Original')
-plt.axis('off')
+ax1.imshow(m)
+ax1.set_title('Original')
+ax1.axis('off')
 
-plt.subplot(1, 2, 2)
-plt.imshow(m >= t, interpolation='nearest')
-plt.title('Local Otsu ($r=%d$)' % radius)
-plt.axis('off')
+ax2.imshow(m >= t, interpolation='nearest')
+ax2.set_title('Local Otsu ($r=%d$)' % radius)
+ax2.axis('off')
 
 """
 
@@ -501,27 +473,24 @@ opening = minimum(maximum(noisy_image, disk(5)), disk(5))
 grad = gradient(noisy_image, disk(5))
 
 # display results
-fig = plt.figure(figsize=[10, 7])
+fig, ax = plt.subplots(2, 2, figsize=[10, 7])
+ax1, ax2, ax3, ax4 = ax.ravel()
 
-plt.subplot(2, 2, 1)
-plt.imshow(noisy_image, cmap=plt.cm.gray)
-plt.title('Original')
-plt.axis('off')
+ax1.imshow(noisy_image, cmap=plt.cm.gray)
+ax1.set_title('Original')
+ax1.axis('off')
 
-plt.subplot(2, 2, 2)
-plt.imshow(closing, cmap=plt.cm.gray)
-plt.title('Gray-level closing')
-plt.axis('off')
+ax2.imshow(closing, cmap=plt.cm.gray)
+ax2.set_title('Gray-level closing')
+ax2.axis('off')
 
-plt.subplot(2, 2, 3)
-plt.imshow(opening, cmap=plt.cm.gray)
-plt.title('Gray-level opening')
-plt.axis('off')
+ax3.imshow(opening, cmap=plt.cm.gray)
+ax3.set_title('Gray-level opening')
+ax3.axis('off')
 
-plt.subplot(2, 2, 4)
-plt.imshow(grad, cmap=plt.cm.gray)
-plt.title('Morphological gradient')
-plt.axis('off')
+ax4.imshow(grad, cmap=plt.cm.gray)
+ax4.set_title('Morphological gradient')
+ax4.axis('off')
 
 """
 
@@ -554,19 +523,15 @@ import matplotlib.pyplot as plt
 
 image = data.camera()
 
-plt.figure(figsize=(10, 4))
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
 
-plt.subplot(1, 2, 1)
-plt.imshow(image, cmap=plt.cm.gray)
-plt.title('Image')
-plt.colorbar()
-plt.axis('off')
+fig.colorbar(ax1.imshow(image, cmap=plt.cm.gray), ax=ax1)
+ax1.set_title('Image')
+ax1.axis('off')
 
-plt.subplot(1, 2, 2)
-plt.imshow(entropy(image, disk(5)), cmap=plt.cm.jet)
-plt.title('Entropy')
-plt.colorbar()
-plt.axis('off')
+fig.colorbar(ax2.imshow(entropy(image, disk(5)), cmap=plt.cm.jet), ax=ax2)
+ax2.set_title('Entropy')
+ax2.axis('off')
 
 """
 
@@ -651,12 +616,12 @@ for r in e_range:
 
 rec = np.asarray(rec)
 
-plt.figure()
-plt.title('Performance with respect to element size')
-plt.ylabel('Time (ms)')
-plt.title('Element radius')
-plt.plot(e_range, rec)
-plt.legend(['filter.rank.maximum', 'morphology.dilate'])
+fig, ax = plt.subplots()
+ax.set_title('Performance with respect to element size')
+ax.set_ylabel('Time (ms)')
+ax.set_xlabel('Element radius')
+ax.plot(e_range, rec)
+ax.legend(['filter.rank.maximum', 'morphology.dilate'])
 
 """
 
@@ -679,12 +644,12 @@ for s in s_range:
 
 rec = np.asarray(rec)
 
-plt.figure()
-plt.title('Performance with respect to image size')
-plt.ylabel('Time (ms)')
-plt.title('Image size')
-plt.plot(s_range, rec)
-plt.legend(['filter.rank.maximum', 'morphology.dilate'])
+fig, ax = plt.subplots()
+ax.set_title('Performance with respect to image size')
+ax.set_ylabel('Time (ms)')
+ax.set_xlabel('Image size')
+ax.plot(s_range, rec)
+ax.legend(['filter.rank.maximum', 'morphology.dilate'])
 
 
 """
@@ -714,13 +679,13 @@ for r in e_range:
 
 rec = np.asarray(rec)
 
-plt.figure()
-plt.title('Performance with respect to element size')
-plt.plot(e_range, rec)
-plt.legend(['filter.rank.median', 'filter.median_filter',
-            'scipy.ndimage.percentile'])
-plt.ylabel('Time (ms)')
-plt.title('Element radius')
+fig, ax = plt.subplots()
+ax.set_title('Performance with respect to element size')
+ax.plot(e_range, rec)
+ax.legend(['filter.rank.median', 'filter.median_filter',
+           'scipy.ndimage.percentile'])
+ax.set_ylabel('Time (ms)')
+ax.set_xlabel('Element radius')
 
 """
 .. image:: PLOT2RST.current_figure
@@ -729,10 +694,10 @@ Comparison of outcome of the three methods:
 
 """
 
-plt.figure()
-plt.imshow(np.hstack((rc, rctmf, rndi)))
-plt.title('filter.rank.median vs filtermedian_filter vs scipy.ndimage.percentile')
-plt.axis('off')
+fig, ax = plt.subplots()
+ax.imshow(np.hstack((rc, rctmf, rndi)))
+ax.set_title('filter.rank.median vs filtermedian_filter vs scipy.ndimage.percentile')
+ax.axis('off')
 
 """
 .. image:: PLOT2RST.current_figure
@@ -755,13 +720,13 @@ for s in s_range:
 
 rec = np.asarray(rec)
 
-plt.figure()
-plt.title('Performance with respect to image size')
-plt.plot(s_range, rec)
-plt.legend(['filter.rank.median', 'filter.median_filter',
-            'scipy.ndimage.percentile'])
-plt.ylabel('Time (ms)')
-plt.title('Image size')
+fig, ax = plt.subplots()
+ax.set_title('Performance with respect to image size')
+ax.plot(s_range, rec)
+ax.legend(['filter.rank.median', 'filter.median_filter',
+           'scipy.ndimage.percentile'])
+ax.set_ylabel('Time (ms)')
+ax.set_xlabel('Image size')
 
 """
 .. image:: PLOT2RST.current_figure
diff --git a/doc/examples/plot_canny.py b/doc/examples/plot_canny.py
index f8269d70..f1caf264 100644
--- a/doc/examples/plot_canny.py
+++ b/doc/examples/plot_canny.py
@@ -35,24 +35,21 @@ edges1 = filter.canny(im)
 edges2 = filter.canny(im, sigma=3)
 
 # display results
-plt.figure(figsize=(8, 3))
+fig, (ax1, ax2, ax3) = plt.subplots(nrows=1, ncols=3, figsize=(8, 3))
 
-plt.subplot(131)
-plt.imshow(im, cmap=plt.cm.jet)
-plt.axis('off')
-plt.title('noisy image', fontsize=20)
+ax1.imshow(im, cmap=plt.cm.jet)
+ax1.axis('off')
+ax1.set_title('noisy image', fontsize=20)
 
-plt.subplot(132)
-plt.imshow(edges1, cmap=plt.cm.gray)
-plt.axis('off')
-plt.title('Canny filter, $\sigma=1$', fontsize=20)
+ax2.imshow(edges1, cmap=plt.cm.gray)
+ax2.axis('off')
+ax2.set_title('Canny filter, $\sigma=1$', fontsize=20)
 
-plt.subplot(133)
-plt.imshow(edges2, cmap=plt.cm.gray)
-plt.axis('off')
-plt.title('Canny filter, $\sigma=3$', fontsize=20)
+ax3.imshow(edges2, cmap=plt.cm.gray)
+ax3.axis('off')
+ax3.set_title('Canny filter, $\sigma=3$', fontsize=20)
 
-plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9,
+fig.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9,
                     bottom=0.02, left=0.02, right=0.98)
 
 plt.show()
diff --git a/doc/examples/plot_contours.py b/doc/examples/plot_contours.py
index d23b2ae5..e686b4ac 100644
--- a/doc/examples/plot_contours.py
+++ b/doc/examples/plot_contours.py
@@ -29,12 +29,13 @@ r = np.sin(np.exp((np.sin(x)**3 + np.cos(y)**2)))
 contours = measure.find_contours(r, 0.8)
 
 # Display the image and plot all contours found
-plt.imshow(r, interpolation='nearest')
+fig, ax = plt.subplots()
+ax.imshow(r, interpolation='nearest', cmap=plt.cm.gray)
 
 for n, contour in enumerate(contours):
-    plt.plot(contour[:, 1], contour[:, 0], linewidth=2)
+    ax.plot(contour[:, 1], contour[:, 0], linewidth=2)
 
-plt.axis('image')
-plt.xticks([])
-plt.yticks([])
+ax.axis('image')
+ax.set_xticks([])
+ax.set_yticks([])
 plt.show()
diff --git a/doc/examples/plot_convex_hull.py b/doc/examples/plot_convex_hull.py
index 31398e6d..e811c1ee 100644
--- a/doc/examples/plot_convex_hull.py
+++ b/doc/examples/plot_convex_hull.py
@@ -40,11 +40,9 @@ image[chull] += 1
 # [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]]
 
 
-fig = plt.subplots(figsize=(10, 6))
-plt.subplot(1, 2, 1)
-plt.title('Original picture')
-plt.imshow(original_image, cmap=plt.cm.gray, interpolation='nearest')
-plt.subplot(1, 2, 2)
-plt.title('Transformed picture')
-plt.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 6))
+ax1.set_title('Original picture')
+ax1.imshow(original_image, cmap=plt.cm.gray, interpolation='nearest')
+ax2.set_title('Transformed picture')
+ax2.imshow(image, cmap=plt.cm.gray, interpolation='nearest')
 plt.show()
diff --git a/doc/examples/plot_corner.py b/doc/examples/plot_corner.py
index a3233188..4e962e17 100644
--- a/doc/examples/plot_corner.py
+++ b/doc/examples/plot_corner.py
@@ -29,9 +29,9 @@ image[230:280, 60:110] = 1
 coords = corner_peaks(corner_harris(image), min_distance=5)
 coords_subpix = corner_subpix(image, coords, window_size=13)
 
-plt.gray()
-plt.imshow(image, interpolation='nearest')
-plt.plot(coords[:, 1], coords[:, 0], '.b', markersize=3)
-plt.plot(coords_subpix[:, 1], coords_subpix[:, 0], '+r', markersize=15)
-plt.axis((0, 350, 350, 0))
+fig, ax = plt.subplots()
+ax.imshow(image, interpolation='nearest', cmap=plt.cm.gray)
+ax.plot(coords[:, 1], coords[:, 0], '.b', markersize=3)
+ax.plot(coords_subpix[:, 1], coords_subpix[:, 0], '+r', markersize=15)
+ax.axis((0, 350, 350, 0))
 plt.show()
diff --git a/doc/examples/plot_daisy.py b/doc/examples/plot_daisy.py
index fbde5c8d..128d6e67 100644
--- a/doc/examples/plot_daisy.py
+++ b/doc/examples/plot_daisy.py
@@ -20,8 +20,9 @@ img = data.camera()
 descs, descs_img = daisy(img, step=180, radius=58, rings=2, histograms=6,
                          orientations=8, visualize=True)
 
-plt.axis('off')
-plt.imshow(descs_img)
+fig, ax = plt.subplots()
+ax.axis('off')
+ax.imshow(descs_img)
 descs_num = descs.shape[0] * descs.shape[1]
-plt.title('%i DAISY descriptors extracted:' % descs_num)
+ax.set_title('%i DAISY descriptors extracted:' % descs_num)
 plt.show()
diff --git a/doc/examples/plot_entropy.py b/doc/examples/plot_entropy.py
index 9f208c73..b92a5bf3 100644
--- a/doc/examples/plot_entropy.py
+++ b/doc/examples/plot_entropy.py
@@ -22,11 +22,11 @@ fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 4))
 img0 = ax0.imshow(image, cmap=plt.cm.gray)
 ax0.set_title('Image')
 ax0.axis('off')
-plt.colorbar(img0, ax=ax0)
+fig.colorbar(img0, ax=ax0)
 
 img1 = ax1.imshow(entropy(image, disk(5)), cmap=plt.cm.jet)
 ax1.set_title('Entropy')
 ax1.axis('off')
-plt.colorbar(img1, ax=ax1)
+fig.colorbar(img1, ax=ax1)
 
 plt.show()
diff --git a/doc/examples/plot_equalize.py b/doc/examples/plot_equalize.py
index 3b3f3788..2d04e22a 100644
--- a/doc/examples/plot_equalize.py
+++ b/doc/examples/plot_equalize.py
@@ -70,7 +70,7 @@ img_eq = exposure.equalize_hist(img)
 img_adapteq = exposure.equalize_adapthist(img, clip_limit=0.03)
 
 # Display results
-f, axes = plt.subplots(nrows=2, ncols=4, figsize=(8, 5))
+fig, axes = plt.subplots(nrows=2, ncols=4, figsize=(8, 5))
 
 ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0])
 ax_img.set_title('Low contrast image')
@@ -92,5 +92,5 @@ ax_cdf.set_ylabel('Fraction of total intensity')
 ax_cdf.set_yticks(np.linspace(0, 1, 5))
 
 # prevent overlap of y-axis labels
-plt.subplots_adjust(wspace=0.4)
+fig.subplots_adjust(wspace=0.4)
 plt.show()
diff --git a/doc/examples/plot_glcm.py b/doc/examples/plot_glcm.py
index 7ee54d57..af9b4168 100644
--- a/doc/examples/plot_glcm.py
+++ b/doc/examples/plot_glcm.py
@@ -53,44 +53,45 @@ for i, patch in enumerate(grass_patches + sky_patches):
     ys.append(greycoprops(glcm, 'correlation')[0, 0])
 
 # create the figure
-plt.figure(figsize=(8, 8))
+fig = plt.figure(figsize=(8, 8))
+
+# display original image with locations of patches
+ax = fig.add_subplot(3, 2, 1)
+ax.imshow(image, cmap=plt.cm.gray, interpolation='nearest',
+          vmin=0, vmax=255)
+for (y, x) in grass_locations:
+    ax.plot(x + PATCH_SIZE / 2, y + PATCH_SIZE / 2, 'gs')
+for (y, x) in sky_locations:
+    ax.plot(x + PATCH_SIZE / 2, y + PATCH_SIZE / 2, 'bs')
+ax.set_xlabel('Original Image')
+ax.set_xticks([])
+ax.set_yticks([])
+ax.axis('image')
+
+# for each patch, plot (dissimilarity, correlation)
+ax = fig.add_subplot(3, 2, 2)
+ax.plot(xs[:len(grass_patches)], ys[:len(grass_patches)], 'go',
+        label='Grass')
+ax.plot(xs[len(grass_patches):], ys[len(grass_patches):], 'bo',
+        label='Sky')
+ax.set_xlabel('GLCM Dissimilarity')
+ax.set_ylabel('GLVM Correlation')
+ax.legend()
 
 # display the image patches
 for i, patch in enumerate(grass_patches):
-    plt.subplot(3, len(grass_patches), len(grass_patches) * 1 + i + 1)
-    plt.imshow(patch, cmap=plt.cm.gray, interpolation='nearest',
-               vmin=0, vmax=255)
-    plt.xlabel('Grass %d' % (i + 1))
+    ax = fig.add_subplot(3, len(grass_patches), len(grass_patches)*1 + i + 1)
+    ax.imshow(patch, cmap=plt.cm.gray, interpolation='nearest',
+              vmin=0, vmax=255)
+    ax.set_xlabel('Grass %d' % (i + 1))
 
 for i, patch in enumerate(sky_patches):
-    plt.subplot(3, len(grass_patches), len(grass_patches) * 2 + i + 1)
-    plt.imshow(patch, cmap=plt.cm.gray, interpolation='nearest',
-               vmin=0, vmax=255)
-    plt.xlabel('Sky %d' % (i + 1))
+    ax = fig.add_subplot(3, len(sky_patches), len(sky_patches)*2 + i + 1)
+    ax.imshow(patch, cmap=plt.cm.gray, interpolation='nearest',
+              vmin=0, vmax=255)
+    ax.set_xlabel('Sky %d' % (i + 1))
 
-# display original image with locations of patches
-plt.subplot(3, 2, 1)
-plt.imshow(image, cmap=plt.cm.gray, interpolation='nearest',
-           vmin=0, vmax=255)
-for (y, x) in grass_locations:
-    plt.plot(x + PATCH_SIZE / 2, y + PATCH_SIZE / 2, 'gs')
-for (y, x) in sky_locations:
-    plt.plot(x + PATCH_SIZE / 2, y + PATCH_SIZE / 2, 'bs')
-plt.xlabel('Original Image')
-plt.xticks([])
-plt.yticks([])
-plt.axis('image')
-
-# for each patch, plot (dissimilarity, correlation)
-plt.subplot(3, 2, 2)
-plt.plot(xs[:len(grass_patches)], ys[:len(grass_patches)], 'go',
-         label='Grass')
-plt.plot(xs[len(grass_patches):], ys[len(grass_patches):], 'bo',
-         label='Sky')
-plt.xlabel('GLCM Dissimilarity')
-plt.ylabel('GLVM Correlation')
-plt.legend()
 
 # display the patches and plot
-plt.suptitle('Grey level co-occurrence matrix features', fontsize=14)
+fig.suptitle('Grey level co-occurrence matrix features', fontsize=14)
 plt.show()
diff --git a/doc/examples/plot_hog.py b/doc/examples/plot_hog.py
index 4a339a1b..4899c590 100644
--- a/doc/examples/plot_hog.py
+++ b/doc/examples/plot_hog.py
@@ -90,16 +90,16 @@ image = color.rgb2gray(data.lena())
 fd, hog_image = hog(image, orientations=8, pixels_per_cell=(16, 16),
                     cells_per_block=(1, 1), visualise=True)
 
-plt.figure(figsize=(8, 4))
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
 
-plt.subplot(121).set_axis_off()
-plt.imshow(image, cmap=plt.cm.gray)
-plt.title('Input image')
+ax1.axis('off')
+ax1.imshow(image, cmap=plt.cm.gray)
+ax1.set_title('Input image')
 
 # Rescale histogram for better display
 hog_image_rescaled = exposure.rescale_intensity(hog_image, in_range=(0, 0.02))
 
-plt.subplot(122).set_axis_off()
-plt.imshow(hog_image_rescaled, cmap=plt.cm.gray)
-plt.title('Histogram of Oriented Gradients')
+ax2.axis('off')
+ax2.imshow(hog_image_rescaled, cmap=plt.cm.gray)
+ax2.set_title('Histogram of Oriented Gradients')
 plt.show()
diff --git a/doc/examples/plot_holes_and_peaks.py b/doc/examples/plot_holes_and_peaks.py
index e40b777e..791c69ee 100644
--- a/doc/examples/plot_holes_and_peaks.py
+++ b/doc/examples/plot_holes_and_peaks.py
@@ -22,13 +22,13 @@ image = data.moon()
 image = rescale_intensity(image, in_range=(50, 200))
 
 # convenience function for plotting images
-def imshow(image, **kwargs):
-    plt.figure(figsize=(5, 4))
-    plt.imshow(image, **kwargs)
-    plt.axis('off')
+def imshow(image, title, **kwargs):
+    fig, ax = plt.subplots(figsize=(5, 4))
+    ax.imshow(image, **kwargs)
+    ax.axis('off')
+    ax.set_title(title)
 
-imshow(image)
-plt.title('original image')
+imshow(image, 'Original image')
 
 """
 .. image:: PLOT2RST.current_figure
@@ -50,8 +50,7 @@ mask = image
 
 filled = reconstruction(seed, mask, method='erosion')
 
-imshow(filled, vmin=image.min(), vmax=image.max())
-plt.title('after filling holes')
+imshow(filled, 'after filling holes',vmin=image.min(), vmax=image.max())
 
 """
 .. image:: PLOT2RST.current_figure
@@ -62,8 +61,8 @@ isolate the dark regions by subtracting the reconstructed image from the
 original image.
 """
 
-imshow(image - filled)
-plt.title('holes')
+imshow(image - filled, 'holes')
+# plt.title('holes')
 
 """
 .. image:: PLOT2RST.current_figure
@@ -78,8 +77,7 @@ intensity instead of the maximum. The remainder of the process is the same.
 seed = np.copy(image)
 seed[1:-1, 1:-1] = image.min()
 rec = reconstruction(seed, mask, method='dilation')
-imshow(image - rec)
-plt.title('peaks')
+imshow(image - rec, 'peaks')
 plt.show()
 
 """
diff --git a/doc/examples/plot_ihc_color_separation.py b/doc/examples/plot_ihc_color_separation.py
index f89f28fe..3e297386 100644
--- a/doc/examples/plot_ihc_color_separation.py
+++ b/doc/examples/plot_ihc_color_separation.py
@@ -61,10 +61,10 @@ h = rescale_intensity(ihc_hed[:, :, 0], out_range=(0, 1))
 d = rescale_intensity(ihc_hed[:, :, 2], out_range=(0, 1))
 zdh = np.dstack((np.zeros_like(h), d, h))
 
-plt.figure()
-plt.imshow(zdh)
-plt.title("Stain separated image (rescaled)")
-plt.axis('off')
+fig, ax = plt.subplots()
+ax.imshow(zdh)
+ax.set_title("Stain separated image (rescaled)")
+ax.axis('off')
 plt.show()
 
 """
diff --git a/doc/examples/plot_join_segmentations.py b/doc/examples/plot_join_segmentations.py
index 8ccc5038..849f087c 100644
--- a/doc/examples/plot_join_segmentations.py
+++ b/doc/examples/plot_join_segmentations.py
@@ -58,5 +58,5 @@ axes[3].set_title('Join')
 
 for ax in axes:
     ax.axis('off')
-plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
+fig.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
 plt.show()
diff --git a/doc/examples/plot_line_hough_transform.py b/doc/examples/plot_line_hough_transform.py
index cd0ae008..bdb05661 100644
--- a/doc/examples/plot_line_hough_transform.py
+++ b/doc/examples/plot_line_hough_transform.py
@@ -77,30 +77,30 @@ image[idx, idx] = 255
 
 h, theta, d = hough_line(image)
 
-plt.figure(figsize=(8, 4))
+fig, ax = plt.subplots(1, 3, figsize=(8, 4))
 
-plt.subplot(131)
-plt.imshow(image, cmap=plt.cm.gray)
-plt.title('Input image')
+ax[0].imshow(image, cmap=plt.cm.gray)
+ax[0].set_title('Input image')
+ax[0].axis('image')
 
-plt.subplot(132)
-plt.imshow(np.log(1 + h),
+ax[1].imshow(np.log(1 + h),
            extent=[np.rad2deg(theta[-1]), np.rad2deg(theta[0]),
                    d[-1], d[0]],
            cmap=plt.cm.gray, aspect=1/1.5)
-plt.title('Hough transform')
-plt.xlabel('Angles (degrees)')
-plt.ylabel('Distance (pixels)')
+ax[1].set_title('Hough transform')
+ax[1].set_xlabel('Angles (degrees)')
+ax[1].set_ylabel('Distance (pixels)')
+ax[1].axis('image')
 
-plt.subplot(133)
-plt.imshow(image, cmap=plt.cm.gray)
+ax[2].imshow(image, cmap=plt.cm.gray)
 rows, cols = image.shape
 for _, angle, dist in zip(*hough_line_peaks(h, theta, d)):
     y0 = (dist - 0 * np.cos(angle)) / np.sin(angle)
     y1 = (dist - cols * np.cos(angle)) / np.sin(angle)
-    plt.plot((0, cols), (y0, y1), '-r')
-plt.axis((0, cols, rows, 0))
-plt.title('Detected lines')
+    ax[2].plot((0, cols), (y0, y1), '-r')
+ax[2].axis((0, cols, rows, 0))
+ax[2].set_title('Detected lines')
+ax[2].axis('image')
 
 # Line finding, using the Probabilistic Hough Transform
 
@@ -108,23 +108,22 @@ image = data.camera()
 edges = canny(image, 2, 1, 25)
 lines = probabilistic_hough_line(edges, threshold=10, line_length=5, line_gap=3)
 
-plt.figure(figsize=(8, 3))
+fig2, ax = plt.subplots(1, 3, figsize=(8, 3))
 
-plt.subplot(131)
-plt.imshow(image, cmap=plt.cm.gray)
-plt.title('Input image')
+ax[0].imshow(image, cmap=plt.cm.gray)
+ax[0].set_title('Input image')
+ax[0].axis('image')
 
-plt.subplot(132)
-plt.imshow(edges, cmap=plt.cm.gray)
-plt.title('Canny edges')
+ax[1].imshow(edges, cmap=plt.cm.gray)
+ax[1].set_title('Canny edges')
+ax[1].axis('image')
 
-plt.subplot(133)
-plt.imshow(edges * 0)
+ax[2].imshow(edges * 0)
 
 for line in lines:
     p0, p1 = line
-    plt.plot((p0[0], p1[0]), (p0[1], p1[1]))
+    ax[2].plot((p0[0], p1[0]), (p0[1], p1[1]))
 
-plt.title('Probabilistic Hough')
-plt.axis('image')
+ax[2].set_title('Probabilistic Hough')
+ax[2].axis('image')
 plt.show()
diff --git a/doc/examples/plot_local_equalize.py b/doc/examples/plot_local_equalize.py
index 309d1e43..b3c11d55 100644
--- a/doc/examples/plot_local_equalize.py
+++ b/doc/examples/plot_local_equalize.py
@@ -72,7 +72,7 @@ img_eq = rank.equalize(img, selem=selem)
 
 
 # Display results
-f, axes = plt.subplots(2, 3, figsize=(8, 5))
+fig, axes = plt.subplots(2, 3, figsize=(8, 5))
 
 ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0])
 ax_img.set_title('Low contrast image')
@@ -87,5 +87,5 @@ ax_cdf.set_ylabel('Fraction of total intensity')
 
 
 # prevent overlap of y-axis labels
-plt.subplots_adjust(wspace=0.4)
+fig.subplots_adjust(wspace=0.4)
 plt.show()
diff --git a/doc/examples/plot_local_otsu.py b/doc/examples/plot_local_otsu.py
index aff7db62..d40af6b6 100644
--- a/doc/examples/plot_local_otsu.py
+++ b/doc/examples/plot_local_otsu.py
@@ -37,28 +37,25 @@ threshold_global_otsu = threshold_otsu(img)
 global_otsu = img >= threshold_global_otsu
 
 
-plt.figure(figsize=(8, 5))
+fig, ax = plt.subplots(2, 2, figsize=(8, 5))
+ax1, ax2, ax3, ax4 = ax.ravel()
 
-plt.subplot(2, 2, 1)
-plt.imshow(img, cmap=plt.cm.gray)
-plt.title('Original')
-plt.colorbar(orientation='horizontal')
-plt.axis('off')
+fig.colorbar(ax1.imshow(img, cmap=plt.cm.gray),
+           ax=ax1, orientation='horizontal')
+ax1.set_title('Original')
+ax1.axis('off')
 
-plt.subplot(2, 2, 2)
-plt.imshow(local_otsu, cmap=plt.cm.gray)
-plt.title('Local Otsu (radius=%d)' % radius)
-plt.colorbar(orientation='horizontal')
-plt.axis('off')
+fig.colorbar(ax2.imshow(local_otsu, cmap=plt.cm.gray),
+           ax=ax2, orientation='horizontal')
+ax2.set_title('Local Otsu (radius=%d)' % radius)
+ax2.axis('off')
 
-plt.subplot(2, 2, 3)
-plt.imshow(img >= local_otsu, cmap=plt.cm.gray)
-plt.title('Original >= Local Otsu' % threshold_global_otsu)
-plt.axis('off')
+ax3.imshow(img >= local_otsu, cmap=plt.cm.gray)
+ax3.set_title('Original >= Local Otsu' % threshold_global_otsu)
+ax3.axis('off')
 
-plt.subplot(2, 2, 4)
-plt.imshow(global_otsu, cmap=plt.cm.gray)
-plt.title('Global Otsu (threshold = %d)' % threshold_global_otsu)
-plt.axis('off')
+ax4.imshow(global_otsu, cmap=plt.cm.gray)
+ax4.set_title('Global Otsu (threshold = %d)' % threshold_global_otsu)
+ax4.axis('off')
 
 plt.show()
diff --git a/doc/examples/plot_marked_watershed.py b/doc/examples/plot_marked_watershed.py
index d7f2c354..0e1f1c58 100644
--- a/doc/examples/plot_marked_watershed.py
+++ b/doc/examples/plot_marked_watershed.py
@@ -51,5 +51,5 @@ ax3.imshow(labels, cmap=plt.cm.spectral, interpolation='nearest', alpha=.7)
 for ax in axes:
     ax.axis('off')
 
-plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
+fig.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
 plt.show()
diff --git a/doc/examples/plot_medial_transform.py b/doc/examples/plot_medial_transform.py
index f0792e50..1a8469b2 100644
--- a/doc/examples/plot_medial_transform.py
+++ b/doc/examples/plot_medial_transform.py
@@ -55,14 +55,12 @@ skel, distance = medial_axis(data, return_distance=True)
 # Distance to the background for pixels of the skeleton
 dist_on_skel = distance * skel
 
-plt.figure(figsize=(8, 4))
-plt.subplot(121)
-plt.imshow(data, cmap=plt.cm.gray, interpolation='nearest')
-plt.axis('off')
-plt.subplot(122)
-plt.imshow(dist_on_skel, cmap=plt.cm.spectral, interpolation='nearest')
-plt.contour(data, [0.5], colors='w')
-plt.axis('off')
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
+ax1.imshow(data, cmap=plt.cm.gray, interpolation='nearest')
+ax1.axis('off')
+ax2.imshow(dist_on_skel, cmap=plt.cm.spectral, interpolation='nearest')
+ax2.contour(data, [0.5], colors='w')
+ax2.axis('off')
 
-plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
+fig.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, right=1)
 plt.show()
diff --git a/doc/examples/plot_otsu.py b/doc/examples/plot_otsu.py
index b3b8e279..059093ab 100644
--- a/doc/examples/plot_otsu.py
+++ b/doc/examples/plot_otsu.py
@@ -28,20 +28,17 @@ image = camera()
 thresh = threshold_otsu(image)
 binary = image > thresh
 
-plt.figure(figsize=(8, 2.5))
-plt.subplot(1, 3, 1)
-plt.imshow(image, cmap=plt.cm.gray)
-plt.title('Original')
-plt.axis('off')
+fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8, 2.5))
+ax1.imshow(image, cmap=plt.cm.gray)
+ax1.set_title('Original')
+ax1.axis('off')
 
-plt.subplot(1, 3, 2, aspect='equal')
-plt.hist(image)
-plt.title('Histogram')
-plt.axvline(thresh, color='r')
+ax2.hist(image)
+ax2.set_title('Histogram')
+ax2.axvline(thresh, color='r')
 
-plt.subplot(1, 3, 3)
-plt.imshow(binary, cmap=plt.cm.gray)
-plt.title('Thresholded')
-plt.axis('off')
+ax3.imshow(binary, cmap=plt.cm.gray)
+ax3.set_title('Thresholded')
+ax3.axis('off')
 
 plt.show()
diff --git a/doc/examples/plot_peak_local_max.py b/doc/examples/plot_peak_local_max.py
index 3a2dccfe..cebd5738 100644
--- a/doc/examples/plot_peak_local_max.py
+++ b/doc/examples/plot_peak_local_max.py
@@ -25,25 +25,23 @@ image_max = ndimage.maximum_filter(im, size=20, mode='constant')
 coordinates = peak_local_max(im, min_distance=20)
 
 # display results
-plt.figure(figsize=(8, 3))
-plt.subplot(131)
-plt.imshow(im, cmap=plt.cm.gray)
-plt.axis('off')
-plt.title('Original')
+fig, ax = plt.subplots(1, 3, figsize=(8, 3))
+ax1, ax2, ax3 = ax.ravel()
+ax1.imshow(im, cmap=plt.cm.gray)
+ax1.axis('off')
+ax1.set_title('Original')
 
-plt.subplot(132)
-plt.imshow(image_max, cmap=plt.cm.gray)
-plt.axis('off')
-plt.title('Maximum filter')
+ax2.imshow(image_max, cmap=plt.cm.gray)
+ax2.axis('off')
+ax2.set_title('Maximum filter')
 
-plt.subplot(133)
-plt.imshow(im, cmap=plt.cm.gray)
-plt.autoscale(False)
-plt.plot([p[1] for p in coordinates], [p[0] for p in coordinates], 'r.')
-plt.axis('off')
-plt.title('Peak local max')
+ax3.imshow(im, cmap=plt.cm.gray)
+ax3.autoscale(False)
+ax3.plot([p[1] for p in coordinates], [p[0] for p in coordinates], 'r.')
+ax3.axis('off')
+ax3.set_title('Peak local max')
 
-plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9,
+fig.subplots_adjust(wspace=0.02, hspace=0.02, top=0.9,
                     bottom=0.02, left=0.02, right=0.98)
 
 plt.show()
diff --git a/doc/examples/plot_phase_unwrap.py b/doc/examples/plot_phase_unwrap.py
index 22e3bc89..ded87ada 100644
--- a/doc/examples/plot_phase_unwrap.py
+++ b/doc/examples/plot_phase_unwrap.py
@@ -26,26 +26,20 @@ image_wrapped = np.angle(np.exp(1j * image))
 # Perform phase unwrapping
 image_unwrapped = unwrap_phase(image_wrapped)
 
-plt.figure()
-plt.subplot(221)
-plt.title('Original')
-plt.imshow(image, cmap='gray', vmin=0, vmax=4 * np.pi)
-plt.colorbar()
+fig, ax = plt.subplots(2, 2)
+ax1, ax2, ax3, ax4 = ax.ravel()
 
-plt.subplot(222)
-plt.title('Wrapped phase')
-plt.imshow(image_wrapped, cmap='gray', vmin=-np.pi, vmax=np.pi)
-plt.colorbar()
+fig.colorbar(ax1.imshow(image, cmap='gray', vmin=0, vmax=4 * np.pi), ax=ax1)
+ax1.set_title('Original')
 
-plt.subplot(223)
-plt.title('After phase unwrapping')
-plt.imshow(image_unwrapped, cmap='gray')
-plt.colorbar()
+fig.colorbar(ax2.imshow(image_wrapped, cmap='gray', vmin=-np.pi, vmax=np.pi), ax=ax2)
+ax2.set_title('Wrapped phase')
 
-plt.subplot(224)
-plt.title('Unwrapped minus original')
-plt.imshow(image_unwrapped - image, cmap='gray')
-plt.colorbar()
+fig.colorbar(ax3.imshow(image_unwrapped, cmap='gray'), ax=ax3)
+ax3.set_title('After phase unwrapping')
+
+fig.colorbar(ax4.imshow(image_unwrapped - image, cmap='gray'), ax=ax4)
+ax4.set_title('Unwrapped minus original')
 
 """
 .. image:: PLOT2RST.current_figure
@@ -70,26 +64,22 @@ image_unwrapped_no_wrap_around = unwrap_phase(image_wrapped,
 image_unwrapped_wrap_around = unwrap_phase(image_wrapped,
                                            wrap_around=(True, False))
 
-plt.figure()
-plt.subplot(221)
-plt.title('Original')
-plt.imshow(np.ma.array(image, mask=mask), cmap='jet')
-plt.colorbar()
+fig, ax = plt.subplots(2, 2)
+ax1, ax2, ax3, ax4 = ax.ravel()
 
-plt.subplot(222)
-plt.title('Wrapped phase')
-plt.imshow(image_wrapped, cmap='jet', vmin=-np.pi, vmax=np.pi)
-plt.colorbar()
+fig.colorbar(ax1.imshow(np.ma.array(image, mask=mask), cmap='jet'), ax=ax1)
+ax1.set_title('Original')
 
-plt.subplot(223)
-plt.title('Unwrapped without wrap_around')
-plt.imshow(image_unwrapped_no_wrap_around, cmap='jet')
-plt.colorbar()
+fig.colorbar(ax2.imshow(image_wrapped, cmap='jet', vmin=-np.pi, vmax=np.pi),
+           ax=ax2)
+ax2.set_title('Wrapped phase')
 
-plt.subplot(224)
-plt.title('Unwrapped with wrap_around')
-plt.imshow(image_unwrapped_wrap_around, cmap='jet')
-plt.colorbar()
+fig.colorbar(ax3.imshow(image_unwrapped_no_wrap_around, cmap='jet'),
+           ax=ax3)
+ax3.set_title('Unwrapped without wrap_around')
+
+fig.colorbar(ax4.imshow(image_unwrapped_wrap_around, cmap='jet'), ax=ax4)
+ax4.set_title('Unwrapped with wrap_around')
 
 plt.show()
 
diff --git a/doc/examples/plot_piecewise_affine.py b/doc/examples/plot_piecewise_affine.py
index a0ad7d8f..b1fb551f 100644
--- a/doc/examples/plot_piecewise_affine.py
+++ b/doc/examples/plot_piecewise_affine.py
@@ -35,7 +35,8 @@ out_rows = image.shape[0] - 1.5 * 50
 out_cols = cols
 out = warp(image, tform, output_shape=(out_rows, out_cols))
 
-plt.imshow(out)
-plt.plot(tform.inverse(src)[:, 0], tform.inverse(src)[:, 1], '.b')
-plt.axis((0, out_cols, out_rows, 0))
+fig, ax = plt.subplots()
+ax.imshow(out)
+ax.plot(tform.inverse(src)[:, 0], tform.inverse(src)[:, 1], '.b')
+ax.axis((0, out_cols, out_rows, 0))
 plt.show()
diff --git a/doc/examples/plot_pyramid.py b/doc/examples/plot_pyramid.py
index c4853c2f..8e309edf 100644
--- a/doc/examples/plot_pyramid.py
+++ b/doc/examples/plot_pyramid.py
@@ -30,5 +30,6 @@ for p in pyramid[1:]:
     composite_image[i_row:i_row + n_rows, cols:cols + n_cols] = p
     i_row += n_rows
 
-plt.imshow(composite_image)
+fig, ax = plt.subplots()
+ax.imshow(composite_image)
 plt.show()
diff --git a/doc/examples/plot_radon_transform.py b/doc/examples/plot_radon_transform.py
index 0dab82e6..04cf7bab 100644
--- a/doc/examples/plot_radon_transform.py
+++ b/doc/examples/plot_radon_transform.py
@@ -60,22 +60,20 @@ from skimage.transform import radon, rescale
 image = imread(data_dir + "/phantom.png", as_grey=True)
 image = rescale(image, scale=0.4)
 
-plt.figure(figsize=(8, 4.5))
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))
 
-plt.subplot(121)
-plt.title("Original")
-plt.imshow(image, cmap=plt.cm.Greys_r)
+ax1.set_title("Original")
+ax1.imshow(image, cmap=plt.cm.Greys_r)
 
 theta = np.linspace(0., 180., max(image.shape), endpoint=True)
 sinogram = radon(image, theta=theta, circle=True)
-plt.subplot(122)
-plt.title("Radon transform\n(Sinogram)")
-plt.xlabel("Projection angle (deg)")
-plt.ylabel("Projection position (pixels)")
-plt.imshow(sinogram, cmap=plt.cm.Greys_r,
+ax2.set_title("Radon transform\n(Sinogram)")
+ax2.set_xlabel("Projection angle (deg)")
+ax2.set_ylabel("Projection position (pixels)")
+ax2.imshow(sinogram, cmap=plt.cm.Greys_r,
            extent=(0, 180, 0, sinogram.shape[0]), aspect='auto')
 
-plt.subplots_adjust(hspace=0.4, wspace=0.5)
+fig.subplots_adjust(hspace=0.4, wspace=0.5)
 plt.show()
 
 """
@@ -103,13 +101,11 @@ error = reconstruction_fbp - image
 print('FBP rms reconstruction error: %.3g' % np.sqrt(np.mean(error**2)))
 
 imkwargs = dict(vmin=-0.2, vmax=0.2)
-plt.figure(figsize=(8, 4.5))
-plt.subplot(121)
-plt.title("Reconstruction\nFiltered back projection")
-plt.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r)
-plt.subplot(122)
-plt.title("Reconstruction error\nFiltered back projection")
-plt.imshow(reconstruction_fbp - image, cmap=plt.cm.Greys_r, **imkwargs)
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))
+ax1.set_title("Reconstruction\nFiltered back projection")
+ax1.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r)
+ax2.set_title("Reconstruction error\nFiltered back projection")
+ax2.imshow(reconstruction_fbp - image, cmap=plt.cm.Greys_r, **imkwargs)
 plt.show()
 
 """
@@ -156,14 +152,12 @@ error = reconstruction_sart - image
 print('SART (1 iteration) rms reconstruction error: %.3g'
       % np.sqrt(np.mean(error**2)))
 
-plt.figure(figsize=(8, 8.5))
-
-plt.subplot(221)
-plt.title("Reconstruction\nSART")
-plt.imshow(reconstruction_sart, cmap=plt.cm.Greys_r)
-plt.subplot(222)
-plt.title("Reconstruction error\nSART")
-plt.imshow(reconstruction_sart - image, cmap=plt.cm.Greys_r, **imkwargs)
+fig, ax = plt.subplots(2, 2, figsize=(8, 8.5))
+ax1, ax2, ax3, ax4 = ax.ravel()
+ax1.set_title("Reconstruction\nSART")
+ax1.imshow(reconstruction_sart, cmap=plt.cm.Greys_r)
+ax2.set_title("Reconstruction error\nSART")
+ax2.imshow(reconstruction_sart - image, cmap=plt.cm.Greys_r, **imkwargs)
 
 # Run a second iteration of SART by supplying the reconstruction
 # from the first iteration as an initial estimate
@@ -173,12 +167,10 @@ error = reconstruction_sart2 - image
 print('SART (2 iterations) rms reconstruction error: %.3g'
       % np.sqrt(np.mean(error**2)))
 
-plt.subplot(223)
-plt.title("Reconstruction\nSART, 2 iterations")
-plt.imshow(reconstruction_sart2, cmap=plt.cm.Greys_r)
-plt.subplot(224)
-plt.title("Reconstruction error\nSART, 2 iterations")
-plt.imshow(reconstruction_sart2 - image, cmap=plt.cm.Greys_r, **imkwargs)
+ax3.set_title("Reconstruction\nSART, 2 iterations")
+ax3.imshow(reconstruction_sart2, cmap=plt.cm.Greys_r)
+ax4.set_title("Reconstruction error\nSART, 2 iterations")
+ax4.imshow(reconstruction_sart2 - image, cmap=plt.cm.Greys_r, **imkwargs)
 plt.show()
 
 """
diff --git a/doc/examples/plot_random_walker_segmentation.py b/doc/examples/plot_random_walker_segmentation.py
index 8c960a8b..fa70e503 100644
--- a/doc/examples/plot_random_walker_segmentation.py
+++ b/doc/examples/plot_random_walker_segmentation.py
@@ -58,20 +58,17 @@ markers[data > 1.3] = 2
 labels = random_walker(data, markers, beta=10, mode='bf')
 
 # Plot results
-plt.figure(figsize=(8, 3.2))
-plt.subplot(131)
-plt.imshow(data, cmap='gray', interpolation='nearest')
-plt.axis('off')
-plt.title('Noisy data')
-plt.subplot(132)
-plt.imshow(markers, cmap='hot', interpolation='nearest')
-plt.axis('off')
-plt.title('Markers')
-plt.subplot(133)
-plt.imshow(labels, cmap='gray', interpolation='nearest')
-plt.axis('off')
-plt.title('Segmentation')
+fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(8, 3.2))
+ax1.imshow(data, cmap='gray', interpolation='nearest')
+ax1.axis('off')
+ax1.set_title('Noisy data')
+ax2.imshow(markers, cmap='hot', interpolation='nearest')
+ax2.axis('off')
+ax2.set_title('Markers')
+ax3.imshow(labels, cmap='gray', interpolation='nearest')
+ax3.axis('off')
+ax3.set_title('Segmentation')
 
-plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
+fig.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
                     right=1)
 plt.show()
diff --git a/doc/examples/plot_ransac.py b/doc/examples/plot_ransac.py
index 2de76619..770e8cd1 100644
--- a/doc/examples/plot_ransac.py
+++ b/doc/examples/plot_ransac.py
@@ -45,11 +45,12 @@ line_x = np.arange(-250, 250)
 line_y = model.predict_y(line_x)
 line_y_robust = model_robust.predict_y(line_x)
 
-plt.plot(data[inliers, 0], data[inliers, 1], '.b', alpha=0.6,
-         label='Inlier data')
-plt.plot(data[outliers, 0], data[outliers, 1], '.r', alpha=0.6,
-         label='Outlier data')
-plt.plot(line_x, line_y, '-k', label='Line model from all data')
-plt.plot(line_x, line_y_robust, '-b', label='Robust line model')
-plt.legend(loc='lower left')
+fig, ax = plt.subplots()
+ax.plot(data[inliers, 0], data[inliers, 1], '.b', alpha=0.6,
+        label='Inlier data')
+ax.plot(data[outliers, 0], data[outliers, 1], '.r', alpha=0.6,
+        label='Outlier data')
+ax.plot(line_x, line_y, '-k', label='Line model from all data')
+ax.plot(line_x, line_y_robust, '-b', label='Robust line model')
+ax.legend(loc='lower left')
 plt.show()
diff --git a/doc/examples/plot_regional_maxima.py b/doc/examples/plot_regional_maxima.py
index f7b125f1..dd676014 100644
--- a/doc/examples/plot_regional_maxima.py
+++ b/doc/examples/plot_regional_maxima.py
@@ -50,7 +50,7 @@ ax3.imshow(image - dilated)
 ax3.set_title('image - dilated')
 ax3.axis('off')
 
-plt.tight_layout()
+fig.tight_layout()
 
 """
 
@@ -98,7 +98,7 @@ ax3.axhline(yslice, color='r', alpha=0.4)
 ax3.set_title('image - dilated')
 ax3.axis('off')
 
-plt.tight_layout()
+fig.tight_layout()
 plt.show()
 
 """
diff --git a/doc/examples/plot_regionprops.py b/doc/examples/plot_regionprops.py
index f675d11c..065bae55 100644
--- a/doc/examples/plot_regionprops.py
+++ b/doc/examples/plot_regionprops.py
@@ -26,7 +26,8 @@ image = rotate(image, angle=15, order=0)
 label_img = label(image)
 regions = regionprops(label_img)
 
-plt.imshow(image)
+fig, ax = plt.subplots()
+ax.imshow(image, cmap=plt.cm.gray)
 
 for props in regions:
     y0, x0 = props.centroid
@@ -36,15 +37,14 @@ for props in regions:
     x2 = x0 - math.sin(orientation) * 0.5 * props.minor_axis_length
     y2 = y0 - math.cos(orientation) * 0.5 * props.minor_axis_length
 
-    plt.plot((x0, x1), (y0, y1), '-r', linewidth=2.5)
-    plt.plot((x0, x2), (y0, y2), '-r', linewidth=2.5)
-    plt.plot(x0, y0, '.g', markersize=15)
+    ax.plot((x0, x1), (y0, y1), '-r', linewidth=2.5)
+    ax.plot((x0, x2), (y0, y2), '-r', linewidth=2.5)
+    ax.plot(x0, y0, '.g', markersize=15)
 
     minr, minc, maxr, maxc = props.bbox
     bx = (minc, maxc, maxc, minc, minc)
     by = (minr, minr, maxr, maxr, minr)
-    plt.plot(bx, by, '-b', linewidth=2.5)
+    ax.plot(bx, by, '-b', linewidth=2.5)
 
-plt.gray()
-plt.axis((0, 600, 600, 0))
+ax.axis((0, 600, 600, 0))
 plt.show()
diff --git a/doc/examples/plot_segmentations.py b/doc/examples/plot_segmentations.py
index 55bb6082..931e92d5 100644
--- a/doc/examples/plot_segmentations.py
+++ b/doc/examples/plot_segmentations.py
@@ -79,7 +79,7 @@ print("Quickshift number of segments: %d" % len(np.unique(segments_quick)))
 
 fig, ax = plt.subplots(1, 3)
 fig.set_size_inches(8, 3, forward=True)
-plt.subplots_adjust(0.05, 0.05, 0.95, 0.95, 0.05, 0.05)
+fig.subplots_adjust(0.05, 0.05, 0.95, 0.95, 0.05, 0.05)
 
 ax[0].imshow(mark_boundaries(img, segments_fz))
 ax[0].set_title("Felzenszwalbs's method")
diff --git a/doc/examples/plot_skeleton.py b/doc/examples/plot_skeleton.py
index a0bfd267..a4a3a60c 100644
--- a/doc/examples/plot_skeleton.py
+++ b/doc/examples/plot_skeleton.py
@@ -47,19 +47,17 @@ image[circle2] = 0
 skeleton = skeletonize(image)
 
 # display results
-plt.figure(figsize=(8, 4.5))
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4.5))
 
-plt.subplot(121)
-plt.imshow(image, cmap=plt.cm.gray)
-plt.axis('off')
-plt.title('original', fontsize=20)
+ax1.imshow(image, cmap=plt.cm.gray)
+ax1.axis('off')
+ax1.set_title('original', fontsize=20)
 
-plt.subplot(122)
-plt.imshow(skeleton, cmap=plt.cm.gray)
-plt.axis('off')
-plt.title('skeleton', fontsize=20)
+ax2.imshow(skeleton, cmap=plt.cm.gray)
+ax2.axis('off')
+ax2.set_title('skeleton', fontsize=20)
 
-plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.98,
+fig.subplots_adjust(wspace=0.02, hspace=0.02, top=0.98,
                     bottom=0.02, left=0.02, right=0.98)
 
 plt.show()
diff --git a/doc/examples/plot_ssim.py b/doc/examples/plot_ssim.py
index 4e657150..319c72fe 100644
--- a/doc/examples/plot_ssim.py
+++ b/doc/examples/plot_ssim.py
@@ -45,7 +45,7 @@ def mse(x, y):
 img_noise = img + noise
 img_const = img + abs(noise)
 
-f, (ax0, ax1, ax2) = plt.subplots(nrows=1, ncols=3, figsize=(8, 4))
+fig, (ax0, ax1, ax2) = plt.subplots(nrows=1, ncols=3, figsize=(8, 4))
 
 mse_none = mse(img, img)
 ssim_none = ssim(img, img, dynamic_range=img.max() - img.min())
diff --git a/doc/examples/plot_swirl.py b/doc/examples/plot_swirl.py
index 98226a20..8f79db94 100644
--- a/doc/examples/plot_swirl.py
+++ b/doc/examples/plot_swirl.py
@@ -74,7 +74,7 @@ from skimage.transform import swirl
 image = data.checkerboard()
 swirled = swirl(image, rotation=0, strength=10, radius=120, order=2)
 
-f, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 3))
+fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 3))
 
 ax0.imshow(image, cmap=plt.cm.gray, interpolation='none')
 ax0.axis('off')
diff --git a/doc/examples/plot_view_as_blocks.py b/doc/examples/plot_view_as_blocks.py
index c58a74a6..2237e50b 100644
--- a/doc/examples/plot_view_as_blocks.py
+++ b/doc/examples/plot_view_as_blocks.py
@@ -60,5 +60,5 @@ ax2.imshow(max_view, cmap=cm.Greys_r)
 ax3.set_title("Block view with\n local median pooling")
 ax3.imshow(median_view, cmap=cm.Greys_r)
 
-plt.subplots_adjust(hspace=0.4, wspace=0.4)
+fig.subplots_adjust(hspace=0.4, wspace=0.4)
 plt.show()
diff --git a/doc/examples/plot_watershed.py b/doc/examples/plot_watershed.py
index f003c122..27f7016e 100644
--- a/doc/examples/plot_watershed.py
+++ b/doc/examples/plot_watershed.py
@@ -61,6 +61,6 @@ ax2.set_title('Separated objects')
 for ax in axes:
     ax.axis('off')
 
-plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
+fig.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
                     right=1)
 plt.show()