diff --git a/doc/examples/plot_harris.py b/doc/examples/plot_harris.py
index 6212ac4a..39c1f29e 100644
--- a/doc/examples/plot_harris.py
+++ b/doc/examples/plot_harris.py
@@ -13,7 +13,7 @@ import numpy as np
 from matplotlib import pyplot as plt
 
 from skimage import data, img_as_float
-from skimage.feature import harris
+from skimage.feature import harris, peak_local_max
 
 
 def plot_harris_points(image, filtered_coords):
@@ -29,12 +29,12 @@ plt.figure(figsize=(8, 3))
 im_lena = img_as_float(data.lena())
 im_text = img_as_float(data.text())
 
-filtered_coords = harris(im_lena, min_distance=4)
+filtered_coords = peak_local_max(harris(im_lena), min_distance=4)
 
 plt.axes([0, 0, 0.3, 0.95])
 plot_harris_points(im_lena, filtered_coords)
 
-filtered_coords = harris(im_text, min_distance=4)
+filtered_coords = peak_local_max(harris(im_text), min_distance=4)
 
 plt.axes([0.2, 0, 0.77, 1])
 plot_harris_points(im_text, filtered_coords)
diff --git a/skimage/feature/interest.py b/skimage/feature/interest.py
index 15b249bc..88280a27 100644
--- a/skimage/feature/interest.py
+++ b/skimage/feature/interest.py
@@ -1,24 +1,31 @@
 import numpy as np
 from scipy import ndimage
+from skimage.color import rgb2grey
 from . import peak
 
 
-def harris(image, eps=1e-6, sigma=1):
+def harris(image, method='k', k=0.05, eps=1e-6, sigma=1):
     """Compute Harris response image.
 
     Parameters
     ----------
     image : ndarray
         Input image.
+    method : {'k', 'eps'}, optional
+        Method to
+    k : float, optional
+        Sensitivity factor to separate corners from edges, typically in range
+        `[0, 0.2]`. Small values of k result in detection of sharp corners.
     eps : float, optional
-        Normalisation factor.
+        Normalisation factor (Noble's corner measure).
     sigma : float, optional
-        Standard deviation used for the Gaussian kernel.
+        Standard deviation used for the Gaussian kernel, which is used as
+        weighting function for the auto-correlation matrix.
 
     Returns
     -------
     response : ndarray
-        Moravec response image.
+        Harris response image.
 
     Examples
     -------
@@ -48,23 +55,24 @@ def harris(image, eps=1e-6, sigma=1):
         image = rgb2grey(image)
 
     # derivatives
-    image = ndimage.gaussian_filter(image, sigma, mode='constant', cval=0)
     imx = ndimage.sobel(image, axis=0, mode='constant', cval=0)
     imy = ndimage.sobel(image, axis=1, mode='constant', cval=0)
 
-    Wxx = ndimage.gaussian_filter(imx * imx, sigma,
+    Axx = ndimage.gaussian_filter(imx * imx, sigma,
                                   mode='constant', cval=0)
-    Wxy = ndimage.gaussian_filter(imx * imy, sigma,
+    Axy = ndimage.gaussian_filter(imx * imy, sigma,
                                   mode='constant', cval=0)
-    Wyy = ndimage.gaussian_filter(imy * imy, sigma,
+    Ayy = ndimage.gaussian_filter(imy * imy, sigma,
                                   mode='constant', cval=0)
 
-    # determinant and trace
-    Wdet = Wxx * Wyy - Wxy**2
-    Wtr = Wxx + Wyy
+    # determinant
+    detA = Axx * Ayy - Axy**2
+    # trace
+    traceA = Axx + Ayy
 
-    # Alternate formula for Harris response.
-    # Alison Noble, "Descriptions of Image Surfaces", PhD thesis (1989)
-    harris = Wdet / (Wtr + eps)
+    if method == 'k':
+        harris = detA - k * traceA**2
+    else:
+        harris = 2 * detA / (traceA + eps)
 
     return harris
diff --git a/skimage/feature/tests/test_interest.py b/skimage/feature/tests/test_interest.py
index 43bf28a3..22de33c3 100644
--- a/skimage/feature/tests/test_interest.py
+++ b/skimage/feature/tests/test_interest.py
@@ -3,13 +3,13 @@ import numpy as np
 from skimage import data
 from skimage import img_as_float
 
-from skimage.feature import harris
+from skimage.feature import moravec, harris, peak_local_max
 
 
 def test_square_image():
     im = np.zeros((50, 50)).astype(float)
     im[:25, :25] = 1.
-    results = harris(im)
+    results = peak_local_max(harris(im))
     assert results.any()
     assert len(results) == 1
 
@@ -18,7 +18,7 @@ def test_noisy_square_image():
     im = np.zeros((50, 50)).astype(float)
     im[:25, :25] = 1.
     im = im + np.random.uniform(size=im.shape) * .5
-    results = harris(im)
+    results = peak_local_max(harris(im))
     assert results.any()
     assert len(results) == 1
 
@@ -27,7 +27,7 @@ def test_squared_dot():
     im = np.zeros((50, 50))
     im[4:8, 4:8] = 1
     im = img_as_float(im)
-    results = harris(im, min_distance=3)
+    results = peak_local_max(harris(im))
     assert (results == np.array([[6, 6]])).all()
 
 
@@ -37,9 +37,9 @@ def test_rotated_lena():
     rotation.
     """
     im = img_as_float(data.lena().mean(axis=2))
-    results = harris(im)
+    results = peak_local_max(harris(im))
     im_rotated = im.T
-    results_rotated = harris(im_rotated)
+    results_rotated = peak_local_max(harris(im_rotated))
     assert (np.sort(results[:, 0]) == np.sort(results_rotated[:, 1])).all()
     assert (np.sort(results[:, 1]) == np.sort(results_rotated[:, 0])).all()