%
% based, in part, on Example 6.3 in Aster & Borchers,  pp. 133-135
%

% Clear everything out.
clear
rand('state',0);
randn('state',0);

addpath c:\local\latest\519\cd_2.1\regutools   % package by P.C. Hansen
addpath c:\local\latest\519\cd_2.1\Lib
addpath c:\local\latest\519\cd_2.1\Examples\chap6\ex_6_3\

img=double(imread('image.png'));  % Load the raw image in.
img = img(1:200,1:200);
sz = size(img);
r = sz(1); 
c = sz(2);

G=blur(r,3,15);  % Build the G matrix.
facG = mean(sum(G));
G = G/facG;

d=G*reshape(img,r*c,1);     % Compute the blurred image.
dn=d + (1.0e-2/facG)*randn(r*c,1);   % Add noise as appropriate.

% Use CGLS (with no explicit regularization) to deblur.
[X,rho,eta]=cgls(G,dn,200);

%
% Plots with no explicit regulariztion.
%

% Plot raw image
figure(1);
subplot(1,2,1);  bookfonts;
imagesc(img);
   colormap(gray);
   set(gca, 'Units', 'Normalized')
   P=get(gca,'Position');
   set(gca, 'Position', [P(1)-P(3)/3 P(2)/2 P(3)*1.4 P(4)])
    H=gca;
     set(H,'XTick',[]);
     set(H,'YTick',[]);
   title('Raw image')
%%print -deps img.eps

% Plot blurred image with noise
subplot(1,2,2);     bookfonts;
imagesc(reshape(dn,r,c));
   colormap(gray);
   set(gca, 'Units', 'Normalized')
   P=get(gca,'Position');
   set(gca, 'Position', [P(1) - P(3)/5 P(2)/2 P(3)*1.4 P(4)])
   H=gca;
     set(H,'XTick',[]);
     set(H,'YTick',[]);
   title('Blurred image with noise')
%print -deps c6fblur.eps

% Plot CGLS solution with 30 iterations 
figure(2)
subplot(1,2,1);   bookfonts;
imagesc(reshape(X(:,30),r,c));
  colormap(gray);
   set(gca, 'Units', 'Normalized')
   P=get(gca,'Position');
   set(gca, 'Position', [P(1)-P(3)/3 P(2)/2 P(3)*1.4 P(4)])
    H=gca;
     set(H,'XTick',[]);
     set(H,'YTick',[]);
  title('CGLS solution after 30 iter')
%print -deps c6fblur30.eps

% Plot CGLS solution with 100 iterations 
subplot(1,2,2);     bookfonts;
imagesc(reshape(X(:,100),r,c));
   colormap(gray);
      set(gca, 'Units', 'Normalized')
   P=get(gca,'Position');
   set(gca, 'Position', [P(1)-P(3)/5 P(2)/2 P(3)*1.4 P(4)])
    H=gca;
     set(H,'XTick',[]);
     set(H,'YTick',[]);
  title('CGLS solution after 100 iter')
%print -deps c6fblur100.eps

% Plot of model norm vs. residual norm for CGLS solution
figure(5)
clf;
bookfonts;
loglog(rho(1:200),eta(1:200),'ko');
  xlabel('|| Gm-d ||');
  ylabel('|| m ||');
  hold on
  loglog(rho(30),eta(30),'*')
  hold off
  title('model norm vs. residual norm for CGLS solution')
%print -deps c6fblurlcurve.eps

%
% Now, use CGLS to solve the explicitly regularized problem for several
% values of alpha.

%q = 10;
%img = zeros(q,q);
%img(:,5)= ones(q,1);


Lmini1 =  [eye(c-1), zeros(c-1,1)];
Lmini2 =  [zeros(c-1,1), eye(c-1)];
Lmini = Lmini1 - Lmini2;

Lrow = kron(speye(r), sparse(Lmini));

Lcol = [speye(c*(r-1)), sparse(c*(r-1),c)];
Lcol2 = [sparse(c*(r-1),c), speye(c*(r-1))];
Lcol = Lcol - Lcol2;
% spy(Lcol)

L = [Lrow; Lcol];

d = reshape(img,r*c,1);

imnot = L*d;
szDd = length(imnot);
imnot = 0*abs(imnot(1:szDd/2)) + 1*abs(imnot((szDd/2+1):end));

imLshow = reshape(imnot, c,r-1);
szL = size(L);
nrowsL = szL(1);

% imagesc(imLshow)
%   colormap(gray); 
 
% =================================
   
alpha = 0.05;           % alpha = 0.001 is still kinda noisy, and alpha=0.01 is nice but very smooth 

Gtemp=[G;  alpha*L];
dtemp=[dn; zeros(nrowsL,1)];
n = length(dtemp);
szG = size(G);  m = szG(1);

Niter = 3;
Lone = zeros(Niter,1);
sqR = speye(n);

X0 = zeros(r*c,1);
tic
for iter = 1:Niter,
   Gtil = sqR*Gtemp;
   dtil = sqR*dtemp;    
   [X,rho,eta] = cglsom(Gtil,dtil,500, X0);
   Xfinal = X(:,end);
     absres = abs(dtemp - Gtemp*Xfinal);      % *not*  dtil, Gtil !!!
     Lone(iter) = sum(absres);
     restaper = max(sqrt(absres), 0.01);
     diagR = 1./restaper;
    %diagR(1:m) = diagR(1:m)*5;
    % diagR((m+1):end) =   diagR((m+1):end)/mean(diagR((m+1):end));          % rescale the weights so they don't kill the G-fit
     sqR = sparse(1:n, 1:n, diagR);
    %  X0 = Xfinal;
 figure
 imagesc(reshape(X(:,end),r,c));
 colormap gray
end
 toc


