Spatial Encoding

To create images of the net magnetization, magnetic field gradients are applied in order to modulate the resonance frequency as a function of position. This section describes the effect of magnetic field gradients during data acquisition to create spatial encoding, including frequency encoding, phase encoding and other k-space trajectories.

Learning Goals

1. Describe how images are formed

   • Describe the effects of gradients on the signal

   • Understand how frequency encoding works

   • Understand how phase encoding works

   • Describe what a k-space trajectory is

Effects of Magnetic Field Gradients

Applied magnetic field gradients change the magnetic field linearly with position:

\[\Delta B_z = \vec{G}\cdot\vec{r}\]

which leads to linear variation in the resonance frequency as a function of position:

\[ f = \bar \gamma \ \vec{G}\cdot\vec{r} \]

This forms the basis for spatial encoding.

For constant gradients, \(\vec{G}(t) = \vec{G}_0\), this variation in frequency can be described by (neglecting relaxation)

\[\begin{split}\begin{align} M_{XY}(\vec{r},t) & = M_{XY}(\vec{r},0) \exp(-t/T_2(\vec{r})) \exp\left( -i \gamma \vec{G}_0 \cdot \vec{r} \ t \right) \\ & \approx M_{XY}(\vec{r}, 0) \exp\left( -i \gamma \vec{G}_0 \cdot \vec{r} \ t \right) \end{align}\end{split}\]

For general time-varying gradients, their effect is to create phase accumulation depending on the cumulative gradient area, or the integral of the gradient. This means we can turn the gradients on and off to start and stop the phase accumulation, and also reverse the gradient polarity to undo any prior phase accumulation.

\[M_{XY}(\vec{r}, t) = M_{XY}(\vec{r}, 0) \exp\left( -i \gamma \int_0^t \vec{G}(\tau) \cdot \vec{r} \ d\tau \right)\]

k-space

Here, we introduce the concept of k-space, which is a simplified representation of the phase accumulation due to magnetic field gradients. It is defined as

\[\vec{k}(t) = [k_X(t), k_Y(t), k_Z(t)]^T = \frac{\gamma}{2\pi} \int_0^t \vec{G}(\tau) \cdot \vec{r} \ d\tau\]

The effect of gradients on the transverse magnetization then becomes

\[M_{XY}(\vec{r}, t) = M_{XY}(\vec{r}, 0) e^{ -i 2 \pi \vec{k}(t) \cdot \vec{r} }\]

The specific motivation for k-space will become apparent soon, when we formulate the signal equation and see that there is a Fourier Transform.

MRI signal

The MRI signal comes from the precession of the transverse magnetization, and thus is proportional to the transverse magnetization. The MRI signal also comes from any precessing magnetization within the sensitive volume of the RF receive coils. In other words, it is also not localized to a single location, but rather is summed over a volume:

\[\begin{split}\begin{align} s(t) & = \int_\mathrm{Volume} M_{XY}(\vec{r},t) \ d\vec{r} \\ & = \int_{\textrm{Volume}} M_{XY}(\vec{r},0) \exp(-i2\pi \vec{k}(t) \cdot \vec{r}) \ d\vec{r} \end{align}\end{split}\]

The amazing result here is that this describes our MRI signal in the form of a Fourier Transform:

\[\mathcal{F} \{ f(\vec{r}) \} = F(\vec{k}) = \int_{-\infty}^\infty f(\vec{r}) \exp(-i 2 \pi \vec{k} \cdot \vec{r}) \ d\vec{r}\]

This is the power of the k-space representation, that it describes how MRI is sampling data in the Fourier Transform domain, or the spatial frequency domain, of the object net magnetization. In other words, MRI signals are a measure of the spatial frequencies of our object.

This result means that, to reconstruct an image we need to put our MRI signals into their k-space location based on the applied gradients, and then use an inverse Fourier Transform. Our signal is the Fourier Transform of the initial transverse magnetization, evaluated at a spatial frequnecy, \(\vec{k}\), that is determined by the k-space trajectory, \(\vec{k}(t)\):

\[ s(t) = \mathcal{F} \{ M_{XY}(\vec{r},0) \} |_{\vec{k} = \vec{k}(t)} \]

From MRI data to images

First we use the following notation for simplification

\[m(\vec{r}) = M_{XY}(\vec{r},0)\]

\[\mathcal{F}\{ m(\vec{r}) \} = M(\vec{k})\]

The flow of the experiment and data is as follows:

1. RF excitation to create transverse magnetization, \(M_{XY}(\vec{r},0)\)

2. Gradients applied as \(\vec{G}(t)\) and data is acquired

3. k-space locations, \(\vec{k}(t)\), are determined based on the applied gradients

4. MR signal acquired represents the Fourier Transform of the transverse magnetization at the k-space location:$\( s(t) = M(\vec{k}(t))\)$

