%-------------------------------------------------------%
% Reconstruction program for complex-valued pMRI datas  %
%                                                       %
%-------------------------------------------------------%


clc
close all
clear all

addpath(genpath('Tools'));
addpath(genpath('Data'));
addpath(genpath('Algo'));
addpath(genpath('Wavelab850'));

%Acceleration factor
R = 5;

%Load sampling matrix Sigma
load('SigmaUNI.mat');

figure
imagesc(Sigma); colormap gray; title('Sampling matrix');
axis image off

%Load datas I, S and Mask
load('Data82.mat');

figure
imagesc(abs(I)); colormap jet; title('Modulus of original image');
axis image off

figure
imagesc(abs(S(:,:,1))); colormap jet; title('Modulus of 4th sensitivity matrix');
axis image off

figure
imagesc(Mask); colormap gray; title('Mask');
axis image off

%Image size (should be squared for wavelet transform!!)
[Nx,Ny] = size(I);

%Number of coils
L = size(S,3);

%Number of samples
D = length(find(Sigma==1));

%noise covariance matrix
s2 = 6e9; Lambda = eye(L)*s2;

%Create noisy data
dclean = zeros(D,L);   
for l = 1:L
        dltemp = fft2(S(:,:,l).*I)/sqrt(Nx*Ny);
        dclean(:,l) = subsamp(dltemp,Sigma);
end
d = addcomplexnoisev(dclean,Lambda);

%wavelet parameters
qmf = MakeONFilter('Symmlet',4);%wavelet basis
jm = 3;%decomposition level

%% RUN 3MG ALGORITHM

%penalty function
%(1) psi(u) = (1-exp(-u.^2./(2*delta^2)));  %NON CONVEX - NON COERCIVE
%(2) psi(u) = (u.^2)./(2*delta^2 + u.^2);   %NON CONVEX - NON COERCIVE
%(3) psi(u) = log(1 + (u.^2)./(delta^2));   %NON CONVEX - COERCIVE
%(4) psi(u) = sqrt(1 + u^2/delta^2)-1;      %CONVEX 
%(5) psi(u) = 1/2 u^2;                      %CONVEX
psi = 4;
%regularization parameters
la = 0.006;
delta = 1700;
%max iteration number
kmax = 500;
%initialization
xinit  =  zeros(Nx,Ny);

%run algorithm
[snr,x,Iout,F,Time] = algoritm3MG(kmax,d,jm,qmf,la,delta,Lambda,S,Sigma,psi,I,xinit,Mask);


%display results

figure
imagesc(abs(Iout),[0 max(abs(I(:)))]); colormap jet; title(['Modulus of restored image with snr = ',num2str(snr(end))]);
axis image off

figure
plot(cumsum(Time),snr); xlabel('Time (s.)'); ylabel('SNR');

    
