%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% A solution to CSE For Your Homework Project 1 %
% Image Deblurring: I Can See Clearly Now %
% Dianne P. O'Leary and James G. Nagy 09/02 %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Read in a blurred image G and matrices A and B. %
% Try to reconstruct F, where g is approximately %
% kron(A,B) * reshape(F,m*n,1) and g = reshape(G,m*n,1) %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load proj1data.mat
figure(1)
imagesc(G)
title('Blurred image')
colormap(gray)
pause
figure(2)
[m, n] = size(G);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Compute the singular value decomposition of A and B. %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Ua, Sa, Va] = svd(A);
[Ub, Sb, Vb] = svd(B);
Ghat = Ub'*G*Ua;
S = diag(Sb)*(diag(Sa))';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Use Tikhonov regularization, with parameters %
% specified in lamtest, to reconstruct the image. %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lamtest = [.05 .01 .005 .0025 .0015 .001 .00095 .00016681 .00005];
mytrue = 0;
if (mytrue)
for lam=lamtest
Fhat = (S.*Ghat) ./ (S.*S+lam^2);
F = Vb*Fhat*Va';
imagesc(F)
colormap(gray)
disp(sprintf('Tikhonov lambda= %f',lam))
disp('Strike any key to continue.')
drawnow
pause
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Use TSVD regularization to reconstruct the image. %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Store the entries of S in a linear array s, then sort this
% set in descending order and construct the permutation
% that locates the ordered set in the unsorted list.
s = reshape(S,n^2,1);
[ss,prm] = sort(s);
prm = flipud(prm);
ss = flipud(ss);
ls = length(s);
iprm = zeros(ls,1);
iprm(prm) = [1:ls]';
% Using the sorted set and permutation constructed above,
% remove all but the p largest singular values and then
% compute the truncated expanded solution. Test using a
% set of values of p.
pset = [100, 250, 500, 1000, 1500, 2000, 2500, 4000, 6000, 8000];
for i=1:length(pset),
p = pset(i);
ssnew = [ss(1:p); zeros(length(ss)-p,1)];
Snew = reshape(ssnew(iprm),n,n);
Fhat = Ghat ./ S;
Fnew = Fhat .* (Snew>0);
F = Vb*Fnew*Va';
imagesc(F)
colormap(gray)
disp(sprintf('TSVD p= %d ', p))
disp('Strike any key to continue.')
drawnow
pause;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% Here's an alternative approach to TSVD regularization,%
% which determines p implicitly by throwing out all %
% singular values less than some value, here %
% determined by the parameters in lamtest. %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fhat = Ghat ./ S;
for lam=lamtest
ind = S > lam;
Fnew = Fhat .* ind;
F = Vb*Fnew*Va';
imagesc(F)
disp(sprintf('TSVD lambda = %f, p = %d',lam,sum(sum(ind))))
disp('Strike any key to continue.')
drawnow
pause
end