How to perform 2D image rotation by FFT based image translation using matlab?

Question

I am studying rotation using 3 translations as described in this scientific paper. According to the authors , it is possible to decompose a rotation into 3 shear translations which will be implemented in the frequency domain by FFT according to time-frequency shifting theorems.

For example the first translation is obtained as in equation (13) by doing the following operation for every pixel considering that the horizental coordinate x is variable and the vertical coordinate y is fixed:

Gx(x,y)=FFT-1(exp(-2*piau*y)*FFT(f(x,y)) where f is the value of the pixel x,y , a is a constant for the translation and FFT-1 is the inverse fourrier transform and u is x mapped into the frequency domain.

To code this transformation, one should think as follow: We have to compute Gx by changing the range of U from -M/(2*M) to (M/2*M) and by iterating y from 1 to N.

Surprisingly, the suggested implementation doesn't follow the procedure above: Can somebody develop further why the implementation is as follow.

close all;
clear all;
I = double(imread('Put image here'))/255;
%imshow(I)
I=rgb2gray(I);
facteur=4;
[M, N, s] = size (I); 


temp=ones(2*M,2*N);
temp(floor(M/2):floor(M/2)+M-1,floor(N/2):floor(N/2)+N-1)=I(:,:);
I=temp; clear temp;
imshow(I)
[M N s]=size(I)
 teta=pi/8;
a=tan(teta/2)
b=-sin(teta)

 Nx = ifftshift(-fix(M/2):ceil(M/2)-1); 
 Ny = ifftshift(-fix(N/2):ceil(N/2)-1);
 Ix=zeros(length(Nx),length(Ny));
 for k=1:M
   Ix(k,:)=(ifft(fft(I(k,:)).*exp(-2*i*pi*(k-floor(M/2))*Ny*a/N)));
 end
 %imshow(Ix);

Answer 1

The following is the completed code that yields the expected result:

close all; clear all; clc;

I = double(imread('Cameraman.tif'))/255;
figure,imshow(I)

[M, N, s] = size (I); 

temp=ones(2*M,2*N);
temp(floor(M/2):floor(M/2)+M-1,floor(N/2):floor(N/2)+N-1)=I;
I=temp; clear temp;

figure,imshow(I)
[M N s]=size(I)
teta = pi/3;
a = tan(teta/2)
b = -sin(teta)

 Nx = ifftshift(-fix(M/2):ceil(M/2)-1); 
 Ny = ifftshift(-fix(N/2):ceil(N/2)-1);
 Ix=zeros(length(Nx),length(Ny));
 for k=1:M
   Ix(k,:)=real(ifft(fft(I(k,:)).*exp(-2*i*pi*(k-floor(M/2))*Ny*a/N)));
 end

 Iy = zeros(size(Ix));
 for k=1:M
   Iy(:,k)=real(ifft(fft(Ix(:,k)).*exp(2*i*pi*(k-floor(M/2))*Ny*b/N)'));
 end

 If = zeros(size(Iy));
 for k=1:M
   If(k,:)=real(ifft(fft(Iy(k,:)).*exp(-2*i*pi*(k-floor(M/2))*Ny*a/N)));
 end

 figure,imshow(Ix);
 figure,imshow(Iy);
 figure,imshow(If);

So this works as expected from a 2D rotation about z-axis. The paper argues that this is a fast and accurate method, but I doubt that, as a simple 2D rotation could easily be implemented using raw trigonometric computations and a bilinear interpolation (or bicubic if a sharper result is desired) other than forward and inverse FFTs.

I can say that the method is slow (on matlab, it can be further vectorized however). On the other hand the perceptual quality of the resulting rotated image is fairly good except some little ringing artefacts towards the edges.

I suggest you the ch-5 of Principles of Computer Graphics 2E by FOLEY for further understanding of the mathematics used in these algorithms and applications.