diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd
index fe43d2a0..088be1ef 100644
--- a/skimage/_shared/interpolation.pxd
+++ b/skimage/_shared/interpolation.pxd
@@ -81,7 +81,9 @@ cdef inline double bilinear_interpolation(double* image, Py_ssize_t rows,
 
 
 cdef inline double quadratic_interpolation(double x, double[3] f):
-    """Quadratic interpolation.
+    """WARNING: Do not use, not implemented correctly.
+
+    Quadratic interpolation.
 
     Parameters
     ----------
@@ -104,7 +106,9 @@ cdef inline double quadratic_interpolation(double x, double[3] f):
 cdef inline double biquadratic_interpolation(double* image, Py_ssize_t rows,
                                              Py_ssize_t cols, double r,
                                              double c, char mode, double cval):
-    """Biquadratic interpolation at a given position in the image.
+    """WARNING: Do not use, not implemented correctly.
+
+    Biquadratic interpolation at a given position in the image.
 
     Parameters
     ----------
@@ -129,7 +133,6 @@ cdef inline double biquadratic_interpolation(double* image, Py_ssize_t rows,
     cdef long r0 = <long>round(r) - 1
     cdef long c0 = <long>round(c) - 1
 
-    # scale position to range [0, 2]
     cdef double xr = r - r0
     cdef double xc = c - c0
 
@@ -139,10 +142,11 @@ cdef inline double biquadratic_interpolation(double* image, Py_ssize_t rows,
     cdef long pr, pc
 
     # row-wise cubic interpolation
-    for pr in range(r0, r0 + 3):
-        for pc in range(c0, c0 + 3):
-            fc[pc - c0] = get_pixel2d(image, rows, cols, pr, pc, mode, cval)
-        fr[pr - r0] = quadratic_interpolation(xc, fc)
+    for pr in range(3):
+        for pc in range(3):
+            fc[pc] = get_pixel2d(image, rows, cols,
+                                 r0 + pr, c0 + pc, mode, cval)
+        fr[pr] = quadratic_interpolation(xc, fc)
 
     # cubic interpolation for interpolated values of each row
     return quadratic_interpolation(xr, fr)
@@ -221,10 +225,10 @@ cdef inline double bicubic_interpolation(double* image, Py_ssize_t rows,
     cdef long pr, pc
 
     # row-wise cubic interpolation
-    for pr in range(r0, r0 + 4):
-        for pc in range(c0, c0 + 4):
-            fc[pc - c0] = get_pixel2d(image, rows, cols, pr, pc, mode, cval)
-        fr[pr - r0] = cubic_interpolation(xc, fc)
+    for pr in range(4):
+        for pc in range(4):
+            fc[pc] = get_pixel2d(image, rows, cols, pr + r0, pc + c0, mode, cval)
+        fr[pr] = cubic_interpolation(xc, fc)
 
     # cubic interpolation for interpolated values of each row
     return cubic_interpolation(xr, fr)
diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 78387cf1..456570f2 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -1134,7 +1134,18 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
 
     out = None
 
-    if order in range(4) and not map_args:
+    if order == 2:
+        # When fixing this issue, make sure to fix the branches further
+        # below in this function
+        warnings.warn("Bi-quadratic interpolation behavior has changed due "
+                      "to a bug in the implementation of scikit-image. "
+                      "The new version now serves as a wrapper "
+                      "around SciPy's interpolation functions, which itself "
+                      "is not verified to be a correct implementation. Until "
+                      "skimage's implementation is fixed, we recommend "
+                      "to use bi-linear or bi-cubic interpolation instead.")
+
+    if order in (0, 1, 3) and not map_args:
         # use fast Cython version for specific interpolation orders and input
 
         matrix = None