
% This file computes the T1 and PD maps of data obtained using the TOF sequence

load akhil_TOF_T1.mat; % Load the data into a 4D array named T1data
   			     %The dimensions are in readout-phase-slice-exam order


slice=71; % Choose a slice
data_1 = T1data(:, :, slice, 1);
data_2 = T1data(:, :, slice, 2);
data_3 = T1data(:, :, slice, 3);
data_4 = T1data(:, :, slice, 4);

theta = [4 10 20 30];
theta = theta*pi/180; % We try four flip angles


% The images are 128x128
% Do the estimation for each pixel
for i = 1:128
    for j = 1:128
        I = [data_1(i, j) data_2(i, j) data_3(i, j) data_4(i,j)];
        y = abs(I./sin(theta));% Do linear regression as in Fram et al. 1987
        x = abs(I./tan(theta));
        if (max(y) > .001) %If the data is not noise
            d = polyfit (x, y, 1); %Find the polynomial fit
            recon_T1(i, j) = d(1); %Take the slope, estimate the exponential term exp(-TR/T1)
        else recon_T1(i,j) = 0;
        end
    end
end

recon_T1 = abs(recon_T1); % There are some negative values for pixels outside the ROI.
T1map = (-log(recon_T1)/20e-3).^-1; %From the exponential term derive the T1 values

%Using the signal equation for SPGR in the absence of MT, pick one of the flip angles and derive the PDmap
theta = 10*pi/180; 
PDmap = data_2./sin(theta).*((1-cos(theta).*recon_T1)./(1-recon_T1));
% Could do least squares for better PDmap.
            
%Generate a synthetic normalization image for the MT fit.
%Using PD and T1, find the signal when flip angle is 10, TR is 32ms, and there is no MT.
S0 = PDmap.*sin(theta).*(1-exp(-32e-3./T1map))./(1-cos(theta).*exp(-32e-3./T1map));

