Speedup 2D warping for affine transformations (scikit-image#2902)

* Speedup 2D warping for affine transformations
* Use exact affine equality check
* Fix syntax error in if condition
* Fix not setting output argument in warping
* Check for valid order argument
* Increase speed for metric warping

skimage/transform/_warps.py

@@ -153,6 +153,11 @@ def resize(image, output_shape, order=1, mode='reflect', cval=0, clip=True,
             tform = AffineTransform()
             tform.estimate(src_corners, dst_corners)
 
+        # Make sure the transform is exactly metric, to ensure fast warping.
+        tform.params[2] = (0, 0, 1)
+        tform.params[0, 1] = 0
+        tform.params[1, 0] = 0
+
         out = warp(image, tform, output_shape=output_shape, order=order,
                    mode=mode, cval=cval, clip=clip,
                    preserve_range=preserve_range)
@@ -373,6 +378,9 @@ def rotate(image, angle, resize=False, center=None, order=1, mode='constant',
         tform4 = SimilarityTransform(translation=translation)
         tform = tform4 + tform
 
+    # Make sure the transform is exactly affine, to ensure fast warping.
+    tform.params[2] = (0, 0, 1)
+
     return warp(image, tform, output_shape=output_shape, order=order,
                 mode=mode, cval=cval, clip=clip, preserve_range=preserve_range)
 

skimage/transform/_warps_cy.pyx

@@ -10,9 +10,9 @@ from .._shared.interpolation cimport (nearest_neighbour_interpolation,
                                       bicubic_interpolation)
 
 
-cdef inline void _matrix_transform(double x, double y, double* H, double *x_,
-                                   double *y_) nogil:
-    """Apply a homography to a coordinate.
+cdef inline void _transform_metric(double x, double y, double* H,
+                                   double *x_, double *y_) nogil:
+    """Apply a metric transformation to a coordinate.
 
     Parameters
     ----------
@@ -24,14 +24,46 @@ cdef inline void _matrix_transform(double x, double y, double* H, double *x_,
         Output coordinate.
 
     """
-    cdef double xx, yy, zz
+    x_[0] = H[0] * x + H[2]
+    y_[0] = H[4] * y + H[5]
 
-    xx = H[0] * x + H[1] * y + H[2]
-    yy = H[3] * x + H[4] * y + H[5]
-    zz =  H[6] * x + H[7] * y + H[8]
 
-    x_[0] = xx / zz
-    y_[0] = yy / zz
+cdef inline void _transform_affine(double x, double y, double* H,
+                                   double *x_, double *y_) nogil:
+    """Apply an affine transformation to a coordinate.
+
+    Parameters
+    ----------
+    x, y : double
+        Input coordinate.
+    H : (3,3) *double
+        Transformation matrix.
+    x_, y_ : *double
+        Output coordinate.
+
+    """
+    x_[0] = H[0] * x + H[1] * y + H[2]
+    y_[0] = H[3] * x + H[4] * y + H[5]
+
+
+cdef inline void _transform_projective(double x, double y, double* H,
+                                       double *x_, double *y_) nogil:
+    """Apply a homography to a coordinate.
+
+    Parameters
+    ----------
+    x, y : double
+        Input coordinate.
+    H : (3,3) *double
+        Transformation matrix.
+    x_, y_ : *double
+        Output coordinate.
+
+    """
+    cdef double z_
+    z_ = H[6] * x + H[7] * y + H[8]
+    x_[0] = (H[0] * x + H[1] * y + H[2]) / z_
+    y_[0] = (H[3] * x + H[4] * y + H[5]) / z_
 
 
 def _warp_fast(cnp.ndarray image, cnp.ndarray H, output_shape=None,
@@ -110,6 +142,15 @@ def _warp_fast(cnp.ndarray image, cnp.ndarray H, output_shape=None,
     cdef Py_ssize_t rows = img.shape[0]
     cdef Py_ssize_t cols = img.shape[1]
 
+    cdef void (*transform_func)(double, double, double*, double*, double*) nogil
+    if M[2, 0] == 0 and M[2, 1] == 0 and M[2, 2] == 1:
+        if M[0, 1] == 0 and M[1, 0] == 0:
+            transform_func = _transform_metric
+        else:
+            transform_func = _transform_affine
+    else:
+        transform_func = _transform_projective
+
     cdef double (*interp_func)(double*, Py_ssize_t, Py_ssize_t, double, double,
                                char, double) nogil
     if order == 0:
@@ -120,11 +161,13 @@ def _warp_fast(cnp.ndarray image, cnp.ndarray H, output_shape=None,
         interp_func = biquadratic_interpolation
     elif order == 3:
         interp_func = bicubic_interpolation
+    else:
+        raise ValueError("Unsupported interpolation order", order)
 
     with nogil:
         for tfr in range(out_r):
             for tfc in range(out_c):
-                _matrix_transform(tfc, tfr, &M[0, 0], &c, &r)
+                transform_func(tfc, tfr, &M[0, 0], &c, &r)
                 out[tfr, tfc] = interp_func(&img[0, 0], rows, cols, r, c,
                                             mode_c, cval)