genintpol.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function Y=GenIntPol(X,x0,y0,dx,dy,alpha,type,M,N);
% GenIntPol  General Interpolation of Images/Patterns, generate textures 
%            Nearest Neighbor or linear interpolation
% Interpolation steps are equal for each direction.
% Interpolation grid may be rotated, and/or translated by an offset.
% If initial image/pattern is too small, it will be periodically extended.
% For images (and matrices) the first axis (x-axis) points downwards and the 
% second axis (y-axis) points to the right. 
% Used as in mainT01 this function generates different textures.
%
% Y=GenIntPol(X,x0,y0,dx,dy,alpha,type,M,N);
%-----------------------------------------------------------------------------------
% arguments:
%   X     - input image/pattern, this is scaled to fill the box given by
%           its four corners; (0,0), (0,1), (1,1) and (1,0).
%   x0    - offset in x-axis for the first point of the grid in Y
%   y0    - offset in y-axis for the first point of the grid in Y
%   dx    - increment in x-axis for the grid in Y
%   dy    - increment in y-axis for the grid in Y
%   alpha - angle (in radians) the grid of Y 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)
%-----------------------------------------------------------------------------------

%-----------------------------------------------------------------------------------
% Copyright (c) 2002.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  [e-mail]   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  10.12.2002  KS: function made
%-----------------------------------------------------------------------------------

Mfile='GenIntPol';
if nargin~=9
    error([Mfile,': wrong number of arguments, see help.']);
end

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)
    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);
    if 1   % try to speed up (it seems to work ok)
        Y1=Y1(:);Y2=Y2(:);
        r1=mod(Y1,1);  r2=mod(Y2,1);
        f1=floor(Y1); f1(find(f1<1))=Mx;
        f2=floor(Y2)-1; f2(find(f2<0))=Nx-1;
        c1=ceil(Y1); c1(find(c1>Mx))=1;
        c2=ceil(Y2)-1; c2(find(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   % this is the slow way of doing this
        for m=1:M
            for n=1:N
                y1=Y1(m,n);     y2=Y2(m,n);
                r1=mod(y1,1);   r2=mod(y2,1);
                f1=floor(y1);if f1<1; f1=Mx; end;
                f2=floor(y2);if f2<1; f2=Nx; end;
                c1=ceil(y1);if c1>Mx; c1=1; end;
                c2=ceil(y2);if c2>Nx; c2=1; end;
                Y(m,n)=r1*r2*X(c1,c2)+r1*(1-r2)*X(c1,f2)+...
                    (1-r1)*r2*X(f1,c2)+(1-r1)*(1-r2)*X(f1,f2);
            end
        end
    end
else       % type==0
    if 0   % try to speed up
        Y1=floor(Y1(:)*size(X,1)+1);
        Y2=floor(Y2(:)*size(X,2));
        X=X(:);
        Y=X(Y2*M+Y1);
        Y=reshape(Y,M,N);
    else
        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
end


return;


% test of function;
clear all
if 1
   K=5;X=[ones(K,K),zeros(K,K);zeros(K,K),ones(K,K)];
   %X=[0,1,0;1,1,1;0,1,0];
   x0=rand(1);y0=rand(1);
   x0=eps;y0=eps;
   dx=0.025;dy=0.025;
   alpha=0*pi/180;
   M=64;N=64;
   type=1;
else
   X=double(GetIm(1));    % 512x512 image Lena
   x0=0.001;y0=0.001;
   dx=0.002;dy=0.002;
   alpha=15*pi/180;
   M=512;N=512;
   type=0;
end
tic;
Y=GenIntPol(X,x0,y0,dx,dy,alpha,type,M,N);
toc;

figure(1);clf;
colormap(gray);
imagesc(Y);


% noen syntetiske teksturer ??
%X=[8,7,6,5; 7,8,7,6; 6,7,8,7; 5,6,7,8]/10;
%X=[X,1-fliplr(X);1-fliplr(X),X];
X1=[1,1,0,1,1,0; 1,1,0,1,1,0; 2,2,0,2,2,2; 1,1,0,1,1,0; 1,1,0,1,1,0; 2,2,2,2,2,0]/2;
X2=[1,0,1,0;2,0,2,2;1,0,1,0;2,2,2,0]/2;
x0=0.01;y0=0.01;
dx=1/(4*pi);dy=1/(4*pi);
alpha=20*pi/180;
M=256;N=256;
type=1;
Y1=GenIntPol(X1,x0,y0,dx,dy,alpha,type,M,N);
Y2=GenIntPol(X2,x0,y0,dx,dy,alpha,type,M,N);
%
Y3=[Y1(:,1:(N/2)),Y2(:,(N/2+1):N)];
figure(1);clf;imagesc(Y3);
colormap(gray);