%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                        %
%  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

imagesc(G)
colormap(gray)
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];

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                        %
%  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)
  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