Rotate an image with bicubic interpolation without imrotate

Question

I have implemented a code for image warping using bilinear interpolation:
Matlab image rotation

I would like to improve the code by using bicubic interpolation to rotate the image WITHOUT using the built-in functions like imrotate or imwarp and interp functions in MATLAB.

score 1 · Answer 1 · edited May 23 '17 at 12:06

I successfully managed to implement a full working example.

Code is based on Anna1994's code: Matlab image rotation
Biqubic code is also based on Java (and C++) implementation posted here: http://www.paulinternet.nl/?page=bicubic

The following code applies image rotation example using biqubic interpolation:

function BicubicInterpolationTest()

close all;
% clear all;

img = 'cameraman.tif';

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

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

th=120*pi/180; %Rotate 120 degrees

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,size(input_image,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);


        %Bilinear interpolation (applies RGB image):
%         x1 = floor(a);
%         y1 = floor(b);
%         x2 = x1 + 1;
%         y2 = y1 + 1;        
%         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

        x1 = floor(a);
        y1 = floor(b);

        %Bicubic interpolation (applies grayscale image)
        if ((x1 >= 2) && (y1 >= 2) && (x1 <= W-2) && (y1 <= H-2))
            %Load 4x4 pixels
            P = input_image(y1-1:y1+2, x1-1:x1+2);

            %Interpolation wieghts
            dx = a - x1;
            dy = b - y1;

            %Bi-bicubic interpolation
            output_image(j, i) = bicubicInterpolate(P, dx, 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', 'cubic');
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)]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%http://www.paulinternet.nl/?page=bicubic
%double cubicInterpolate (double p[4], double x) {
%   return p[1] + 0.5 * x*(p[2] - p[0] + x*(2.0*p[0] - 5.0*p[1] + 4.0*p[2] - p[3] + x*(3.0*(p[1] - p[2]) + p[3] - p[0])));
%}
function q = cubicInterpolate(p, x)
q = p(2) + 0.5 * x*(p(3) - p(1) + x*(2.0*p(1) - 5.0*p(2) + 4.0*p(3) - p(4) + x*(3.0*(p(2) - p(3)) + p(4) - p(1))));


%http://www.paulinternet.nl/?page=bicubic
% double bicubicInterpolate (double p[4][4], double x, double y) {
%   double arr[4];
%   arr[0] = cubicInterpolate(p[0], y);
%   arr[1] = cubicInterpolate(p[1], y);
%   arr[2] = cubicInterpolate(p[2], y);
%   arr[3] = cubicInterpolate(p[3], y);
%   return cubicInterpolate(arr, x);
% }
function q = bicubicInterpolate(p, x, y)
q1 = cubicInterpolate(p(1,:), x);
q2 = cubicInterpolate(p(2,:), x);
q3 = cubicInterpolate(p(3,:), x);
q4 = cubicInterpolate(p(4,:), x);
q = cubicInterpolate([q1, q2, q3, q4], y);

I verified implementation by comparing to Matalb build in function imwarp

Result:

The following example uses the "CachedBicubicInterpolator" code version, and also supports RGB image:

function BicubicInterpolationTest2()

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=120*pi/180; %Rotate 120 degrees

