I am trying to implement Bicubic interpolation for image scaling in matlab. The issue is that it properly works for grayscale images, however, for colored images, the result is again in grayscale. Please help me find out what the problem is. Thank you.
For bicubic interpoaltion i have used the matrix containing the gradients. The matix can be found at https://en.wikipedia.org/wiki/Bicubic_interpolation.
Here is my code.
input_image = im2double(imread('peppers.png'));
x_res = 700;
y_res = 700;
imshow(input_image, []);
M_inv = [
1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0;
-3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0;
2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0;
0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0;
0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0;
-3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0;
0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0;
9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1;
-6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1;
2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0;
0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0;
-6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1;
4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1
];
I = input_image;
[j k c] = size(I);
%{
if c > 1
I = double(rgb2gray(I));
end
%}
x_new = x_res;
y_new = y_res;
x_scale = x_new./(j-1);
y_scale = y_new./(k-1);
temp_image = zeros(x_new,y_new);
Ix = double(zeros(j,k));
for count1 = 1:j
for count2 = 1:k
if( (count2==1) || (count2==k) )
Ix(count1,count2)=0;
else
Ix(count1,count2)=(0.5).*(I(count1,count2+1)-I(count1,count2-1));
end
end
end
Iy = double(zeros(j,k));
for count1 = 1:j
for count2 = 1:k
if( (count1==1) || (count1==j) )
Iy(count1,count2)=0;
else
Iy(count1,count2)=(0.5).*(I(count1+1,count2)-I(count1-1,count2));
end
end
end
Ixy = double(zeros(j,k));
for count1 = 1:j
for count2 = 1:k
if( (count1==1) || (count1==j) || (count2==1) || (count2==k) )
Ixy(count1,count2)=0;
else
Ixy(count1,count2)=(0.25).*((I(count1+1,count2+1)+I(count1-1,count2-1)) - (I(count1+1,count2-1)+I(count1-1,count2+1)));
end
end
end
for count1 = 0:x_new-1
for count2 = 0:y_new-1
W = -(((count1./x_scale)-floor(count1./x_scale))-1);
H = -(((count2./y_scale)-floor(count2./y_scale))-1);
I11_index = [1+floor(count1./x_scale),1+floor(count2./y_scale)];
I21_index = [1+floor(count1./x_scale),1+ceil(count2./y_scale)];
I12_index = [1+ceil(count1./x_scale),1+floor(count2./y_scale)];
I22_index = [1+ceil(count1./x_scale),1+ceil(count2./y_scale)];
I11 = I(I11_index(1),I11_index(2));
I21 = I(I21_index(1),I21_index(2));
I12 = I(I12_index(1),I12_index(2));
I22 = I(I22_index(1),I22_index(2));
Ix11 = Ix(I11_index(1),I11_index(2));
Ix21 = Ix(I21_index(1),I21_index(2));
Ix12 = Ix(I12_index(1),I12_index(2));
Ix22 = Ix(I22_index(1),I22_index(2));
Iy11 = Iy(I11_index(1),I11_index(2));
Iy21 = Iy(I21_index(1),I21_index(2));
Iy12 = Iy(I12_index(1),I12_index(2));
Iy22 = Iy(I22_index(1),I22_index(2));
Ixy11 = Ixy(I11_index(1),I11_index(2));
Ixy21 = Ixy(I21_index(1),I21_index(2));
Ixy12 = Ixy(I12_index(1),I12_index(2));
Ixy22 = Ixy(I22_index(1),I22_index(2));
beta = [I11 I21 I12 I22 Ix11 Ix21 Ix12 Ix22 Iy11 Iy21 Iy12 Iy22 Ixy11 Ixy21 Ixy12 Ixy22];
alpha = M_inv*beta';
temp_p=0;
for count3 = 1:16
w_temp = floor((count3-1)/4);
h_temp = mod(count3-1,4);
temp_p = temp_p + alpha(count3).*((1-W)^(w_temp)).*((1-H)^(h_temp));
end
temp_image(count1+1,count2+1)=temp_p;
end
end
output_image = temp_image;
figure;
imshow(output_image);
You need to adjust your code so that it performs the interpolation on each channel individually. Simply put, you need to make I each channel, create the output_image so that it has multiple channels and loop through each channel individually.
With these changes, we thus have:
input_image = im2double(imread('peppers.png'));
x_res = 700;
y_res = 700;
imshow(input_image, []);
% input_image - an image on which to perform bicubic interpolation
% x_res - the new horizontal dimensions (in pixels)
% y_res - the new vertical dimensions (in pixels)
%Define the inverted weighting matrix, M^(-1), no need to recalculate it
%ever again
M_inv = [1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0;
-3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0;
2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0;
0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0;
0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0;
-3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0;
0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0;
9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1;
-6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1;
2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0;
0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0;
-6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1;
4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1
];
%Make a copy of the input image
%I = input_image; % Ray: Not needed here
%Determine the dimensions of the source image
%Note that we will have three values - width, height, and the number
%of color vectors, 3
[j k c] = size(input_image); % Ray: Change
%Specify the new image dimensions we want for our larger output image
x_new = x_res;
y_new = y_res;
%Determine the ratio of the old dimensions compared to the new dimensions
%Referred to as S1 and S2 in my tutorial
x_scale = x_new./(j-1);
y_scale = y_new./(k-1);
% Change by Ray - Declare new output image here with c channels
output_image = zeros(x_new, y_new, c);
% Change by Ray - Now loop through each channel
for z = 1 : c
%Declare and initialize an output image buffer
temp_image = zeros(x_new,y_new);
% New - Change by Ray. Access the right channel
I = input_image(:,:,z);
Ix = double(zeros(j,k));
for count1 = 1:j
for count2 = 1:k
if( (count2==1) || (count2==k) )
Ix(count1,count2)=0;
else
Ix(count1,count2)=(0.5).*(I(count1,count2+1)-I(count1,count2-1));
end
end
end
Iy = double(zeros(j,k));
for count1 = 1:j
for count2 = 1:k
if( (count1==1) || (count1==j) )
Iy(count1,count2)=0;
else
Iy(count1,count2)=(0.5).*(I(count1+1,count2)-I(count1-1,count2));
end
end
end
Ixy = double(zeros(j,k));
for count1 = 1:j
for count2 = 1:k
if( (count1==1) || (count1==j) || (count2==1) || (count2==k) )
Ixy(count1,count2)=0;
else
Ixy(count1,count2)=(0.25).*((I(count1+1,count2+1)+I(count1-1,count2-1)) - (I(count1+1,count2-1)+I(count1-1,count2+1)));
end
end
end
for count1 = 0:x_new-1
for count2 = 0:y_new-1
%Calculate the normalized distance constants, h and w
W = -(((count1./x_scale)-floor(count1./x_scale))-1);
H = -(((count2./y_scale)-floor(count2./y_scale))-1);
%Determine the indexes/address of the 4 neighbouring pixels from
%the source data/image
I11_index = [1+floor(count1./x_scale),1+floor(count2./y_scale)];
I21_index = [1+floor(count1./x_scale),1+ceil(count2./y_scale)];
I12_index = [1+ceil(count1./x_scale),1+floor(count2./y_scale)];
I22_index = [1+ceil(count1./x_scale),1+ceil(count2./y_scale)];
%Calculate the four nearest function values
I11 = I(I11_index(1),I11_index(2));
I21 = I(I21_index(1),I21_index(2));
I12 = I(I12_index(1),I12_index(2));
I22 = I(I22_index(1),I22_index(2));
%Calculate the four nearest horizontal derivatives
Ix11 = Ix(I11_index(1),I11_index(2));
Ix21 = Ix(I21_index(1),I21_index(2));
Ix12 = Ix(I12_index(1),I12_index(2));
Ix22 = Ix(I22_index(1),I22_index(2));
%Calculate the four nearest vertical derivatives
Iy11 = Iy(I11_index(1),I11_index(2));
Iy21 = Iy(I21_index(1),I21_index(2));
Iy12 = Iy(I12_index(1),I12_index(2));
Iy22 = Iy(I22_index(1),I22_index(2));
%Calculate the four nearest cross derivatives
Ixy11 = Ixy(I11_index(1),I11_index(2));
Ixy21 = Ixy(I21_index(1),I21_index(2));
Ixy12 = Ixy(I12_index(1),I12_index(2));
Ixy22 = Ixy(I22_index(1),I22_index(2));
%Create our beta-vector
beta = [I11 I21 I12 I22 Ix11 Ix21 Ix12 Ix22 Iy11 Iy21 Iy12 Iy22 Ixy11 Ixy21 Ixy12 Ixy22];
alpha = M_inv*beta';
temp_p=0;
for count3 = 1:16
w_temp = floor((count3-1)/4);
h_temp = mod(count3-1,4);
temp_p = temp_p + alpha(count3).*((1-W)^(w_temp)).*((1-H)^(h_temp));
end
temp_image(count1+1,count2+1)=temp_p;
end
end
% New - Change by Ray - Assign to output channel
output_image(:,:,z) = temp_image;
end
figure;
imshow(output_image);