function Irot =rotation(I,theta); 

[m,n] = size(I); 
xc = (n-1)/2; 
yc = (m-1)/2; 
D=sqrt(m^2+n^2)/2; 
MM = ceil(D-m/2); 
NN = ceil(D-n/2); 
mm = m + 2*MM; nn = n+2*NN;   
Irot = zeros(mm,nn);  
for j=1:mm 
    for k=1:nn 
        x = - NN + k -1; 
        y = - MM + j-1;
        X = xc + cos(theta)*(x-xc) - sin(theta)*(y-yc); 
        Y = yc  + sin(theta)*(x-xc) + cos(theta)*(y-yc);
        if max(X*(X-(n-1)),Y*(Y-(m-1)))>0
                Irot(j,k) = 0; 
        else 
            t = X - floor(X);  
            s = Y - floor(Y); 
            indexXinf = floor(X)+1; 
            indexXsup = ceil(X)+1; 
            indexYinf = floor(Y)+1; 
            indexYsup = ceil(Y) +1; 
            Iaux1 = (1-t)*I(indexYinf,indexXinf) + t*I(indexYinf,indexXsup); 
            Iaux2 = (1-t)*I(indexYsup,indexXinf) + t*I(indexYsup,indexXsup); 
            Irot(j,k) = (1-s)*Iaux1 + s*Iaux2;
        end        
    end
end     
end

