Image rotation and scaling the frequency domain?

Question

I'm writing some code to recover the rotation, scaling and translation of a test image relative to a template using phase correlation, a la Reddy & Chatterji 1996. I take the FFT of my original test image in order to find the scale factor and angle of rotation, but I then need the FFT of the rotated and scaled test image in order to get the translation.

Now I could apply rotation and scaling in the spatial domain then take the FFT, but that seems a bit inefficient - is it possible to obtain the Fourier coefficients of the rotated/scaled image directly in the frequency domain?

Edit 1: OK, I had a play around following user1816548's suggestion. I can get vaguely sensible-looking rotations for angles that are multiples of 90o, albeit with odd changes in the polarity of the image. Angles that aren't multiples of 90o give me pretty zany results.

Edit 2: I've applied zero padding to the image, and I'm wrapping the edges of the FFT when I rotate it. I'm quite certain that I'm rotating about the DC component of the FFT, but I still get weird results for angles that aren't multiples of 90o.

Example output:

---

Executable Numpy/Scipy code:

---

import numpy as np
from scipy.misc import lena
from scipy.ndimage.interpolation import rotate,zoom
from scipy.fftpack import fft2,ifft2,fftshift,ifftshift
from matplotlib.pyplot import subplots,cm

def testFourierRotation(angle):

    M = lena()
    newshape = [2*dim for dim in M.shape]
    M = procrustes(M,newshape)

    # rotate, then take the FFT
    rM = rotate(M,angle,reshape=False)
    FrM = fftshift(fft2(rM))

    # take the FFT, then rotate
    FM = fftshift(fft2(M))
    rFM = rotatecomplex(FM,angle,reshape=False)
    IrFM = ifft2(ifftshift(rFM))

    fig,[[ax1,ax2,ax3],[ax4,ax5,ax6]] = subplots(2,3)

    ax1.imshow(M,interpolation='nearest',cmap=cm.gray)
    ax1.set_title('Original')
    ax2.imshow(rM,interpolation='nearest',cmap=cm.gray)
    ax2.set_title('Rotated in spatial domain')
    ax3.imshow(abs(IrFM),interpolation='nearest',cmap=cm.gray)
    ax3.set_title('Rotated in Fourier domain')
    ax4.imshow(np.log(abs(FM)),interpolation='nearest',cmap=cm.gray)
    ax4.set_title('FFT')
    ax5.imshow(np.log(abs(FrM)),interpolation='nearest',cmap=cm.gray)
    ax5.set_title('FFT of spatially rotated image')
    ax6.imshow(np.log(abs(rFM)),interpolation='nearest',cmap=cm.gray)
    ax6.set_title('Rotated FFT')
    fig.tight_layout()

    pass

def rotatecomplex(a,angle,reshape=True):
    r = rotate(a.real,angle,reshape=reshape,mode='wrap')
    i = rotate(a.imag,angle,reshape=reshape,mode='wrap')
    return r+1j*i

def procrustes(a,target,padval=0):
    b = np.ones(target,a.dtype)*padval
    aind = [slice(None,None)]*a.ndim
    bind = [slice(None,None)]*a.ndim
    for dd in xrange(a.ndim):
        if a.shape[dd] > target[dd]:
            diff = (a.shape[dd]-target[dd])/2.
            aind[dd] = slice(np.floor(diff),a.shape[dd]-np.ceil(diff))
        elif a.shape[dd] < target[dd]:
            diff = (target[dd]-a.shape[dd])/2.
            bind[dd] = slice(np.floor(diff),target[dd]-np.ceil(diff))
    b[bind] = a[aind]
    return b

Answer 1

I am not sure if this has been resolved yet or not, but I believe I have the solution to your problem regarding the observed effect in your third figure:

This weird effect you observe is due to the origin from which you actually calculate the FFT. Essentially, FFT starts in at the very first pixel of the array at M[0][0] . However, you define your rotation around M[size/2+1,size/2+1] , which is the right way, but wrong :) . The Fourier domain has been calculated from M[0][0]! If you now rotate in Fourier domain, you are rotating around M[0][0] and not around M[size/2+1,size/2+1]. I can't fully explain what's really happening here, but you get same effect I used to get, too. In order to rotate the original image in Fourier domain you must first apply the 2D fftShift to the original image M, then calculate the FFT, rotate, IFFT and then apply ifftShift . This way your rotation center of the image and the center of the Fourier domain come into sync.

AFAI remember we were also rotating real and imaginary components in two separate arrays and merged them afterwards. We also tested various interpolation algorithms on the complex numbers with not much effect :) . It's in our package pytom.

However, this may be super los less, but with the two additional shifts not really fast, unless you specify some funky array index arithmetic.

Answer 2

Well, the rotated and scaled image results in a rotated and scaled (with inverse scale) fourier transform.

Also note that rotation and scaling are both linear in the number of pixels, whereas FFT is O(w*logw*h*logh) so it's actually not that expensive in the end.