s0 = 0.8;
s1 = 0.8;

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,size(input_image,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);

        x1 = floor(a);
        y1 = floor(b);

        %Bicubic interpolation (applies grayscale image)
        if ((x1 >= 2) && (y1 >= 2) && (x1 <= W-2) && (y1 <= H-2))
            %Load 4x4 pixels
            P = input_image(y1-1:y1+2, x1-1:x1+2, :);

            %Interpolation wieghts
            dx = a - x1;
            dy = b - y1;

            %Bi-bicubic interpolation
            output_image(j, i, :) = bicubicInterpolate(P, dx, 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', 'cubic');
figure;imshow(ref_image)

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



function [p0, p1, p2, p3] = list4(P)
P = squeeze(P);
p0 = P(1, :);
p1 = P(2, :);
p2 = P(3, :);
p3 = P(4, :);



%http://www.paulinternet.nl/?page=bicubic
% public void updateCoefficients (double[][] p) {
%     a00 = p[1][1];
%     a01 = -.5*p[1][0] + .5*p[1][2];
%     a02 = p[1][0] - 2.5*p[1][1] + 2*p[1][2] - .5*p[1][3];
%     a03 = -.5*p[1][0] + 1.5*p[1][1] - 1.5*p[1][2] + .5*p[1][3];
%     a10 = -.5*p[0][1] + .5*p[2][1];
%     a11 = .25*p[0][0] - .25*p[0][2] - .25*p[2][0] + .25*p[2][2];
%     a12 = -.5*p[0][0] + 1.25*p[0][1] - p[0][2] + .25*p[0][3] + .5*p[2][0] - 1.25*p[2][1] + p[2][2] - .25*p[2][3];
%     a13 = .25*p[0][0] - .75*p[0][1] + .75*p[0][2] - .25*p[0][3] - .25*p[2][0] + .75*p[2][1] - .75*p[2][2] + .25*p[2][3];
%     a20 = p[0][1] - 2.5*p[1][1] + 2*p[2][1] - .5*p[3][1];
%     a21 = -.5*p[0][0] + .5*p[0][2] + 1.25*p[1][0] - 1.25*p[1][2] - p[2][0] + p[2][2] + .25*p[3][0] - .25*p[3][2];
%     a22 = p[0][0] - 2.5*p[0][1] + 2*p[0][2] - .5*p[0][3] - 2.5*p[1][0] + 6.25*p[1][1] - 5*p[1][2] + 1.25*p[1][3] + 2*p[2][0] - 5*p[2][1] + 4*p[2][2] - p[2][3] - .5*p[3][0] + 1.25*p[3][1] - p[3][2] + .25*p[3][3];
%     a23 = -.5*p[0][0] + 1.5*p[0][1] - 1.5*p[0][2] + .5*p[0][3] + 1.25*p[1][0] - 3.75*p[1][1] + 3.75*p[1][2] - 1.25*p[1][3] - p[2][0] + 3*p[2][1] - 3*p[2][2] + p[2][3] + .25*p[3][0] - .75*p[3][1] + .75*p[3][2] - .25*p[3][3];
%     a30 = -.5*p[0][1] + 1.5*p[1][1] - 1.5*p[2][1] + .5*p[3][1];
%     a31 = .25*p[0][0] - .25*p[0][2] - .75*p[1][0] + .75*p[1][2] + .75*p[2][0] - .75*p[2][2] - .25*p[3][0] + .25*p[3][2];
%     a32 = -.5*p[0][0] + 1.25*p[0][1] - p[0][2] + .25*p[0][3] + 1.5*p[1][0] - 3.75*p[1][1] + 3*p[1][2] - .75*p[1][3] - 1.5*p[2][0] + 3.75*p[2][1] - 3*p[2][2] + .75*p[2][3] + .5*p[3][0] - 1.25*p[3][1] + p[3][2] - .25*p[3][3];
%     a33 = .25*p[0][0] - .75*p[0][1] + .75*p[0][2] - .25*p[0][3] - .75*p[1][0] + 2.25*p[1][1] - 2.25*p[1][2] + .75*p[1][3] + .75*p[2][0] - 2.25*p[2][1] + 2.25*p[2][2] - .75*p[2][3] - .25*p[3][0] + .75*p[3][1] - .75*p[3][2] + .25*p[3][3];
% }

% public double getValue (double x, double y) {
%     double x2 = x * x;
%     double x3 = x2 * x;
%     double y2 = y * y;
%     double y3 = y2 * y;
% 
%     return (a00 + a01 * y + a02 * y2 + a03 * y3) +
%            (a10 + a11 * y + a12 * y2 + a13 * y3) * x +
%            (a20 + a21 * y + a22 * y2 + a23 * y3) * x2 +
%            (a30 + a31 * y + a32 * y2 + a33 * y3) * x3;
% }
function q = bicubicInterpolate(P, x, y)
[p00, p01, p02, p03] = list4(P(1, :, :));
[p10, p11, p12, p13] = list4(P(2, :, :));
[p20, p21, p22, p23] = list4(P(3, :, :));
[p30, p31, p32, p33] = list4(P(4, :, :));

a00 = p11;
a01 = -.5*p10 + .5*p12;
a02 = p10 - 2.5*p11 + 2*p12 - .5*p13;
a03 = -.5*p10 + 1.5*p11 - 1.5*p12 + .5*p13;
a10 = -.5*p01 + .5*p21;
a11 = .25*p00 - .25*p02 - .25*p20 + .25*p22;
a12 = -.5*p00 + 1.25*p01 - p02 + .25*p03 + .5*p20 - 1.25*p21 + p22 - .25*p23;
a13 = .25*p00 - .75*p01 + .75*p02 - .25*p03 - .25*p20 + .75*p21 - .75*p22 + .25*p23;
a20 = p01 - 2.5*p11 + 2*p21 - .5*p31;
a21 = -.5*p00 + .5*p02 + 1.25*p10 - 1.25*p12 - p20 + p22 + .25*p30 - .25*p32;
a22 = p00 - 2.5*p01 + 2*p02 - .5*p03 - 2.5*p10 + 6.25*p11 - 5*p12 + 1.25*p13 + 2*p20 - 5*p21 + 4*p22 - p23 - .5*p30 + 1.25*p31 - p32 + .25*p33;
a23 = -.5*p00 + 1.5*p01 - 1.5*p02 + .5*p03 + 1.25*p10 - 3.75*p11 + 3.75*p12 - 1.25*p13 - p20 + 3*p21 - 3*p22 + p23 + .25*p30 - .75*p31 + .75*p32 - .25*p33;
a30 = -.5*p01 + 1.5*p11 - 1.5*p21 + .5*p31;
a31 = .25*p00 - .25*p02 - .75*p10 + .75*p12 + .75*p20 - .75*p22 - .25*p30 + .25*p32;
a32 = -.5*p00 + 1.25*p01 - p02 + .25*p03 + 1.5*p10 - 3.75*p11 + 3*p12 - .75*p13 - 1.5*p20 + 3.75*p21 - 3*p22 + .75*p23 + .5*p30 - 1.25*p31 + p32 - .25*p33;
a33 = .25*p00 - .75*p01 + .75*p02 - .25*p03 - .75*p10 + 2.25*p11 - 2.25*p12 + .75*p13 + .75*p20 - 2.25*p21 + 2.25*p22 - .75*p23 - .25*p30 + .75*p31 - .75*p32 + .25*p33;


x2 = x * x;
x3 = x2 * x;
y2 = y * y;
y3 = y2 * y;

% q = (a00 + a01 * y + a02 * y2 + a03 * y3) +...
%     (a10 + a11 * y + a12 * y2 + a13 * y3) * x +...
%     (a20 + a21 * y + a22 * y2 + a23 * y3) * x2 +...
%     (a30 + a31 * y + a32 * y2 + a33 * y3) * x3;

q = (a00 + a01 * x + a02 * x2 + a03 * x3) +...
    (a10 + a11 * x + a12 * x2 + a13 * x3) * y +...
    (a20 + a21 * x + a22 * x2 + a23 * x3) * y2 +...
    (a30 + a31 * x + a32 * x2 + a33 * x3) * y3;

Result:

Typically you can flag the deleted question and state that the user deleted after you answered. A mod will open the question back up for you. — Suever, Oct 30 '16 at 22:49
Hi! The reason the question was deleted (if it was from Anna) its most likely because it was a homework question, and you solved it so they quicly try to delete any proof of getting help in the internet. I know that because Anna is very very likely to be one of my students, as this exact exercise, with almost the same code, is given by us to the students to do this task. My point: Try not to answer obvious homework questions with full working code. They will delete it, and it doesn't really help them learn. Come to the MATLAB chatroom to talk more about it, if you wish. — Ander Biguri, Nov 01 '16 at 09:53
@AnderBiguri No it was not from Anna. I guessed it was a homework question, but I think the subject is relevant for many other users. After posting bilinear solution, I really wanted to post bicubic solution, and tempted to solve the student's homework. I deleted the remark. — Rotem, Nov 01 '16 at 11:31

Rotate an image with bicubic interpolation without imrotate

1 Answers1

Linked