% 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);
specIm(5:11,(Nx/4):(3*Nx/4))=1;
%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
%
%Shear 1
%
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;
%
%
% Hmat is Nx \times (Nl*Nx) matrix
%
% apply t to each spectral channel (punch)
%
H=kron(eye(Nl),diag(t))*H;
%
% 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;
clear Hs;
HpreS=reshape(H*ones(Nl*Nx,1),Nx,Nl)';
%
% presmash data
%
subplot(4,1,1);imagesc(reshape(H*reshape(f',Nl*Nx,1),Nx,Nl)');
title('(a)');
%
% 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,2);plot([t; y'+2]','-k');axis([0 256 0 10]);
title('(b)');

% 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,3 );imagesc(x_twist);
title('(c)'); axis off;
%
% alternative local reconstruction
%
%
xl=zeros(Nl,Nx);
np=20;
for ct=(np/2):(Nx-np/2)
    xl(:,ct)=pinv(HpreS(:,(ct-(np/2-1)):(ct+np/2))')*y((ct-(np/2-1)):(ct+np/2));
end
subplot(4,1,4 );imagesc(xl);
title('(d)'); axis off;
%