% Optical Imaging and Spectroscopy
%
% David J. Brady
% Duke University
% www.opticalimaging.org
%
% Figure 10.54
%
% modification of http://www.lx.it.pt/~bioucas/TwIST/TwIST.htm
%
%function demo_l2_TV
%
%
% 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
%
%
%
% cassi component added by D. Brady, July 2008
%

% Measurement noise standard deviation
% generate phantom
%f = phantom(n);
rand('seed',333);
Nl=16;Nx=256;
t=zeros(1,Nx);
t(16:16:Nx)=1;
%t=round(rand(1,Nx));
spec=zeros(1,Nx);spec((Nx/4):(3*Nx/4))=1;
ls=zeros(1,Nl);ls(5)=1;ls(11)=1;
specIm=spec'*ls;specIm=specIm';
specIm=zeros(Nl,Nx);
%uncomment the next line for fig. 10.50
specIm(5:11,(Nx/4):(3*Nx/4))=1;
% uncomment the next two lines for fig. 10.54
%specIm(5,(Nx/2+10):(5*Nx/6))=1;
%specIm(12,(Nx/6-5):(Nx/2-10))=1;
f=specIm;
f=f/max(max(f));
%
% plot data
%
figure(1);set(gcf,'color','white');colormap gray;
subplot(4,1,1);imagesc(f);
title('(a)');
xlabel('j');ylabel('n');
%
% measurement matrix
%
%
%
% Hmat is Nx \times (Nl*Nx) matrix
%
% apply t to each spectral channel (punch)
%
H=kron(eye(Nl),diag(t));
%
% shift (shear)
%
Hs=zeros(Nl*Nx,Nl*Nx);
for ct=1:Nl;
Hs((Nx*(ct-1)+(1:Nx)),(Nx*(ct-1)+(1:Nx)))=circshift(eye(Nx),ct);
end
H=Hs*H;
subplot(4,1,2);imagesc(reshape(H*reshape(f',Nl*Nx,1),Nx,Nl)');
title('(b)');
clear Hs;
%
% presmash data
%
subplot(4,1,2);imagesc(reshape(H*reshape(f',Nl*Nx,1),Nx,Nl)');
title('(b)');
%
% smash
%
% smash operator is Nx \times (Nl*Nx) takes 1 component from each spectrum
%
H=kron(ones(1,Nl),eye(Nx))*H;
ph=pinv(H);
% create handle for observation operator and transpose
hR = @(x) H*reshape(x',Nl*Nx,1);
hRT = @(x) reshape(ph*x,Nx,Nl)';

% vector of observations
%y_noiseless = cassiF(f,t);
%y = y_noiseless + sigma*randn(size(y_noiseless));
y=H*reshape(f',Nl*Nx,1);
subplot(4,1,3);plot([t; y'+2]','-k');axis([0 256 0 3.2]);
title('(c)');

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

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

% regularization parameters (empirical)
tau = 0.1;
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(4,1,4 );imagesc(x_twist);
title('(d)'); axis off;