FFT-based translation, rotation, and scale-invariant image registration

fft_registration.py

"""
B. S. Reddy and B. N. Chatterji, “An FFT-based technique for translation, rotation, and scale-invariant image registration,” IEEE Transactions on Image Processing, vol. 5, no. 8, pp. 1266–1271, Aug. 1996, doi: 10.1109/83.506761.

Implemented in Numpy/OpenCV and Numpy/Scikit-Image

The OpenCV version is about 4 times faster than the Scikit-Image version.
"""
# Numpy/OpenCV implementation
import cv2
from scipy import signal
import scipy.fft as sfft
import numpy as np


# OpenCV implementation
def find_rotation_cv(img_ref: np.ndarray, img_moved: np.ndarray) -> float:
    h, w = img_ref.shape[:2]
    hamming_w = signal.windows.hamming(w)
    hamming_h = signal.windows.hamming(h)
    hamming = np.outer(hamming_h, hamming_w)

    F_ref = np.log(np.abs(sfft.fftshift(sfft.fft2(img_ref * hamming))))
    F_moved = np.log(np.abs(sfft.fftshift(sfft.fft2(img_moved * hamming))))

    center_x = w // 2
    center_y = h // 2
    radius = min(center_x, center_y)

    # Define the desired size of the output polar image
    polar_width = radius
    polar_height = int(np.ceil(radius * np.pi / 2))

    # Perform the polar transformation
    F_ref_warpped = cv2.warpPolar(
        F_ref, 
        (polar_width, polar_height), 
        (center_x, center_y), 
        radius, 
        cv2.WARP_POLAR_LOG + cv2.INTER_CUBIC + cv2.WARP_FILL_OUTLIERS
    )
    F_moved_warpped = cv2.warpPolar(
        F_moved, 
        (polar_width, polar_height), 
        (center_x, center_y), 
        radius, 
        cv2.WARP_POLAR_LOG + cv2.INTER_CUBIC + cv2.WARP_FILL_OUTLIERS
    )

    ret = cv2.phaseCorrelate(F_ref_warpped[:180], F_moved_warpped[:180])
    theta_shift = 360 / polar_height * -ret[0][1]
    return theta_shift



from skimage.transform import warp_polar
from skimage.registration import phase_cross_correlation


# scikit-image implementation
def find_rotation(img_ref: np.ndarray, img_moved: np.ndarray) -> float:
    h, w = img_ref.shape[:2]
    hamming_w = signal.windows.hamming(w)
    hamming_h = signal.windows.hamming(h)
    hamming = np.outer(hamming_h, hamming_w)

    F_ref = np.log(np.abs(sfft.fftshift(sfft.fft2(img_ref * hamming))))
    F_moved = np.log(np.abs(sfft.fftshift(sfft.fft2(img_moved * hamming))))

    F_ref_warpped = warp_polar(F_ref, scaling="log")
    F_moved_warpped = warp_polar(F_moved, scaling="log")

    ret = phase_cross_correlation(F_ref_warpped[:180, ], F_moved_warpped[:180, ])
    return ret[0][0]

Translation is trivially computed through the cross correlation (and therefore Fourier/Phase correlation). Using OpenCV's cv2.phaseCorrelation. Note that phaseCorrelate takes two arrays of type float32 or float64, so you might need to convert the type

import cv2

img1 = cv2.imread("...")
img2 = cv2.imread("...")

# if images are originally uint8 [0-255], convert to double, normalized between [0-1]
img1 /= 255.
img2 /= 255.

# translation is a 2-tuple (x, y)
translation, phase_shift = cv2.phaseCorrelate(img1, img2)