Tsai and Huang Method Code

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tsai and Huang Multiframe Image Enhancement
% -------------------------------------------
% AUTHOR: Stephen Rose, Maher Khoury
%    DATE: March 12, 1999
% PURPOSE:
%         Computes a SR frame using the method developed by
%         Tsai and Huang
%
%
% Variables:
%   -p              = Number of frames
%   -deltaK         = Matrix containing the x,y shift components
%   -F              = Fourier Transform of the SR image
%   -image          = SR image
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%% Read in and store images %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
p = 2;

for num = 1:p,
temp = pgmRead(strcat('books/i.seq',num2str(num),'.pgm'));
frame(:,:,num) = temp(129:256,129:256);
end

M = size(frame(:,:,1),1);
N = size(frame(:,:,1),2);
[m,n] = meshgrid(1:M,1:N);

%%%%%%%%% Find Translation Parameters %%%%%%%%%%%%%%%%%%%%%
K(1,:) = zeros(1,4);

for num= 2:p,
K(num,:) = affine(frame(:,:,num),frame(:,:,1));
end

deltaK = K(:,3:4);

%%%%%%%%%% Initialization %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Lu = 1;
Lv = 1;
Tx = 2;
Ty = 2;

wsx = 1./Tx;
wsy = 1./Ty;

F1 = fft2(frame(:,:,1));
F2 = fft2(frame(:,:,2));

F = zeros(2*Lu*M,2*Lv*N);

%%%%%%%%%%% Calculating values in matrix F %%%%%%%%%%%%%%%%%%
for n = 1:N,
disp(n);
for m = 1:M,
Gmn = [F1(m,n);F2(m,n)];

    for k=1:p,
      for r = 1:4
PHImn(k,r) = exp(j*2*pi.*(deltaK(k,1)*(m./M*Tx)-(Lu/Tx)+deltaK(k,2)*(n./N*Ty)-(Lv/Ty)))*...
     exp(j*2*pi.*((deltaK(k,1)./Tx)*rem(r-1,2*Lu)+(deltaK(k,2)*floor((r-1)./2*Lv))));
      end
    end

    Fmn = pinv(PHImn)*Gmn;

    F(m,n) = Fmn(4,1);
    F(m+M,n) = Fmn(2,1);
    F(m,n+N) = Fmn(3,1);
    F(m+M,n+N) = Fmn(1,1);

end
end

%%%%%%%%%%% Get f from F and display %%%%%%%%%%%%%%%%%%%%%%%
image = abs(ifft2(fftshift(F)));
figure;
showIm(image);