diff --git a/skimage/filter/harris.py b/skimage/filter/harris.py
index 303431dd..2d48f7ec 100644
--- a/skimage/filter/harris.py
+++ b/skimage/filter/harris.py
@@ -22,7 +22,7 @@ def _compute_harris_response(image, eps=1e-6):
 
     Returns
     --------
-    features: (M, 2) ndarray
+    image: (M, N) ndarray
         Harris image response
     """
     if len(image.shape) == 3:
@@ -85,12 +85,10 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6):
 
     # get coordinates of candidates
     candidates = harrisim_t.nonzero()
-    coords = np.concatenate((candidates[0].reshape((len(candidates[0]), 1)),
-                             candidates[1].reshape((len(candidates[0]), 1))),
-                            axis=1)
+    coords = np.transpose(candidates)
 
     # ...and their values
-    candidate_values = [harrisim[c[0]][c[1]] for c in coords]
+    candidate_values = harrisim[candidates]
 
     # sort candidates
     index = np.argsort(candidate_values)
@@ -103,7 +101,7 @@ def harris(image, min_distance=10, threshold=0.1, eps=1e-6):
     # select the best points taking min_distance into account
     filtered_coords = []
     for i in index:
-        if allowed_locations[coords[i][0]][coords[i][1]] == 1:
+        if allowed_locations[tuple(coords[i])] == 1:
             filtered_coords.append(coords[i])
             allowed_locations[
               (coords[i][0] - min_distance):(coords[i][0] + min_distance),
diff --git a/skimage/filter/tests/test_harris.py b/skimage/filter/tests/test_harris.py
index c78d3e26..853a00a0 100644
--- a/skimage/filter/tests/test_harris.py
+++ b/skimage/filter/tests/test_harris.py
@@ -1,8 +1,10 @@
 import numpy as np
 
-from skimage.filter import harris
+from skimage import data
 from skimage import img_as_float
 
+from skimage.filter import harris
+
 
 class TestHarris():
     def test_square_image(self):
@@ -24,5 +26,17 @@ class TestHarris():
         im = np.zeros((50, 50))
         im[4:8, 4:8] = 1
         im = img_as_float(im)
+        results = harris(im, min_distance=3)
+        print results
+        assert (results == np.array([[6, 6]])).all()
+
+    def test_rotated_lena(self):
+        """
+        The harris filter should yield the same results with an image and it's
+        rotation.
+        """
+        im = img_as_float(data.lena().mean(axis=2))
         results = harris(im)
-        assert results == np.array([6, 6])
+        im_rotated = im.T
+        results_rotated = harris(im_rotated)
+        assert (results[:, 0] == results_rotated[:, 1]).all()