% Module 5:  Distortion.  (And some other ugliness of GE EPI)

% Download new version of EPI code.  No substantial changes, just bugfix.
% Download new pineapple data (sagittalPineapple.mat)

% Problem 1.  T2* and echo time.  Since we've been running field maps, we
% also have T2* measurements from the magnitude images.  So this pineapple
% has, instead of beautiful long T2, realistic short T2*.  This is one of
% the reasons that the pineapple makes such a miserable EPI phantom - we'll
% demonstrate here.
load sagittalPineapple

TR = 1000; TE = 30; alpha = 65; % Typical EPI experiment
FOV = .064; res = 32;
figure; set(gcf,'Name','Long TE, true T2*')
[image kdata] = EPI(FOV,res,pd,T1,T2star,TR,TE,alpha,[],freqMap,[],1);
% zoikes!  There's nothing there.
% Fix it by making TE as short as possible
TE = 3;
figure; set(gcf,'Name','Short TE, true T2*')
[image kdata] = EPI(FOV,res,pd,T1,T2star,TR,TE,alpha,[],freqMap,[],1);
% as you watch the code run between these first three examples, you can see
% that there's just no signal from most of the pineapple with the long echo
% time acquisition.  Short TE recovers a good chunk of the signal.

% So, for the rest of the problem set, we'll artifically inflate T2star to
% keep the signal decay fro messing with our ability to see what else is
% going on.
T2star = 20*T2star;

% Now we'll run a short echo time - watch the evolution of the frequency
% map with time, during the readout.  The lines get wiggly, but are still
% lines.
TE = 3;
figure; set(gcf,'Name','Short TE, long T2*')
[image kdata] = EPI(FOV,res,pd,T1,T2star,TR,TE,alpha,[],freqMap,[],1);
% Now run a long echo time.  BW and resolution haven't changed, so the
% read-out duration is the same, so we're waiting 25ms or so before we
% start acquiring data.  This is the dumbest possbile thing to do, from a
% SNR stand-point.  And you can see the disorganization of the phase
% pattern across the image, which leads to plenty of degradation.
TE = 30;
figure; set(gcf,'Name','Long TE, long T2*')
[image kdata] = EPI(FOV,res,pd,T1,T2star,TR,TE,alpha,[],freqMap,[],1);
% So we'll go back to a moderate TE for the rest of the problem set.
TE = 6;
% We're also going to smooth out the off-resonance map to improve signal
% quality a bit - get rid of this unwanted dephasing that's making the
% pineapple unrecognizable.
freqMap = medfilt2(freqMap);

% Problem 2.  Effect of bandwidth on distortion.  
figure; set(gcf,'Name','Short TE, high BW')
[image kdata] = EPI(FOV,res,pd,T1,T2star,TR,TE,alpha,[],freqMap,[],1);
% If you look carefully, you can see that the top of the pineapple is
% misshapen.  We're going to exaggerate this by increasing the read-out
% time (the time during which errors accumulate in the PE direction), which
% we accomplish by lowering the bandwidth.  Bonus:  higher SNR.
BW = 50000;
figure; set(gcf,'Name','Short TE, low BW')
[image kdata] = EPI(FOV,res,pd,T1,T2star,TR,TE,alpha,[],freqMap,BW,1);

% Problem 3.  Effect of phase polarity on distortion.
% For this problem, you need to edit the EPI code you downloaded.  On line 116,
% there's a variable named phasePolarity that's currently set equal to 1.
% Set this to -1, and then re-run the above example.
figure; set(gcf,'Name','Short TE, high BW, negative phase polarity')
[image kdata] = EPI(FOV,res,pd,T1,T2star,TR,TE,alpha,[],freqMap,BW,1);
% Now, instead of being shifted up, the pineapple is shifted down. The
% phase polarity is the polarity of the phase encode gradient.  When it was
% set to 1, then the gradient was "positive", meaning that spins with
% higher frequency were farther from the isocenter in the positive z
% direction (since this is a sagittal image of the pineapple, the
% slice-select gradient is the x, the read-out gradient is y, and the
% phase-encode gradient is the z.  You can see in the frequency map that
% the field in the bottom of the pineapple (the pineapple's upside down) is
% higher than it should be.  Therefore, with postive phase-encode gradient
% polarity, these spins are displaced away from the isocenter of the
% magnet.  With negative phase-encode polarity, these spins are
% reconstructed at locations closer to the isocenter of the magnet.

% Problem 4.  Calculation of voxel displacement map.
% Quantitative prediction of the displacement of spins is not difficult.
% All you need to know is your phase polarity and your total read-out time
% (time from excitation to the acquisition of the last data point).
%  The total read-out time is determined by the bandwidth and the sampling
%  rate - for an EPI image, it's the dwellTime (1/BW) x Nro x Npe, where
%  Nro and Npe are the number of read-out and phase-encode steps (matrix
%  size).

% To turn in:

% 4A) What is the dwell time of the last simulation (50kHz BW)?
%
% 4B) What is the pixel bandwidth in the PE direction? Remember pixel
% bandwidth can be calculated 2 ways - as the inverse of the total read-out
% time, or as the BW divided by (Npe*Nro).
%
% 4C) Calculate a voxel displacement map (VDM) for the simulated acquisition.
% Starting with the frequency map (the one in memory has been filtered; you
% can use this or re-load the raw data), scale the map (in Hz) by the pixel
% bandwidth, so you have a map of how much each pixel will be displaced in
% the phase encode direction.  Remember to check that you've got the phase
% polarity right.  Turn in a set of images showing the proton density map,
% the simulated image, and the VDM, so a side-by-side comparison can show 
% how well things worked..