% ICE4:  Gradients and k-space
%
% This script is designed to simulate and illustrate the imaging process.
% In the first example, the sample is a cube, 192cm on a side, i.e. a
% uniform proton density throughout the field of view.
object = ones(192); % 1mm resolution

% Imaging parameters
FOV = .192; % a 19.2cm FOV will fit most heads, and provide even voxel dimensions
Gro = 15e-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
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 = 64; % resolution in the read-out direction
Npe = 64; % resolution in the phase-encode direction

% First, visualize what is happening to the relative phase of the spins in
% the sample when a gradient is applied. 

% Immediately after excitation, before any gradient is applied, the
% magnetization vectors are all aligned on the positive x-axis (for the
% sake of argument; depending on the phase of the RF pulse, newly excited
% spins can have any orientation in the transverse plane).
phase0 = zeros(size(object)); 

% Gradients are designed so that, in the exact center of the magnet, the
% field is not modified.  So when a single gradient is applied, a line
% through the middle of the sample is not affected; spins to one side
% accumulate a positive relative phase, while spins on the other side
% accumulate a negative relative phase.  For now we are ignoring T2 or T2*
% dephasing, so as long as no gradient is applied, all the spins will
% precess together at the reference frequency, and the relative phase will
% remain zero.  When a gradient is applied, phase (relative to a vector
% precessing at the reference frequency) accumulates according to the
% equation: theta = omega*t.  If omega (frequency, = gamma*B0) is measured
% in radians per second, then angles are measured in radians.  If, however,
% omega is a frequency in Hz, then angles are measured in fractions of a
% cycle.  A very non-standard unit of angle, but certainly convenient, so
% we'll go with gamma = 42.58MHz/T, and measure our phases in fractions of
% a cycle.

% Make a time-series of images showing the accumulated phase (relative to
% precession of reference frequency at isocenter) as a function of time.
% Choose time increments equal to the dwelltime in the acquisition process,
% and calculate 128 steps.
sampleTimes = 0:dwellTime:127*dwellTime;
xcoord = -.0955:.001:.0955; % measured in m
for it = 1:length(sampleTimes)
    % the offset frequency at any point is equal to the gradient (Hz/m) times
    % distance (m)
    phase(:,:,it) = phase0 + mod(0.5 + repmat(Gro*xcoord*sampleTimes(it),[192 1]),1) - 0.5;
end

% The definitions (limitations) of field of view (FOV) and spatial
% resolution (voxel size) are inherent in the physics of the process just
% simulated.  The extent of the FOV is the width of a 2pi phase wrap at the
% first time point.  This is the lowest spatial frequency that can be
% resolved - one cycle per image.  Anything larger than that (spins outside
% the field of view, resonating at a frequency higher than the acquisition
% bandwidth) will be aliased back into the field of view.  Another way of
% thinking of this is to think of the image as a frequency map.  When the
% read-out gradient is applied, spins on the very left side of the image
% are resonating at -BW/2 relative to isocenter; spins on the very right
% side of the image are resonating at BW/2 relative to isocenter.  

% In a similar fashion, the voxel size is determined by how fast phase is
% wrapping when the last data point is collected.  The accumulation of
% phase is determined by the product of the gradient and time. In
% the example above, we stopped after 128 timepoints.  At that point, 128
% 2*pi phase wraps had been superimposed on the image.  Therefore, the
% voxel size of the reconstructed image would have been 192/128 = 1.5mm.
% For the rest of the simulation, we'll stop after 64 time-points (I'm off
% by a factor of 2 here somehow ...)

% From this description it might be obvious why "Hz/pixel" is a common way
% of describing the frequency resolution of an image.  As we'll see below,
% the frequency resolution is rarely the same in both the read-out and
% phase-encode directions.  But in the read-out direction, the full image
% has a bandwidth of 100kHz, so an image with 64 pixels can be described as
% having a pixel bandwidth of 800 Hz/pixel.  This number is useful when
% calculating distortion in images - a frequency offset near a
% susceptibility interface will generate a field distortion, which can be
% converted to Hz using the gyromagnetic ratio.  Then distortions (in
% numbers of pixels) can be predicted from the frequency offset.

