% Module 3:  EPI vs. FLASH
%
% Problem 1.  Nyquist ghost.  EPI images are prone to a particular kind of
% ghost - the half field-of-view, or Nyquist.  This is because every other
% line is acquired with a read gradient with the opposite polarity - if
% there's even the slightes imbalance in the gradients or timing, then a
% consistent every-other-line modulation is introduced into the data.  
% Given the Fourier relationship between the data matrix and the image
% matrix, an every-other-line modulation in k-space is equivalent to a half
% FOV modulation in the image.  First we'll try this out with data
% generated for show-and-tell.
%
% Problem 1a) Nyquist ghost from alternate line modulation

[x y] = meshgrid(-31:32,-31:32); % make 2 matrices, with x and y coordinates
r = sqrt(x.^2 + y.^2); % a matrix where each value is the distance from center
circle = r < 24; % if anyone knows an easier way to make a circle, let me know!

Circle = fftshift(fft2(fftshift(circle))); % often we use just capital letters 
% to distinguish between a function and its Fourier transform

subplot(2,2,1); imagesc(circle); axis image off; title('object')
subplot(2,2,2); imagesc(abs(Circle)); axis image off; title('FFT')

% now just attenuate every other line a bit
modCircle = Circle;
modCircle(2:2:end,:) = .5*modCircle(2:2:end,:);

ghostCircle = fftshift(ifft2(fftshift(modCircle)));
subplot(2,2,3); imagesc(abs(ghostCircle)); axis image off; title('object')
subplot(2,2,4); imagesc(abs(modCircle)); axis image off; title('FFT')
% if you're curious, you can look at the real and imaginary parts of the
% transformed circle, and see how they behave
colormap(gray)

% Problem 1b) Nyquist ghost from data shift
% now just attenuate every other line a bit
modCircle = Circle;
modCircle(3:end,2:2:end) = modCircle(1:(end-2),2:2:end);

ghostCircle = fftshift(ifft2(fftshift(modCircle)));
subplot(2,2,3); imagesc(abs(ghostCircle)); axis image off; title('object')
subplot(2,2,4); imagesc(abs(modCircle)); axis image off; title('FFT')
% The shifting added a phase roll to the ghost.
% if you're curious, you can look at the real and imaginary parts of the
% transformed circle, and see how they behave
colormap(jet)

% Problem 2.  Chemical shift artifact.  Because FLASH sequences have no
% systematic differnces between phase encode lines, and any frequency-shift
% errors that accumulate between excitation and echo are the same for every
% phase-encode line, FLASH sequences are for the most part immune to
% susceptibility-induced artifacts in the phase-encode direction.  Any
% sensitivity they have to frequency errors is in the read direction, but
% with a 5ms echo time, it has to be a big error to show up.  EPI images,
% on the other hand, are run with a very high bandwidth, so the time
% between adjacent read-out points is very small (as is the time from start
% to finish on a read-out line). Therefore, the read-out direction is not
% so susceptible to artifacts.  The phase encode direction, however, takes
% a long time to cover - the same as saying that it has a very low
% effective bandwidth.  Therefore, all the big problems in EPI images show
% up in the phase-encode direction.  We're going to simulate a new phantom
% data set with a (fake) fat band around the edge of the pineapple.  Fat
% protons resonate at a frequency 3.5 ppm different from water protons - at
% 3T (127.7MHz) this means that they're 447 Hz off resonance.  The
% following simulations investigate the effects of this chemical shift in
% FLASH and EPI sequences with various bandwidths.

load phantom_fat
FOV = 0.096; TR = 1500; TE = 30; alpha = 70;
% Problem 2a) Bad problem:  high res, low bandwidth
BW = 100000; res = 64;
[res64bw100k kData] = EPI(FOV,res,pd,T1,T2,TR,TE,alpha,[],freqMap,BW);
% Problem 2b) Option 1:  higher bandwidth
BW = 150000; res = 64;
[res64bw150k kData] = EPI(FOV,res,pd,T1,T2,TR,TE,alpha,[],freqMap,BW);
% Problem 2c) Option 2: lower resolution
BW = 100000; res = 32;
[res32bw100k kData] = EPI(FOV,res,pd,T1,T2,TR,TE,alpha,[],freqMap,BW);
% Problem 2d) Option 2: lower resolution
BW = 100000; res = 64; TR = 10; TE = 5; alpha = 7;
[flash64bw100k kData] = FLASH(FOV,res,pd,T1,T2,TR,TE,alpha,[],freqMap,BW);
% In this code, the phase-encode direction is always the vertical dimension
% on the image.  Hard to tell at first glance if the fat is shifted ...
% Problem 2e) 
BW = 10000; res = 64;
[flash64bw50k kData] = FLASH(FOV,res,pd,T1,T2,TR,TE,alpha,[],freqMap,BW);

% Problem 3.  Calculation of magnitude of voxel displacement in FLASH and
% EPI images.  There are various ways to approach this question.  My
% preference is to go back to thinking of the image as a frequency map - my
% bandwidth (sampling rate) determines the range of frequencies represented
% in the image.  With a bandwidth of 100kHz and field of view of 9.6 cm,
% water protons 4.8 cm to the left of isocenter should be resonating at
% -50kHz when the read gradient is turned on, and protons 4.8 cm to the
% right of isocenter should be resonating at +50kHz (relative to the
% resonant frequency at the center of the magnet.  If my matrix size is 64,
% then the center of each voxel is separated by 1.5mm or 1.56kHz.
% Therefore, a chemical shift of 1.56kHz is required before the signal will
% be displaced a voxel in my image.  
%
% Problem 3a:  Verify that the displacement of the "fat" in Problem 2e (3
% voxels) is appropriate for a BW of 10kHz and resolution of 64.
%
% Problem 3b:  EPI is trickier.  Because phase encode lines are read
% sequentially after one excitation, the effective bandwidth in the phase
% encode direction is 1/Npe times the bandwidth in the read-out direction.
% So estimating voxel displacement maps requires dividing the bandwidth by
% the phase resolution before dividing by the read resolution.  A better
% way to think about this is that the time it takes you to get from one
% side of k-space to the other in the phase-encode direction is Npe times
% longer than the time it takes to cover in the read-out direction.  Look
% again at the images in Problem 2 and verify, for each, that the observed
% and predicted shifts match.