diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd
index 40df6524..ef880109 100644
--- a/skimage/_shared/interpolation.pxd
+++ b/skimage/_shared/interpolation.pxd
@@ -8,8 +8,12 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols,
                                           double r, double c, char mode,
                                           double cval)
 
-cdef inline double cubic_interpolation(double x, double[4] f)
+cdef inline double quadratic_interpolation(double x, double[3] f)
+cdef inline double biquadratic_interpolation(double* image, int rows, int cols,
+                                             double r, double c, char mode,
+                                             double cval)
 
+cdef inline double cubic_interpolation(double x, double[4] f)
 cdef inline double bicubic_interpolation(double* image, int rows, int cols,
                                          double r, double c, char mode,
                                          double cval)
diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx
index ad8aae4b..3150f41b 100644
--- a/skimage/_shared/interpolation.pyx
+++ b/skimage/_shared/interpolation.pyx
@@ -75,8 +75,80 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols,
     return (1 - dr) * top + dr * bottom
 
 
+cdef inline double quadratic_interpolation(double x, double[3] f):
+    """Quadratic interpolation.
+
+    Parameters
+    ----------
+    x : double
+        Position in the interval [-1, 1].
+    f : double[4]
+        Function values at positions [-1, 0, 1].
+
+    Returns
+    -------
+    value : double
+        Interpolated value.
+
+    """
+    return f[1] - 0.25 * (f[0] - f[2]) * x
+
+
+cdef inline double biquadratic_interpolation(double* image, int rows, int cols,
+                                             double r, double c, char mode,
+                                             double cval):
+    """Biquadratic interpolation at a given position in the image.
+
+    Parameters
+    ----------
+    image : double array
+        Input image.
+    rows, cols : int
+        Shape of image.
+    r, c : int
+        Position at which to interpolate.
+    mode : {'C', 'W', 'R', 'N'}
+        Wrapping mode. Constant, Wrap, Reflect or Nearest.
+    cval : double
+        Constant value to use for constant mode.
+
+    Returns
+    -------
+    value : double
+        Interpolated value.
+
+    """
+
+    cdef int r0 = <int>round(r)
+    cdef int c0 = <int>round(c)
+    if r < 0:
+        r0 -= 1
+    if c < 0:
+        c0 -= 1
+    # scale position to range [-1, 1]
+    cdef double xr = (r - r0) - 1
+    cdef double xc = (c - c0) - 1
+    if r == r0:
+        xr += 1
+    if c == c0:
+        xc += 1
+
+    cdef double fc[3], fr[3]
+
+    cdef int pr, pc
+
+    # row-wise cubic interpolation
+    for pr in range(r0, r0 + 3):
+        for pc in range(c0, c0 + 3):
+            fc[pc - c0] = get_pixel(image, rows, cols, pr, pc, mode, cval)
+        fr[pr - r0] = quadratic_interpolation(xc, fc)
+
+    # cubic interpolation for interpolated values of each row
+    return quadratic_interpolation(xr, fr)
+
+
 cdef inline double cubic_interpolation(double x, double[4] f):
-    """Ccubic interpolation.
+    """Cubic interpolation.
 
     Parameters
     ----------
diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 02b6c161..270198aa 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -833,7 +833,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
     bands = ishape[2]
 
     # use fast Cython version for specific interpolation orders
-    if order in (0, 1, 3) and not map_args:
+    if order in range(4) and not map_args:
         matrix = None
         if isinstance(inverse_map, HOMOGRAPHY_TRANSFORMS):
             matrix = inverse_map._matrix
diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx
index bedbfa61..ce400ed6 100644
--- a/skimage/transform/_warps_cy.pyx
+++ b/skimage/transform/_warps_cy.pyx
@@ -7,6 +7,7 @@ cimport numpy as np
 import numpy as np
 from skimage._shared.interpolation cimport (nearest_neighbour_interpolation,
                                             bilinear_interpolation,
+                                            biquadratic_interpolation,
                                             bicubic_interpolation)
 
 
@@ -72,6 +73,7 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
         Order of interpolation::
         * 0: Nearest-neighbour interpolation.
         * 1: Bilinear interpolation (default).
+        * 2: Biquadratic interpolation (default).
         * 3: Bicubic interpolation.
     mode : {'constant', 'reflect', 'wrap'}
         How to handle values outside the image borders.
@@ -113,6 +115,8 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
         interp_func = nearest_neighbour_interpolation
     elif order == 1:
         interp_func = bilinear_interpolation
+    elif order == 2:
+        interp_func = biquadratic_interpolation
     elif order == 3:
         interp_func = bicubic_interpolation