% make an informal movie of this for show-and-tell
       for it=1:size(phase,3)
          imagesc(phase(:,:,it)); axis image off
          pause(.1)
       end   

% Now repeat the above process with a pre-defocusing gradient, to show the
% formation of a gradient echo
sampleTimes = 0:dwellTime:159*dwellTime;
for it = 1:32
    phase(:,:,it) = phase0 + repmat(Gro*xcoord*sampleTimes(it),[192 1]);
end
% Now reverse the polarity of the gradient and let it run out to 64 points
for it = 1:64
    phase(:,:,it+32) = phase(:,:,32) - repmat(Gro*xcoord*sampleTimes(it),[192 1]);
end 
% Now take it back again
for it = 1:64
    phase(:,:,it+96) = phase(:,:,96) + repmat(Gro*xcoord*sampleTimes(it),[192 1]);
end

% Show a couple gradient echoes for a uniform proton density
       for it=1:size(phase,3)
          imagesc(mod(0.5+phase(1:2:end,1:2:end,it),1)-0.5,[-.5 .5]); axis image off
          pause(0.1);
       end

% This process could be repeated forever IF it weren't for T2 and T2*
% decay.  We'll add those in later, but first we need to invent a good
% way for showing the net signal in the sample.  The signal is going to be
% the vector sum of all the magnetization vectors, weigted by the local
% proton density.
clf; h = gcf; set(h,'Color',[1 1 1]); colormap(hsv)
object = ones(192); % 1mm resolution
for it = 1:size(phase,3)
    subplot(1,2,1)
    imagesc(object.*mod(0.5+phase(:,:,it),1)-0.5,[-.5 .5]); axis image off
    subplot(1,2,2)
    % each pixel of the image is proton density * phase
    weightedVectors = object.*exp(2*pi*i*phase(:,:,it));
    % adding complex numbers gives us the vector sum
    S(it) = sum(weightedVectors(:));
    plot(sampleTimes(1:it),abs(S(1:it)),'k'); hold on
    plot(sampleTimes(1:it),real(S(1:it)),'r')
    plot(sampleTimes(1:it),imag(S(1:it)),'g'); hold off
    xlabel('time (seconds)'); ylabel('signal'); axis([ 0 .0016 -1e4 4e4])
    M(it) = getframe(h);
end
movie2avi(M,'uniformSample')
% Unfortunately, this first case of uniform signal throughout the imaging
% field of view is informative but trivial.  What we've simulated is
% reading out the Ky=0 line of k-space (no phase encode gradient, one
% read-out line).  It is true that the temporal profile of the read-out
% is the 1D Fourier transform of the object (the FT of the object projected
% onto one dimension).  However, the FT of a uniform signal is non-zero
% only at the center of k-space (Kx=0, Ky=0), the DC component.  So when
% the echo forms, the signal will be non-zero.  Elsewhere, it should be
% zero.

% Make a slightly more interesting object - an ellipse with a major axis
% in the y-direction twice as long as its minor axis in the x-direction.
[X Y] = meshgrid(xcoord,xcoord);
object = sqrt(2*X.^2 + Y.^2) < .08; imagesc(object); axis image off;

% Now repeat the acquisition of the center line of k-space (Ky = 0);
clf; clear M; h = gcf; set(h,'Color',[1 1 1]); 
colormap(hsv); map = colormap; map(33,:) = [0 0 0]; colormap(map)
for it = 1:96
    subplot(1,2,1)
    if (it == 1) | (it == 64)
        imagesc(0.05*object,[-.5 .5]); axis image off
    else
        imagesc(object.*(mod(0.5+phase(:,:,it),1)-0.5),[-.5 .5]); axis image off
    end
    subplot(1,2,2)
    weightedVectors = object.*exp(2*pi*i*phase(:,:,it));
    S(it) = sum(weightedVectors(:));
    plot(sampleTimes(1:it),abs(S(1:it)),'k'); hold on
    plot(sampleTimes(1:it),real(S(1:it)),'r')
    plot(sampleTimes(1:it),imag(S(1:it)),'g'); hold off
    xlabel('time (seconds)'); ylabel('signal'); axis([ 0 .00096 -1e4 2e4])
    title('Center line of k-space, ellipse')
    M(it) = getframe(h);
