%% load ROI mask tissue field map

addpath('utils/')

load ../Phase' Processing'/msk.mat
load ../Phase' Processing'/nfm_final_200iter.mat          % load foreground field map


%% L1 or L2 regularized QSM


N = size(nfm_final);
FOV = N .* [240/256, 240/256, 2.5];  % (in milimeters)
center = N/2 + 1;
figure(1),  imagesc( nfm_final(:,:,32), [-.05,.05] ), axis image off, colormap(gray), drawnow


% regularized QSM with L1 / L2 penalty

param = init;

% Regularization type selection:
% pNorm = 1 for L1, pNorm = 2 for L2 norm penalty
param.pNorm = 1;

param.data = nfm_final;
param.kernel = chi_kernel(FOV,N,center);     
param.padsize = 2*N;
param.TV = TVOP;
param.Itnlim = 10;

% Regularization parameter selection using L-curve:
%   param.TVWeight : 2e-4 for L1, 1.5e-2 for L2 regularization
param.TVWeight = 2e-4;  
%param.TVWeight = 1.5e-2;  
chi_res = zeros(N);

tic
for n = 1:10
    chi_res = fnlCg_susc_padding(chi_res,param); 
end
toc

chi_res = chi_res .* msk;      
for t = 1:size(chi_res,3)
    figure(1),  imagesc( chi_res(:,:,t), [-.2,.2] ), axis image off, colormap(gray), pause(.2)
end


