I’d like to apply a homography matrix to transform a very large (50kx50k pixels). On the small scale, both cv2 and skimage apply a matrix in a predictable way. However, SciPy (which has large image support in Dask Image) doesn’t transform the image in a predictable way and I can’t work out why. If I could get SciPy to work, I can easily affine transform my large image with Dask.

Strangely, SciPy seems to correctly rotate the image but not scale. If the affine matrix is inverted, the scaling is correct but the rotation is inverted. Here’s example code to replicate:

import scipy
import cv2
import skimage
import numpy as np
import matplotlib.pyplot as plt

# Create the homography matrix
translation = np.array([[1, 0, 40],  # Translate by (0.2, 0.1)
                        [0, 1, -48],
                        [0, 0, 1]])

# Rotate by 45 degrees
rotation_angle = np.deg2rad(45)  
rotation = np.array([[np.cos(rotation_angle), -np.sin(rotation_angle), 0],
                     [np.sin(rotation_angle), np.cos(rotation_angle), 1],
                     [0, 0, 1]])

# Scale by 2x
scaling = np.array([[2, 0, 0],  
                    [0, 2, 0],
                    [0, 0, 1]])

# Combine transformations
affine_matrix = translation @ rotation @ scaling

# Apply transformations
image = skimage.data.astronaut()
skimage_result = skimage.transform.warp(
        image,
        np.linalg.inv(affine_matrix),
    )
cv2_result = cv2.warpAffine(image, 
                        affine_matrix[:2],
                        dsize=(image.shape[1], image.shape[0])
                        )
scipy_result = []
scipy_result_inverted = []
for ch in range(image.shape[2]):
     # Apply to each RGB channel independently
    loop_result = scipy.ndimage.affine_transform(
            image[..., ch],
            affine_matrix,
            order=0,
        )
    loop_result_inverted = scipy.ndimage.affine_transform(
            image[..., ch],
            np.linalg.inv(affine_matrix),
            order=0,
        )
    scipy_result.append(loop_result)
    scipy_result_inverted.append(loop_result_inverted)

scipy_result = np.stack(scipy_result).transpose(1, 2, 0)
scipy_result_inverted = np.stack(scipy_result_inverted).transpose(1, 2, 0)

fig, ax = plt.subplots(1, 5, figsize=(20, 10))
ax[0].imshow(image)
ax[1].imshow(skimage_result)
ax[2].imshow(cv2_result)
ax[3].imshow(scipy_result)
ax[4].imshow(scipy_result_inverted)

ax[1].set_title("skimage.transform.warp")
ax[2].set_title("cv2.warpAffine")
ax[3].set_title("scipy.ndimage.affine_transform")
ax[4].set_title("scipy.ndimage.affine_transformnInverted matrix")

Example transformations plot