end
movie2avi(M,'ellipseReadOutK0')

% Just to check that this is indeed FT of the the x-projection of the
% ellipse:
plot(abs(fftshift(ifft(fftshift(S(33:96)))))/512); hold on
plot(mean(object(2:3:end,2:3:end),1),'r'); hold off
%% Looks good except for some scaling issues ... 

% To simulate the entire imaging process, we'll need to generate 64 times
% as many phase images ... which will make one big matrix.  So we'll reduce
% the resolution of our object to match our final imaging resolution.
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; imagesc(object); axis image off;

% First, do the thought experiment of applying simultaneous gradients on
% two axes.
Gpe = Gro;
phase0 = zeros(size(object)); clear phase
for it = 1:32
    phase(:,:,it) = phase0 + Gro*X*sampleTimes(it) + Gpe*Y*sampleTimes(it);
end
for it = 1:32
    phase(:,:,it+32) = phase(:,:,32) - Gro*X*sampleTimes(it) - Gpe*Y*sampleTimes(it);
end
% Now repeat the acquisition of the center line of k-space (Ky = 0);
clf; clear M; h = gcf; set(h,'Color',[1 1 1]); 
colormap(hsv); map = colormap; map(33,:) = [0 0 0]; colormap(map)
for it = 1:64
    subplot(1,2,1)
    if (it == 1) | (it == 64)
        imagesc(0.05*object,[-.5 .5]); axis image off
    else
        imagesc(object.*(mod(0.5+phase(:,:,it),1)-0.5),[-.5 .5]); axis image off
    end
    subplot(1,2,2)
    weightedVectors = object.*exp(2*pi*i*phase(:,:,it));
    S(it) = sum(weightedVectors(:));
    plot(sampleTimes(1:it),abs(S(1:it)),'k'); hold on
    plot(sampleTimes(1:it),real(S(1:it)),'r')
    plot(sampleTimes(1:it),imag(S(1:it)),'g'); hold off
    xlabel('time (seconds)'); ylabel('signal'); axis([ 0 .00064 -1e4 2e4])
    title('Center line of k-space, ellipse')
    M(it) = getframe(h);
end
movie2avi(M,'diagonalGradient')


% Check it out ... we've just projected the image on a diagonal axis.
% Simultaneous gradients simply rotate the coordinate system - we can't get
% something for nothing.

% Do a FLASH image
% First, calculate the effect of hte read pre-dephasing gradient, since
% this is going to be the same for each excitation/phase-encode line.
Nro = 32; Npe = 32;
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; imagesc(object); axis image off;

phase0 = zeros(size(object)); % phase right after excitation
phaseReadDephased = phase0 - Gro*X*Nro/2*dwellTime;
clear phase M S
S = sum(object(:)/2)*ones(size(object));
% Now, step through all 64 phase encode steps
sampleTimes = 0:dwellTime:(Nro-1)*dwellTime;
PEtable = -Npe/2:(Npe/2-1);
w=waitbar(0,'Looping through PEtable ...');
clf; clear M; h = gcf; set(h,'Color',[1 1 1]); 
colormap(hsv); map = colormap; map(33,:) = [0 0 0]; colormap(map)
for iPE = 1:length(PEtable)
    waitbar(iPE/length(PEtable),w)
    phaseDephased = phaseReadDephased + Gpe*PEtable(iPE)*Y*dwellTime;
    for it = 1:Nro
        phase(:,:,iPE,it) =  phaseDephased + Gro*X*sampleTimes(it);
        subplot(1,2,1)
        if (it == (Nro/2+1)) & (iPE == (Npe/2+1))
            imagesc(0.05*object,[-.5 .5]); axis image off
        else
            imagesc(object.*(mod(0.5+phase(:,:,iPE,it),1)-0.5),[-.5 .5]); axis image off
        end
        weighted = object.*exp(2*pi*i*phase(:,:,iPE,it));
        S(iPE,it) = sum(weighted(:));
        subplot(1,2,2)
        imagesc(abs(S),[0 sum(object(:))]); axis image off
        M((iPE-1)*Nro + it) = getframe(h);
    end