5. MR signal over time is stored in a data matrix with known k-space locations to create \(M(\vec{k})\)

6. Inverse Fourier Transform is applied to reconstruct an image of the transverse magnetization$\(\mathcal{F}^{-1}\{ M(\vec{k} )\} = m(\vec{r})\)$

We have now have the incredible result that we used magnetic field gradients, the k-space framework, and appropriate acquisition and ordering of the MR signal to create an IMAGE!!

Cartesian K-space (Frequency and Phase encoding, or FT Imaging)

By far the most common way to perform spatial encoding is using a Cartesian k-space trajectory, which the name refers to the fact that data is sampled on a regularly spaced grid in k-space. This is also known FT imaging, since we reconstruct images with a 2D Fourier Transform.

It consists of 2 types of encoding:

1. Frequency encoding - one dimension is encoded using a constant gradient applied during data acquisition

2. Phase encoding - additional dimensions are encoded using gradient applied before data acquisition. This gradient is incremented to provide complete spatial encoding. Phase encoding is applied in 1 dimension for 2D imaging, and 2 dimensions for 3D imaging

Frequency Encoding

Typically 1 dimension of the object is encoded using “frequency encoding”. This means that, after RF excitation, a magnetic field gradient is turned on and the signal is read out. The frequencies present in the signal correspond to given spatial locations. By convention, this is applied in the x-direction (but in practice can be rotated to any direction);

\[f = \bar \gamma G_{xr} x\]

where \(G_{xr}\) is the readout gradient amplitude. Therefore, the position is proportional to the frequency:

\[ x = \frac{f}{\bar \gamma G_{xr}}\]

From frequency encoding data alone, a 1D image can be reconstructed with an inverse Fourier Transform

% simulate frequency encoding
N = 8;
Mxy= ones(N+1,N+1);

x = [-N/2:N/2];
[x,y] = meshgrid(x,x);
Splot = 0.5;

kFE = 1/2;
dt = 0.1;
Tfe = 2;

GAMMA = 42.58;


for tfe = linspace(0,1,6)*Tfe
    phase_x = 2*pi*kFE*x *tfe/Tfe;
        
    Mxy_FE = Mxy .* exp(i*phase_x);
        
    figure
    subplot(121), plot([0:dt:Tfe+dt],  [0, ones(1,Tfe/dt)*kFE/(Tfe/dt), 0], tfe*ones(1,2), [-0.05 0.05]), ylabel('G_X')
    title('Frequency Encoding')
    subplot(122)
    quiver(x-real(Mxy_FE)*Splot/2,y-imag(Mxy_FE)*Splot/2,real(Mxy_FE), imag(Mxy_FE), Splot)
    xlabel('x'), ylabel('y')
    xlim([-N/2-1, N/2+1])
    ylim([-N/2-1, N/2+1])
    title('M_{XY}')
    drawnow
end

  Cell In[1], line 5
    x = [-N/2:N/2];
             ^
SyntaxError: invalid syntax

Phase Encoding

Typically the 2nd (and optionally 3rd) dimensions of the object are encoded using “phase encoding”. This means that, after RF excitation but before the frequency encoding gradient, a pulsed gradient is applied such that the location is encoded in the phase of the next magnetization:

\[ \Phi(n) = \gamma (-G_{yp} + (n-1) G_{yi} ) t_y y\]

This measurement is repeated for \(n = 1, \ldots, N_{PE}\). \(G_{yp}\) is the maximum phase encoding gradient strength, \(G_{yi}\) is the phase encoding gradient amplitude increment, and \(t_y\) is the phase encoding gradient duration. Note that \(2 G_{yp} = (N_{PE} - 1) G_{yi}\).

These additional dimensions are fully encoded by repeating this pulsed gradient with different amplitudes. This is equivalent to taking different samples of a frequency encoding gradient.

% 2D object
% simulate frequencies + phase encoding

N = 8;
Mxy= ones(N+1,N+1);

x = [-N/2:N/2];
[x,y] = meshgrid(x,x);
Splot = 0.5;

kPE = [-N/2:N/2]/N; %
dt = 0.1;
Tpe = 1;

GAMMA = 42.58;

for Ipe = 1:length(kPE)
    
    for tpe = Tpe %[0:dt:Tpe]
        phase_y = 2*pi*kPE(Ipe)*y *tpe/Tpe;
        
        Mxy_PE = Mxy .* exp(i*phase_y);
        
        figure(Ipe)
        
        subplot(121)
        plot([0:dt:Tpe+dt], [0, ones(1,Tpe/dt)*kPE(Ipe)/(Tpe/dt), 0]), ylim([-0.05 0.05]), ylabel('G_Y') 
        title(['Phase encoding, Step ' int2str(Ipe)])
        subplot(122)
        quiver(x-real(Mxy_PE)*Splot/2,y-imag(Mxy_PE)*Splot/2,real(Mxy_PE), imag(Mxy_PE), Splot)
        xlabel('x'), ylabel('y')
        xlim([-N/2-1, N/2+1])
        ylim([-N/2-1, N/2+1])
        title('M_{XY}')
        drawnow
    end
