diff --git a/doc/examples/plot_matching.py b/doc/examples/plot_matching.py
index bae6d6e2..9fa0b41c 100644
--- a/doc/examples/plot_matching.py
+++ b/doc/examples/plot_matching.py
@@ -27,7 +27,8 @@ from matplotlib import pyplot as plt
 
 from skimage import data
 from skimage.util import img_as_float
-from skimage.feature import corner_harris, corner_subpix, corner_peaks
+from skimage.feature import (corner_harris, corner_subpix, corner_peaks,
+                             plot_matches)
 from skimage.transform import warp, AffineTransform
 from skimage.exposure import rescale_intensity
 from skimage.color import rgb2gray
@@ -117,28 +118,21 @@ print(tform.scale, tform.translation, tform.rotation)
 print(model.scale, model.translation, model.rotation)
 print(model_robust.scale, model_robust.translation, model_robust.rotation)
 
-
-# visualize correspondences
-img_combined = np.concatenate((img_orig_gray, img_warped_gray), axis=1)
-
+# visualize correspondence
 fig, ax = plt.subplots(nrows=2, ncols=1)
+
 plt.gray()
 
-ax[0].imshow(img_combined, interpolation='nearest')
+inlier_idxs = np.nonzero(inliers)
+plot_matches(ax[0], img_orig_gray, img_warped_gray, src, dst,
+             inlier_idxs, inlier_idxs, matches_color='b')
 ax[0].axis('off')
-ax[0].axis((0, 400, 200, 0))
 ax[0].set_title('Correct correspondences')
-ax[1].imshow(img_combined, interpolation='nearest')
+
+outlier_idxs = np.nonzero(outliers)
+plot_matches(ax[1], img_orig_gray, img_warped_gray, src, dst,
+             outlier_idxs, outlier_idxs, matches_color='r')
 ax[1].axis('off')
-ax[1].axis((0, 400, 200, 0))
 ax[1].set_title('Faulty correspondences')
 
-
-for ax_idx, (m, color) in enumerate(((inliers, 'g'), (outliers, 'r'))):
-    ax[ax_idx].plot((src[m, 1], dst[m, 1] + 200), (src[m, 0], dst[m, 0]), '-',
-                    color=color)
-    ax[ax_idx].plot(src[m, 1], src[m, 0], '.', markersize=10, color=color)
-    ax[ax_idx].plot(dst[m, 1] + 200, dst[m, 0], '.', markersize=10,
-                    color=color)
-
 plt.show()
diff --git a/skimage/feature/util.py b/skimage/feature/util.py
index 2d33261f..a27ab1ba 100644
--- a/skimage/feature/util.py
+++ b/skimage/feature/util.py
@@ -36,7 +36,6 @@ class DescriptorExtractor(object):
 def plot_matches(ax, image1, image2, keypoints1, keypoints2,
                  indices1, indices2, keypoints_color='k', matches_color=None,
                  only_matches=False):
-
     """Plot matched features.
 
     Parameters
@@ -70,6 +69,9 @@ def plot_matches(ax, image1, image2, keypoints1, keypoints2,
     image1 = img_as_float(image1)
     image2 = img_as_float(image2)
 
+    indices1 = np.squeeze(indices1)
+    indices2 = np.squeeze(indices2)
+
     new_shape1 = image1.shape
     new_shape2 = image2.shape