%
% Optical Imaging and Spectroscopy
%
% David J. Brady
% Duke University
% www.opticalimaging.org
%
% Figure 8.17 use of TWIST algorithm with truncated SVD
%
% The Two-step iterative shrinkage/thresholding algorithm was developed by
% Jose Biocaus-Dias and Mario Figueiredo. Code is posted at
% http://www.lx.it.pt/~bioucas/TwIST/TwIST.htm
%
% This file is a modificaiton of the demo_l2_tv.m from the original
% distribution
%
% Shift Code demo
%
% This demo illustrates the twIST algorithm computing the solution  of  
%
%     xe = arg min 0.5*||A x-y||^2 + tau TV(x),
%             x
%  where x is an image, Ax is the Fourier Transform of x  computed in a 
%  "small" set of frequencies, y is a noisy observation and TV(x) is
%  the total variation of x.
% 
% The proximal operator, i.e., the solution of
%
%     xe = arg min 0.5 ||x-y||^2 + tau TV(x)
%             x
% is  given by the Chambolle algorithm
% 
% A. Chambolle, “An algorithm for total variation minimization and
% applications,” Journal of Mathematical Imaging and Vision, vol. 20,
% pp. 89-97, 2004.
%
%
% For further details about the TwIST algorithm, see the paper:
%
% J. Bioucas-Dias and M. Figueiredo, "A New TwIST: Two-Step
% Iterative Shrinkage/Thresholding Algorithms for Image 
% Restoration",  IEEE Transactions on Image processing, 2007.
%
%
% (available at   http://www.lx.it.pt/~bioucas/publications.html)
%
%
% Authors: Jose Bioucas-Dias and Mario Figueiredo, 
% Instituto Superior Técnico, October, 2007
%
%
%
% shift code component added by D. Brady, May 2008
%
close all;

N=4;
sigma = 0.01;
lambda=0.3;
% Measurement noise standard deviation
% generate phantom
%f = phantom(n);
f=abs(phantom(256));
% load flowers;
% f=dog;clear dog;
f=f/max(max(f));
% f(1:20,:)=0;
% f(:,237:256)=0;
% f(:,1:20)=0;
% f(237:256,:)=0;
f=double(f);
%
%
% generate measurement operator
Hl=toeplitz([1 zeros(1,255)],[ones(1, N),zeros(1,256-N)])/N;

[al bl cl]=svd(Hl);
phl=pinv(Hl);
rl=125;
hl=al(:,1:rl)*bl(1:rl,1:rl)*cl(:,1:rl)';
hlpi=cl(:,1:rl)*pinv(bl(1:rl,1:rl))*al(:,1:rl)';
% create handle for observation operator and transpose
hR = @(x)  hl*x*hl';
hRT = @(x) hlpi*x*hlpi';

% vector of observations
y_noiseless = Hl*f*Hl';
y = y_noiseless + sigma*randn(size(y_noiseless));

figure(10); set(gcf,'color','white');
% subplot(3,2,1);imagesc(f)
colormap(gray)
% axis image
% axis off
% title('Original')
% drawnow

% ty=hRT(y);
% subplot(3,2,3);imagesc(ty,[0 1])
% axis image
% axis off
% title(['truncated SVD (a) mse=' num2str(mean(mean((ty-f).^2)),'%5.2e')]);
% drawnow
% % tikhonov regularized image
% dNoise=pinv(Hl'*Hl+lambda^2*eye(256))*Hl'*y*Hl*pinv(Hl*Hl'+lambda^2*eye(256));
% subplot(3,2,2);imagesc(dNoise,[0 1])
% axis image
% axis off
% title(['Tikhonov (a) mse=' num2str(mean(mean((dNoise-f).^2)),'%5.2e')]);
% drawnow
%
% project the observations onto the data space

% set the denoising function
Psi = @(x,th) mex_vartotale(x,th,'itmax',10);

% set the penalty function, to compute the objective
Phi = @(x) TVnorm(x);

% regularization parameters (empirical)
tau = sqrt(3)*sigma;
tau = 0.001;
y=al(:,1:rl)*al(:,1:rl)'*y*al(:,1:rl)*al(:,1:rl)';
tolA = 1e-4;
% -- TwIST ---------------------------
% stop criterium:  the relative change in the objective function 
% falls below 'ToleranceA'         
 [x_twist,dummy,obj_twist,...
    times_twist,dummy,mse_twist]= ...
         TwIST(y,hR,tau,...
         'Lambda', 1e-3, ...
         'AT', hRT, ...
         'Psi', Psi, ...
         'Phi',Phi, ...
         'True_x', f, ...
         'Monotone',1,...
         'MaxiterA', 10000, ...
         'Initialization',0,...
         'StopCriterion',1,...
       	 'ToleranceA',tolA,...
         'Verbose', 1);
   
% [x_twist,dummy,obj_twist,...
%     times_twist,dummy,mse_twist]= ...
%          SpaRSA(y,hR,tau,...
%          'AT', hRT, ...
%          'Psi', Psi, ...
%          'Phi',Phi, ...
%          'True_x', f, ...
%          'BB_factor', 0.8, ...
%          'Monotone',1,...
%          'Initialization',0,...
%          'StopCriterion',4,...
%        	 'ToleranceA',obj_twist(end),...
%          'Verbose', 1);

subplot(2,2,1);imagesc(x_twist,[0 1]);
axis image
title(['(a) mse=' num2str(mean(mean((x_twist-f).^2)),'%5.2e')]); axis off;
subplot(2,2,2);imagesc(x_twist(100:150,100:150));axis off; axis image;

% subplot(3,2,2);plot(diag(bl),'-k'); title('singular value spectra');
% 
% hold on;

Hl=toeplitz([1 zeros(1,255)],[1 1 1 0 0 1 0 0,zeros(1,256-8)])/4;
[al bl cl]=svd(Hl);
phl=pinv(Hl);
hl=al(:,1:rl)*bl(1:rl,1:rl)*cl(:,1:rl)';
hlpi=cl(:,1:rl)*pinv(bl(1:rl,1:rl))*al(:,1:rl)';
% create handle for observation operator and transpose
hR = @(x)  hl*x*hl';
hRT = @(x) hlpi*x*hlpi';

% vector of observations
y_noiseless = hR(f);
y = y_noiseless + sigma*randn(size(y_noiseless));

% ty=hRT(y);
% subplot(3,2,5);imagesc(ty,[0 1])
% axis image
% axis off
% title(['truncated SVD (b) mse=' num2str(mean(mean((ty-f).^2)),'%5.2e')])
% drawnow

% set the denoising function
Psi = @(x,th) mex_vartotale(x,th,'itmax',10);

% set the penalty function, to compute the objective
Phi = @(x) TVnorm(x);

% regularization parameters (empirical)
tau = sqrt(3)*sigma;
tau = 0.001;
y=al(:,1:rl)*al(:,1:rl)'*y*al(:,1:rl)*al(:,1:rl)';
tolA = 1e-4;
% -- TwIST ---------------------------
% stop criterium:  the relative change in the objective function 
% falls below 'ToleranceA'         
 [x_twist,dummy,obj_twist,...
    times_twist,dummy,mse_twist]= ...
         TwIST(y,hR,tau,...
         'Lambda', 1e-3, ...
         'AT', hRT, ...
         'Psi', Psi, ...
         'Phi',Phi, ...
         'True_x', f, ...
         'Monotone',1,...
         'MaxiterA', 10000, ...
         'Initialization',0,...
         'StopCriterion',1,...
       	 'ToleranceA',tolA,...
         'Verbose', 1);
     
     
subplot(2,2,3);imagesc(x_twist,[0 1]);
axis image
title(['(b) mse=' num2str(mean(mean((x_twist-f).^2)),'%5.2e')]); axis off;
subplot(2,2,4);imagesc(x_twist(100:150,100:150));axis off; axis image;
% subplot(3,2,2);plot(diag(bl),'-k');