diff --git a/doc/examples/plot_edge_filter.py b/doc/examples/plot_edge_filter.py
index 7eba57df..0bb6c567 100644
--- a/doc/examples/plot_edge_filter.py
+++ b/doc/examples/plot_edge_filter.py
@@ -8,10 +8,11 @@ They are discrete differentiation operators, computing an approximation of the
 gradient of the image intensity function.
 
 """
+import numpy as np
 import matplotlib.pyplot as plt
 
 from skimage.data import camera
-from skimage.filters import roberts, sobel
+from skimage.filters import roberts, sobel, scharr
 
 
 image = camera()
@@ -28,4 +29,55 @@ ax1.imshow(edge_sobel, cmap=plt.cm.gray)
 ax1.set_title('Sobel Edge Detection')
 ax1.axis('off')
 
+plt.tight_layout()
+
+"""
+.. image:: PLOT2RST.current_figure
+
+Different operators compute different finite-difference approximations of the
+gradient. For example, the Scharr filter results in a better rotational
+variance than other filters such as the Sobel filter [1]_ [2]_. The difference
+between the two filters is illustrated below on an image that is the
+discretization of a rotation-invariant continuous function. The discrepancy
+between the two filters is stronger for regions of the image where the
+direction of the gradient is close to diagonal, and for regions with high
+spatial frequencies.
+
+.. [1] http://en.wikipedia.org/wiki/Sobel_operator#Alternative_operators
+
+.. [2] B. Jaehne, H. Scharr, and S. Koerkel. Principles of filter design. In
+       Handbook of Computer Vision and Applications. Academic Press, 1999.
+"""
+
+x, y = np.ogrid[:100, :100]
+# Rotation-invariant image with different spatial frequencies
+img = np.exp(1j * np.hypot(x, y)**1.3 / 20.).real
+
+edge_sobel = sobel(img)
+edge_scharr = scharr(img)
+
+fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(nrows=2, ncols=2)
+
+ax0.imshow(edge_sobel, cmap=plt.cm.gray)
+ax0.set_title('Sobel Edge Detection')
+ax0.axis('off')
+
+ax1.imshow(edge_scharr, cmap=plt.cm.gray)
+ax1.set_title('Scharr Edge Detection')
+ax1.axis('off')
+
+ax2.imshow(img, cmap=plt.cm.gray)
+ax2.set_title('Original image')
+ax2.axis('off')
+
+ax3.imshow(edge_scharr - edge_sobel, cmap=plt.cm.jet)
+ax3.set_title('difference (Scharr - Sobel)')
+ax3.axis('off')
+
+plt.tight_layout()
+
 plt.show()
+
+"""
+.. image:: PLOT2RST.current_figure
+"""
diff --git a/skimage/filters/edges.py b/skimage/filters/edges.py
index 7fb1fb21..e382bc95 100644
--- a/skimage/filters/edges.py
+++ b/skimage/filters/edges.py
@@ -72,14 +72,32 @@ def sobel(image, mask=None):
     output : 2-D array
         The Sobel edge map.
 
+    See also
+    --------
+    scharr, prewitt, roberts, feature.canny
+
     Notes
     -----
     Take the square root of the sum of the squares of the horizontal and
     vertical Sobels to get a magnitude that's somewhat insensitive to
     direction.
 
+    The 3x3 convolution kernel used in the horizontal and vertical Sobels is
+    an approximation of the gradient of the image (with some slight blurring
+    since 9 pixels are used to compute the gradient at a given pixel). As an
+    approximation of the gradient, the Sobel operator is not completely
+    rotation-invariant. The Scharr operator should be used for a better
+    rotation invariance.
+
     Note that ``scipy.ndimage.sobel`` returns a directional Sobel which
     has to be further processed to perform edge detection.
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> camera = data.camera()
+    >>> from skimage import filters
+    >>> edges = filters.sobel(camera)
     """
     assert_nD(image, 2)
     out = np.sqrt(sobel_h(image, mask)**2 + sobel_v(image, mask)**2)
@@ -230,17 +248,30 @@ def scharr(image, mask=None):
     output : 2-D array
         The Scharr edge map.
 
+    See also
+    --------
+    sobel, prewitt, canny
+
     Notes
     -----
     Take the square root of the sum of the squares of the horizontal and
-    vertical Scharrs to get a magnitude that's somewhat insensitive to
-    direction.
+    vertical Scharrs to get a magnitude that is somewhat insensitive to
+    direction. The Scharr operator has a better rotation invariance than
+    other edge filters such as the Sobel or the Prewitt operators.
 
     References
     ----------
     .. [1] D. Kroon, 2009, Short Paper University Twente, Numerical
            Optimization of Kernel Based Image Derivatives.
 
+    .. [2] http://en.wikipedia.org/wiki/Sobel_operator#Alternative_operators
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> camera = data.camera()
+    >>> from skimage import filters
+    >>> edges = filters.scharr(camera)
     """
     out = np.sqrt(scharr_h(image, mask)**2 + scharr_v(image, mask)**2)
     out /= np.sqrt(2)
@@ -410,10 +441,26 @@ def prewitt(image, mask=None):
     output : 2-D array
         The Prewitt edge map.
 
+    See also
+    --------
+    sobel, scharr
+
     Notes
     -----
     Return the square root of the sum of squares of the horizontal
-    and vertical Prewitt transforms.
+    and vertical Prewitt transforms. The edge magnitude depends slightly
+    on edge directions, since the approximation of the gradient operator by
+    the Prewitt operator is not completely rotation invariant. For a better
+    rotation invariance, the Scharr operator should be used. The Sobel operator
+    has a better rotation invariance than the Prewitt operator, but a worse
+    rotation invariance than the Scharr operator.
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> camera = data.camera()
+    >>> from skimage import filters
+    >>> edges = filters.prewitt(camera)
     """
     assert_nD(image, 2)
     out = np.sqrt(prewitt_h(image, mask)**2 + prewitt_v(image, mask)**2)
@@ -563,6 +610,18 @@ def roberts(image, mask=None):
     -------
     output : 2-D array
         The Roberts' Cross edge map.
+
+    See also
+    --------
+    sobel, scharr, prewitt, feature.canny
+
+    Examples
+    --------
+    >>> from skimage import data
+    >>> camera = data.camera()
+    >>> from skimage import filters
+    >>> edges = filters.roberts(camera)
+
     """
     assert_nD(image, 2)
     out = np.sqrt(roberts_pos_diag(image, mask)**2 +