end ; close(w)
movie2avi(M,'FLASHimage')
% FLASH images are esentially error-free in the phase-encode direction -
% the exact opposite of EPI images.  This is because every phase-encode
% line has exactly the same phase history - errors accumulate during
% read-out time, but all of these errors are the same in each phase-encode
% step, so no distortions or displacement of signal occur in that
% direction.

%% make movie of FLASH sequence

% If magic happened, then we should recover the image by doing the inverse
% Fourier transform of the S(t):
clf
image = fftshift(ifft2(fftshift(S)));
subplot(2,2,1); imagesc(angle(image)); axis image; title('phase')
subplot(2,2,2); imagesc(abs(image)); axis image; title('magnitude')
subplot(2,2,3); imagesc(imag(image)); axis image; title('imaginary'); colorbar
subplot(2,2,4); imagesc(real(image)); axis image; title('real'); colorbar
colormap(gray)


% The object we're imaging is real-valued (no phase in the proton density
% map; only magnitude).  Therefore the Fourier transform that we generated
% with all those gradients has Hermitian symmetry.  If a point is reflected
% through the origin (Kx=Ky=0), the symmetry point will be the complex
% conjugate of the reflected point (same real portion, same imaginary
% component but with opposite sign).

% Now make an EPI image.
Nro = 16; Npe = 16; 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; 

phase0 = zeros(size(object)); % phase right after excitation
% Apply refocusing (defocusing) gradients in both read and phase
phaseDephased = phase0 - (Gro*X*Nro/2 + Gpe*Y*Npe/2)*dwellTime;
phase = zeros([Nro Npe Npe Nro]);
% Set up figure and colormap for movie-making
clf; clear M; h = gcf; set(h,'Color',[1 1 1]); 
colormap(hsv); map = colormap; map(33,:) = [0 0 0]; colormap(map)

% Now, zig-zag through k-space.  
sampleTimes = 0:dwellTime:Nro*dwellTime;
blip = Gpe*dwellTime; 
readPolarity = ones([1 Npe]); readPolarity(2:2:end) = -1;
S = sum(object(:)/2)*ones(size(object)); % for display purposes
w=waitbar(0,'Looping through PEtable ...');
for iPE = 1:Npe
    waitbar(iPE/length(PEtable),w)
    for it = 1:Nro
        if readPolarity(iPE) > 0
          phase(:,:,iPE,it) = phaseDephased + ...
            readPolarity(iPE)*Gro*X*sampleTimes(it);
          weighted = object.*exp(2*pi*i*phase(:,:,iPE,it));
          S(iPE,it) = sum(weighted(:));
        else
          phase(:,:,iPE,it) = phaseDephased + ...
            readPolarity(iPE)*Gro*X*sampleTimes(it);
          weighted = object.*exp(2*pi*i*phase(:,:,iPE,it));
          % to generate Nyquist ghost, add phase shift to every-other line
          %  weighted = object.*exp(pi/4+2*pi*i*phase(:,:,iPE,it));
          S(iPE,Nro+1-it) = sum(weighted(:));
        end
        subplot(1,2,1)
        if (it == (Nro/2+1)) & (iPE == (Npe/2+1))
            imagesc(0.05*object,[-.5 .5]); axis image off
        else
            imagesc(object.*(mod(0.5+phase(:,:,iPE,it),1)-0.5),[-.5 .5]); axis image off
        end
        subplot(1,2,2)
        imagesc(abs(S),[0 sum(object(:))]); axis image off
        M((iPE-1)*Nro + it) = getframe(h);
    end
    phaseDephased = phase(:,:,iPE,it) + blip*Y;
