diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py
index c1e7f845..ec81849f 100644
--- a/skimage/feature/__init__.py
+++ b/skimage/feature/__init__.py
@@ -2,6 +2,6 @@ from ._hog import hog
 from .texture import greycomatrix, greycoprops, local_binary_pattern
 from .peak import peak_local_max
 from .corner import (corner_kitchen_rosenfeld, corner_harris, corner_shi_tomasi,
-                     corner_foerstner)
+                     corner_foerstner, corner_subpix)
 from .corner_cy import corner_moravec
 from .template import match_template
diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py
index 4c1dd89c..9c7a2fdf 100644
--- a/skimage/feature/corner.py
+++ b/skimage/feature/corner.py
@@ -1,6 +1,8 @@
 import numpy as np
 from scipy import ndimage
+from scipy import stats
 from skimage.color import rgb2grey
+from skimage.util import img_as_float
 from . import peak
 
 
@@ -19,8 +21,8 @@ def _compute_derivatives(image):
 
     """
 
-    imx = ndimage.sobel(image, axis=0, mode='constant', cval=0)
-    imy = ndimage.sobel(image, axis=1, mode='constant', cval=0)
+    imy = ndimage.sobel(image, axis=0, mode='constant', cval=0)
+    imx = ndimage.sobel(image, axis=1, mode='constant', cval=0)
 
     return imx, imy
 
@@ -44,7 +46,7 @@ def _compute_auto_correlation(image, sigma):
     """
 
     if image.ndim == 3:
-        image = rgb2grey(image)
+        image = img_as_float(rgb2grey(image))
 
     imx, imy = _compute_derivatives(image)
 
@@ -279,3 +281,123 @@ def corner_foerstner(image, sigma=1):
     q = 4 * detA / traceA**2
 
     return w, q
+
+
+def corner_subpix(image, corners, window_size=11, alpha=0.99):
+    """Determine subpixel position of corners.
+
+    Parameters
+    ----------
+    image : ndarray
+        Input image.
+    corners : (N, 2) ndarray
+        Corner coordinates `(row, cols)`.
+    window_size : int, optional
+        Search window size for subpixel estimation.
+    alpha : float, optional
+        Significance level for point classification.
+
+    Returns
+    -------
+    positions : (N, 2) ndarray
+        Subpixel corner positions. NaN for "not classified" corners.
+
+    """
+
+    # window extent in one direction
+    wext = (window_size - 1) / 2
+
+    # normal equation arrays
+    N_dot = np.zeros((2, 2), dtype=np.double)
+    N_edge = np.zeros((2, 2), dtype=np.double)
+    b_dot = np.zeros((2, ), dtype=np.double)
+    b_edge = np.zeros((2, ), dtype=np.double)
+
+    # critical statistical test values
+    redundancy = window_size**2 - 2
+    t_crit_dot = stats.f.isf(1 - alpha, redundancy, redundancy)
+    t_crit_edge = stats.f.isf(alpha, redundancy, redundancy)
+
+    # coordinates of pixels within window
+    y, x = np.mgrid[- wext:wext + 1, - wext:wext + 1]
+
+    # output arrays
+    corners_subpix = np.zeros_like(corners, dtype=np.double)
+    corners_class = np.zeros(corners.shape[0], dtype=np.int8)
+
+    for i, (y0, x0) in enumerate(corners):
+
+        # crop window around corner + border for sobel operator
+        miny = y0 - wext - 1
+        maxy = y0 + wext + 2
+        minx = x0 - wext - 1
+        maxx = x0 + wext + 2
+        window = image[miny:maxy, minx:maxx]
+
+        winx, winy = _compute_derivatives(window)
+
+        # compute gradient suares and remove border
+        winx_winx = (winx * winx)[1:-1, 1:-1]
+        winx_winy = (winx * winy)[1:-1, 1:-1]
+        winy_winy = (winy * winy)[1:-1, 1:-1]
+
+        # sum of squared differences (mean instead of gaussian filter)
+        Axx = np.sum(winx_winx)
+        Axy = np.sum(winx_winy)
+        Ayy = np.sum(winy_winy)
+
+        # sum of squared differences weighted with coordinates
+        # (mean instead of gaussian filter)
+        bxx_x = np.sum(winx_winx * x)
+        bxx_y = np.sum(winx_winx * y)
+        bxy_x = np.sum(winx_winy * x)
+        bxy_y = np.sum(winx_winy * y)
+        byy_x = np.sum(winy_winy * x)
+        byy_y = np.sum(winy_winy * y)
+
+        # normal equations for subpixel position
+        N_dot[0, 0] = Axx
+        N_dot[0, 1] = N_dot[1, 0] = - Axy
+        N_dot[1, 1] = Ayy
+
+        N_edge[0, 0] = Ayy
+        N_edge[0, 1] = N_edge[1, 0] = Axy
+        N_edge[1, 1] = Axx
+
+        b_dot[:] = bxx_y - bxy_x, byy_x - bxy_y
+        b_edge[:] = byy_y + bxy_x, bxx_x + bxy_y
+
+        # estimated positions
+        est_dot = np.linalg.solve(N_dot, b_dot)
+        est_edge = np.linalg.solve(N_edge, b_edge)
+
+        # residuals
+        ry_dot = y - est_dot[0]
+        rx_dot = x - est_dot[1]
+        ry_edge = y - est_edge[0]
+        rx_edge = x - est_edge[1]
+        # squared residuals
+        rxx_dot = rx_dot * rx_dot
+        rxy_dot = rx_dot * ry_dot
+        ryy_dot = ry_dot * ry_dot
+        rxx_edge = rx_edge * rx_edge
+        rxy_edge = rx_edge * ry_edge
+        ryy_edge = ry_edge * ry_edge
+
+        # determine corner class (dot or edge)
+        w_dot = np.sum(winx_winx * ryy_dot - 2 * winx_winy * rxy_dot \
+                       + winy_winy * rxx_dot)
+        w_edge = np.sum(winy_winy * ryy_edge + 2 * winx_winy * rxy_edge \
+                        + winx_winx * rxx_edge)
+        t = w_edge / w_dot
+        # 1 for edge, -1 for dot, 0 for "not classified"
+        corner_class = (t < t_crit_edge) - (t > t_crit_dot)
+
+        if corner_class == - 1:
+            corners_subpix[i, :] = y0 + est_dot[0], x0 + est_dot[1]
+        elif corner_class == 0:
+            corners_subpix[i, :] = np.nan, np.nan
+        elif corner_class == 1:
+            corners_subpix[i, :] = y0 + est_edge[0], x0 + est_edge[1]
+
+    return corners_subpix