Matlab image rotation

Question

I am new to image processing, I have implemented a code for image warping and it works perfectly. I would like to improve the code by using linear interpolation to rotate the image WITHOUT using the built-in function (interp). Here is my code:

close all;
clear all;

img = 'woods.jpg';

input_image =double(imread(img))./255;

H=size(input_image,1); 
W=size(input_image,2);

th=pi/4;

s0 = 2;
s1 = 2;

x0 = -W/2;
x1 = -H/2;

T=[1 0 x0 ; ...
   0 1 x1 ; ...
   0 0 1];

RST = [ (s0*cos(th))   (-s1*sin(th)) ((s0*x0*cos(th))-(s1*x1*sin(th))); ...
        (s0*sin(th))   (s1*cos(th))   ((s0*x0*sin(th))+(s1*x1*cos(th))); ...
        0   0   1];

M=inv(T)*R;
N = inv(M);


output_image=zeros(H,W,3);

for i=1:W
    for j=1:H

        x = [i ; j ; 1];
        y = N * x;

        a = y(1)/y(3);
        b = y(2)/y(3);

        a = round(a);
        b = round(b);

        if (a>0 && a<=W && b>0 && b<=H) 
           output_image(j,i,:)=input_image(b,a,:);
        end

    end
end

imgshow(output_image);

Answer 1

Check the following solution:

I verified implementation by comparing to Matalb build in function [imwarp][1].

close all;
clear all;

img = 'peppers.png';

input_image =double(imread(img))./255;

H=size(input_image,1);  % height
W=size(input_image,2);  % width

th=pi/4;

s0 = 2;
s1 = 2;

x0 = -W/2;
x1 = -H/2;

T=[1 0 x0 ; ...
   0 1 x1 ; ...
   0 0 1];

RST = [ (s0*cos(th))   (-s1*sin(th)) ((s0*x0*cos(th))-(s1*x1*sin(th))); ...
        (s0*sin(th))   (s1*cos(th))   ((s0*x0*sin(th))+(s1*x1*cos(th))); ...
        0   0   1];

M=inv(T)*RST;
N = inv(M);

output_image=zeros(H,W,3);

for i=1:W
    for j=1:H

        x = [i ; j ; 1];
        y = N * x;

        a = y(1)/y(3);
        b = y(2)/y(3);

        %Nearest neighbor
        %a = round(a);
        %b = round(b);


        x1 = floor(a);
        y1 = floor(b);
        x2 = x1 + 1;
        y2 = y1 + 1;

        %Bi-linear interpolation ilsutration:
        %Image coordinates style (horizontal index first)
        %
        %(x1,y1)    |           (x2,y1)
        %           | 1-dy
        %    1-dx   |       dx
        %   ------(a,b)------------
        %           |
        %           |
        %           |
        %           | dy
        %           |
        %           |
        %(x1,y2)    |           (x2,y2)

        if ((x1 >= 1) && (y1 >= 1) && (x2 <= W) && (y2 <= H))
            %Load 2x2 pixels
            i11 = input_image(y1, x1, :); %Top left pixel
            i21 = input_image(y2, x1, :); %Bottom left pixel
            i12 = input_image(y1, x2, :); %Top right pixel
            i22 = input_image(y2, x2, :); %Bottom right pixel

            %Interpolation wieghts
            dx = x2 - a;
            dy = y2 - b;

            %Bi-lienar interpolation
            output_image(j, i, :) = i11*dx*dy + i21*dx*(1-dy) + i12*(1-dx)*dy + i22*(1-dx)*(1-dy);
        end
    end
end

imshow(output_image);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Verify implementation by comparing with Matalb build in function imwarp:
tform = affine2d(M');
ref_image = imwarp(input_image, tform, 'OutputView', imref2d(size(input_image)), 'Interp', 'linear');
figure;imshow(ref_image)

figure;imshow(output_image - ref_image)
max_diff = max(abs(output_image(:) - ref_image(:)));
disp(['Maximum difference from imwarp = ', num2str(max_diff)]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Result:

---

Remark:
I missed the following execelnt post that resizes an image using bilinear interpolation (bilinear interpolation is explained better):
Resize an image with bilinear interpolation without imresize