function Y = FTCM_interpolate(X, x0, y0, dx, dy, alpha, type, M, N)
% FTCM_interpolate Make an image by interpolating on a (floor of) tile(s)
%                  as in figure 2 in the EURASIP JASP 2006 FTCM paper 
%
%  use: Y = FTCM_interpolate(X, x0, y0, dx, dy, alpha, type, M, N);
% ----------------------------------------------------------------------------
%   Y     - the resulting image after sampling the tile (the tiled floor), 
%           size is MxN
%   X     - the tile, the corners are: (0,0), (0,1), (1,1) and (1,0)
%   x0    - offset in x-axis for the first point of the grid giving Y
%   y0    - offset in y-axis for the first point of the grid giving Y
%   dx    - increment in x-axis for the grid 
%   dy    - increment in y-axis for the grid 
%   alpha - angle (in radians) the grid is rotated, center of rotation is (x0,y0)
%           negative angle of the grid gives 'positive' rotation of image
%   type  - 0 for nearest neighbor, or 1 for linear interpolation
%   M,N   - Y will be of size MxN (M in x-direction, and N in y direction)
% ----------------------------------------------------------------------------
% ex:
% X = [ 1,0,1,0; 2,0,2,2; 1,0,1,0; 2,2,2,0 ]/2;
% Y = FTCM_interpolate(X, rand(1), rand(1), 0.047, 0.049, 0, 1, 256, 256);
% imagesc(Y); colormap(gray);

%----------------------------------------------------------------------
% Copyright (c) 2002-2013.  Karl Skretting.  All rights reserved.
% University of Stavanger, Signal Processing Group
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.uis.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  04.02.2013  Made function based on older functions 
%                       (makeextexture.m, GenIntPol.m, InterpolateImage.m)
%----------------------------------------------------------------------

c = cos(alpha); s = sin(alpha);
Y1 = (c*dx*(0:(M-1))')*ones(1,N)-ones(M,1)*(s*dy*(0:(N-1)));
Y2 = (s*dx*(0:(M-1))')*ones(1,N)+ones(M,1)*(c*dy*(0:(N-1)));
Y1 = mod(Y1+x0+eps,1);       % range 0-1
Y2 = mod(Y2+y0+eps,1);       % range 0-1
if type    % type~=0  (==1)  Linear interpolation
    Y1 = Y1*size(X,1)+0.5;        % range 0.5-size(X,1)+0.5
    Y2 = Y2*size(X,2)+0.5;        % range 0.5-size(X,2)+0.5
    Mx = size(X,1);
    Nx = size(X,2);
    Y1 = Y1(:);       Y2 = Y2(:);
    r1 = mod(Y1,1);   r2 = mod(Y2,1);
    f1 = floor(Y1);   f1(f1<1) = Mx;
    f2 = floor(Y2)-1; f2(f2<0) = Nx-1;
    c1 = ceil(Y1);    c1(c1>Mx) = 1;
    c2 = ceil(Y2)-1;  c2(c2>(Nx-1)) = 0;
    X = X(:);
    Y = r1.*r2.*X(c1+c2*Mx)+(1-r1).*r2.*X(f1+c2*Mx)+...
        r1.*(1-r2).*X(c1+f2*Mx)+(1-r1).*(1-r2).*X(f1+f2*Mx);
    Y = reshape(Y,M,N);
else       % type==0   nearest neighbor interpolation
    Y1 = floor(Y1*size(X,1)+1);
    Y2 = floor(Y2*size(X,2)+1);
    Y = zeros(M,N);
    for m=1:M
        for n=1:N
            Y(m,n) = X(Y1(m,n),Y2(m,n));
        end
    end
end

return
