% Module 2:  relaxation, TR and TE
%
% This script loads phantom.mat (a Matlab variable that contains proton
% density, T1, T2 & shim information for a single slice) and calls FLASH.m
% to simulate image acquisition.
%
%
% Assignment:
%
%
% load phantom
% 1) use FLASH.m to simulate images with the following parameters:
%    Common parameters:  FOV = 64 mm, resolution = 64; alpha = 10; BW = 100kHz
%    a) TE = 5ms; TR = 10ms; TI = 50 ms;
%    b) TE = 5ms; TR = 10ms; TI = 100 ms;
%    c) TE = 5ms; TR = 10ms; TI = 400 ms;
%    d) TE = 5ms; TR = 10ms; TI = 800 ms;
%    e) TE = 5ms; TR = 10ms; TI = 1600 ms;
%    f) TE = 5ms; TR = 10ms; TI = 3200 ms;
%    g) TE = 5ms; TR = 10ms; TI = 4800 ms;
%
%    Collect the final images from each and use these to estimate a T1 map
%    for the sample.
%
% 2) Use FLASH.m to simulate images with the following parameters
%    Common parameters:  FOV = 64 mm, resolution = 64; alpha = 10; BW = 100kHz
%    a) TE = 15ms; TR = 10ms; TI = [];
%    b) TE = 30ms; TR = 10ms; TI = [];
%    c) TE = 60ms; TR = 10ms; TI = [];
%    d) TE = 90ms; TR = 10ms; TI = [];
%
% 3) Using the "data" acquired in problem 2 (or lab), calculate a T2 map
%    for the pineapple slice.
%
% 4) Optional.  Using the data from the inversion recovery experiment,
% calculate a T1 map for the slice.  Note:  the fact that we reconstruct
% magnitude images makes this difficult - an easy way around is to just use
% the data acquired before any of the tissues reach their null point, and
% calculate the map just as in the T2 map.

load phantom

% These are parameters that will be constant throughout
alpha = 10; %degrees
FOV = .064; % ms
res = 64; % ms
BW = 100e3; % Hz (equivalent to 10 microseconds to acquire each point)

% Problem 1
TE = 5; TR = 10; TIs = [50 100 200 400 800];
for n = 1:length(TIs)
    TI = TIs(n);
    figure(1)
    [TIstack(:,:,n),kData,Mss,M0] = FLASH(FOV,32,pd,T1,T2,freqMap,BW,TE,TR,alpha,TI);
    figure(2); colormap(gray)
    subplot(1,length(TIs),n); imagesc(abs(TIstack(:,:,n)),[0 2000]); axis image off
end

% Problem 2 
TI = []; TR = 100; TEs = [15 30 60 90];
for n = 1:length(TEs)
    TE = TEs(n);
	figure(1)
    [TEstack(:,:,n),kData,Mss,M0] = FLASH(FOV,32,pd,T1,T2,freqMap,BW,TE,TR,alpha,TI);
    figure(2); colormap(gray)
    subplot(1,length(TEs),n); imagesc(abs(TEstack(:,:,n)),[0 20]); axis image
end

% Problem 3
w=waitbar(0,'T2 the hard way ...');
for nX = 1:size(TEstack,1)
    waitbar(nX/size(TEstack,1),w)
    for nY = 1:size(TEstack,2)
        p = polyfit(TEs',log(squeeze(abs(TEstack(nX,nY,:)))),1);
        T2est(nX,nY) = -1/p(1);
    end
end; close(w)
figure(3); 
mask = abs(TEstack(:,:,1)) > 5;
imagesc(T2est.*mask,[0 60]); axis image; colorbar; colormap(gray)