function [B] = ResizeImageMatrix(varargin)
% ResizeImageMatrix : This program can be used to resize an image matrix or
% resize a 2D matrix with various interpolation methods
% Interpolations Implemented : Nearest Neighbour
% : Bilinear
% : Bicubic 2x2
% : Bicubic 4x4
% : Bicubic 6x6
% : Bicubic 8x8
%
% THIS IMPLEMENTATION IS FOR EDUCATIONAL PERPOUSES ONLY . IT IS NOT
% OPTAMIZED AT ALL AND NOT EFFICIENT IN TERMS OF SPEED AND MEMORY.THE
% ONLY THING YOU WILL LEARN FROM THIS IS HOW INTERPOLATION IS DONE
%
% B = ResizeImageMatrix(A,ResizeRatioY,ResizeRatioX,InterpolationMethod)
% B = ResizeImageMatrix(A,ResizeRatio,InterpolationMethod)
%
% Input : A = Image Matrix
% Interpolation Method : 'near' = nearest neighbour
% : 'bili' = bilinear
% : 'bic2' = bicubic 2x2
% : 'bic4' = bicubic 4x4
% : 'bic6' = bicubic 6x6
% : 'bic8' = bicubic 8x8
% ResizeRatio = resize ratio, if the same factor of resize is desired
% along both axis
% ResizeRatioY = resize ratio along Y - axis
% ResizeRatioX = resize ratio along X - axis
%
% Example 1:
%
% I=imread('cameraman.tif');
% I = double(I)/255;
% B = ResizeImageMatrix(I,1.2,'bic4');
% imtool(B);
%
% Example 2:
%
% A=zeros(11,11);
% A(6,6)=1;
% B = ResizeImageMatrix(A,22,24,'bic6');
% imtool(B);
%
% Example 3
%
% A = uint8(round(255-255*rand(5,5,3)));
% B = ResizeImageMatrix(A,100,100,'bic4');
% imtool(B);
if length(varargin)==4
A = varargin{1};
ResizeRatioY = varargin{2};
ResizeRatioX = varargin{3};
if ResizeRatioY<=0 | ResizeRatioX<=0
str ='Resizeratio should be positive';
error(str,'1');
end
InterpolationMethod = varargin{4};
elseif length(varargin)==3
A = varargin{1};
ResizeRatioY = varargin{2};
ResizeRatioX = varargin{2};
if ResizeRatioY<=0
str ='Resizeratio should be positive';
error(str,'1');
end
InterpolationMethod = varargin{3};
else
str='Number of arguments should be 3 or 4 \n'
str =[str 'If 3 then args should be Image,Ratio,Interp \nIf 4 then args should be Image,RatioY,RatioX,Interp'];
error(str,'1');
end
ClassFlag=0;
if strcmp(class(A),'uint8');
A = double(A);
A = A/255;
ClassFlag=1;
elseif ~strcmp(class(A),'double');
error('only double or uint8 classes supported for A');
end
% m = No of Rows in the image matrix A
% n = No of Columns in the image matrix A
% o = No of Layer for colour Image
[m n o] = size(A);
if o == 3
B(:,:,1) = ResizeImageMatrix(A(:,:,1),ResizeRatioY,ResizeRatioX,InterpolationMethod);
B(:,:,2) = ResizeImageMatrix(A(:,:,2),ResizeRatioY,ResizeRatioX,InterpolationMethod);
B(:,:,3) = ResizeImageMatrix(A(:,:,3),ResizeRatioY,ResizeRatioX,InterpolationMethod);
if ClassFlag
B = uint8(B*255);
end
return;
elseif o~=1
error('only M*N*1 or M*N*3 matrices supported');
B = [];
end
% Calculation of size of new image matrix rounded downwards
NewRows = floor(m*ResizeRatioY);
NewColumns = floor(n*ResizeRatioX);
% Recalculation of resize ratios
% if ResizeRatioY = 1.4 & m = 6 then NewRows = round(8.4) = 8 this creats
% problems in interpolation as mapping of pixels can be from 6 to 8 but and
% not 8.4. So the ResizeRatio is recalculated
% Therefore ResizeRatioY = NewRows/m; = 8/6 = 1.33
ResizeRatioY = NewRows/m;
ResizeRatioX = NewColumns/n;
if ResizeRatioY==0 | ResizeRatioX==0
str ='Resizeratio too small';
error(str,'1');
end
% This is the all important transformation matrix
a = [ResizeRatioY,0,0;0,ResizeRatioX,0;0.5*(1-ResizeRatioY),0.5*(1-ResizeRatioX),1];
% Nearest Neighbour Interpolation
if strcmp(InterpolationMethod,'near')
B(1:NewRows,1:NewColumns)=0;
for i = 1:NewColumns
for j =1:NewRows
% Reverse Mapping of coordinates to find out what pixel in the
% the matrix depended on which pixel in the old image matrix
temp = [j i 1]*inv(a);
% As the name suggests each new pixel depends only on one old
% pixel
B(j,i)=A(round(temp(1)),round(temp(2)));
end
end
elseif strcmp(InterpolationMethod,'bili')
% In case any of ResizeRatioY or ResizeRatioX are < 1 an Anti
% Aliasing filter is applied .
% Anti Aliasing Filter Starts
if ResizeRatioX<1
A = ReSample(A,ResizeRatioX);
end
if ResizeRatioY<1
A = ReSample(A',ResizeRatioY);
A =A';
end
% Anti Aliasing Filter Starts
% Resizeing starts here
for i = 1:NewColumns
for j =1:NewRows
% Reverse Mapping of coordinates to find out what pixel in the
% new image matrix depended on which pixel in the old image
% matrix.
temp = [j i 1]*inv(a);
% Now Comes the fun part.
% Linear Intepolation works as follows: if i have a function
% F(x)=y; such the y=[1,4] for x=[1,2] we assume that F(x) is a
% straing line so if i want to find the value of say F(1.4)
% i assume di = 1 - 0.4 = 0.6 or di = 1 - fractional part of x
% Hence F(floor(t) + fractional(t)) = y(t)*dt + y(t+1)*(1-dt);
% So for t=1.4 F(1.4) = 0.6*1 + 0.4*4 = 2.2;
% This technique when applied to each dimention individually is
% called Bilinear Interpolation.
% di dj are fractional values in y & x direction respectively
if floor(temp(1))==temp(1)
di=1;
else
di=temp(1)-floor(temp(1));
di = 1 - di;
end
if floor(temp(2))==temp(2)
dj=1;
else
dj=temp(2)-floor(temp(2));
dj = 1 - dj;
end
% y & x are reversed mapped coordinate of pixel i,j in final image
y = temp(1);
x = temp(2);
% B(i,j) = function of (A(Fy,Fx),A(Cy,Fx),A(Fy,Cx),A(Cy,Cx),di,dj)
Fx=floor(x);
Fy=floor(y);
Cx=ceil(x);
Cy=ceil(y);
if Fx>0 && Cx<(n+1) && Fy>0 && Cy<(m+1) % regular case if pixel is in interior of image
B(j,i)=di*dj*A(Fy,Fx)+(1-di)*dj*A(Cy,Fx) + di*(1-dj)*A(Fy,Cx) + (1-di)*(1-dj)*A(Cy,Cx);
elseif (Fx==0&&Fy==0) % if pixel is top left corner of image
B(j,i)=A(1,1);
elseif (Cx==(n+1)&&Cy==(m+1)) % if pixel is bottom right corner of image
B(j,i)=A(m,n);
elseif (Cx==(n+1)&&Fy==0) % if pixel is top right corner of image
B(j,i)=A(1,n);
elseif (Fx==0&&Cy==(m+1)) % if pixel is bottom left corner of image
B(j,i)=A(m,1);
elseif Fx==0 % if pixel is on leftmost boundary of image
B(j,i)=di*A(Fy,Cx) + (1-di)*A(Cy,Cx);
elseif Fy==0 % if pixel is on top boundary of image
B(j,i)=dj*A(Cy,Fx) + (1-dj)*A(Cy,Cx);
elseif Cx==(n+1) % if pixel is on rightmost boundary of image
B(j,i)=di*A(Fy,Fx)+(1-di)*A(Cy,Fx);
elseif Cy==(m+1) % if pixel is on bottom boundary of image
B(j,i)=dj*A(Fy,Fx)+ (1-dj)*A(Fy,Cx);
end
end
end
elseif strcmp(InterpolationMethod,'bic4')
%Anti Alising Filter
if ResizeRatioX<1
A = ReSample(A,ResizeRatioX);
end
if ResizeRatioY<1
A = ReSample(A',ResizeRatioY);
A =A';
end
% The diagrams below illustrate the pixels involved in one-dimensional
% interpolation. Point s0 is the interpolation kernel key position.
% xfrac and yfrac, indicated by the dots, represent the point of
% interpolation between two pixels
% Horizontal Vertical
%
% s_ s0 . s1 s2 s_
% ^
% xfrac s0
% .< yfrac
% s1
%
% s2
%
% Point s00 is reference for xfract and yfract
%
% s__ s_0 s_1 s_2
%
%
%
% s0_ s00 s01 s02
%
% . < yfract
%
% s1_ s10 s11 s12
%
%
%
% s2_ s20 s21 s22
% ^
% xfrac
%
%
% in this case dj = xfract & di = yfract
afact = -0.5;
for i=1:NewRows
for j=1:NewColumns
temp = [i j 1]*inv(a);
di=temp(1)-floor(temp(1)); % di = yfract
dj=temp(2)-floor(temp(2)); % dj = xfract
y = floor(temp(1));
x = floor(temp(2));
% Since Bicubic interpolation requires 4x4 points that means
% with s(i,j) as my refference pixel I require 16 pixels
% from s(i-1:i+2,j-1:j+2) the following code makes sure that
% if the pixel lies on or outside a boundary then the value
% such as s(0,0) replicates s(1,1) as s(0,0) does not exist
% Similarly Bicubic 2x2,6x6,8x8 depend on 4,36 & 64 pixels
% respectively . Even for their implementationa the replicate
% method is used
aX=max(x-1,1);
bX=max(x,1);
cX=min(x+1,n);
dX=min(x+2,n);
aY=max(y-1,1);
bY=max(y,1);
cY=min(y+1,m);
dY=min(y+2,m);
% first calling for each of the 4 rows (X direction
% interpolation)
temp1=bicubic4x4([A(aY,aX),A(aY,bX),A(aY,cX),A(aY,dX)],dj,afact);
temp2=bicubic4x4([A(bY,aX),A(bY,bX),A(bY,cX),A(bY,dX)],dj,afact);
temp3=bicubic4x4([A(cY,aX),A(cY,bX),A(cY,cX),A(cY,dX)],dj,afact);
temp4=bicubic4x4([A(dY,aX),A(dY,bX),A(dY,cX),A(dY,dX)],dj,afact);
% Then calling bicubic for these 4 interpolated values (Y
% direction interpolation)
B(i,j)=bicubic4x4([temp1,temp2,temp3,temp4],di,afact);
end
end
elseif strcmp(InterpolationMethod,'bic2')
%Anti Alising Filter
if ResizeRatioX<1
A = ReSample(A,ResizeRatioX);
end
if ResizeRatioY<1
A = ReSample(A',ResizeRatioY);
A =A';
end
for i=1:NewRows
for j=1:NewColumns
temp = [i j 1]*inv(a);
di=temp(1)-floor(temp(1));
dj=temp(2)-floor(temp(2));
y = floor(temp(1));
x = floor(temp(2));
aX=max(x,1);
bX=min(x+1,n);
aY=max(y,1);
bY=min(y+1,m);
% first calling for each of the 2 rows (X direction
% interpolation)
temp1=bicubic2x2([A(aY,aX),A(aY,bX)],dj);
temp2=bicubic2x2([A(bY,aX),A(bY,bX)],dj);
% Then calling bicubic for these 2 interpolated values (Y
% direction interpolation)
B(i,j)=bicubic2x2([temp1,temp2],di);
end
end
elseif strcmp(InterpolationMethod,'bic6')
%Anti Alising Filter
if ResizeRatioX<1
A = ReSample(A,ResizeRatioX);
end
if ResizeRatioY<1
A = ReSample(A',ResizeRatioY);
A =A';
end
for i=1:NewRows
for j=1:NewColumns
temp = [i j 1]*inv(a);
di=temp(1)-floor(temp(1));
dj=temp(2)-floor(temp(2));
y = floor(temp(1));
x = floor(temp(2));
aX=max(x-2,1);
bX=max(x-1,1);
cX=max(x,1);
dX=min(x+1,n);
eX=min(x+2,n);
fX=min(x+3,n);
aY=max(y-2,1);
bY=max(y-1,1);
cY=max(y,1);
dY=min(y+1,m);
eY=min(y+2,m);
fY=min(y+3,m);
% first calling for each of the 6 rows (X direction
% interpolation)
temp1=bicubic6x6([A(aY,aX),A(aY,bX),A(aY,cX),A(aY,dX),A(aY,eX),A(aY,fX)],dj);
temp2=bicubic6x6([A(bY,aX),A(bY,bX),A(bY,cX),A(bY,dX),A(bY,eX),A(bY,fX)],dj);
temp3=bicubic6x6([A(cY,aX),A(cY,bX),A(cY,cX),A(cY,dX),A(cY,eX),A(cY,fX)],dj);
temp4=bicubic6x6([A(dY,aX),A(dY,bX),A(dY,cX),A(dY,dX),A(dY,eX),A(dY,fX)],dj);
temp5=bicubic6x6([A(eY,aX),A(eY,bX),A(eY,cX),A(eY,dX),A(eY,eX),A(eY,fX)],dj);
temp6=bicubic6x6([A(fY,aX),A(fY,bX),A(fY,cX),A(fY,dX),A(fY,eX),A(fY,fX)],dj);
% Then calling bicubic for these 6 interpolated values (Y
% direction interpolation)s
B(i,j)=bicubic6x6([temp1,temp2,temp3,temp4,temp5,temp6],di);
end
end
elseif strcmp(InterpolationMethod,'bic8')
%Anti Alising Filter
if ResizeRatioX<1
A = ReSample(A,ResizeRatioX);
end
if ResizeRatioY<1
A = ReSample(A',ResizeRatioY);
A =A';
end
for i=1:NewRows
for j=1:NewColumns
temp = [i j 1]*inv(a);
di=temp(1)-floor(temp(1));
dj=temp(2)-floor(temp(2));
y = floor(temp(1));
x = floor(temp(2));
aX=max(x-3,1);
bX=max(x-2,1);
cX=max(x-1,1);
dX=max(x,1);
eX=min(x+1,n);
fX=min(x+2,n);
gX=min(x+3,n);
hX=min(x+4,n);
aY=max(y-3,1);
bY=max(y-2,1);
cY=max(y-1,1);
dY=max(y,1);
eY=min(y+1,m);
fY=min(y+2,m);
gY=min(y+3,m);
hY=min(y+4,m);
% first calling for each of the 8 rows (X direction
% interpolation)
temp1=bicubic8x8([A(aY,aX),A(aY,bX),A(aY,cX),A(aY,dX),A(aY,eX),A(aY,fX),A(aY,gX),A(aY,hX)],dj);
temp2=bicubic8x8([A(bY,aX),A(bY,bX),A(bY,cX),A(bY,dX),A(bY,eX),A(bY,fX),A(bY,gX),A(bY,hX)],dj);
temp3=bicubic8x8([A(cY,aX),A(cY,bX),A(cY,cX),A(cY,dX),A(cY,eX),A(cY,fX),A(cY,gX),A(cY,hX)],dj);
temp4=bicubic8x8([A(dY,aX),A(dY,bX),A(dY,cX),A(dY,dX),A(dY,eX),A(dY,fX),A(dY,gX),A(dY,hX)],dj);
temp5=bicubic8x8([A(eY,aX),A(eY,bX),A(eY,cX),A(eY,dX),A(eY,eX),A(eY,fX),A(eY,gX),A(eY,hX)],dj);
temp6=bicubic8x8([A(fY,aX),A(fY,bX),A(fY,cX),A(fY,dX),A(fY,eX),A(fY,fX),A(fY,gX),A(fY,hX)],dj);
temp7=bicubic8x8([A(gY,aX),A(gY,bX),A(gY,cX),A(gY,dX),A(gY,eX),A(gY,fX),A(gY,gX),A(gY,hX)],dj);
temp8=bicubic8x8([A(hY,aX),A(hY,bX),A(hY,cX),A(hY,dX),A(hY,eX),A(hY,fX),A(hY,gX),A(hY,hX)],dj);
% Then calling bicubic for these 8 interpolated values (Y
% direction interpolation)
B(i,j)=bicubic8x8([temp1,temp2,temp3,temp4,temp5,temp6,temp7,temp8],di);
end
end
end
if ClassFlag
B = uint8(B*255);
end
return;
% ------------------------------------------------ %
% This is the final bicubic interpolation polynomial for 4x4 bicubic
% interpolation
function [X] = bicubic4x4(v,fact,afact)
X = 0;
if fact==0
X = v(2);
return
end
Y = [1 0 1 2];
Z = [1 1 -1 -1];
% r(x) = (a + 2)|x|^3 - (a + 3)|x|^2 + 1 , 0 <= |x| < 1
% r(x) = a|x|^3 - 5a|x|^2 + 8a|x| - 4a , 1 <= |x| < 2
% where x is {fact,fact+1,abs(fact-1),abs(fact-2)}
% then r(x) will contain the interpolation ratio which will be multiplied
% by the point value.
% The sum of r(fact+1)*v1+r(fact)*v2+r(abs(fact-1))*v3+r(abs(fact-2))*v4
% will give the interpolated value.
% The value a is taken as -0.5 as per Matlab. Play
% around with the value of a and you get some interesting results.
% ALL THE CUBIC EQUATION HAVE BEEN PICKED FROM RESEARCH PAPERS I DO NOT
% KNOW THEIR DERIVATION OR WHY THE FACTOR a EXISTS IN 4x4 but not in the
% others. The only reason I can think of is that for 6x6 & 8x8 bicubic
% equations are continous till the second differential that
% is : F(x-) = F(x+), F'(x-) = F'(x+) & is F''(x-) = F''(x+)
% while 4x4 bicubic is only : F(x-) = F(x+), F'(x-) = F'(x+)
% a = -0.5;
a = afact;
A=(2+a);
B=-(3+a);
C=1;
D=a;
E=-5*a;
F=8*a;
G=-4*a;
for i = 1:4
if v(i)~=0
Fi = Y(i)+fact*Z(i);
if Fi<1
X = X + v(i)*(1 + Fi*Fi*(B + Fi*A));
elseif Fi<2
X = X + v(i)*(G + Fi*(F + Fi*(E + Fi*D)));
end
end
end
% ------------------------------------------------ %
function [X] = bicubic2x2(v,fact)
X = 0;
if fact==0
X = v(1);
return
end
% 2x2 bicubic polynomial
% r(x) = 2|x|^3 - 3|x|^2 + 1 , 0 <= |x| < 1
% The sum of r(fact)*v1+r(1-fact)*v2
A = 2;
B = -3;
C = 1;
Fi = fact;
if v(1) ~=0
X = X +v(1)*(C + Fi*Fi*(B+Fi*A));
end
Fi = 1-fact;
if v(2) ~=0
X = X +v(2)*(C + Fi*Fi*(B+Fi*A));
end
return;
% ------------------------------------------------ %
function [X] = bicubic6x6(v,fact)
X = 0;
if fact==0
X = v(3);
return
end
Y = [2 1 0 1 2 3];
Z = [1 1 1 -1 -1 -1];
% 6x6 bicubic polynomial
% r(x) = (6/5)|x|^3 - (11/5)|x|^2 + 1 , 0 <= |x| < 1
% r(x) = -(3/5)|x|^3 + (16/5)|x|^2 -(27/5)|x| + 14/5 , 1 <= |x| < 2
% r(x) = (1/5)|x|^3 - (8/5)|x|^2 +(21/5)|x| - 18/5 , 2 <= |x| < 3
%
% The sum of r(2+fact)*v1 + r(1+fact)*v2 + r(fact)*v3 + r(1-fact)*v4 +
% r(2-fact)*v5 + r(3-fact)*v6 represents the interpolated value
A = 6/5;
B = -11/5;
C = 1;
D = -3/5;
E = 16/5;
F = -27/5;
G = 14/5;
H = 1/5;
I = -8/5;
J = 21/5;
K = -18/5;
for i = 1:6
if v(i)~=0
Fi = Y(i)+fact*Z(i);
if Fi<1
X = X + v(i)*(1 + Fi*Fi*(B + Fi*A));
elseif Fi<2
X = X + v(i)*(G + Fi*(F + Fi*(E + Fi*D)));
elseif Fi<3
X = X + v(i)*(K + Fi*(J + Fi*(I + Fi*H)));
end
end
end
return;
% ------------------------------------------------ %
function [X] = bicubic8x8(v,fact)
X = 0;
if fact==0
X = v(4);
return
end
Y = [3 2 1 0 1 2 3 4];
Z = [1 1 1 1 -1 -1 -1 -1];
% 8x8 bicubic polynomial
% r(x) = (67/56)|x|^3 - (123/56)|x|^2 + 1 , 0 <= |x| < 1
% r(x) = -(33/56)|x|^3 + (177/56)|x|^2 - (75/14)|x| + (39/14), 1 <= |x| < 2
% r(x) = (9/56)|x|^3 - (75/56)|x|^2 + (51/14)|x| - (45/14), 2 <= |x| < 3
% r(x) = -(3/56)|x|^3 + (33/56)|x|^2 - (15/7)|x| + (18/7) , 3 <= |x| < 4
% The sum of r(3+fact)*v1 + r(2+fact)*v2 + r(1+fact)*v3 + r(fact)*v4 +
% r(1-fact)*v5 + r(2-fact)*v6 + r(3-fact)*v7 + r(4-fact)*v8
A = 67/56;
B = -123/56;
C = 1;
D = -33/56;
E = 177/56;
F = -75/14;
G = 39/14;
H = 9/56;
I = -75/56;
J = 51/14;
K = -45/14;
L = -3/56;
M = 33/56;
N = -15/7;
O = 18/7;
for i = 1:8
if v(i)~=0
Fi = Y(i)+fact*Z(i);
if Fi<1
X = X + v(i)*(1 + Fi*Fi*(B + Fi*A));
elseif Fi<2
X = X + v(i)*(G + Fi*(F + Fi*(E + Fi*D)));
elseif Fi<3
X = X + v(i)*(K + Fi*(J + Fi*(I + Fi*H)));
elseif Fi<4
X = X + v(i)*(O + Fi*(N + Fi*(M + Fi*L)));
end
end
end
return;
% ------------------------------------------------ %
% lowpass filter default size 11
% Anti Aliasing Filter (GOOGLE IT)
function [A] =ReSample(A,ResizeRatio)
[m n] =size(A);
N=11;
Wn = ResizeRatio;
vec = 1:5;
vec2 = vec*pi;
wind = zeros(1,5);
by =zeros(1,6);
for i = 1:5
wind(i) = .54-.46*cos(2*pi*(vec(i)-1)/(N-1));
end
for i = 1:5
by(1,6-i) = (sin(Wn*vec2(i))/vec2(i))*wind(6-i);
by(1,6+i) = by(1,6-i);
end
by(1,6) = ResizeRatio;
by = by/abs(sum(by));
temp = ones(1,5);
for i=1:m
Anew(i,:) = [temp*(A(i,1)) A(i,:) temp*(A(i,m))];
end
A=Anew;
Anew = []; % free memory
B(1:m,1:n)=0;
for i = 1:m
for j = 6:n+5
for k =1:11
B(i,j-5) = B(i,j-5)+by(k)*A(i,j-6+k);
end
end
end
A=B;
return;