Irani and Peleg Method Matlab Code

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IraniPeleg
% ---------------
% AUTHOR: Stephen Rose, Maher Khoury
%    DATE: March 1, 1999
% PURPOSE:
%         Computes a SR frame using the method developped by
%         Irani and Peleg
%
%
% Notes:
%   -The initial guess for the SR frame is done by interpolating
%    the first frame
%   -Assuming a gaussian PSF for h and p
%   -Replaced the calculation function for e by RMSE
%   -Replaced while loop that minimizes e with finite number of
%    iterations (it was converging but way too slowly)
%
% Variables:
%   -ref            = LR reference frame
%   -upref          = HR reference frame
%   -NumberOfFrames = Number of pixel frames to consider
%   -h              = Blurring operator (Gaussian PSF)
%   -p              = "backprojection" kernel
%   -k              = affine parameters for forward projection
%   -inv_k          = affine parameters for backward projection
%   -e              = error function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%% Initialization %%%%%%
NumberOfFrames = 5;
ref = pgmRead('frame00.pgm');

%%% Getting an initial guess %%%
[x,y] = meshgrid(1:size(ref,2),1:size(ref,1));
[X,Y] = meshgrid(1:2.*size(ref,2),1:2.*size(ref,1));
upref = interp2(x,y,ref,X./2,Y./2,'linear');
upref(isnan(upref)) = 0;

figure;
showIm(upref);
drawnow;

%%% Defining the filter kernels and normalizing them%%%
h = mkGaussian(5);
h = h./sum(sum(h));
p = mkGaussian(5,.5);
p = p./sum(sum(p));

disp('Done Initializing');

%%% Calculating the affine parameters %%%
disp('Calculating the affine parameters');
k = zeros(NumberOfFrames,4);
inv_k = zeros(NumberOfFrames,4);

for num = 1:NumberOfFrames,

if (num<10)
frame(:,:,num) = pgmRead(strcat('frame0',num2str(num),'.pgm'));
else
frame(:,:,num) = pgmRead(strcat('frame',num2str(num),'.pgm'));
end

k(num,:) = affine(frame(:,:,num),ref);
inv_k(num,:) = affine(ref,frame(:,:,num));

end
disp('Done');

%%% Iterating until error is minimal
for count = 1:10,

cumWarped = zeros(size(upref));
cumErr = 0;

%%% Computing the LR model %%%
g = conv2(warp(upref,k),h);
g = g(3:size(g,1)-2,3:size(g,2)-2);
g = corrDn(g,namedFilter('binom5'),'reflect1',[2 2]);

for fr = 1:NumberOfFrames,
    f = upConv((frame(:,:,fr)-g),namedFilter('binom5'),'reflect1',[2 2]);
    f = conv2(f,p);
    f = f(3:size(f,1)-2,3:size(f,2)-2);

cumWarped = cumWarped + warp(f,inv_k);
end

%%% Computing the SR frame %%%
upref = upref + (1./NumberOfFrames).*cumWarped;
upref(upref>255) = 255;
upref(upref<0) = 0;

%%% Calculate the cummulative error %%%
for fr = 1:NumberOfFrames,
cumErr = cumErr + mse(g,frame(:,:,fr));
end

e = sqrt((1./NumberOfFrames).*cumErr)

end

pgmWrite(upref,'SRframe.pgm');