%
% demo_RecPF
%
close all; clear;
addpath('utilities');
addpath('solver');

load Data\Phantom22LsNoNs.mat F m n picks B WT W U0; 
% loaded data:
%   F is the original Shepp-Logan phantom
%   [m,n] = size(F);
%   Let FB = fft2(F)/sqrt(m*n) be the 2D Fourier transform of F
%   then, B = FB(picks).
%   WT = []; W = [];
%   U0 is the back-project solution from B.



[m, n]=size(F);
FB2 = fft2(F)/sqrt(m*n);
B2= FB2(picks);
Diff=B2-B;

F = double(imread('pics/3DMR_Chest.jpg'));
F=imresize(F, [256 256]);
F=im256(F);
[m,n] = size(F);
FB = fft2(F)/sqrt(m*n);
% FB = fft2(F);
tmp=randperm(m*n); k=256*256/1.5;
% picks=zeros(m, n); picks(tmp(1:k))=1;
% picks=logical(picks);
B = FB(picks);



aTV = 1e-10;   % total variation penalty parameter of TV(u)
               % 1e-4 is good for noisy (sigma = .01) under-sampled phantom

%%---------------------
aL1 = 0;       % penalty parameter of |PsiT*U|_1
% For wavelets, download Rice Wavelet Tools version 2.4, and unzip at ./rwt
% addpath('rwt');
wav = daubcqf(2);
W = @(x) midwt(x,wav);
WT = @(x) mdwt(x,wav);
%%---------------------

opts = [];
opts.maxItr = 50;        % max # of iterations
opts.gamma = 1.6;        % noiseless choice = 1.6
opts.beta = 8;           % noiseless choice = 200
opts.relchg_tol = 1e-5;  % stopping tolerance based on relative change
opts.real_sol = true;    % return solution of real type (i.e., U = real(U))

opts.normalize = true;  % New parameter/data normalization was added to make 
                        % parameters rather independent of the image size, 
                        % pixel intensity range, and number of CS measurements.


pick=false(m,n);pick(picks)=true; % convert "picks" to a logical array "pick"

disp('RecPF is running ...');
tic
[U,Out_RecPF] = RecPF(m,n,aTV,aL1,pick,B,2,opts,WT,W,range(F(:)),F);
toc
disp('done!');


% plot image
figure(1);
subplot(131); imshow(F,[]); title('Original');
subplot(132); imshow(U0,[]); title('Back projection');
subplot(133); imshow(U,[]); title(sprintf('Recon. RelErr=%4.2f%% SNR=%4.1f',norm(U(:)-F(:))*100/norm(F(:)),snr(U,F)));