end

Frequency and Phase Encoding

The following simulation of the net magnetizations shows how first the phase encoding gradient (\(G_Y\)) creates some phase variation in \(y\), and then during the frequency encoding gradient (\(G_X\)) the net magnetizations rotate at varying frequencies depending on their \(x\) position:

Instead of viewing the net magnetizations, we can also visualize this encoding as a map of the phase of the transverse magnetization:

(See spatial_encoding_Mxy_illustration.m for code generating this movie)

K-space Trajectories

K-space is a very general method for capturing the effect of spatial encoding gradients, and the “k-space trajectory” is defined as the pattern created over time by the gradients:

\[\vec{k}(t) = \frac{\gamma}{2\pi} \int_0^t \vec{G}(\tau) d\tau\]

Note that k-space trajectories always start at the center of k-space, \(\vec{k}(0) = 0\) and are only defined once there is transverse magnetization (e.g. after an RF pulse).

The following simulation of the net magnetizations shows how rotations and k-space trajectory during a typical Cartesian (or 2D FT) gradient pulse sequence, which is differs from the simulation above in that an initial dephasing gradient in the frequency encoding direction is applied to sample both positive and negative spatial frequencies in k-space:

(See spatial_encoding_Mxy_illustration.m for code generating this movie)

K-space trajectories

The k-space pattern during a MRI experiment is referred to as the k-space “trajectory”. The most common are Cartesian trajectories, in which parallel lines of k-space are covered to sample a 2D (or 3D) grid. K-space trajectories with other patterns, such as radial lines, spirals, rastered lines (echo-planar trajectories), or blades can also be used.

%  plot 2D k-space trajectories

% Cartesian
N = 8;

k = [-N/2:N/2]/N;
[ky,kx] = meshgrid(k,k);

figure
plot(kx, ky, 'LineWidth',10), xlim([-.6 .6]), ylim([-.6 .6])
xlabel('k_x'),ylabel('k_y')
title('Cartesian trajectory')

% echo-planar
kx_ep = kx; kx_ep(:,2:2:end) = kx_ep(end:-1:1,2:2:end);
kx_ep = kx_ep(:);

ky_ep = ky(:);

figure
plot(kx_ep, ky_ep, 'LineWidth',10), xlim([-.6 .6]), ylim([-.6 .6])
xlabel('k_x'),ylabel('k_y')
title('Echo-Planar trajectory')

% radial

k_theta = exp(i*2*pi*[1:N]/(2*N));
k_radial = k.' * k_theta;

figure
plot(real(k_radial), imag(k_radial), 'LineWidth',10), xlim([-.6 .6]), ylim([-.6 .6])
xlabel('k_x'),ylabel('k_y')
title('Radial trajectory')


% spiral
n = linspace(0,1,201);
Nturns = N/2;
k_spiral = 1/2*n.*exp(i*2*pi*Nturns*n);

figure
plot(real(k_spiral), imag(k_spiral), 'LineWidth',10), xlim([-.6 .6]), ylim([-.6 .6])
xlabel('k_x'),ylabel('k_y')
title('Spiral trajectory')

%  plot 3D k-space trajectories

% Cartesian
N = 8;

k = [-N/2:N/2]/N;
[ky,kx,kz] = meshgrid(k,k,k);


figure
plot3(reshape(kx, [N+1 (N+1)^2]), reshape(ky, [N+1 (N+1)^2]), reshape(kz, [N+1 (N+1)^2]), 'LineWidth',10)
xlim([-.6 .6]), ylim([-.6 .6]), zlim([-.6 .6])
xlabel('k_x'),ylabel('k_y'), zlabel('k_z')
title('Cartesian trajectory')

% radial

N_radial = N^3;

n = 1:N_radial;

kz = (2*n - N_radial - 1) / N_radial;
kx = cos(sqrt(N_radial*pi) .* asin(kz)) .* sqrt(1-kz.^2);
ky = sin(sqrt(N_radial*pi) .* asin(kz)) .* sqrt(1-kz.^2);

figure
plot3([zeros(1,N_radial); kx], [zeros(1,N_radial); ky], [zeros(1,N_radial); kz])
xlim([-.6 .6]), ylim([-.6 .6]), zlim([-.6 .6])
xlabel('k_x'),ylabel('k_y'), zlabel('k_z')
title('Radial trajectory')

These movies illustrate the phase accumulation during non-Cartesian trajectories

(See spatial_encoding_Mxy_illustration.m for code generating this movie)