lambda1=0.003; %weight of similarity between slices 
lambda2=0.001; %weight of similarity of average to thick slice
lambda3=0.005;%weight of sparsity in wavelet

%% Load data
load data 
%Normalize ground truth, high SNR, thin slices data
ground_truth_slice1=ground_truth_slice1./max(ground_truth_slice1(:));
ground_truth_slice2=ground_truth_slice2./max(ground_truth_slice2(:));
%Normalize acquired data (2 high resolution, low SNR thin slices and a
%single high SNR thick slice
high_res_low_snr_slice1=high_res_low_snr_slice1./max(high_res_low_snr_slice1(:)); 
high_res_low_snr_slice2=high_res_low_snr_slice2./max(high_res_low_snr_slice2(:)); 
low_res_high_snr_slice=low_res_high_snr_slice./max(low_res_high_snr_slice(:)); 



%% Use acquired data in k-space domain and solve the minimization problem:  
%% ||A(Y-FX)||_2+lambda1||X_1-X_2||_1+lamda2||0.5(X1+X2)-X3||_1+lambda3||WX||_1 
%% where X=[X1 X2 X3], X1,X2 are the low SNR thin slices and X3 is the high SNR thick slice, 
%% A is a weighting matrix that gives more weight to high SNR samples
A=[ones(size(high_res_low_snr_slice1)); ones(size(high_res_low_snr_slice2)); ones(size(low_res_high_snr_slice))];
k_space_thin_slice_1=fft2c(high_res_low_snr_slice1);
k_space_thin_slice_2=fft2c(high_res_low_snr_slice2);
k_space_thick_slice=fft2c(low_res_high_snr_slice);
[res_im_cs, res_no_cs]=FISTA_based_solver(k_space_thin_slice_1,k_space_thin_slice_2,k_space_thick_slice,lambda1,lambda2,lambda3,A,10);

%% Show results for a single thin slice
 figure; imagesc(flipud(high_res_low_snr_slice1'),[0 0.75]); colormap(gray); brighten(0.5); title('Original noisy slice');
 figure; imagesc(flipud(res_im_cs(1:512,:)'),[0 0.75]); colormap(gray); brighten(0.5);title('Results using reference based MRI');
 figure; imagesc(flipud(ground_truth_slice1'),[0 0.75]); colormap(gray); brighten(0.5);title('Ground truth');
    