%% This script calculates the k and f maps using the MT TOF sequence data

load akhil_TOF_MT.mat; %load the data into a 4D array named MTdata
   			     %The dimensions are in readout-phase-slice-exam order
load mask.mat; %load the brain mask created by BET into an array named brainMask

slice = 71; %Pick the slice number

%Mask the MT data, and normalize it with the synthetic image created S0 (mentioned in the paper by Yarnykh)
data_1 = MTdata(:, :, slice, 1).*brainMask(:,:,slice)./S0; 
data_2 = MTdata(:, :, slice, 2).*brainMask(:,:,slice)./S0;
data_3 = MTdata(:, :, slice, 3).*brainMask(:,:,slice)./S0;
data_4 = MTdata(:, :, slice, 4).*brainMask(:,:,slice)./S0;

x0 = [2, 7]'; %this is our initial guess for k and f respectively

% Set certain pulse parameters 
t_m = 8e-3;
t_s = 5e-3;
t_r = 19e-3;
T2_B = 11e-6;
w1rms = 190;
T1 = T1map;%Load the T1map values into T1
R1 = 1./T1;

delta = [3 6 9 12]'*1e3; % offset frequencies for MT pulses

for i = 1:4 %Compute the W_B for each exam through the super-lorentzian line shape
    W_B(i) = lorentzian (delta(i), T2_B);
end;
    
W_F = (w1rms./(2*pi*delta)).^2/.055; % Find W_F from the formula in the paper

options = optimset('Largescale','off','Display','off'); %Set optimization options largescale can also be turned on
recon_MT = zeros(128,128,2); %Set the k and f map matrices to all 0 initially

%For each pixel do the MT fit
for i = 1:128
    i %Print the rows, to get an idea of progress
    for j = 1:128
        if (brainMask(i,j,slice)==1 && T1map(i,j)<=2) % If the data is not zeroed out, and if it is not CSF
        m = [data_1(i, j) data_2(i, j) data_3(i, j) data_4(i,j)]';
        recon_MT (i, j, :) = lsqnonlin (@mt_fit, x0, [], [], options, m, W_B, W_F, delta, R1(i, j), t_m, t_s, t_r); % This function does the fit
        end
    end
end