end    ; close(w) 
%movie2avi(M,'EPIimage')

% If magic happened, then we should recover the image by doing the inverse
% Fourier transform of the S(t):
clf
image = fftshift(ifft2(fftshift(S)));
subplot(2,2,1); imagesc(angle(image)); axis image off; title('phase')
subplot(2,2,2); imagesc(abs(image)); axis image off; title('magnitude')
subplot(2,2,3); imagesc(imag(image)); axis image off; title('imaginary'); colorbar
subplot(2,2,4); imagesc(real(image)); axis image off; title('real'); colorbar
colormap(gray)

% In this case, however, the effective bandwidth in the phase-encode 
% direction is 1/Nro times the bandwidth in the read-out direction.  This
% is because the time between phase-encode points is Nro times longer than
% the time bewteen read-out points - therefore, phase errors have Nro times
% longer to accumulate.  To demonstrate this difference bewteen FLASH and
% EPI read-outs, we're going to introduce a ring of fat around the brain.
% Fat resonates at 3.5ppm relative to water, so if water has a resonant
% frequency of 128MHz, fat will have a resonant frequency of 128.9996 MHz
% Imaging parameters
FOV = .192; % a 19.2cm FOV will fit most heads, and provide even voxel dimensions
Gro = 10e-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 = 32; Npe = 32; 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; 
fat = bwperim(object); imagesc(object + 2*fat); axis image off
object = object + 1.5*fat; % proton density map
%% now make a map of resonant frequencies (in Hz, relative to ref)
offResonance = zeros(size(object)) - fat*3.5e-6*128e6;
%% or just make a global shift in B0
%offResonance = 500*ones(size(object)); 

phase0 = zeros(size(object)); % phase right after excitation
phaseDephased = phase0 - (Gro*X*Nro/2 + Gpe*Y*Npe/2)*dwellTime;
phase = zeros([Nro Npe Npe Nro]);
% Set up figure and colormap for movie-making
clf; clear M; h = gcf; set(h,'Color',[1 1 1]); 
colormap(hsv); map = colormap; map(33,:) = [0 0 0]; colormap(map)

% Make EPI image  
sampleTimes = 0:dwellTime:Nro*dwellTime;
blip = Gpe*dwellTime; 
readPolarity = ones([1 Npe]); readPolarity(2:2:end) = -1;
S = sum(object(:))/2*ones(size(object)); % for display purposes
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(:,:,iPE,it) = phaseDephased + ...
            readPolarity(iPE)*Gro*X*sampleTimes(it);
          % add in the evolution of spins due to off-resonance effects
          evolvedPhase = exp(2*pi*i*offResonance*timeSinceExcitation);
          weighted = object.*evolvedPhase.*exp(2*pi*i*phase(:,:,iPE,it));
          S(iPE,it) = sum(weighted(:));
          timeSinceExcitation = timeSinceExcitation + dwellTime;
        else
          phase(:,:,iPE,it) = phaseDephased + ...
            readPolarity(iPE)*Gro*X*sampleTimes(it);
          evolvedPhase = exp(2*pi*i*offResonance*timeSinceExcitation);
          weighted = object.*evolvedPhase.*exp(2*pi*i*phase(:,:,iPE,it));
          % to generate Nyquist ghost, add phase shift to every-other line
          %  weighted = object.*exp(pi/4+2*pi*i*phase(:,:,iPE,it));
          S(iPE,Nro+1-it) = sum(weighted(:));
          timeSinceExcitation = timeSinceExcitation + dwellTime;
        end
        subplot(1,2,1)
        if (it == (Nro/2+1)) & (iPE == (Npe/2+1))
            imagesc(0.05*object,[-.5 .5]); axis image off
        else
            imagesc(object.*(mod(0.5+phase(:,:,iPE,it),1)-0.5),[-.5 .5]); axis image off
        end
        subplot(1,2,2)
        imagesc(abs(S),[0 sum(object(:))]); axis image off
