Image rotation by Matlab without using imrotate

Question

I am trying to rotate an image with Matlab without using imrotate function. I actually made it by using transformation matrix.But it is not good enough.The problem is, the rotated image is "sliding".Let me tell you with pictures.

This is my image which I want to rotate:

But when I rotate it ,for example 45 degrees, it becomes this:

I am asking why this is happening.Here is my code,is there any mathematical or programming mistakes about it?

image=torso;

%image padding
[Rows, Cols] = size(image); 
Diagonal = sqrt(Rows^2 + Cols^2); 
RowPad = ceil(Diagonal - Rows) + 2;
ColPad = ceil(Diagonal - Cols) + 2;
imagepad = zeros(Rows+RowPad, Cols+ColPad);
imagepad(ceil(RowPad/2):(ceil(RowPad/2)+Rows-1),ceil(ColPad/2):(ceil(ColPad/2)+Cols-1)) = image;

degree=45;

%midpoints
midx=ceil((size(imagepad,1)+1)/2);
midy=ceil((size(imagepad,2)+1)/2);

imagerot=zeros(size(imagepad));

%rotation
for i=1:size(imagepad,1)
    for j=1:size(imagepad,2)

         x=(i-midx)*cos(degree)-(j-midy)*sin(degree);
         y=(i-midx)*sin(degree)+(j-midy)*cos(degree);
         x=round(x)+midx;
         y=round(y)+midy;

         if (x>=1 && y>=1)
              imagerot(x,y)=imagepad(i,j); % k degrees rotated image         
         end

    end
end

 figure,imagesc(imagerot);
 colormap(gray(256));

Answer 1

The reason you have holes in your image is because you are computing the location in imagerot of each pixel in imagepad. You need to do the computation the other way around. That is, for each pixel in imagerot interpolate in imagepad. To do this, you just need to apply the inverse transform, which in the case of a rotation matrix is just the transpose of the matrix (just change the sign on each sin and translate the other way).

Loop over pixels in imagerot:

imagerot=zeros(size(imagepad)); % midx and midy same for both

for i=1:size(imagerot,1)
    for j=1:size(imagerot,2)

         x= (i-midx)*cos(rads)+(j-midy)*sin(rads);
         y=-(i-midx)*sin(rads)+(j-midy)*cos(rads);
         x=round(x)+midx;
         y=round(y)+midy;

         if (x>=1 && y>=1 && x<=size(imagepad,2) && y<=size(imagepad,1))
              imagerot(i,j)=imagepad(x,y); % k degrees rotated image         
         end

    end
end

Also note that your midx and midy need to be calculated with size(imagepad,2) and size(imagepad,1) respectively, since the first dimension refers to the number of rows (height) and the second to width.

NOTE: The same approach applies when you decide to adopt an interpolation scheme other than nearest neighbor, as in Rody's example with linear interpolation.

EDIT: I'm assuming you are using a loop for demonstrative purposes, but in practice there is no need for loops. Here's an example of nearest neighbor interpolation (what you are using), keeping the same size image, but you can modify this to produce a larger image that includes the whole source image:

imagepad = imread('peppers.png');
[nrows ncols nslices] = size(imagepad);
midx=ceil((ncols+1)/2);
midy=ceil((nrows+1)/2);

Mr = [cos(pi/4) sin(pi/4); -sin(pi/4) cos(pi/4)]; % e.g. 45 degree rotation

% rotate about center
[X Y] = meshgrid(1:ncols,1:nrows);
XYt = [X(:)-midx Y(:)-midy]*Mr;
XYt = bsxfun(@plus,XYt,[midx midy]);

xout = round(XYt(:,1)); yout = round(XYt(:,2)); % nearest neighbor!
outbound = yout<1 | yout>nrows | xout<1 | xout>ncols;
zout=repmat(cat(3,1,2,3),nrows,ncols,1); zout=zout(:);
xout(xout<1) = 1; xout(xout>ncols) = ncols;
yout(yout<1) = 1; yout(yout>nrows) = nrows;
xout = repmat(xout,[3 1]); yout = repmat(yout,[3 1]);
imagerot = imagepad(sub2ind(size(imagepad),yout,xout,zout(:))); % lookup
imagerot = reshape(imagerot,size(imagepad));
imagerot(repmat(outbound,[1 1 3])) = 0; % set background value to [0 0 0] (black)

To modify the above to linear interpolation, compute the 4 neighboring pixels to each coordinate in XYt and perform a weighted sum using the fractional components product as the weights. I'll leave that as an exercise, since it would only serve to bloat my answer further beyond the scope of your question. :)

Answer 2

The method you are using (rotate by sampling) is the fastest and simplest, but also the least accurate.

Rotation by area mapping, as given below (this is a good reference), is much better at preserving color.

But: note that this will only work on greyscale/RGB images, but NOT on colormapped images like the one you seem to be using.

image = imread('peppers.png');

figure(1), clf, hold on
subplot(1,2,1)
imshow(image);

degree = 45;

switch mod(degree, 360)
    % Special cases
    case 0
        imagerot = image;
    case 90
        imagerot = rot90(image);
    case 180
        imagerot = image(end:-1:1, end:-1:1);
    case 270
        imagerot = rot90(image(end:-1:1, end:-1:1));

    % General rotations
    otherwise

        % Convert to radians and create transformation matrix
        a = degree*pi/180;
        R = [+cos(a) +sin(a); -sin(a) +cos(a)];

        % Figure out the size of the transformed image
        [m,n,p] = size(image);
        dest = round( [1 1; 1 n; m 1; m n]*R );
        dest = bsxfun(@minus, dest, min(dest)) + 1;
        imagerot = zeros([max(dest) p],class(image));

        % Map all pixels of the transformed image to the original image
        for ii = 1:size(imagerot,1)
            for jj = 1:size(imagerot,2)
                source = ([ii jj]-dest(1,:))*R.';
                if all(source >= 1) && all(source <= [m n])

                    % Get all 4 surrounding pixels
                    C = ceil(source);
                    F = floor(source);

                    % Compute the relative areas
                    A = [...
                        ((C(2)-source(2))*(C(1)-source(1))),...
                        ((source(2)-F(2))*(source(1)-F(1)));
                        ((C(2)-source(2))*(source(1)-F(1))),...
                        ((source(2)-F(2))*(C(1)-source(1)))];

                    % Extract colors and re-scale them relative to area
                    cols = bsxfun(@times, A, double(image(F(1):C(1),F(2):C(2),:)));

                    % Assign                     
                    imagerot(ii,jj,:) = sum(sum(cols),2);

                end
            end
        end        
end

subplot(1,2,2)
imshow(imagerot);

Output: