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

% load Data\Phantom22Ls001Noise.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)+noise.
%   WT = []; W = [];
%   U0 is the back-project solution from B.

F = double(imread('pics/3DMR_Chest.jpg'));
F=imresize(F, [256 256]);
F=im256(F);

F2=ifft2(fft2(F));
figure; subplot(1,2,1); imshow(F, []);

[m,n] = size(F);
FB = fft2(F)/sqrt(m*n);
% FB = fft2(F);
tmp=randperm(m*n); k=256*256/4;
picks=zeros(m, n); picks(tmp(1:k))=1;
picks=logical(picks);
B = FB(picks)+0.0000*randn(k, 1);

aTV = 1e-3;   % 1e-4 is good for noisy (sigma = .01) under-sampled phantom

%%---------------------
aL1 = 0.035;       % 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 = 100; 
opts.gamma = 1.0;
opts.beta = 10;
opts.relchg_tol = 5e-4;
opts.real_sol = true;    % return solution of real type (i.e., U = real(U))

opts.normalize = false;  % 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)));
