%% load ROI mask and unwrapped phase

clear all
addpath('utils/')

load msk
load uphase


%%

uphase = uphase .* msk;

for t = 1:size(msk,3)
    figure(1), imagesc( uphase(:,:,t) ), colormap(gray), axis image off, title(num2str(t)), pause(.1) 
end


%% generate k-space kernel and apply padding 

N = size(uphase);
FOV = N .* [240/256, 240/256, 2.5];  % (in milimeters)

B0 = 1.5;   % Tesla
gyro = 2*pi*42.58;
TE = 40e-3;  % second

nfm = -uphase / (B0 * gyro * TE);   % normalized field map, unitless, in ppm

% pad the image to twice the size and use kernel corresponding to this FOV
nfm_pad = padarray(nfm,N/2);
msk_pad = padarray(msk,N/2);

D_spatial = ifft3c(get_aniso_kernel( FOV, N, N/2+1 ));
D = fft3c(padarray(D_spatial,N/2)) * sqrt(8);     % scale by sqrt(8) due to zero padding
clear D_spatial; clear nfm; clear msk; clear uphase


%% Conjugate gradient dipole fitting

x = zeros(2*N);
num_iter = 200;

tic
x = CG_dipole(nfm_pad, x, D, msk_pad, num_iter);
toc


x_dipole = x .* (1-msk_pad);
nfm_dipole = ifft3c(D .* fft3c(x_dipole));

nfm_final = msk_pad.*(nfm_pad - nfm_dipole);
nfm_final = nfm_final(257-128:256+128, 257-128:256+128, 63-31:62+31);

% nfm_final is the foregound field map to be used with regularized qsm
for t = 1:62
    figure(1), imagesc( nfm_final(:,:,t), [-.1,.1] ), axis image off, colormap(gray), pause(.2)
end