%        M((iPE-1)*Nro + it) = getframe(h);
    end
    phaseDephased = phase(:,:,iPE,it) + blip*Y;
end    ; close(w) 
%movie2avi(M,'EPIimage')

% If magic happened, then we should recover the image by doing the inverse
% Fourier transform of the S(t):
clf
image = fftshift(ifft2(fftshift(S)));
subplot(1,3,1); imagesc(object); axis image off; title('brain + fat')
subplot(1,3,2); imagesc(offResonance); axis image off; colorbar; title('off resonance map')
subplot(1,3,3); imagesc(abs(image)); axis image off; title('EPI image')
colormap(gray)

% And, last but not least, show that the susceptibility to shift artifacts
% is in the read-out direction for FLASH, since the time elapsed since
% excitation is the same for each column (PE direction) and only varies as
% you go down the row in a FLASH read-out.

% Imaging parameters
FOV = .192; % a 19.2cm FOV will fit most heads, and provide even voxel dimensions
Gro = .5e-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 = 32; Npe = 32; 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; 
fat = bwperim(object); imagesc(object + 2*fat); axis image off
object = object + 1.5*fat; % proton density map
%% now make a map of resonant frequencies (in Hz, relative to ref)
offResonance = zeros(size(object)) - fat*3.5e-6*128e6;
%% or just make a global shift in B0
%offResonance = 500*ones(size(object)); 

phase = zeros([Nro Npe]);
phase0 = zeros(size(object)); % phase right after excitation
phaseReadDephased = phase0 - Gro*X*Nro/2*dwellTime;
S = sum(object(:)/2)*ones(size(object));
% Now, step through all 64 phase encode steps
sampleTimes = 0:dwellTime:(Nro-1)*dwellTime;
PEtable = -Npe/2:(Npe/2-1);
w=waitbar(0,'Looping through PEtable ...');
clf; clear M; h = gcf; set(h,'Color',[1 1 1]); 
colormap(hsv); map = colormap; map(33,:) = [0 0 0]; colormap(map)
for iPE = 1:length(PEtable)
    waitbar(iPE/length(PEtable),w)
    phaseDephased = phaseReadDephased + Gpe*PEtable(iPE)*Y*dwellTime;
    timeSinceExcitation = 0; %% this is the differnce b/w EPI and FLASH
    for it = 1:Nro
        phase =  phaseDephased + Gro*X*sampleTimes(it);
        subplot(1,2,1)
        if (it == (Nro/2+1)) & (iPE == (Npe/2+1))
            imagesc(0.05*object,[-.5 .5]); axis image off
        else
            imagesc(object.*(mod(0.5+phase,1)-0.5),[-.5 .5]); axis image off
        end
        evolvedPhase = exp(2*pi*i*offResonance*timeSinceExcitation);
        weighted = object.*evolvedPhase.*exp(2*pi*i*phase);
        S(iPE,it) = sum(weighted(:));
        timeSinceExcitation = timeSinceExcitation + dwellTime;
        subplot(1,2,2)
        imagesc(abs(S),[0 sum(object(:))]); axis image off
       % M((iPE-1)*Nro + it) = getframe(h);
    end
end ; close(w)
%movie2avi(M,'FLASHimage')

clf
image = fftshift(ifft2(fftshift(S)));
subplot(1,3,1); imagesc(object); axis image off; title('brain + fat')
subplot(1,3,2); imagesc(offResonance); axis image off; colorbar; title('off resonance map')
subplot(1,3,3); imagesc(abs(image)); axis image off; title('FLASH image')
colormap(gray)