% ICE5:  Distortion

% Past simulations used an oval for a head.  But now we'll load some real
% images and some real fieldmaps, and simulate the distortion that
% originates from field perturbations in EPI images.

load filesForUnwarping
whos
%brain = imresize(brain,2);
%fieldMap = imresize(fieldMap,2);

% There's "brain", which was taken with a short TE FLASH sequence, which is
% about as close to "ground truth" as we get.
% There's fieldMap, which was taken with the same sequence, which was
% actually a dual-echo sequence, and this was constructed from the
% difference between the phases of the images of the first and second echo
% times.  In more anterior slices, there's some phase wrapping near that
% auditory canals - that means that the frequency offset down there was so
% large that the spins got more than 2pi out of phase between the 2 echos,
% which were 4.89 and 7.35 ms.  But this slice was selected far enough back
% that phase wrapping wasn't a problem
dTE = (7.35 - 4.89)*.001; % time for phase evolution
phasePolarity = -1;
freqMap = fieldMap/4096/dTE; % in Hz
freqMap = freqMap - mean(freqMap(:));
freqMap = freqMap(:,end:-1:1)';
imagesc(freqMap); axis image off; colormap(gray); colorbar


% Now we'll set up some reasonable imaging parameters
FOV = .192; % a 19.2cm FOV will fit most heads, and provide even voxel dimensions
Gro = 40e-3; % 12 mT/m, or 12 G/cm, would be a moderate gradient strength
% for simplicity in the following calculations, convert the gradient from
% T/m to Hz/m
Gro = Gro*42.58e6; % gyromagnetic ratio for the proton is 42.58 MHz/T
Gpe = Gro; % keep things square for simplicity
BW = FOV*Gro; % sampling rate and bandwidth determine field of view 
dwellTime = 1/BW; % this is the time it takes to read each data point

Nro = size(brain,2); Npe = size(brain,1); xres = FOV/Nro;
xcoord = (-FOV/2 + xres/2):xres:(FOV/2 - xres/2);
[X Y] = meshgrid(xcoord,xcoord); object = sqrt(2*X.^2 + Y.^2) < .08; 
% the object will be the brain
object = brain(:,end:-1:1)'; % rotated to match field map
imagesc(object); axis image off


phase0 = zeros(size(object)); % phase right after excitation
phaseDephased = phase0 - (Gro*X*Nro/2 + phasePolarity*Gpe*Y*Npe/2)*dwellTime;

% Make EPI image  
sampleTimes = 0:dwellTime:Nro*dwellTime;
blip = phasePolarity*Gpe*dwellTime; 
readPolarity = ones([1 Npe]); readPolarity(2:2:end) = -1;
S = zeros(size(object));
timeSinceExcitation = 0;
w=waitbar(0,'Stepping through PE  ...');
for iPE = 1:Npe
    waitbar(iPE/Npe,w)
    for it = 1:Nro
        if readPolarity(iPE) > 0
          phase = phaseDephased + readPolarity(iPE)*Gro*X*sampleTimes(it);
          % add in the evolution of spins due to off-resonance effects
          evolvedPhase = exp(2*pi*i*freqMap*timeSinceExcitation);
          weighted = object.*evolvedPhase.*exp(2*pi*i*phase);
          S(Npe + 1 - iPE,Nro+1-it) = sum(weighted(:));
          timeSinceExcitation = timeSinceExcitation + dwellTime;
        else
          phase = phaseDephased + readPolarity(iPE)*Gro*X*sampleTimes(it);
          evolvedPhase = exp(2*pi*i*freqMap*timeSinceExcitation);
          weighted = object.*evolvedPhase.*exp(2*pi*i*phase);
          % to generate Nyquist ghost, add phase shift to every-other line
          %  weighted = object.*exp(pi/4+2*pi*i*phase(:,:,iPE,it));
          S(Npe + 1 - iPE,it) = sum(weighted(:));
          timeSinceExcitation = timeSinceExcitation + dwellTime;
        end
    end
    phaseDephased = phase + blip*Y;
end    ; close(w) 

% If magic happened, then we should recover the image by doing the inverse
% Fourier transform of the S(t):
clf
if phasePolarity == 1
    image = fftshift(ifft2(fftshift(S)));
else
    image = fftshift(ifft2(fftshift(S(end:-1:1,:))));
end
subplot(1,3,1); imagesc(object); axis image off; title('brain ')
subplot(1,3,2); imagesc(freqMap); axis image off; colorbar; title('off resonance map')
subplot(1,3,3); imagesc(abs(image)); axis image off; title(['T_R_O = ' num2str(1000*timeSinceExcitation,2) 'ms'])
colormap(gray)

