function im2 = rotateImage(im, theta, method)
%ROTATEIMAGE Image Rotation
% for a square image of even size n, this function creates a rotated 
% version of that image such that the rows are aligned with the direction 
% (cos(theta),sin(theta)). It is therefore a rotation of angle (-theta).
% 
% Assumes that the image is supported inside the circle enclosed into the 
% square image.
% INPUT ARGS:  im     : image to be rotated
%              theta  : clockwise rotation angle
%              method : rotation method, see possible values below
%
% OUTPUT ARGS: im2    : rotated image of the same size as im
%
% possible values for 'method' variable:
% - '4pt': 4-point interpolation
% - 'spectral': based on shearings and resamplings, this method has
% spectral accuracy

% Author: Francois Monard
% Date: 06/21/2011

theta = mod(theta,2*pi);

if theta<=pi/4,
    im2 = rotateImageElem(im,theta,method);
elseif theta<=pi/2,
    im2 = flipud(rotateImageElem(flipud(rot90(im,-1)),pi/2-theta,method));
elseif theta<=3*pi/4,
    im2 = rotateImageElem(rot90(im,-1),theta-pi/2,method);
elseif theta<=pi,
    im2 = flipud(rotateImageElem(flipud(rot90(im,-2)),pi-theta,method));
elseif theta<=5*pi/4
    im2 = rotateImageElem(rot90(im,-2),theta-pi,method);
elseif theta<=3*pi/2
    im2 = flipud(rotateImageElem(flipud(rot90(im,-3)),3*pi/2-theta,method));
elseif theta<=7*pi/4
    im2 = rotateImageElem(rot90(im,-3),theta-3*pi/2,method);
else
    im2 = flipud(rotateImageElem(flipud(im),2*pi-theta,method));
end

end  

function im2 = rotateImageElem(im, theta, method)

if theta==0
    im2 = im;
elseif theta<=pi/4
    switch method
        case 'spectral'
            im2 = slantResample(im,tan(theta),1/cos(theta));
            im2 = rot90(slantResample(rot90(im2),sin(theta),cos(theta)),-1);
        case '4pt'
            im2 = rotateImage4pt(im,theta);
        otherwise
            error('Invalid rotation method.')
    end
else
    error(['theta = ',num2str(theta),' is not a valid angle.']);
end

end

% elementary rotation based on 4-point interpolation
function im2 = rotateImage4pt(im, theta)

im2 = zeros(size(im));

if norm(im,'fro')==0
    return
else
    n = size(im,1);
    xs = (-1+1/n):2/n:(1.-1/n);
    ys = xs';
    cth = cos(theta); sth = sin(theta);
    
    for a=1:n,
        for b=1:n,
            if xs(b)^2 + ys(a)^2 < 1.
                x_rs = n/2. * (xs(b) * cth + ys(a) * sth + 1.-1./n);
                y_rs = n/2. * (- xs(b) * sth + ys(a) * cth + 1.-1./n);
                i0 = max(min(floor(y_rs), n-2), 0.);
                j0 = max(min(floor(x_rs), n-2), 0.);
                xbar = x_rs - j0;
                ybar = y_rs - i0;
                i0 = i0+1; j0=j0+1;
                im2(a,b) = im(i0, j0) * (1.-xbar) * (1.-ybar) ...
                    + im(i0+1, j0) * ybar * (1.-xbar) ...
                    + im(i0, j0+1) * xbar * (1.-ybar) ...
                    + im(i0+1, j0+1) * xbar * ybar;
                else
                    im2(a,b) = 0.;
            end
        end
    end

    % vectorized version to be fixed
%     [X,Y] = meshgrid(xs,ys);
%     X_rs = n/2. * (X * cth + Y * sth + 1.-1./n);
%     Y_rs = n/2. * (- X * sth + Y * cth + 1.-1./n);
%     i0 = zeros(n,n); j0 = zeros(n,n);
%     for a = 1:n,
%         for b =1:n
%             i0(a,b) = max(min(floor(Y_rs(a,b)), n-2), 0.)+1;
%             j0(a,b) = max(min(floor(X_rs(a,b)), n-2), 0.)+1;
%         end
%     end
%     X_bar = X_rs - j0; Y_bar = Y_rs - i0;
%     
%     im2 = im(i0, j0) .* (1.-X_bar) .* (1.-Y_bar) ...
%         + im(i0+1, j0) .* Y_bar .* (1.-X_bar) ...
%         + im(i0, j0+1) .* X_bar .* (1.-Y_bar) ...
%         + im(i0+1, j0+1) .* X_bar .* Y_bar;
end

end

% slant and resample altogether
function res = slantResample(im,t,h)

res = zeros(size(im));

if norm(im,'fro')==0
    return
else
    [nr,nc] = size(im);
  
    shifts = (-1+1/nc:2/nc:1-1/nc)*t/2;
    N = 2*nr;
    %--- slant params
    temp1 = exp(-i*(pi/N)*(-N+1)*(0:N-1)');
    % temp2 = conj(temp1);
    
    %--- resampling params (precomputations for the FRFT of coeff -h/N
    z = exp(-i*pi*(0:N).^2*h/N); z = [z,z(N:-1:2)]; Fz = fft(z).';
    factors = exp(i*pi*h/N*(0:N-1).^2).'; 
    ns = nr-1/2 - h*(nr/2-1/2);
    
    % preproc = exp(-i*(pi/N)*(-N+1)*(0:N-1).');
    % temp2.*preproc = 1
    midproc = factors.*exp(i*(pi/N)*((-N+1):2:(N-1)).'*ns);
    postproc = factors(1:nr).*exp(i*(pi/N)*(-N+1)*h*(0:(N/2-1)).');
  
    phases = -i*pi/2*((-N+1):2:(N-1))';
    
    for c=1:nc,
        % slant
        colTempN = [zeros(nr/2,1); im(:,c); zeros(nr/2,1)].*temp1;
        colTempN = fft(colTempN).*exp(shifts(c)*phases);
%         colTemp = ifft(colTempN).*temp2;
        
        % resampling
               
%         colTemp = ifft(fft(fft(colTempN).*midproc,2*N).*Fz);
        colTemp2N = ifft(fft(colTempN.*midproc,2*N).*Fz);
        res(:,c) = 1/N*postproc.*colTemp2N(1:nr);
    end
    
    if all(all(isreal(im)))
        res = real(res);
    end
    
    %% TODO : try out the decimation scheme in C
end   
    
end