| SOCR ≫ | BPAD1 Website ≫ | BPAD GitHub ≫ |
Biomedical Physics with Applications to Disease (BPAD)
Chapter 4: Magnetic Resonance Imaging
Magdalena I. Ivanova, Ivo D. Dinov (BPAD Team + SOCR)
August 05, 2026
Chapter 4: Magnetic Resonance Imaging
MRI is the only modality in this book in which the contrast and the position of every voxel are set by two completely independent levers. Relaxation and pulse timing decide what is bright. Gradients and the Fourier transform decide its location. Because those levers are separable, one machine and one population of protons can be re-tasked — by software alone — to answer questions that would otherwise require entirely different instruments.
How to Use This Chapter
Magnetic resonance imaging produces exquisite soft-tissue contrast, in any plane, without ionizing radiation, and it can be tuned to report on anatomy, blood flow, water diffusion, brain function, perfusion, and tissue chemistry. All of it rests on one phenomenon: nuclear spins — overwhelmingly the hydrogen protons of tissue water and fat — placed in a strong magnetic field, nudged with radiofrequency pulses, and listened to as they relax.
The chapter follows the natural sequence of an MRI experiment:
\[\underbrace{\text{align}}_{B_0,\ M_0} \rightarrow \underbrace{\text{excite}}_{B_1,\ \alpha} \rightarrow \underbrace{\text{relax}}_{T_1,\ T_2,\ T_2^*} \rightarrow \underbrace{\text{encode}}_{\text{gradients},\ k\text{-space}} \rightarrow \underbrace{\text{reconstruct}}_{\text{inverse FFT}} \rightarrow \underbrace{\text{measure}}_{\text{parameter maps}} \rightarrow \underbrace{\text{apply}}_{\text{disease}}\]
More than any other modality in this book, MRI is applied Chapter 1 mathematics. The Bloch equations are a system of ordinary differential equations; relaxation is exponential recovery and decay; the precessing magnetization is a rotating complex phasor detected in quadrature; k-space and the image are a two-dimensional Fourier-transform pair; field of view and image artifacts are governed by the Nyquist sampling criterion; magnitude-image noise follows the Rician model; and the diffusion tensor of DTI is analyzed by the eigen-decomposition that Chapter 1 introduced. Table 1 maps each connection to the section that uses it.
Each section closes with a Section summary suitable for a lecture slide and a Checkpoint question suitable for discussion. Three end-of-chapter laboratories build the signature computational skills of MRI: k-space and Fourier reconstruction, artifact simulation, and diffusion-tensor analysis.
Chapter-level clinical puzzle. A patient arrives with sudden weakness — is it an acute stroke, and is there salvageable brain? A neurosurgeon must remove a tumour without cutting the motor tracts or the language cortex. A hepatologist needs to characterize a liver lesion. And a fourth patient with a cardiac implant needs any of these safely. One scanner answers the first three by changing only the pulse timing and the reconstruction; the fourth is answered by a completely different branch of the physics. Which property does each measurement read out, and which one is not about contrast at all?
Spin story for the chapter. Follow one proton. In the bore of the magnet it precesses about \(B_0\) like a tiny gyroscope at the Larmor frequency, and it is one of only about ten in every million that contribute any net signal at all. A resonant RF pulse tips it into the transverse plane, where it induces a faint voltage in a coil. It begins to recover (\(T_1\)) and to fall out of step with its neighbours (\(T_2\)). Gradients stamp its location onto its frequency and phase, so that when millions of such signals are summed, a Fourier transform can sort them back into a picture. Change the timing, and the same proton instead reports how freely its water diffuses, or whether the blood beside it just delivered oxygen to firing neurons.
Key idea. MRI contrast is created by tissue relaxation and chosen by pulse timing (TR, TE, TI, flip angle); spatial position is encoded by gradients into k-space and recovered by the Fourier transform. Contrast and position are independent levers, and that separation is the whole power of the method.
Learning Objectives
| # | Objective | Section |
|---|---|---|
| 1 | Explain why only certain nuclei are MR-active and why \(^1\)H dominates clinical imaging. | 4.1 |
| 2 | Compute the Boltzmann polarization and the equilibrium magnetization \(M_0\), and explain quantitatively why signal grows with field. | 4.1 |
| 3 | Relate magnetic moment, angular momentum, and gyromagnetic ratio, and compute the Larmor frequency. | 4.1 |
| 4 | Use the rotating frame to derive the flip-angle relation \(\alpha = \gamma B_1 \tau\). | 4.1 |
| 5 | Write the full Bloch equations and solve them for \(T_1\) recovery and \(T_2\) decay from arbitrary initial conditions. | 4.2 |
| 6 | Distinguish \(T_1\), \(T_2\), and \(T_2^*\), and state which dephasing a spin echo reverses. | 4.2 |
| 7 | Explain how \(T_1\) and \(T_2\) change with field strength and why 1.5 T contrast recipes must be retuned at 3 T. | 4.2 |
| 8 | Use the spin-echo signal equation to predict T1-, T2-, and PD-weighting from TR and TE. | 4.3 |
| 9 | Compute inversion times for FLAIR and STIR, including the finite-TR correction. | 4.3 |
| 10 | Explain slice selection and frequency/phase encoding, and compute slice thickness and readout bandwidth from gradient amplitudes. | 4.4 |
| 11 | Describe k-space as the Fourier transform of the image and relate \(\Delta k\) and \(k_{\max}\) to FOV and resolution. | 4.4 |
| 12 | Compare spin-echo and gradient-echo sequences and derive the Ernst angle from the steady-state signal. | 4.5 |
| 13 | Relate resolution, FOV, voxel size, bandwidth, and SNR correctly, and explain acceleration methods. | 4.6 |
| 14 | Explain why magnitude MRI noise is Rician and quantify the resulting bias at low SNR. | 4.6 |
| 15 | Identify major MRI artifacts and compute the chemical-shift displacement in pixels. | 4.7 |
| 16 | State the mechanisms and limits of MRI safety: SAR, \(dB/dt\), projectiles, acoustic noise, implants, and contrast agents. | 4.8 |
| 17 | Explain DWI/ADC, DTI/FA, BOLD fMRI, perfusion, and spectroscopy, and connect each to disease. | 4.10 |
| 18 | Quantify how the Rician noise floor biases ADC at high b-value. | 4.10 |
| 19 | Implement Fourier reconstruction, artifact simulation, and diffusion-tensor FA in R. |
4.14 |
Notation
| Symbol | Meaning | Units |
|---|---|---|
| \(B_0,\ B_1\) | static field; RF (transmit) field | T |
| \(\gamma,\ \bar\gamma = \gamma/2\pi\) | gyromagnetic ratio | rad s\(^{-1}\)T\(^{-1}\); Hz T\(^{-1}\) |
| \(\omega_0,\ f_0\) | Larmor angular frequency; frequency | rad s\(^{-1}\); Hz |
| \(\vec\mu,\ \vec L,\ I\) | magnetic moment; angular momentum; nuclear spin | J T\(^{-1}\); J s; — |
| \(\vec M,\ M_0,\ M_z,\ M_{xy}\) | magnetization vector; equilibrium, longitudinal, transverse | A m\(^{-1}\) |
| \(\alpha,\ \alpha_E\) | flip angle; Ernst angle | rad or degrees |
| \(\tau_p\) | RF pulse duration | s |
| \(T_1,\ T_2,\ T_2',\ T_2^*\) | relaxation times | ms |
| \(\rho\) | proton (spin) density | — |
| TR, TE, TI | repetition, echo, inversion time | ms |
| \(G_x, G_y, G_z\) | gradient amplitudes | T m\(^{-1}\) (mT m\(^{-1}\)) |
| \(k_x, k_y,\ \Delta k,\ k_{\max}\) | spatial frequency; sampling step; extent | m\(^{-1}\) |
| FOV, \(\Delta x\), \(N\) | field of view; pixel size; matrix size | m, m, — |
| BW, BW\(_{\text{px}}\) | readout bandwidth; bandwidth per pixel | Hz |
| NSA | number of signal averages | — |
| \(R,\ g\) | parallel-imaging acceleration; g-factor | — |
| \(b\), ADC | diffusion weighting; apparent diffusion coefficient | s mm\(^{-2}\); mm\(^2\)s\(^{-1}\) |
| \(\mathbf D,\ \lambda_i\) | diffusion tensor; eigenvalues | mm\(^2\)s\(^{-1}\) |
| MD, FA | mean diffusivity; fractional anisotropy | mm\(^2\)s\(^{-1}\); — |
| SAR | specific absorption rate | W kg\(^{-1}\) |
| \(\sigma\) | noise standard deviation per channel | — |
Bridges from Earlier Chapters
| Earlier result | Chapter | Role here | Section |
|---|---|---|---|
| Systems of ODEs | 1 | the Bloch equations | 4.2 |
| Exponential growth and decay | 1, 2, 3 | \(T_1\) recovery, \(T_2\) decay, diffusion attenuation | 4.2, 4.10 |
| Complex exponentials and phasors | 1 | precessing \(M_{xy}\); quadrature detection | 4.1, 4.2 |
| Rotating coordinate frames | 1 | flip angle from \(B_1\) | 4.1 |
| Boltzmann factor and statistical populations | 1, 2 | polarization and the size of \(M_0\) | 4.1 |
| Two-dimensional Fourier transform | 1 | k-space and reconstruction | 4.4, Lab 1 |
| Sampling, Nyquist, aliasing | 1, 3 | FOV, wraparound, and undersampling | 4.4, 4.6 |
| Windowing and convolution | 1 | Gibbs ringing from truncated k-space | 4.7, Lab 2 |
| Rician and Gaussian noise models | 1, 8 | magnitude-image noise and ADC bias | 4.6, 4.10 |
| Eigenvalues and eigenvectors | 1 | diffusion tensor, FA, MD, tractography | 4.10, Lab 3 |
| Least squares and conditioning | 1 | tensor fitting; parameter mapping | 4.10 |
| Optical absorption of haemoglobin | 2 | BOLD and fNIRS measure the same haemodynamics | 4.10 |
| Fourier pair: decay \(\leftrightarrow\) spectrum | 2, 3 | FID \(\rightarrow\) MRS spectrum, as FTIR interferogram \(\rightarrow\) IR spectrum | 4.2, 4.10 |
4.1 Spin Physics: The Origin of the MR Signal
4.1.1 Spin-active nuclei
Only nuclei with non-zero nuclear spin (\(I \neq 0\)) possess a magnetic moment and can participate in magnetic resonance. Hydrogen (\(^1\)H, a single proton, \(I = \tfrac12\)) is the workhorse of clinical MRI because it is enormously abundant — water and fat fill the body — and because it has the largest gyromagnetic ratio of the common nuclei, giving the strongest signal per nucleus. Other MR-active nuclei such as \(^{13}\)C, \(^{19}\)F, \(^{23}\)Na, and \(^{31}\)P are used in research and spectroscopy, while nuclei with paired spins (zero net moment) are MR-silent.
4.1.2 Energy splitting, polarization, and the size of the signal
A nucleus with spin behaves like a tiny magnet whose magnetic moment \(\vec\mu\) is proportional to its angular momentum \(\vec L\):
\[\begin{equation} \vec\mu = \gamma\,\vec L . \tag{1} \end{equation}\]
In a field \(B_0\) the spin energy is \(E = -\vec\mu\cdot\vec B\), so the two states of a spin-\(\tfrac12\) nucleus are split by
\[\begin{equation} \Delta E = \gamma\hbar B_0 = \hbar\omega_0 . \tag{2} \end{equation}\]
At thermal equilibrium the lower-energy state is more populated, and the Boltzmann factor of Chapter 1 gives the fractional excess directly. Because \(\Delta E \ll k_B T\) by an enormous margin, the exponential can be expanded to first order:
\[\begin{equation} \frac{N_\uparrow - N_\downarrow}{N_\uparrow + N_\downarrow} = \tanh\!\left(\frac{\gamma\hbar B_0}{2k_BT}\right) \approx \frac{\gamma\hbar B_0}{2k_BT} . \tag{3} \end{equation}\]
Summing the moments of that tiny excess over a spin density \(N\) gives the equilibrium magnetization
\[\begin{equation} M_0 = \frac{N\gamma^{2}\hbar^{2}B_0}{4k_BT} \qquad\left(\text{spin-}\tfrac12;\ \text{generally } \frac{N\gamma^2\hbar^2 I(I+1)B_0}{3k_BT}\right). \tag{4} \end{equation}\]
Equation (4) is the quantitative reason higher-field scanners give more signal: \(M_0 \propto B_0\), so doubling the field doubles the available magnetization.
B0 <- seq(0.1, 7, by = 0.02)
pol <- gamma_H*hbar*B0/(2*k_B*T_body)
p_pol <- ggplot(data.frame(B0 = B0, ppm = pol*1e6), aes(B0, ppm)) +
geom_line(linewidth = 0.95, colour = bpad_pal[1]) +
geom_point(data = data.frame(B0 = c(1.5, 3, 7),
ppm = gamma_H*hbar*c(1.5,3,7)/(2*k_B*T_body)*1e6),
size = 2.4, colour = bpad_pal[2]) +
geom_text(data = data.frame(B0 = c(1.5, 3, 7),
ppm = gamma_H*hbar*c(1.5,3,7)/(2*k_B*T_body)*1e6),
aes(label = sprintf("%.1f ppm", ppm)), vjust = -0.8, size = 3) +
labs(x = expression("Field strength "*B[0]*" (T)"),
y = "Polarization (parts per million)",
title = "Boltzmann polarization at body temperature")
M0rel <- B0/1.5
f0 <- gbar_H*B0/1e6
p_m0 <- ggplot(data.frame(B0 = B0, M0 = M0rel, f0 = f0), aes(B0)) +
geom_line(aes(y = M0, colour = "M0 (relative to 1.5 T)"), linewidth = 0.95) +
geom_line(aes(y = f0/63.87, colour = "Larmor frequency (relative)"),
linewidth = 0.95, linetype = "dashed") +
scale_colour_manual(values = bpad_pal[c(1,2)]) +
labs(x = expression(B[0]*" (T)"), y = "Relative to 1.5 T",
title = "Both magnetization and frequency scale with B0")
bpad_grid(p_pol, p_m0, ncol = 2)
Figure 1: How little of the spin population MRI actually uses. Left: the Boltzmann polarization, Eq. (polarization), is of order ten parts per million at clinical field strengths, so roughly one proton in 100,000 contributes net signal at 3 T; MRI works only because a voxel contains of order 10 to the 21 protons. Right: equilibrium magnetization and Larmor frequency both rise linearly with field, Eqs. (m0) and (larmor), which is the physical basis of the move from 1.5 T to 3 T and beyond.
knitr::kable(data.frame(
B0_T = c(0.55, 1.5, 3, 7),
polarization_ppm = round(gamma_H*hbar*c(0.55,1.5,3,7)/(2*k_B*T_body)*1e6, 2),
excess_per_million = round(gamma_H*hbar*c(0.55,1.5,3,7)/(2*k_B*T_body)*1e6, 1),
f0_MHz = round(gbar_H*c(0.55,1.5,3,7)/1e6, 1),
M0_rel_to_1p5T = round(c(0.55,1.5,3,7)/1.5, 2)),
col.names = c("B0 (T)","Polarization (ppm)","Excess spins per million",
"1H frequency (MHz)","M0 relative to 1.5 T"),
caption = "Polarization, resonance frequency, and equilibrium magnetization across clinical field strengths. Even at 7 T fewer than 25 protons per million contribute.")| B0 (T) | Polarization (ppm) | Excess spins per million | 1H frequency (MHz) | M0 relative to 1.5 T |
|---|---|---|---|---|
| 0.55 | 1.81 | 1.8 | 23.4 | 0.37 |
| 1.50 | 4.94 | 4.9 | 63.9 | 1.00 |
| 3.00 | 9.89 | 9.9 | 127.7 | 2.00 |
| 7.00 | 23.07 | 23.1 | 298.0 | 4.67 |
The number that should astonish you. At 3 T only about 10 protons in every million are in excess in the lower-energy state. MRI images that tiny residue. It is viable only because a \(1\times1\times3\) mm voxel of water contains of order \(10^{20}\) protons, so even a \(10^{-5}\) fractional excess leaves \(10^{15}\) spins contributing coherently. This is also why MRI is a low-sensitivity technique in absolute terms compared with the nuclear methods of Chapter 6, where every decaying nucleus is counted.
Common misconception. “Spins flip up or down like little bar magnets.”
Correction. The classical two-state picture is a useful shorthand for the energy argument of Eq. (2), but it misleads about dynamics. Almost all of MRI is described by the behaviour of a continuous net magnetization vector \(\vec M\) — the statistical sum over an enormous spin ensemble — which precesses and relaxes smoothly and can point in any direction, including entirely within the transverse plane. A single spin never “points” transversely in the classical sense; the ensemble does.
4.1.3 Torque, precession, and the Larmor frequency
A magnetic moment tilted relative to \(B_0\) feels a torque \(\vec\tau = \vec\mu\times\vec B_0\) that changes its direction but not its magnitude, so the magnetization precesses about \(B_0\) like a spinning top:
\[\begin{equation} \frac{d\vec M}{dt} = \gamma\,\vec M\times\vec B . \tag{5} \end{equation}\]
For \(\vec B = B_0\hat z\) the transverse components obey \(\dot M_x = \gamma B_0 M_y\), \(\dot M_y = -\gamma B_0 M_x\), a linear ODE system whose solution rotates at the Larmor frequency
\[\begin{equation} \omega_0 = \gamma B_0, \qquad f_0 = \frac{\gamma}{2\pi}B_0 = \bar\gamma B_0 . \tag{6} \end{equation}\]
Connection to Chapter 1. Writing \(M_{xy} = M_x + iM_y\) collapses the two real equations into one complex equation \(\dot M_{xy} = -i\omega_0 M_{xy}\), whose solution \(M_{xy}(t) = M_\perp e^{-i\omega_0 t}\) is exactly the rotating phasor of the Chapter 1 treatment of complex numbers. MRI receivers detect this complex signal in quadrature (two channels \(90^\circ\) apart), recovering both magnitude and phase — which is why phase can be used to encode position (Section 4.4), velocity (phase-contrast angiography), and susceptibility.
gamma_bar <- c("1H" = 42.577, "19F" = 40.078, "31P" = 17.235,
"23Na" = 11.262, "13C" = 10.708)
B0g <- seq(0, 7.5, by = 0.05)
clin <- data.frame(B0 = c(0.55, 1.5, 3, 7),
f0 = gamma_bar["1H"]*c(0.55, 1.5, 3, 7))
p_f <- ggplot(data.frame(B0 = B0g, f0 = gamma_bar["1H"]*B0g), aes(B0, f0)) +
geom_line(colour = bpad_pal[1], linewidth = 0.95) +
geom_point(data = clin, colour = bpad_pal[2], size = 2.4) +
geom_text(data = clin, aes(label = sprintf("%.2g T\n%.0f MHz", B0, f0)),
vjust = -0.35, hjust = 1.1, size = 2.7) +
labs(x = expression("Field strength "*B[0]*" (T)"),
y = expression(""^1*"H Larmor frequency (MHz)"),
title = "Proton resonance versus field strength")
df_g <- data.frame(nucleus = names(gamma_bar), gamma = as.numeric(gamma_bar))
df_g$nucleus <- factor(df_g$nucleus, levels = df_g$nucleus[order(-df_g$gamma)])
p_g <- ggplot(df_g, aes(nucleus, gamma)) +
geom_col(fill = bpad_pal[1]) +
geom_text(aes(label = sprintf("%.2f", gamma)), vjust = -0.35, size = 3) +
coord_cartesian(ylim = c(0, 50)) +
labs(x = NULL, y = expression(gamma/2*pi*" (MHz/T)"),
title = "Gyromagnetic ratios")
bpad_grid(p_f, p_g, ncol = 2)
Figure 2: Larmor frequencies. Left: the proton resonance frequency rises linearly with field strength, Eq. (larmor); labelled points mark clinical field strengths, and 3 T lands at 128 MHz, squarely in the FM broadcast band, which is why scanners sit inside a radiofrequency-shielded room. Right: gyromagnetic ratios of MR-active nuclei; 1H has both the largest ratio and by far the greatest abundance, which is why it dominates clinical imaging.
4.1.4 The rotating frame and the flip angle
At equilibrium the magnetization is purely longitudinal (\(M_z = M_0\), \(M_{xy} = 0\)) and induces no signal. To create signal we must tip it, and the cleanest way to see how is to change coordinates.
Transform to a frame rotating about \(\hat z\) at the RF carrier frequency \(\omega_{\rm RF}\). In that rotating frame the apparent static field is reduced to the off-resonance residual \(\Delta B = B_0 - \omega_{\rm RF}/\gamma\), so on resonance the huge \(B_0\) field disappears entirely and only the small applied RF field \(B_1\) remains:
\[\begin{equation} \vec B_{\rm eff} = \left(B_0 - \frac{\omega_{\rm RF}}{\gamma}\right)\hat z + B_1\hat x' \;\xrightarrow{\ \omega_{\rm RF}=\omega_0\ }\; B_1\hat x' . \tag{7} \end{equation}\]
The magnetization then simply precesses about \(B_1\) at rate \(\gamma B_1\), so a pulse of duration \(\tau_p\) tips it through the flip angle
\[\begin{equation} \alpha = \gamma\, B_1\, \tau_p \qquad\text{(for a rectangular pulse)}, \qquad \alpha = \gamma\!\int_0^{\tau_p}\! B_1(t)\,dt \ \ \text{(in general)} . \tag{8} \end{equation}\]
Resonance matters because only a field oscillating at \(\omega_0\) appears static in the rotating frame; an off-resonance field spirals and produces no net tip. This is the magnetic analogue of pushing a swing in time with its natural rhythm.
B1_uT <- 10
alpha_deg <- seq(0, 180, by = 1)
tau_ms <- (alpha_deg*pi/180)/(gamma_H*B1_uT*1e-6)*1e3
p_tau <- ggplot(data.frame(a = alpha_deg, tau = tau_ms), aes(a, tau)) +
geom_line(linewidth = 0.95, colour = bpad_pal[1]) +
geom_point(data = data.frame(a = c(90, 180),
tau = (c(90,180)*pi/180)/(gamma_H*B1_uT*1e-6)*1e3),
size = 2.4, colour = bpad_pal[2]) +
geom_text(data = data.frame(a = c(90, 180),
tau = (c(90,180)*pi/180)/(gamma_H*B1_uT*1e-6)*1e3),
aes(label = sprintf("%.0f deg: %.2f ms", a, tau)),
hjust = 1.05, vjust = -0.6, size = 3) +
labs(x = "Flip angle (degrees)", y = "Pulse duration (ms)",
title = expression("Pulse duration for "*B[1]*" = 10 "*mu*"T"))
df_a <- rbind(
data.frame(a = alpha_deg, v = sin(alpha_deg*pi/180), q = "transverse ~ sin(alpha)"),
data.frame(a = alpha_deg, v = cos(alpha_deg*pi/180), q = "longitudinal ~ cos(alpha)"))
p_sc <- ggplot(df_a, aes(a, v, colour = q)) +
geom_hline(yintercept = 0, colour = "grey70") +
geom_line(linewidth = 0.95) +
geom_vline(xintercept = 90, linetype = "dotted", colour = "grey45") +
scale_colour_manual(values = bpad_pal[c(1,2)]) +
labs(x = "Flip angle (degrees)", y = "Relative magnetization",
title = "What a flip angle buys and costs")
bpad_grid(p_tau, p_sc, ncol = 2)
Figure 3: The flip angle in the rotating frame, Eq. (flip-angle). Left: with a 10 microtesla B1 field, the pulse duration required grows linearly with the desired flip angle; a 90-degree pulse takes about 0.6 ms and a 180-degree pulse twice that. Right: transverse signal is proportional to sin(alpha) while residual longitudinal magnetization goes as cos(alpha), so a 90-degree pulse maximizes immediate signal but leaves nothing to recover, which is exactly the trade the Ernst angle of Section 4.5 resolves.
cat(sprintf("With B1 = %d uT: 90-degree pulse takes %.3f ms, 180-degree takes %.3f ms.\n",
B1_uT, (pi/2)/(gamma_H*B1_uT*1e-6)*1e3, pi/(gamma_H*B1_uT*1e-6)*1e3))## With B1 = 10 uT: 90-degree pulse takes 0.587 ms, 180-degree takes 1.174 ms.cat(sprintf("B1 is ~%.0e times weaker than B0 at 3 T, yet it does all the work,\n",
3/(B1_uT*1e-6)))## B1 is ~3e+05 times weaker than B0 at 3 T, yet it does all the work,## because in the rotating frame on resonance B0 has vanished.Worked Example 4.1 (a resonance frequency and a pulse). For \(^1\)H at 3 T, \(f_0 = 42.577 \times 3 = 127.7\) MHz. A rectangular \(B_1 = 10\) µT pulse produces \(\alpha = \gamma B_1\tau_p\), so a \(90^\circ\) tip (\(\alpha = \pi/2\)) needs
\[\tau_p = \frac{\pi/2}{\gamma B_1} = \frac{1.5708}{2.675\times10^{8}\times10^{-5}} = 5.87\times10^{-4}\ \text{s} = 0.59\ \text{ms}.\]
Note the scale separation: \(B_1\) is about 300,000 times weaker than \(B_0\), yet it tips the magnetization completely. That is only possible because, on resonance in the rotating frame, \(B_0\) has been transformed away and \(B_1\) is the only field acting.
Section 4.1 summary.
- MR-active nuclei have non-zero spin; \(^1\)H dominates through abundance and the largest \(\gamma\).
- The Zeeman splitting is \(\Delta E = \gamma\hbar B_0\); the Boltzmann polarization is only about 10 ppm at 3 T, and \(M_0 \propto N\gamma^2\hbar^2B_0/4k_BT\).
- Magnetization precesses at \(\omega_0 = \gamma B_0\); in complex form it is a rotating phasor detected in quadrature.
- In the rotating frame on resonance, \(B_0\) vanishes and \(\alpha = \gamma B_1\tau_p\).
4.2 Relaxation and Tissue Contrast
4.2.1 The full Bloch equations
Precession alone would leave the magnetization spinning forever. Bloch added phenomeno- logical relaxation terms — one restoring \(M_z\) toward \(M_0\), one destroying \(M_x\) and \(M_y\) — giving the complete system that governs essentially all of MRI
\[\begin{equation} \frac{d\vec M}{dt} = \underbrace{\gamma\,\vec M\times\vec B}_{\text{precession}} - \underbrace{\frac{M_x\hat x + M_y\hat y}{T_2}}_{\text{transverse decay}} + \underbrace{\frac{(M_0 - M_z)\hat z}{T_1}}_{\text{longitudinal recovery}} . \tag{9} \end{equation}\]
Component by component, in the laboratory frame with \(\vec B = B_0\hat z\):
\[\begin{equation} \dot M_x = \gamma B_0 M_y - \frac{M_x}{T_2},\qquad \dot M_y = -\gamma B_0 M_x - \frac{M_y}{T_2},\qquad \dot M_z = \frac{M_0 - M_z}{T_1} . \tag{10} \end{equation}\]
This is a linear ODE system of exactly the type solved in Chapter 1. Two features are worth noticing before solving it. The longitudinal and transverse equations are decoupled in the absence of RF, so \(T_1\) and \(T_2\) can be measured independently. And the transverse pair is the real form of a single complex equation, which is why the solution below is a damped phasor.
4.2.2 \(T_1\): longitudinal recovery
The \(M_z\) equation in (10) is a first-order linear ODE with the general solution
\[\begin{equation} M_z(t) = M_0 + \big[M_z(0) - M_0\big]e^{-t/T_1} . \tag{11} \end{equation}\]
Equation (11) is the form to remember, because the three sequences that matter differ only in \(M_z(0)\):
| Preparation | \(M_z(0)\) | Recovery |
|---|---|---|
| Saturation (\(90^\circ\) pulse) | \(0\) | \(M_z(t) = M_0(1 - e^{-t/T_1})\) |
| Inversion (\(180^\circ\) pulse) | \(-M_0\) | \(M_z(t) = M_0(1 - 2e^{-t/T_1})\) |
| Partial (\(\alpha\) pulse) | \(M_0\cos\alpha\) | \(M_z(t) = M_0[1 - (1-\cos\alpha)e^{-t/T_1}]\) |
\(T_1\) (spin-lattice relaxation) is the time constant for spins to return energy to their molecular surroundings. It is long for free water (CSF) and short for fat and protein-rich tissue, because efficient relaxation requires molecular tumbling at rates near the Larmor frequency — and free water tumbles far too fast.
4.2.3 \(T_2\): transverse dephasing
Meanwhile the transverse magnetization decays as spins lose phase coherence through random spin–spin interactions:
\[\begin{equation} M_{xy}(t) = M_{xy}(0)\, e^{-t/T_2}, \qquad M_{xy}(t) = M_{xy}(0)\,e^{-t/T_2}e^{-i\omega_0 t}\ \text{(complex form)} . \tag{12} \end{equation}\]
In the complex form of Eq. (12) the phasor of Eq. (6) reappears under an exponential envelope. \(T_2\) (spin–spin relaxation) is always \(\le T_1\), because any process that returns longitudinal magnetization also destroys transverse coherence, while the converse is not true. Crucially, the spins are not gone: they have fanned out in phase so that their vector sum shrinks.
Connection to Chapter 1. Equations (11)–(12) are the same exponential growth-and-decay processes met in Chapter 1, as Beer–Lambert attenuation in Chapter 2, and as acoustic attenuation in Chapter 3. The complex form of Eq. (12) is a damped phasor, and it is literally the free-induction decay plotted below.
4.2.4 \(T_2^*\) and field inhomogeneity
The observed transverse decay is faster than true \(T_2\), because static inhomogeneities in \(B_0\) make spins at different locations precess at slightly different rates. Writing the inhomogeneous contribution as \(T_2'\),
\[\begin{equation} \frac{1}{T_2^{*}} = \frac{1}{T_2} + \frac{1}{T_2'}, \qquad \frac{1}{T_2'} \approx \gamma\,\Delta B_0 , \qquad T_2^{*} \le T_2 \le T_1 . \tag{13} \end{equation}\]
Shimming improves field uniformity and lengthens \(T_2^*\). The vital distinction, developed in Section 4.5, is that the static part of this dephasing can be reversed by a \(180^\circ\) refocusing pulse, recovering true \(T_2\), whereas the random spin–spin part cannot.
4.2.5 Relaxation depends on field strength
Relaxation times are not material constants: \(T_1\) rises substantially with field while \(T_2\) is roughly flat or falls slightly. This is not a footnote. A T1-weighted protocol optimized at 1.5 T produces visibly different contrast at 3 T unless TR is lengthened, and the FLAIR inversion time must be re-derived entirely.
t1g <- seq(0, 5000, by = 10)
df_t1 <- do.call(rbind, lapply(seq_len(nrow(tissues)), function(i)
data.frame(t = t1g, M = 1 - exp(-t1g/tissues$T1[i]), tissue = tissues$tissue[i])))
df_t1$tissue <- factor(df_t1$tissue, levels = tissues$tissue)
p1 <- ggplot(df_t1, aes(t, M, colour = tissue)) +
geom_line(linewidth = 0.9) +
geom_hline(yintercept = 1 - exp(-1), linetype = "dotted", colour = "grey55") +
scale_colour_manual(values = cols) +
labs(x = "Time after 90-degree pulse (ms)", y = expression(M[z]/M[0]),
title = "T1 recovery (longitudinal)")
t2g <- seq(0, 400, by = 1)
df_t2 <- do.call(rbind, lapply(seq_len(nrow(tissues)), function(i)
data.frame(t = t2g, M = exp(-t2g/tissues$T2[i]), tissue = tissues$tissue[i])))
df_t2$tissue <- factor(df_t2$tissue, levels = tissues$tissue)
p2 <- ggplot(df_t2, aes(t, M, colour = tissue)) +
geom_line(linewidth = 0.9) +
geom_hline(yintercept = exp(-1), linetype = "dotted", colour = "grey55") +
scale_colour_manual(values = cols) +
labs(x = "Time after 90-degree pulse (ms)", y = expression(M[xy]/M[xy](0)),
title = "T2 decay (transverse)")
fld <- rbind(
data.frame(tissue = relax$tissue, ratio = relax$T1_3/relax$T1_1p5, q = "T1 (3T / 1.5T)"),
data.frame(tissue = relax$tissue, ratio = relax$T2_3/relax$T2_1p5, q = "T2 (3T / 1.5T)"))
p3 <- ggplot(fld, aes(reorder(tissue, ratio), ratio, fill = q)) +
geom_hline(yintercept = 1, colour = "grey55") +
geom_col(position = "dodge") + coord_flip() +
scale_fill_manual(values = bpad_pal[c(1,2)]) +
labs(x = NULL, y = "Ratio 3 T / 1.5 T", title = "Field dependence of relaxation")
bpad_grid(p1, p2, p3, ncol = 3)
Figure 4: Relaxation generates contrast, and the contrast depends on field. Left: T1 recovery, Eq. (t1-general), fast for fat and very slow for CSF. Centre: T2 decay, Eq. (t2-decay); tissues separate at intermediate echo times, which is exactly where contrast is created. Right: T1 rises by 6 to 60 percent from 1.5 T to 3 T while T2 barely moves, which is why 1.5 T timing recipes must be retuned rather than transplanted.
knitr::kable(relax, col.names = c("Tissue","T1 at 1.5 T (ms)","T1 at 3 T (ms)",
"T2 at 1.5 T (ms)","T2 at 3 T (ms)"),
caption = "Representative relaxation times. Values vary between studies and with tissue state; they are used here for teaching, not for protocol design.")| Tissue | T1 at 1.5 T (ms) | T1 at 3 T (ms) | T2 at 1.5 T (ms) | T2 at 3 T (ms) |
|---|---|---|---|---|
| Fat | 260 | 370 | 80 | 130 |
| White matter | 780 | 830 | 90 | 80 |
| Gray matter | 920 | 1330 | 100 | 110 |
| CSF | 3500 | 4300 | 1800 | 1800 |
| Muscle | 880 | 1410 | 47 | 50 |
| Liver | 570 | 810 | 46 | 42 |
cat(sprintf("Mean T1 increase from 1.5 T to 3 T: %.0f%%; mean T2 change: %+.0f%%\n",
100*(mean(relax$T1_3/relax$T1_1p5) - 1),
100*(mean(relax$T2_3/relax$T2_1p5) - 1)))## Mean T1 increase from 1.5 T to 3 T: 36%; mean T2 change: +10%fs <- 8000; tt <- seq(0, 2, by = 1/fs)
f_off <- 60; T2star <- 0.04
sig_c <- exp(-tt/T2star)*exp(1i*2*pi*f_off*tt) # complex FID (quadrature)
df_fid <- rbind(
data.frame(t = tt*1000, v = Re(sig_c), ch = "real channel"),
data.frame(t = tt*1000, v = Im(sig_c), ch = "imaginary channel"),
data.frame(t = tt*1000, v = exp(-tt/T2star), ch = "envelope exp(-t/T2*)"))
p_fid <- ggplot(subset(df_fid, t <= 200), aes(t, v, colour = ch)) +
geom_line(linewidth = 0.7) +
scale_colour_manual(values = c("real channel" = bpad_pal[1],
"imaginary channel" = bpad_pal[6],
"envelope exp(-t/T2*)" = bpad_pal[2])) +
labs(x = "Time (ms)", y = "Signal", title = "Free-induction decay (complex)")
N <- length(sig_c)
Z <- fft(sig_c)/N
fr <- (0:(N - 1))*fs/N
fr <- ifelse(fr >= fs/2, fr - fs, fr)
o <- order(fr); frs <- fr[o]
absn <- Re(Z[o]); absn <- absn/max(absn) # absorption mode (phased)
magn <- Mod(Z[o]); magn <- magn/max(magn) # magnitude mode
fwhm_of <- function(f, y) {
p <- which.max(y)
lo <- approx(y[1:p], f[1:p], 0.5)$y
hi <- approx(rev(y[p:length(y)]), rev(f[p:length(y)]), 0.5)$y
hi - lo
}
k <- frs > 0 & frs < 200
spec <- rbind(data.frame(f = frs[k], S = absn[k], mode = "absorption (Re, phased)"),
data.frame(f = frs[k], S = magn[k], mode = "magnitude |Z|"))
p_spec <- ggplot(spec, aes(f, S, colour = mode)) +
geom_hline(yintercept = 0.5, linetype = "dotted", colour = "grey45") +
geom_line(linewidth = 0.9) +
coord_cartesian(xlim = c(20, 100)) +
scale_colour_manual(values = bpad_pal[c(3,2)]) +
labs(x = "Frequency (Hz)", y = "Normalized amplitude",
title = "Lorentzian lineshape")
bpad_grid(p_fid, p_spec, ncol = 2)
Figure 5: The free-induction decay and its spectrum. Left: the complex FID that quadrature detection actually records, shown as its real and imaginary channels inside the exp(-t/T2) envelope. Right: the Fourier transform, displayed two ways. The absorption-mode spectrum (the real part of a correctly phased signal) is a true Lorentzian of full width 1/(pi T2); the magnitude spectrum, which discards phase, is a square-root Lorentzian and is wider by exactly the square root of 3. Both widths are verified numerically below. This is why MR spectroscopy is phased and displayed in absorption mode rather than as a magnitude.
cat(sprintf("absorption mode: theory 1/(pi T2*) = %.3f Hz, measured %.3f Hz\n",
1/(pi*T2star), fwhm_of(frs[k], absn[k])))## absorption mode: theory 1/(pi T2*) = 7.958 Hz, measured 7.972 Hzcat(sprintf("magnitude mode: theory sqrt(3)/(pi T2*) = %.3f Hz, measured %.3f Hz\n",
sqrt(3)/(pi*T2star), fwhm_of(frs[k], magn[k])))## magnitude mode: theory sqrt(3)/(pi T2*) = 13.783 Hz, measured 13.791 Hzcat(sprintf("ratio = %.4f = sqrt(3); discarding phase costs %.0f%% in linewidth.\n",
sqrt(3), 100*(sqrt(3) - 1)))## ratio = 1.7321 = sqrt(3); discarding phase costs 73% in linewidth.## A poorly shimmed voxel with T2* = 10 ms would broaden the absorption linecat(sprintf("to %.1f Hz, which at 3 T is %.2f ppm and would merge adjacent metabolites.\n",
1/(pi*0.010), (1/(pi*0.010))/(gbar_H*3)*1e6))## to 31.8 Hz, which at 3 T is 0.25 ppm and would merge adjacent metabolites.Section 4.2 summary.
- The Bloch equations, Eq. (9), combine precession with two relaxation terms and decouple into independent \(T_1\) and \(T_2\) behaviour without RF.
- Recovery follows \(M_z(t) = M_0 + [M_z(0)-M_0]e^{-t/T_1}\); saturation, inversion, and partial excitation differ only in \(M_z(0)\).
- \(T_2 \le T_1\) always; \(T_2^* \le T_2\) adds static inhomogeneity, and only the static part is refocusable.
- \(T_1\) rises markedly with field while \(T_2\) does not, so protocols must be retuned when moving from 1.5 T to 3 T.
- The FID and the spectrum are a Fourier pair with Lorentzian linewidth \(1/(\pi T_2^*)\).
Checkpoint 4.2. Why can a \(180^\circ\) pulse recover signal lost to \(T_2^*\) but not signal lost to true \(T_2\)? Frame your answer in terms of which dephasing is a deterministic function of position and which is random in time.
4.3 Image Contrast: TR, TE, and Weighting
4.3.1 The spin-echo signal equation
Image contrast is chosen by two timing parameters: the repetition time TR (between successive excitations) and the echo time TE (between excitation and readout). For a spin-echo sequence with \(90^\circ\) excitation, the signal from a tissue of proton density \(\rho\) is
\[\begin{equation} S \;\propto\; \rho\left(1 - e^{-TR/T_1}\right)e^{-TE/T_2} . \tag{14} \end{equation}\]
The first factor carries \(T_1\) information (set by TR) and the second carries \(T_2\) information (set by TE). Choosing TR and TE therefore selects the weighting:
| Weighting | TR | TE | Bright | Clinical use |
|---|---|---|---|---|
| T1 | short (400–700 ms) | short (10–20 ms) | fat, gadolinium; CSF dark | anatomy, post-contrast |
| T2 | long (>2000 ms) | long (80–120 ms) | fluid, oedema | pathology |
| Proton density | long (>2000 ms) | short (10–20 ms) | high \(\rho\) | subtle grey/white contrast |
From the equation to the image. Equation (14) is a recipe. Want CSF bright, to find a fluid collection? Make TE long so its very long \(T_2\) dominates. Want fat bright and crisp anatomy? Make TR short so short-\(T_1\) tissues stand out. The same tissue is bright or dark depending only on timing — which is why an MRI report must always state the sequence, and why “bright on MRI” is a meaningless phrase on its own.
sig_eq <- function(rho, T1, T2, TR, TE) rho*(1 - exp(-TR/T1))*exp(-TE/T2)
TE_g <- seq(0, 160, by = 1)
df_c <- do.call(rbind, lapply(seq_len(nrow(tissues)), function(i)
data.frame(TE = TE_g, S = sig_eq(1, tissues$T1[i], tissues$T2[i], 2500, TE_g),
tissue = tissues$tissue[i])))
df_c$tissue <- factor(df_c$tissue, levels = tissues$tissue)
p_te <- ggplot(df_c, aes(TE, S, colour = tissue)) +
geom_line(linewidth = 0.9) +
geom_vline(xintercept = c(15, 100), linetype = "dotted", colour = "grey45") +
annotate("text", x = 15, y = 0.98, label = "T1w TE", size = 2.8, hjust = -0.1) +
annotate("text", x = 100, y = 0.98, label = "T2w TE", size = 2.8, hjust = -0.1) +
scale_colour_manual(values = cols) +
labs(x = "Echo time TE (ms)", y = "Relative signal",
title = "Signal versus TE (TR = 2500 ms)")
TRv <- seq(100, 3000, length.out = 120)
TEv <- seq(2, 150, length.out = 120)
gm <- tissues[tissues$tissue == "Gray matter", ]
wm <- tissues[tissues$tissue == "White matter", ]
grid <- expand.grid(TR = TRv, TE = TEv)
grid$C <- sig_eq(1, gm$T1, gm$T2, grid$TR, grid$TE) -
sig_eq(1, wm$T1, wm$T2, grid$TR, grid$TE)
p_map <- ggplot(grid, aes(TR, TE, fill = C)) +
geom_raster() +
geom_contour(data = grid, aes(TR, TE, z = C), breaks = 0,
colour = "black", linewidth = 0.5, inherit.aes = FALSE) +
scale_fill_gradient2(low = bpad_pal[1], mid = "white", high = bpad_pal[2],
midpoint = 0) +
theme(legend.position = "right", legend.title = element_text(size = 8)) +
labs(x = "TR (ms)", y = "TE (ms)", fill = "GM - WM",
title = "Grey-white contrast over the TR-TE plane")
bpad_grid(p_te, p_map, ncol = 2)
Figure 6: Contrast is engineered, not fixed. Left: signal versus TE at long TR; at short TE the tissues cluster, while at long TE CSF overtakes everything. Right: the contrast between grey and white matter across the full TR-TE plane, with red showing grey matter brighter and blue showing white matter brighter. The zero-contrast contour running through the plane is the set of timings at which the two tissues are indistinguishable, and no amount of receiver gain can recover a difference that the physics did not create.
cat(sprintf("At TR=500, TE=15 (T1w): GM %.3f, WM %.3f -> WM brighter by %.3f\n",
sig_eq(1, gm$T1, gm$T2, 500, 15), sig_eq(1, wm$T1, wm$T2, 500, 15),
sig_eq(1, wm$T1, wm$T2, 500, 15) - sig_eq(1, gm$T1, gm$T2, 500, 15)))## At TR=500, TE=15 (T1w): GM 0.361, WM 0.401 -> WM brighter by 0.040cat(sprintf("At TR=2500, TE=100 (T2w): GM %.3f, WM %.3f -> GM brighter by %.3f\n",
sig_eq(1, gm$T1, gm$T2, 2500, 100), sig_eq(1, wm$T1, wm$T2, 2500, 100),
sig_eq(1, gm$T1, gm$T2, 2500, 100) - sig_eq(1, wm$T1, wm$T2, 2500, 100)))## At TR=2500, TE=100 (T2w): GM 0.344, WM 0.316 -> GM brighter by 0.028## The grey-white contrast literally reverses sign between the two protocols.4.3.2 Inversion recovery: FLAIR and STIR
Adding a \(180^\circ\) inversion pulse before excitation provides a third contrast lever through the inversion time TI. From the inversion row of the table in Section 4.2.2, \(M_z(TI) = M_0(1 - 2e^{-TI/T_1})\), which passes through zero when
\[\begin{equation} TI_{\text{null}} = T_1\ln 2 \approx 0.693\,T_1 \qquad (\text{TR} \gg T_1) . \tag{15} \end{equation}\]
Equation (15) is the formula every textbook quotes, and for fat suppression it works. For FLAIR it does not, because CSF has such a long \(T_1\) that TR is not \(\gg T_1\) and the magnetization never fully recovers between repetitions. Solving the steady state of the inversion-recovery sequence gives the corrected null time
\[\begin{equation} TI_{\text{null}} = T_1\ln\!\left(\frac{2}{1 + e^{-TR/T_1}}\right) . \tag{16} \end{equation}\]
ti_simple <- function(T1) T1*log(2)
ti_finite <- function(T1, TR) T1*log(2/(1 + exp(-TR/T1)))
tg <- seq(0, 5000, by = 5)
inv <- rbind(
data.frame(t = tg, M = 1 - 2*exp(-tg/370), tissue = "Fat (T1 = 370 ms, 3 T)"),
data.frame(t = tg, M = 1 - 2*exp(-tg/4300), tissue = "CSF (T1 = 4300 ms, 3 T)"))
p_inv <- ggplot(inv, aes(t, M, colour = tissue)) +
geom_hline(yintercept = 0, colour = "grey55") +
geom_line(linewidth = 0.95) +
geom_vline(xintercept = c(ti_simple(370), ti_simple(4300)),
linetype = "dotted", colour = "grey45") +
scale_colour_manual(values = bpad_pal[c(5,3)]) +
labs(x = "Inversion time TI (ms)", y = expression(M[z]/M[0]),
title = "Inversion recovery: nulling a tissue")
TRg <- seq(1000, 12000, by = 50)
nul <- rbind(
data.frame(TR = TRg, TI = ti_finite(4300, TRg), tissue = "CSF, T1 = 4300 ms"),
data.frame(TR = TRg, TI = ti_finite(370, TRg), tissue = "Fat, T1 = 370 ms"))
p_nul <- ggplot(nul, aes(TR, TI, colour = tissue)) +
geom_line(linewidth = 0.95) +
geom_hline(yintercept = c(ti_simple(4300), ti_simple(370)),
linetype = "dashed", colour = "grey45") +
annotate("text", x = 11500, y = ti_simple(4300) + 130, hjust = 1, size = 2.8,
colour = "grey30", label = "T1 ln2 asymptote") +
scale_colour_manual(values = bpad_pal[c(3,5)]) +
labs(x = "Repetition time TR (ms)", y = "Null time TI (ms)",
title = "Finite TR shortens the null time")
bpad_grid(p_inv, p_nul, ncol = 2)
Figure 7: Inversion recovery and the finite-TR correction. Left: inverted longitudinal magnetization crossing zero, with the null time marked for fat and CSF. Right: null time versus TR from Eq. (ti-null-finite), approaching the textbook T1 ln 2 only as TR becomes much longer than T1. For fat the correction is negligible, but for CSF at 3 T the simple formula overestimates the FLAIR inversion time by about 500 ms, and the corrected value is what clinical protocols actually use.
knitr::kable(data.frame(
application = c("FLAIR 1.5 T", "FLAIR 3 T", "STIR 1.5 T", "STIR 3 T"),
tissue = c("CSF","CSF","Fat","Fat"),
T1 = c(3500, 4300, 260, 370), TR = c(9000, 9000, 2500, 2500),
TI_simple = round(ti_simple(c(3500,4300,260,370))),
TI_finite = round(ti_finite(c(3500,4300,260,370), c(9000,9000,2500,2500))),
difference = round(ti_simple(c(3500,4300,260,370)) -
ti_finite(c(3500,4300,260,370), c(9000,9000,2500,2500)))),
col.names = c("Application","Nulled tissue","T1 (ms)","TR (ms)",
"TI from T1 ln2","TI corrected","Difference (ms)"),
caption = "Inversion times for the two standard nulling sequences. Clinical 3 T FLAIR uses TI near 2500 ms, which matches the corrected formula and not the textbook one.")| Application | Nulled tissue | T1 (ms) | TR (ms) | TI from T1 ln2 | TI corrected | Difference (ms) |
|---|---|---|---|---|---|---|
| FLAIR 1.5 T | CSF | 3500 | 9000 | 2426 | 2168 | 258 |
| FLAIR 3 T | CSF | 4300 | 9000 | 2981 | 2481 | 500 |
| STIR 1.5 T | Fat | 260 | 2500 | 180 | 180 | 0 |
| STIR 3 T | Fat | 370 | 2500 | 256 | 256 | 0 |
Why the textbook formula fails for FLAIR. \(TI = T_1\ln 2\) assumes the magnetization has fully recovered to \(+M_0\) before each inversion. CSF at 3 T has \(T_1 \approx 4300\) ms while FLAIR uses \(TR \approx 9000\) ms, so recovery reaches only \(1 - e^{-9000/4300} = 88\%\). The inversion therefore starts from \(-0.88M_0\) rather than \(-M_0\), and the zero crossing arrives earlier. Using the uncorrected value would leave residual CSF signal and defeat the entire purpose of the sequence.
Section 4.3 summary.
- \(S \propto \rho(1-e^{-TR/T_1})e^{-TE/T_2}\): TR selects \(T_1\) weighting and TE selects \(T_2\) weighting.
- Grey–white contrast reverses sign between T1- and T2-weighted timings, and there is a zero-contrast contour where the tissues are indistinguishable.
- Inversion recovery nulls a chosen tissue at \(T_1\ln2\) only when \(TR \gg T_1\); FLAIR requires the finite-TR correction, Eq. (16).
Checkpoint 4.3. Using the TR–TE contrast map, find a pair of timings at which grey and white matter are equally bright. What clinical consequence follows if a protocol accidentally lands near that contour, and which weighting would you switch to?
4.4 Spatial Encoding and k-Space
4.4.1 Slice selection
Gradients are weak, spatially varying magnetic fields added to \(B_0\). A slice-select gradient \(G_z\) makes the field, and therefore the Larmor frequency, depend on position:
\[\begin{equation} \omega(z) = \gamma\big(B_0 + G_z z\big) . \tag{17} \end{equation}\]
An RF pulse of centre frequency \(f_c\) and bandwidth \(\mathrm{BW}_{\rm RF}\) then excites only the slab in resonance, of thickness
\[\begin{equation} \Delta z = \frac{\mathrm{BW}_{\rm RF}}{\bar\gamma\, G_z} , \tag{18} \end{equation}\]
and shifting \(f_c\) moves the slice without moving the patient. Frequency has become an address.
4.4.2 Frequency and phase encoding
Within the slice, two more gradients encode in-plane position. Frequency encoding \(G_x\), applied during readout, makes precession frequency depend on \(x\). Phase encoding \(G_y\), applied briefly before readout, makes accumulated phase depend on \(y\):
\[\begin{equation} \omega(x) = \gamma(B_0 + G_x x), \qquad \phi(y) = \gamma\, G_y\, y\, \tau_{\rm pe} . \tag{19} \end{equation}\]
Repeating the experiment with stepped \(G_y\) builds the second dimension line by line — which is why the phase-encode axis is sampled slowly, over the whole scan, and why motion artifacts appear along it (Section 4.7).
4.4.3 k-space and the Fourier relationship
Collecting these effects, the measured signal as a function of accumulated gradient area is
\[\begin{equation} S(k_x, k_y) = \iint \rho(x,y)\, e^{-i2\pi(k_x x + k_y y)}\, dx\,dy, \qquad k_{x}(t) = \bar\gamma\!\int_0^t\! G_x(t')\,dt' , \tag{20} \end{equation}\]
which is precisely the two-dimensional Fourier transform of the image. The raw data live in k-space (spatial-frequency space); the image is recovered by the inverse transform:
\[\begin{equation} I(x,y) = \mathcal{F}^{-1}\big\{S(k_x,k_y)\big\} . \tag{21} \end{equation}\]
Because k-space and the image are a Fourier pair, their geometries are reciprocal:
\[\begin{equation} \mathrm{FOV} = \frac{1}{\Delta k}, \qquad \Delta x = \frac{1}{2k_{\max}} = \frac{\mathrm{FOV}}{N} . \tag{22} \end{equation}\]
The sampling step \(\Delta k\) sets the field of view; the extent \(k_{\max}\) sets the resolution. The centre of k-space (low spatial frequencies) carries gross contrast and brightness; the periphery carries edges and fine detail. Every k-space point contributes to every image pixel.
Connection to Chapter 1 (the central one). MRI does not measure pixels; it measures the Fourier transform of the image, point by point, as gradients trace a path through k-space. Reconstruction is the inverse 2D Fourier transform, in practice the 2D FFT. This is why so much of MRI — resolution, field of view, aliasing, ringing, motion ghosting — is most clearly understood in the Fourier domain. Equation (22) is the Nyquist criterion in disguise: sample k-space too coarsely and anatomy beyond the FOV aliases back into the image. Lab 1 makes all of this concrete.
Gz <- seq(2, 40, by = 0.2) # mT/m
sl <- do.call(rbind, lapply(c(500, 1000, 2000, 4000), function(BW)
data.frame(Gz = Gz, dz = BW/(gbar_H*Gz*1e-3)*1e3,
BW = factor(sprintf("%d Hz", BW),
levels = sprintf("%d Hz", c(500,1000,2000,4000))))))
p_sl <- ggplot(sl, aes(Gz, dz, colour = BW)) +
geom_line(linewidth = 0.9) +
geom_hline(yintercept = 1, linetype = "dotted", colour = "grey45") +
annotate("text", x = 38, y = 1.35, hjust = 1, size = 2.9, colour = "grey30",
label = "1 mm slice") +
scale_y_log10() + scale_colour_manual(values = bpad_pal[1:4]) +
labs(x = expression(G[z]*" (mT/m)"), y = "Slice thickness (mm, log scale)",
title = "Slice thickness versus gradient and RF bandwidth")
FOVv <- c(0.20, 0.25, 0.30, 0.40)
Nmat <- c(64, 128, 192, 256, 384, 512)
recip <- expand.grid(FOV = FOVv, N = Nmat)
recip$dx_mm <- recip$FOV/recip$N*1e3
recip$kmax <- 1/(2*recip$FOV/recip$N)
p_rc <- ggplot(recip, aes(N, dx_mm, colour = factor(FOV*100))) +
geom_line(linewidth = 0.9) + geom_point(size = 1.6) +
scale_x_log10(breaks = Nmat) + scale_y_log10() +
scale_colour_manual(values = bpad_pal[1:4], name = "FOV (cm)") +
labs(x = "Matrix size N (log scale)", y = "Pixel size (mm, log scale)",
title = "Pixel size = FOV / N")
bpad_grid(p_sl, p_rc, ncol = 2)
Figure 8: Encoding in physical units. Left: slice thickness from Eq. (slice-thickness); a stronger gradient or narrower RF bandwidth gives a thinner slice, and reaching 1 mm requires either a strong gradient or a long, narrow-band pulse. Right: the reciprocal geometry of Eq. (fov-resolution); the k-space sampling step sets the field of view while the k-space extent sets the pixel size, so a scan that halves the pixel size must double k_max and therefore double the readout duration or the gradient strength.
knitr::kable(data.frame(
Gz_mT_m = c(5, 10, 20, 40),
dz_at_1kHz_mm = round(1000/(gbar_H*c(5,10,20,40)*1e-3)*1e3, 2),
dz_at_2kHz_mm = round(2000/(gbar_H*c(5,10,20,40)*1e-3)*1e3, 2),
readout_BW_kHz_25cmFOV = round(gbar_H*c(5,10,20,40)*1e-3*0.25/1e3, 1)),
col.names = c("Gz (mT/m)","Slice at 1 kHz BW (mm)","Slice at 2 kHz BW (mm)",
"Readout BW over 25 cm FOV (kHz)"),
caption = "Gradient amplitudes translated into the quantities a protocol actually specifies.")| Gz (mT/m) | Slice at 1 kHz BW (mm) | Slice at 2 kHz BW (mm) | Readout BW over 25 cm FOV (kHz) |
|---|---|---|---|
| 5 | 4.70 | 9.39 | 53.2 |
| 10 | 2.35 | 4.70 | 106.4 |
| 20 | 1.17 | 2.35 | 212.9 |
| 40 | 0.59 | 1.17 | 425.8 |
Worked Example 4.2 (from gradient to slice). With \(G_z = 10\) mT/m and \(\mathrm{BW}_{\rm RF} = 2\) kHz,
\[\bar\gamma G_z = (42.577\times10^{6}\ \text{Hz T}^{-1})(0.010\ \text{T m}^{-1}) = 4.258\times10^{5}\ \text{Hz m}^{-1},\]
\[\Delta z = \frac{2000}{4.258\times10^{5}} = 4.70\times10^{-3}\ \text{m} = 4.7\ \text{mm}.\]
To halve the slice to 2.35 mm, either double \(G_z\) (limited by gradient hardware and by peripheral nerve stimulation, Section 4.8) or halve the RF bandwidth — which doubles the pulse duration and therefore lengthens TE.
Section 4.4 summary.
- Gradients turn position into frequency (\(z\), \(x\)) and phase (\(y\)).
- Slice thickness is \(\mathrm{BW}_{\rm RF}/(\bar\gamma G_z)\).
- The acquired signal is the 2D Fourier transform of the image; \(\mathrm{FOV} = 1/\Delta k\) and \(\Delta x = 1/(2k_{\max})\).
- k-space centre carries contrast, periphery carries detail; reconstruction is the inverse FFT.
Checkpoint 4.4. If you acquired only the central quarter of k-space in each direction, what would the image look like, what would happen to the pixel size, and what would happen to the field of view? Answer using Eq. (22) rather than intuition.
4.5 Pulse Sequences
4.5.1 Spin echo
A spin-echo (SE) sequence applies a \(90^\circ\) excitation, waits \(\tau\), applies a \(180^\circ\) refocusing pulse, and reads the echo at \(TE = 2\tau\):
\[\begin{equation} 90^\circ \rightarrow \tau \rightarrow 180^\circ \rightarrow \tau \rightarrow \text{echo}, \qquad TE = 2\tau . \tag{23} \end{equation}\]
The \(180^\circ\) pulse inverts the accumulated phases so that dephasing which is a deterministic function of position — the static \(T_2'\) part — unwinds and re-forms an echo at TE. Random spin–spin dephasing does not unwind, because a spin that has already diffused into a different local field cannot retrace its history. Spin echo therefore measures true \(T_2\), free of \(T_2^*\), at the cost of longer scans and much higher RF deposition (Section 4.8).
4.5.2 Gradient echo and the Ernst angle
A gradient-echo (GE) sequence omits the \(180^\circ\) pulse and instead reverses a gradient to rephase the signal. It is far faster and uses small flip angles \(\alpha\), but because nothing refocuses static dephasing it is inherently \(T_2^*\)-weighted and sensitive to susceptibility.
For repeated excitations at interval TR, the spoiled steady-state signal is
\[\begin{equation} S(\alpha) = M_0 \sin\alpha \, \frac{1 - E_1}{1 - E_1\cos\alpha}\, e^{-TE/T_2^{*}}, \qquad E_1 = e^{-TR/T_1} . \tag{24} \end{equation}\]
Differentiating with respect to \(\alpha\) and setting the derivative to zero gives the Ernst angle, the flip angle that maximizes steady-state signal:
\[\begin{equation} \cos\alpha_E = E_1 = e^{-TR/T_1} . \tag{25} \end{equation}\]
The derivation is short enough to be worth seeing. Writing \(S \propto \sin\alpha/(1 - E_1\cos\alpha)\) and differentiating, \[\frac{dS}{d\alpha} \propto \frac{\cos\alpha(1 - E_1\cos\alpha) - \sin\alpha(E_1\sin\alpha)}{(1-E_1\cos\alpha)^2} = \frac{\cos\alpha - E_1}{(1-E_1\cos\alpha)^2},\] using \(\sin^2 + \cos^2 = 1\). The numerator vanishes exactly when \(\cos\alpha = E_1\).
gre <- function(alpha, TR, T1) {
E1 <- exp(-TR/T1)
sin(alpha)*(1 - E1)/(1 - E1*cos(alpha))
}
a_g <- seq(0, pi/2, length.out = 400)
T1ref <- 800
TRs <- c(5, 20, 100, 500)
df_e <- do.call(rbind, lapply(TRs, function(TR)
data.frame(a = a_g*180/pi, S = gre(a_g, TR, T1ref),
TR = factor(sprintf("TR = %d ms", TR),
levels = sprintf("TR = %d ms", TRs)))))
opt <- data.frame(TR = factor(sprintf("TR = %d ms", TRs),
levels = sprintf("TR = %d ms", TRs)),
a = acos(exp(-TRs/T1ref))*180/pi,
S = gre(acos(exp(-TRs/T1ref)), TRs, T1ref))
p_e1 <- ggplot(df_e, aes(a, S, colour = TR)) +
geom_line(linewidth = 0.9) +
geom_point(data = opt, size = 2.4) +
scale_colour_manual(values = bpad_pal[1:4]) +
labs(x = "Flip angle (degrees)", y = "Steady-state signal",
title = expression("Spoiled gradient echo, "*T[1]*" = 800 ms"))
rat <- 10^seq(-3, 1, length.out = 400)
p_e2 <- ggplot(data.frame(r = rat, a = acos(exp(-rat))*180/pi), aes(r, a)) +
geom_line(linewidth = 0.95, colour = bpad_pal[1]) +
scale_x_log10() +
labs(x = expression("TR / "*T[1]*" (log scale)"), y = "Ernst angle (degrees)",
title = "Ernst angle versus TR/T1")
bpad_grid(p_e1, p_e2, ncol = 2)
Figure 9: The gradient-echo steady state and the Ernst angle. Left: signal versus flip angle at four repetition times, Eq. (gre-signal), with the optimum marked; short TR forces small flip angles, which is why fast gradient-echo imaging looks nothing like a 90-degree spin echo. Right: the Ernst angle as a function of TR/T1, Eq. (ernst). Note that the peak is broad, so being a few degrees off the optimum costs very little signal, and that tissues with different T1 have different optima, which is itself a source of contrast.
knitr::kable(data.frame(
TR_ms = TRs, E1 = round(exp(-TRs/T1ref), 4),
ernst_deg = round(acos(exp(-TRs/T1ref))*180/pi, 1),
signal_at_ernst = round(gre(acos(exp(-TRs/T1ref)), TRs, T1ref), 4),
signal_at_90 = round(gre(pi/2, TRs, T1ref), 4)),
col.names = c("TR (ms)","E1","Ernst angle (deg)","Signal at Ernst","Signal at 90 deg"),
caption = "At short TR the Ernst angle is small and substantially outperforms a 90-degree pulse, because a large flip leaves no longitudinal magnetization to recover before the next excitation.")| TR (ms) | E1 | Ernst angle (deg) | Signal at Ernst | Signal at 90 deg |
|---|---|---|---|---|
| 5 | 0.9938 | 6.4 | 0.0559 | 0.0062 |
| 20 | 0.9753 | 12.8 | 0.1118 | 0.0247 |
| 100 | 0.8825 | 28.1 | 0.2498 | 0.1175 |
| 500 | 0.5353 | 57.6 | 0.5502 | 0.4647 |
Spin echo versus gradient echo. SE uses \(180^\circ\) refocusing, so it measures true \(T_2\), is robust to field inhomogeneity, and is slow and RF-intensive. GE uses gradient reversal and a small flip angle, so it is fast, low-SAR, and \(T_2^*\)-weighted — which is a feature for detecting haemorrhage and for BOLD fMRI, and a bug near air, bone, and metal. The choice is not about image quality in the abstract; it is about whether susceptibility is the signal you want or the noise you must avoid.
Section 4.5 summary.
- Spin echo refocuses position-dependent dephasing and measures true \(T_2\).
- Gradient echo is fast and \(T_2^*\)-weighted, with steady-state signal Eq. (24).
- The Ernst angle \(\cos\alpha_E = e^{-TR/T_1}\) follows in three lines from that equation, and is small whenever \(TR \ll T_1\).
Checkpoint 4.5. Why is gradient echo the sequence of choice for detecting a small old haemorrhage (haemosiderin)? Which term in Eq. (24) carries the effect, and what would a spin echo do to it?
4.6 Reconstruction, Resolution, and SNR
4.6.1 The SNR relation, stated carefully
SNR in MRI is governed by how much magnetization a voxel contains and how long it is observed:
\[\begin{equation} \mathrm{SNR} \;\propto\; \underbrace{\Delta x\,\Delta y\,\Delta z}_{\text{voxel volume}} \;\sqrt{\frac{N_x N_y N_z \cdot \mathrm{NSA}}{\mathrm{BW}}} \;=\; V_{\rm voxel}\sqrt{T_{\rm acq}} . \tag{26} \end{equation}\]
A trap that catches almost everyone. It is often said that “halving each voxel dimension reduces SNR by a factor of 8.” That is true only if the number of acquired samples is held fixed, which is not what happens on a scanner. In the realistic case — fixed field of view, finer matrix — halving the in-plane voxel dimensions doubles \(N_x\) and \(N_y\), so \(N_{\rm samples}\) rises by 4 and Eq. (26) gives
\[\mathrm{SNR} \propto \frac{V}{4}\times\sqrt{4} = \frac{V}{2}:\]
SNR falls by 2, not 8. Always state which quantities are held constant before quoting an SNR scaling. The table below works through the three scenarios that actually arise.
scen <- data.frame(
scenario = c("fixed FOV, finer matrix (in-plane only)",
"fixed FOV, finer matrix in all 3 dims",
"fixed samples (FOV shrinks with voxel)"),
V_rel = c(1/4, 1/8, 1/8),
N_rel = c(4, 8, 1))
scen$SNR_rel <- scen$V_rel*sqrt(scen$N_rel)
scen$SNR_drop <- 1/scen$SNR_rel
p_s1 <- ggplot(scen, aes(reorder(scenario, -SNR_drop), SNR_drop)) +
geom_col(fill = bpad_pal[1]) +
geom_text(aes(label = sprintf("%.1fx", SNR_drop)), hjust = -0.15, size = 3.2) +
coord_flip(ylim = c(0, 10)) +
labs(x = NULL, y = "Factor by which SNR falls",
title = "Halving voxel dimensions: three scenarios")
nsa <- 1:16
p_s2 <- ggplot(data.frame(n = nsa, s = sqrt(nsa)), aes(n, s)) +
geom_line(linewidth = 0.95, colour = bpad_pal[1]) +
geom_point(size = 1.8, colour = bpad_pal[1]) +
geom_segment(aes(x = 4, xend = 4, y = 0, yend = 2), linetype = "dotted",
colour = bpad_pal[2]) +
geom_segment(aes(x = 0, xend = 4, y = 2, yend = 2), linetype = "dotted",
colour = bpad_pal[2]) +
annotate("text", x = 4.4, y = 1.5, hjust = 0, size = 3, colour = bpad_pal[2],
label = "4x the time\nfor 2x the SNR") +
labs(x = "Number of signal averages (NSA)", y = "Relative SNR",
title = "Averaging obeys a square-root law")
bpad_grid(p_s1, p_s2, ncol = 2)
Figure 10: SNR scaling depends entirely on what is held fixed. Left: three ways to improve resolution, each with a different SNR cost; the naive fixed-sample scenario is the harshest and the least physically realistic. Right: the square-root law for signal averaging, showing why doubling SNR by averaging costs a fourfold increase in scan time and why NSA is usually the last resort rather than the first.
knitr::kable(data.frame(
scenario = scen$scenario, voxel_volume = c("1/4","1/8","1/8"),
samples = c("x4","x8","x1"), SNR_relative = round(scen$SNR_rel, 3),
SNR_falls_by = sprintf("%.2fx", scen$SNR_drop)),
col.names = c("Scenario","Voxel volume","Samples","Relative SNR","SNR falls by"),
caption = "The same phrase 'halve the voxel dimensions' produces three different answers. The middle row is the usual clinical situation for a 3D acquisition; the last row is the textbook claim and requires the field of view to shrink with the voxel.")| Scenario | Voxel volume | Samples | Relative SNR | SNR falls by |
|---|---|---|---|---|
| fixed FOV, finer matrix (in-plane only) | 1/4 | x4 | 0.500 | 2.00x |
| fixed FOV, finer matrix in all 3 dims | 1/8 | x8 | 0.354 | 2.83x |
| fixed samples (FOV shrinks with voxel) | 1/8 | x1 | 0.125 | 8.00x |
4.6.2 Magnitude images are Rician, not Gaussian
The scanner records a complex signal with independent Gaussian noise of standard deviation \(\sigma\) in each channel. The displayed image is usually the magnitude, and the magnitude of a complex Gaussian is Rician:
\[\begin{equation} p(m \mid A,\sigma) = \frac{m}{\sigma^{2}} \exp\!\left(-\frac{m^{2}+A^{2}}{2\sigma^{2}}\right) I_0\!\left(\frac{mA}{\sigma^{2}}\right), \tag{27} \end{equation}\]
with \(A\) the true signal and \(I_0\) the modified Bessel function. Two consequences matter. At high SNR the distribution is nearly Gaussian with mean \(\approx A\). At low SNR it is skewed and biased upward, and in pure noise (\(A = 0\)) it reduces to a Rayleigh distribution with mean \(\sigma\sqrt{\pi/2}\) — a noise floor that no amount of averaging removes.
set.seed(1)
sg <- 1; nsim <- 120000
lv <- c(0, 1, 2, 5)
ric <- do.call(rbind, lapply(lv, function(A)
data.frame(m = Mod(complex(real = A + rnorm(nsim, 0, sg),
imaginary = rnorm(nsim, 0, sg))),
A = factor(sprintf("true A = %d", A),
levels = sprintf("true A = %d", lv)))))
p_r1 <- ggplot(ric, aes(m, fill = A)) +
geom_density(alpha = 0.45, colour = NA) +
geom_vline(xintercept = sg*sqrt(pi/2), linetype = "dotted", colour = "grey35") +
scale_fill_manual(values = bpad_pal[c(2,5,3,1)]) +
coord_cartesian(xlim = c(0, 9)) +
labs(x = "Measured magnitude", y = "Density",
title = "Rician distribution")
Ag <- seq(0, 8, by = 0.1)
bias <- vapply(Ag, function(A) {
m <- Mod(complex(real = A + rnorm(20000, 0, sg), imaginary = rnorm(20000, 0, sg)))
mean(m) - A }, numeric(1))
p_r2 <- ggplot(data.frame(A = Ag, b = bias), aes(A, b)) +
geom_hline(yintercept = 0, colour = "grey60") +
geom_line(linewidth = 0.95, colour = bpad_pal[1]) +
geom_hline(yintercept = sg*sqrt(pi/2), linetype = "dotted", colour = "grey40") +
annotate("text", x = 7.8, y = sg*sqrt(pi/2) + 0.06, hjust = 1, size = 2.8,
colour = "grey30", label = "noise floor at A = 0") +
labs(x = expression("True signal A (in units of "*sigma*")"),
y = "Bias in the magnitude", title = "The bias never vanishes")
S0 <- 100; ADCtrue <- 0.8e-3; sg2 <- 5
bvals <- seq(0, 4000, by = 100)
set.seed(3)
meas <- vapply(bvals, function(b) {
St <- S0*exp(-b*ADCtrue)
mean(Mod(complex(real = St + rnorm(40000, 0, sg2), imaginary = rnorm(40000, 0, sg2))))
}, numeric(1))
p_r3 <- ggplot(rbind(
data.frame(b = bvals, S = S0*exp(-bvals*ADCtrue), q = "true signal"),
data.frame(b = bvals, S = meas, q = "measured magnitude")),
aes(b, S, colour = q)) +
geom_line(linewidth = 0.9) +
geom_hline(yintercept = sg2*sqrt(pi/2), linetype = "dotted", colour = "grey40") +
scale_y_log10() + scale_colour_manual(values = bpad_pal[c(2,1)]) +
labs(x = expression("b-value (s/mm"^2*")"), y = "Signal (log scale)",
title = "Noise floor lifts high-b signal")
bpad_grid(p_r1, p_r2, p_r3, ncol = 3)
Figure 11: Magnitude MRI noise is Rician. Left: the distribution at four true signal levels; at high SNR it is Gaussian and unbiased, but as the true signal approaches zero it becomes Rayleigh with a strictly positive mean. Centre: the resulting bias, which exceeds 20 percent of the noise level below an SNR of about 2 and never reaches zero. Right: the practical consequence, developed further in Section 4.10 – the noise floor makes high-b diffusion signal appear larger than it is, biasing the fitted ADC downward.
cat(sprintf("Rayleigh mean at A = 0: theory sigma*sqrt(pi/2) = %.3f, simulated %.3f\n",
sg*sqrt(pi/2), mean(ric$m[ric$A == "true A = 0"])))## Rayleigh mean at A = 0: theory sigma*sqrt(pi/2) = 1.253, simulated 1.254adc_est <- vapply(c(1000, 2000, 3000, 4000), function(b)
-log(meas[bvals == b]/S0)/b, numeric(1))
knitr::kable(data.frame(
b = c(1000, 2000, 3000, 4000),
true_S = round(S0*exp(-c(1000,2000,3000,4000)*ADCtrue), 2),
measured_S = round(meas[bvals %in% c(1000,2000,3000,4000)], 2),
ADC_estimate = signif(adc_est, 3),
bias_pct = sprintf("%+.1f%%", 100*(adc_est/ADCtrue - 1))),
col.names = c("b (s/mm2)","True signal","Measured magnitude","ADC estimate","ADC bias"),
caption = "Rician bias propagating into a fitted ADC at a b = 0 SNR of 20. The error is negligible at b = 1000 but reaches nearly 20 percent at b = 4000, which is why high-b diffusion protocols require noise-floor correction.")| b (s/mm2) | True signal | Measured magnitude | ADC estimate | ADC bias |
|---|---|---|---|---|
| 1000 | 44.93 | 45.20 | 0.000794 | -0.7% |
| 2000 | 20.19 | 20.83 | 0.000784 | -1.9% |
| 3000 | 9.07 | 10.60 | 0.000748 | -6.5% |
| 4000 | 4.08 | 7.25 | 0.000656 | -18.0% |
4.6.3 Acceleration
Filling all of k-space is slow, so modern MRI acquires less and reconstructs the rest.
| Method | Principle | Main trade-off |
|---|---|---|
| Partial Fourier | approximate Hermitian symmetry of k-space | mild blurring, noise amplification |
| Parallel imaging (SENSE, GRAPPA) | multi-coil sensitivity encoding | \(\sqrt{R}\) SNR loss and a spatially varying g-factor |
| Simultaneous multi-slice | RF-encoded slice separation | inter-slice leakage |
| Compressed sensing | image sparsity plus incoherent sampling | reconstruction assumptions, compute cost |
| Deep-learning reconstruction | learned image priors | generalization and hallucination risk (Ch. 8) |
For parallel imaging with acceleration \(R\) and geometry factor \(g\), Eq. (26) acquires an extra penalty:
\[\begin{equation} \mathrm{SNR}_{R} = \frac{\mathrm{SNR}_{R=1}}{g\sqrt{R}} . \tag{28} \end{equation}\]
Section 4.6 summary.
- \(\mathrm{SNR} \propto V_{\rm voxel}\sqrt{T_{\rm acq}}\); always state what is held fixed before quoting a scaling.
- Averaging obeys a square-root law: doubling SNR costs four times the scan time.
- Magnitude images are Rician, biased upward at low SNR, with a noise floor \(\sigma\sqrt{\pi/2}\) that biases high-b ADC.
- Acceleration trades SNR or assumptions for speed; parallel imaging costs \(g\sqrt{R}\).
Checkpoint 4.6. A scan aliases: anatomy wraps into the image. In k-space terms, which quantity was sampled incorrectly, what is the fix, and what does the fix cost?
4.7 Artifacts and Corrections
4.7.1 Why k-space makes MRI artifact-prone
Because k-space is filled over seconds to minutes, MRI is uniquely vulnerable to anything that changes during acquisition. If anatomy moves while different k-space lines are collected, those lines describe different object positions, and the Fourier transform smears the inconsistency across the image — usually as ghosting along the phase-encode direction, because that axis is sampled slowly, line by line, over the whole scan. Table 9 organizes the major artifacts by the assumption each violates.
4.7.2 A taxonomy of artifacts
| Artifact | Cause | Fourier-domain view | Correction |
|---|---|---|---|
| Motion / ghosting | movement during acquisition | inconsistent phase across k-space lines | gating, navigators, radial sampling, registration |
| Gibbs / truncation ringing | finite k-space extent | sinc point spread from windowing | more k-space, apodization |
| Chemical shift | fat–water offset (3.5 ppm) | misplacement along the frequency axis | higher bandwidth, fat suppression |
| Susceptibility | field distortion near air, bone, metal | local \(T_2^*\) dephasing and geometric warp | spin echo, shimming, shorter TE, higher BW |
| Aliasing / wraparound | FOV smaller than the object | \(\Delta k\) too large (Nyquist violation) | larger FOV, oversampling |
| Zipper / RF leakage | external RF entering the room | a stripe of corrupted k-space | shielding, door seal |
4.7.3 Chemical shift, quantified
Fat and water protons resonate about 3.5 ppm apart. Because the frequency-encode axis interprets frequency as position, that offset displaces fat relative to water by
\[\begin{equation} \Delta x_{\rm fat} = \frac{\Delta f}{\mathrm{BW}_{\rm px}}\ \text{pixels}, \qquad \Delta f = 3.5\times10^{-6}\,\bar\gamma\,B_0 . \tag{29} \end{equation}\]
The offset grows linearly with field, which is why chemical-shift artifact is far more conspicuous at 3 T than at 1.5 T, and why high-field protocols use wider receiver bandwidths despite the SNR cost.
B0c <- seq(0.2, 7.5, by = 0.05)
p_cs1 <- ggplot(data.frame(B0 = B0c, df = ppm_fat*gbar_H*B0c), aes(B0, df)) +
geom_line(linewidth = 0.95, colour = bpad_pal[1]) +
geom_point(data = data.frame(B0 = c(1.5, 3, 7), df = ppm_fat*gbar_H*c(1.5,3,7)),
size = 2.4, colour = bpad_pal[2]) +
geom_text(data = data.frame(B0 = c(1.5, 3, 7), df = ppm_fat*gbar_H*c(1.5,3,7)),
aes(label = sprintf("%.0f Hz", df)), vjust = -0.8, size = 3) +
labs(x = expression(B[0]*" (T)"), y = "Fat-water offset (Hz)",
title = "Chemical shift grows with field")
bwpp <- seq(50, 800, by = 5)
cs <- do.call(rbind, lapply(c(1.5, 3, 7), function(B)
data.frame(bw = bwpp, px = ppm_fat*gbar_H*B/bwpp,
B0 = factor(sprintf("%.1f T", B), levels = sprintf("%.1f T", c(1.5,3,7))))))
p_cs2 <- ggplot(cs, aes(bw, px, colour = B0)) +
geom_line(linewidth = 0.9) +
geom_hline(yintercept = 1, linetype = "dotted", colour = "grey45") +
annotate("text", x = 780, y = 1.25, hjust = 1, size = 2.9, colour = "grey30",
label = "1 pixel") +
scale_y_log10() + scale_colour_manual(values = bpad_pal[1:3]) +
labs(x = "Receiver bandwidth per pixel (Hz)", y = "Fat displacement (pixels, log)",
title = "Misregistration versus bandwidth")
bpad_grid(p_cs1, p_cs2, ncol = 2)
Figure 12: Chemical-shift artifact, quantified by Eq. (chemical-shift). Left: the fat-water frequency offset grows linearly with field strength, reaching about 447 Hz at 3 T. Right: the resulting spatial misregistration in pixels; at a typical 200 Hz per pixel bandwidth the fat image is displaced by more than two pixels at 3 T, which is enough to produce the characteristic bright and dark rims at fat-water boundaries. Raising the bandwidth suppresses the artifact but costs SNR as one over the square root of bandwidth.
knitr::kable(data.frame(
B0_T = c(1.5, 3, 7),
offset_Hz = round(ppm_fat*gbar_H*c(1.5,3,7)),
px_at_100 = round(ppm_fat*gbar_H*c(1.5,3,7)/100, 2),
px_at_200 = round(ppm_fat*gbar_H*c(1.5,3,7)/200, 2),
px_at_400 = round(ppm_fat*gbar_H*c(1.5,3,7)/400, 2)),
col.names = c("B0 (T)","Fat-water offset (Hz)","Shift at 100 Hz/px",
"at 200 Hz/px","at 400 Hz/px"),
caption = "Chemical-shift displacement in pixels. The artifact is a frequency-encode phenomenon only; it never appears along the phase-encode axis, which is a useful diagnostic clue.")| B0 (T) | Fat-water offset (Hz) | Shift at 100 Hz/px | at 200 Hz/px | at 400 Hz/px |
|---|---|---|---|---|
| 1.5 | 224 | 2.24 | 1.12 | 0.56 |
| 3.0 | 447 | 4.47 | 2.24 | 1.12 |
| 7.0 | 1043 | 10.43 | 5.22 | 2.61 |
Connection to Chapter 1. Every artifact here is a Fourier or sampling phenomenon. Gibbs ringing is the sinc point-spread function produced by truncating (windowing) k-space — the convolution theorem of Chapter 1. Aliasing is a direct Nyquist violation. Motion ghosting is a phase inconsistency in the Fourier data. Lab 2 reproduces ringing and ghosting by manipulating k-space directly, with no simulation of anatomy at all.
Checkpoint 4.7. Why do motion ghosts propagate along the phase-encode direction rather than the readout direction, and why does chemical-shift artifact do exactly the opposite?
4.8 Safety: SAR, Gradients, Projectiles, and Contrast Agents
MRI uses no ionizing radiation, but “non-ionizing” is not “without hazard”. MRI has four distinct hazard mechanisms, three of them from the three fields the scanner applies, and one from what is injected. Each has its own physics, its own limit, and its own failure mode — and the most dangerous of them has nothing to do with imaging at all.
4.8.1 RF deposition: the specific absorption rate
The \(B_1\) pulses of Section 4.1.4 induce currents in conductive tissue, and those currents dissipate as heat. The regulated quantity is the specific absorption rate, and it scales as
\[\begin{equation} \mathrm{SAR} \;\propto\; \sigma_{\rm tissue}\, B_0^{2}\, \alpha^{2}\, D\, , \tag{30} \end{equation}\]
with \(\alpha\) the flip angle and \(D\) the duty cycle. The \(B_0^2\) dependence is the central fact of high-field MRI: everything else being equal, moving from 1.5 T to 3 T quadruples RF heating, and a train of \(180^\circ\) refocusing pulses costs four times as much again as \(90^\circ\) pulses.
From the equation to the protocol. Equation (30) explains a design decision that otherwise looks arbitrary. Fast spin-echo sequences at 3 T use reduced refocusing flip angles (often \(120^\circ\) or less rather than \(180^\circ\)) and longer echo spacings. This is not a compromise on image quality for its own sake: a full \(180^\circ\) echo train at 3 T carries \(16\times\) the SAR of the same train at 1.5 T with \(90^\circ\) pulses, and would exceed the regulatory limit. The physics of Eq. (30) is written directly into the pulse sequence.
Bs <- seq(0.5, 7, by = 0.05)
sar_rel <- function(B, a) (B/1.5)^2*(a/90)^2
df_sar <- do.call(rbind, lapply(c(90, 120, 180), function(a)
data.frame(B = Bs, S = sar_rel(Bs, a),
alpha = factor(sprintf("%d deg", a), levels = sprintf("%d deg", c(90,120,180))))))
p_sar1 <- ggplot(df_sar, aes(B, S, colour = alpha)) +
geom_line(linewidth = 0.95) +
geom_hline(yintercept = 4, linetype = "dashed", colour = "grey45") +
scale_y_log10() + scale_colour_manual(values = bpad_pal[1:3]) +
labs(x = expression(B[0]*" (T)"), y = "Relative SAR (log scale)",
title = "SAR relative to a 90-degree pulse at 1.5 T")
gsar <- expand.grid(B = seq(0.5, 7, length.out = 90),
a = seq(20, 180, length.out = 90))
gsar$S <- sar_rel(gsar$B, gsar$a)
p_sar2 <- ggplot(gsar, aes(B, a, fill = log10(S))) +
geom_raster() +
geom_contour(data = gsar, aes(B, a, z = log10(S)), breaks = log10(c(1,4,16)),
colour = "white", linewidth = 0.5, inherit.aes = FALSE) +
scale_fill_viridis_c() +
theme(legend.position = "right", legend.title = element_text(size = 8)) +
labs(x = expression(B[0]*" (T)"), y = "Flip angle (degrees)",
fill = "log10 SAR", title = "Contours at 1x, 4x, and 16x")
bpad_grid(p_sar1, p_sar2, ncol = 2)
Figure 13: Radiofrequency heating scales as the square of both field strength and flip angle, Eq. (sar). Left: relative SAR across clinical field strengths for three refocusing flip angles, normalized to a 90-degree pulse at 1.5 T; the dashed line marks a four-fold increase. Right: the same information as a heat map over the field and flip-angle plane. This is why 7 T protocols are SAR-limited rather than time-limited, and why reduced-flip-angle refocusing is standard at 3 T.
knitr::kable(data.frame(
configuration = c("1.5 T, 90 deg","1.5 T, 180 deg","3 T, 90 deg",
"3 T, 180 deg","7 T, 90 deg","7 T, 180 deg"),
relative_SAR = c(sar_rel(1.5,90), sar_rel(1.5,180), sar_rel(3,90),
sar_rel(3,180), round(sar_rel(7,90),1), round(sar_rel(7,180),1))),
col.names = c("Configuration","Relative SAR"),
caption = "IEC limits are 2 W/kg whole-body and 3.2 W/kg head-averaged in normal operating mode. A 180-degree echo train at 3 T is 16 times the reference condition, which is why reduced refocusing angles are used.")| Configuration | Relative SAR |
|---|---|
| 1.5 T, 90 deg | 1.0 |
| 1.5 T, 180 deg | 4.0 |
| 3 T, 90 deg | 4.0 |
| 3 T, 180 deg | 16.0 |
| 7 T, 90 deg | 21.8 |
| 7 T, 180 deg | 87.1 |
4.8.2 Gradient switching: \(dB/dt\) and acoustic noise
Rapidly switched gradients induce electric fields in tissue by Faraday’s law. Above a threshold in \(dB/dt\) this causes peripheral nerve stimulation — an unpleasant twitching sensation — which is the practical limit on gradient slew rate in fast sequences such as echo-planar imaging. The same switching drives Lorentz forces on the gradient coils inside \(B_0\), producing the acoustic noise that can exceed 100 dB and requires hearing protection for every patient.
4.8.3 The static field: the projectile hazard
The most dangerous aspect of MRI is not electromagnetic exposure but mechanical. The static field is always on, even when no scan is running and the scanner appears inert. Ferromagnetic objects are accelerated into the bore with lethal force, and the attractive force rises steeply as an object approaches. Deaths have occurred. Zone-based access control, ferromagnetic screening, and MR-conditional equipment are not bureaucracy; they are the primary safety system.
The static field also interacts with implants: pacemakers and defibrillators (many modern devices are MR-conditional under specified conditions), cochlear implants, neurostimulators, aneurysm clips, and retained metallic foreign bodies, especially intraocular fragments. Screening every patient and every accompanying person is mandatory.
4.8.4 Contrast-agent risk
Gadolinium chelates (Section 4.11) carry two distinct concerns. Nephrogenic systemic fibrosis is a rare but serious fibrosing condition associated with linear gadolinium agents in severe renal impairment, and the risk has fallen sharply with macrocyclic agents and renal screening. Separately, gadolinium retention in brain and bone has been documented after repeated administration; the clinical significance is still uncertain, but it has driven a general principle of using the lowest effective dose and preferring macrocyclic agents.
Four hazards, four physics. The RF field heats tissue and scales as \(B_0^2\alpha^2\). The gradient field induces nerve stimulation through \(dB/dt\) and generates acoustic noise through Lorentz forces. The static field is a mechanical projectile hazard and an implant hazard, and it never turns off. The contrast agent is a pharmacological risk governed by renal function and chelate chemistry. Only the first two are limited by software; the third is limited by access control and screening, and it is the one that has killed people.
Section 4.8 summary.
- \(\mathrm{SAR} \propto B_0^2\alpha^2 D\); the \(B_0^2\) term is why 3 T and 7 T protocols are SAR-limited and use reduced refocusing angles.
- Gradient \(dB/dt\) limits slew rate through peripheral nerve stimulation; gradient Lorentz forces produce >100 dB acoustic noise.
- The static field is always on and is a projectile and implant hazard: the leading cause of serious MRI incidents.
- Gadolinium carries NSF risk in severe renal impairment and documented tissue retention.
Checkpoint 4.8. A protocol validated at 1.5 T with a \(180^\circ\) fast-spin-echo train is transferred unchanged to a 3 T scanner and immediately exceeds the SAR limit. Using Eq. (30), compute the factor by which SAR increased, and give two distinct changes that would bring it back within limits, stating the image-quality cost of each.
4.9 From Images to Measurements
Reconstruction yields an image, but quantitative science needs measurements. The standard pipeline, developed in detail in Chapter 7, proceeds
\[\text{image} \rightarrow \text{registration} \rightarrow \text{segmentation} \rightarrow \text{quantitative maps}.\]
Registration aligns images across time, sequences, or subjects into a common coordinate system; segmentation delineates tissues or lesions; and quantitative MRI fits the signal models of this chapter voxel by voxel to produce parameter maps — \(T_1\), \(T_2\), \(T_2^*\), proton density, ADC, perfusion, and magnetic susceptibility.
The distinction matters more in MRI than in any other modality. A weighted image is qualitative: its brightness depends on coil sensitivity, receiver gain, sequence timing, field strength, and vendor, none of which are recorded in the pixel value. A parameter map aims to be a reproducible physical measurement in physical units. Only the second kind is a legitimate input to the multi-site predictive models of Chapter 8, and Chapter 8’s discussion of acquisition shift is largely a discussion of what happens when this distinction is ignored.
4.10 Advanced and Functional MRI
The sequences so far image anatomy. By sensitizing the signal to motion of water, to blood oxygenation, to perfusion, or to chemistry, MRI becomes a functional and molecular tool.
4.10.1 Diffusion MRI: DWI and ADC
Diffusion-weighted imaging (DWI) adds a pair of strong gradient pulses that dephase and then rephase spins. Stationary water is perfectly refocused; water that has diffused between the two pulses is not, and loses signal. The attenuation is
\[\begin{equation} \frac{S}{S_0} = e^{-b\,\mathrm{ADC}}, \qquad b = (\gamma G \delta)^{2}\left(\Delta - \frac{\delta}{3}\right), \tag{31} \end{equation}\]
where ADC is the apparent diffusion coefficient, \(G\) and \(\delta\) are the amplitude and duration of the diffusion gradients, and \(\Delta\) their separation. Where diffusion is restricted — as in acute ischaemic stroke, where cytotoxic oedema swells cells and traps water — DWI stays bright and the ADC map is dark. DWI detects stroke within minutes, long before CT or conventional MRI.
Two traps in ADC measurement. First, “restricted diffusion appears bright on DWI” is only partly a diffusion statement: a long-\(T_2\) tissue also stays bright at long TE, an effect called \(T_2\) shine-through, and it is precisely why the ADC map — which divides out \(S_0\) — must always be read alongside the DWI image. Second, at high b-value the signal approaches the Rician noise floor of Section 4.6.2, which lifts the measured magnitude and biases the fitted ADC downward. The demonstration in Section 4.6.2 found a 19% ADC error at \(b = 4000\) s mm\(^{-2}\) with a b = 0 SNR of 20.
bg <- seq(0, 3000, by = 10)
adc <- c("Restricted (acute stroke)" = 0.5e-3, "Normal brain" = 0.8e-3,
"Free water (CSF)" = 3.0e-3)
df_d <- do.call(rbind, lapply(names(adc), function(nm)
data.frame(b = bg, S = exp(-bg*adc[nm]),
tissue = factor(nm, levels = names(adc)))))
p_d1 <- ggplot(df_d, aes(b, S, colour = tissue)) +
geom_line(linewidth = 0.9) +
scale_y_log10() + scale_colour_manual(values = bpad_pal[c(2,1,3)]) +
labs(x = expression("b-value (s/mm"^2*")"), y = "Relative signal (log scale)",
title = "DWI attenuation: slope = -ADC")
## T2 shine-through: same ADC, different T2
TE_dwi <- 90
shine <- data.frame(
tissue = c("Lesion: ADC low, T2 normal", "Lesion: ADC normal, T2 very long"),
ADC = c(0.5e-3, 0.8e-3), T2 = c(90, 250))
bb <- seq(0, 1500, by = 10)
df_s <- do.call(rbind, lapply(seq_len(nrow(shine)), function(i)
data.frame(b = bb, S = exp(-TE_dwi/shine$T2[i])*exp(-bb*shine$ADC[i]),
tissue = shine$tissue[i])))
p_d2 <- ggplot(df_s, aes(b, S, colour = tissue)) +
geom_line(linewidth = 0.9) +
geom_vline(xintercept = 1000, linetype = "dotted", colour = "grey45") +
annotate("text", x = 1020, y = 0.55, hjust = 0, size = 2.9, colour = "grey30",
label = "typical clinical b") +
scale_colour_manual(values = bpad_pal[c(2,5)]) +
labs(x = expression("b-value (s/mm"^2*")"), y = "DWI signal (includes T2 weighting)",
title = "T2 shine-through mimics restriction")
bpad_grid(p_d1, p_d2, ncol = 2)
Figure 14: Diffusion-weighted signal. Left: attenuation versus b-value, Eq. (dwi); the slope on this logarithmic axis is the negative ADC, so restricted water in acute stroke decays slowly and stays bright while free CSF decays fastest. Right: T2 shine-through, the reason a bright DWI image alone does not establish restricted diffusion. Two tissues with identical ADC but different T2 produce different DWI brightness at a typical echo time, and only the ADC map separates them.
A_dwi <- exp(-TE_dwi/90)*exp(-1000*0.5e-3) # RESTRICTED diffusion, normal T2
B_dwi <- exp(-TE_dwi/250)*exp(-1000*0.8e-3) # NORMAL diffusion, long T2
cat(sprintf("At b = 1000 and TE = %d ms:\n", TE_dwi))## At b = 1000 and TE = 90 ms:## lesion A (ADC 0.5e-3, RESTRICTED, T2 = 90 ms): DWI = 0.2231## lesion B (ADC 0.8e-3, NORMAL, T2 = 250 ms): DWI = 0.3135cat(sprintf("Lesion B is %.0f%% BRIGHTER despite having 60%% HIGHER (normal) diffusivity.\n",
100*(B_dwi/A_dwi - 1)))## Lesion B is 40% BRIGHTER despite having 60% HIGHER (normal) diffusivity.## A reader looking only at the DWI image would call B the acute infarct and miss A.## The ADC map, which divides out S0 and hence the T2 weighting, reverses the ranking.4.10.2 Diffusion-tensor imaging
In white matter, diffusion is anisotropic: water moves more easily along axons than across them. A single scalar ADC is then insufficient, and diffusion is described by a \(3\times3\) symmetric positive-definite diffusion tensor \(\mathbf D\), so that
\[\begin{equation} \frac{S(\hat g)}{S_0} = e^{-b\,\hat g^{\mathsf T}\mathbf D\,\hat g} \tag{32} \end{equation}\]
for a diffusion direction \(\hat g\). Because \(\mathbf D\) is symmetric it has six independent entries, which is why DTI requires DWI measurements along at least six non-collinear directions plus one \(b = 0\) image. Equation (32) is fitted by least squares across directions, and the eigen-decomposition then gives eigenvalues \(\lambda_1 \ge \lambda_2 \ge \lambda_3\) (diffusivities along principal axes) and eigenvectors (their directions). Two rotation-invariant scalars summarize the tensor:
\[\begin{equation} \mathrm{MD} = \frac{\lambda_1+\lambda_2+\lambda_3}{3}, \qquad \mathrm{FA} = \sqrt{\tfrac{3}{2}}\, \frac{\sqrt{\sum_i(\lambda_i - \bar\lambda)^2}}{\sqrt{\sum_i\lambda_i^2}} . \tag{33} \end{equation}\]
Fractional anisotropy runs from 0 (isotropic, as in CSF) to 1 (perfectly linear diffusion). The principal eigenvector points along the local fibre direction, enabling tractography.
Connection to Chapter 1. DTI is the eigenvalue/eigenvector problem of Chapter 1, applied per voxel to a real symmetric tensor field. FA and MD are rotation-invariant functions of the eigenvalues, which is exactly why they are reportable quantities: they do not depend on how the patient was oriented in the scanner. Lab 3 builds an FA map by eigen-decomposition and verifies the rotation invariance explicitly.
4.10.3 Functional MRI: the BOLD effect
Functional MRI exploits a fortunate coincidence: deoxyhaemoglobin is paramagnetic and shortens \(T_2^*\), whereas oxyhaemoglobin is essentially diamagnetic and does not. When neurons fire, local blood flow over-compensates, flushing in oxygenated blood and lowering deoxyhaemoglobin, so \(T_2^*\) lengthens and the \(T_2^*\)-weighted gradient-echo signal rises.
The effect is small: typically 1–3% at 1.5 T and 2–5% at 3 T, against physiological and thermal noise of comparable size. This is why fMRI requires many repetitions and careful statistical modelling, and why it is one of the clearest examples in this book of a measurement that exists only after inference (Chapter 8).
tg <- seq(0, 30, by = 0.1)
hrf <- function(t) (t^5*exp(-t)/factorial(5)) - 0.35*(t^9*exp(-t)/factorial(9))
h <- hrf(tg); h <- h/max(h)
p_h <- ggplot(data.frame(t = tg, B = h), aes(t, B)) +
geom_hline(yintercept = 0, colour = "grey70") +
geom_line(colour = bpad_pal[2], linewidth = 1) +
geom_vline(xintercept = tg[which.max(h)], linetype = "dotted", colour = "grey45") +
annotate("text", x = tg[which.max(h)] + 0.6, y = 0.85, hjust = 0, size = 3,
colour = "grey30", label = sprintf("peak at %.1f s", tg[which.max(h)])) +
labs(x = "Time after neural event (s)", y = "Relative BOLD signal",
title = "Canonical haemodynamic response")
set.seed(9)
TRf <- 2; nvol <- 200; tv <- (0:(nvol-1))*TRf
block <- as.numeric((tv %% 60) < 30)
kern <- hrf(seq(0, 20, by = TRf)); kern <- kern/sum(kern)
conv <- as.numeric(stats::filter(block, kern, sides = 1))
conv[is.na(conv)] <- 0
conv <- conv/max(conv)
effect <- 0.02
sig_b <- 100*(1 + effect*conv) + rnorm(nvol, 0, 100*0.01) +
0.5*sin(2*pi*tv/90)
df_b <- rbind(data.frame(t = tv, y = sig_b, q = "measured voxel time course"),
data.frame(t = tv, y = 100*(1 + effect*conv), q = "true BOLD response"))
p_b <- ggplot(df_b, aes(t, y, colour = q)) +
geom_line(linewidth = 0.7) +
scale_colour_manual(values = c("measured voxel time course" = bpad_pal[7],
"true BOLD response" = bpad_pal[2])) +
labs(x = "Time (s)", y = "Signal (arbitrary units)",
title = "A 2% effect buried in noise")
bpad_grid(p_h, p_b, ncol = 2)
Figure 15: The BOLD response and why fMRI is a statistics problem. Left: the canonical haemodynamic response to a brief neural event, rising to a peak near 5 seconds and followed by a post-stimulus undershoot; note that the response is an order of magnitude slower than the neural activity that caused it. Right: a simulated block-design time course at a realistic 2 percent effect size with physiological noise. The signal is invisible in a single voxel time course and emerges only after averaging across repetitions, which is the entire methodological basis of the field.
fit <- lm(sig_b ~ conv)
cat(sprintf("single-voxel regression: estimated effect %.3f%% (true 2.00%%), p = %.4f\n",
100*coef(fit)[2]/mean(sig_b), summary(fit)$coefficients[2,4]))## single-voxel regression: estimated effect 2.140% (true 2.00%), p = 0.0000cat(sprintf("noise sd is %.1f%% of the mean, so the effect-to-noise ratio per volume is %.2f\n",
100*sd(residuals(fit))/mean(sig_b), effect/(sd(residuals(fit))/mean(sig_b))))## noise sd is 1.0% of the mean, so the effect-to-noise ratio per volume is 1.95Connection to Chapter 2. BOLD fMRI measures the same neurovascular haemodynamic response as functional NIRS in Chapter 2 — a rise in oxy- and fall in deoxyhaemoglobin after activation — but reads it out magnetically (through \(T_2^*\)) rather than optically (through absorption). The two methods cross-validate each other, trading fMRI’s whole-brain coverage and depth against fNIRS’s portability and motion tolerance.
4.10.4 Perfusion and spectroscopy
Perfusion MRI measures tissue blood supply. Dynamic susceptibility contrast (DSC) tracks a gadolinium bolus through its transient \(T_2^*\) drop; dynamic contrast-enhanced (DCE) imaging models the \(T_1\)-shortening leakage of gadolinium to estimate vascular permeability (\(K^{\text{trans}}\)); and arterial spin labelling (ASL) magnetically tags inflowing blood as an endogenous tracer, requiring no injection at all.
MR spectroscopy (MRS) forgoes imaging to resolve the chemical-shift spectrum within a voxel, quantifying metabolites through the chemical shift of Eq. (29): N-acetylaspartate (a neuronal marker), creatine, choline (membrane turnover), and lactate (anaerobic metabolism). A tumour classically shows rising choline and falling NAA.
Connection to Chapters 2 and 3. MRS recovers its spectrum exactly as FTIR did in Chapter 2: the time-domain free-induction decay and the chemical spectrum are a Fourier-transform pair. The metabolite peaks of MRS are the magnetic analogue of the vibrational peaks of infrared spectroscopy, and the linewidth analysis of Section 4.2.5 is why shimming determines whether two metabolites can be resolved at all.
Section 4.10 summary.
- DWI: \(S/S_0 = e^{-b\,\mathrm{ADC}}\); restricted diffusion is bright on DWI and dark on ADC, but \(T_2\) shine-through means the ADC map must always be read too.
- The Rician noise floor biases high-b ADC downward, by nearly 20% at \(b = 4000\) and moderate SNR.
- DTI needs \(\ge6\) directions because \(\mathbf D\) has six independent entries; FA and MD are rotation-invariant eigenvalue functions.
- BOLD is a 1–5% \(T_2^*\) effect and is detectable only through statistical modelling.
- Perfusion (DSC, DCE, ASL) maps blood supply; MRS resolves chemistry by Fourier transform.
Checkpoint 4.9. Acute stroke is bright on DWI and dark on the ADC map. Explain both observations from Eq. (31), then explain how you would distinguish true restriction from \(T_2\) shine-through using only the two images available.
4.11 Contrast Agents and Angiography
Gadolinium chelates are paramagnetic agents that strongly shorten \(T_1\) wherever they accumulate, brightening those regions on T1-weighted images. Enhancement marks increased vascularity or a disrupted blood–brain barrier, highlighting many tumours, active inflammation, and infection. The relaxivity relation is
\[\begin{equation} \frac{1}{T_1^{\rm obs}} = \frac{1}{T_1^{0}} + r_1\,[\mathrm{Gd}] , \tag{34} \end{equation}\]
with \(r_1\) the agent’s longitudinal relaxivity, typically 3–5 mM\(^{-1}\)s\(^{-1}\) at clinical field strengths. Section 4.8.4 covers the associated risks.
Vessels can be imaged with or without contrast. Time-of-flight (TOF) MRA relies on inflow: stationary tissue is saturated by rapid repeated RF pulses, while fresh blood entering the slice carries full magnetization and appears bright — no injection required. Phase-contrast MRA encodes velocity into signal phase and can quantify flow, using the same complex-phase information introduced in Section 4.1.3. Contrast-enhanced MRA uses a gadolinium bolus for fast, high-quality vascular maps.
4.12 Clinical and Disease Applications
4.12.1 Stroke and cerebrovascular disease
Diffusion-weighted imaging is the front-line test for acute ischaemic stroke: cytotoxic oedema restricts water diffusion within minutes, lighting up on DWI and darkening the ADC map before any change is visible on CT. Combining DWI with perfusion imaging reveals the diffusion–perfusion mismatch: tissue that is underperfused but not yet infarcted, the penumbra, which guides thrombolysis and thrombectomy decisions.
Clinical example. A bright DWI lesion with a dark ADC and a larger surrounding perfusion deficit signals salvageable penumbra — the imaging rationale for emergent reperfusion therapy. The ADC map is essential here rather than optional: a bright DWI with a normal ADC indicates \(T_2\) shine-through, not acute infarction, and would lead to exactly the wrong decision.
4.12.2 Neurological and neuro-oncological disease
Multiple sclerosis lesions are revealed by FLAIR with CSF nulled (Section 4.3.2); neurodegeneration is tracked by volumetric atrophy measurement; epilepsy evaluation uses high-resolution structural imaging to find subtle malformations. Brain tumours are characterized multiparametrically: T1 plus gadolinium for enhancement, T2/FLAIR for oedema, DWI for cellularity, perfusion for vascularity, and MRS for the choline/NAA signature. DTI tractography maps eloquent white-matter tracts and fMRI localizes motor and language cortex, both used to plan maximal safe resection.
4.12.3 Oncology and body imaging
Beyond the brain, DCE perfusion characterizes breast, prostate, and liver lesions through contrast kinetics; whole-body DWI detects and stages metastatic disease; and multiparametric prostate MRI has reshaped prostate-cancer diagnosis. In cardiology, late gadolinium enhancement maps myocardial scar, and phase-contrast imaging quantifies flow and shunts.
4.13 Choosing a Sequence
From the clinical question to the acquisition.
| Question | Sequence | Physical property | Binding limit |
|---|---|---|---|
| What is the anatomy? | T1-weighted SE/GRE | \(T_1\) | contrast reverses with timing |
| Where is the pathology? | T2-weighted / FLAIR | \(T_2\), \(T_1\) (nulling) | FLAIR TI needs the finite-TR correction |
| Is there acute infarction? | DWI + ADC map | water mobility | \(T_2\) shine-through; Rician floor at high b |
| Is there salvageable brain? | DWI + perfusion | mismatch | perfusion needs contrast or ASL |
| Where do the tracts run? | DTI, \(\ge6\) directions | tensor eigenvectors | crossing fibres; low SNR |
| Which cortex is eloquent? | BOLD fMRI | \(T_2^*\) | 1–5% effect; needs statistics |
| Is there old haemorrhage? | GRE / susceptibility-weighted | \(T_2^*\) | same physics that causes artifact |
| Is the barrier broken? | T1 + gadolinium | \(r_1\) relaxivity | renal function, retention |
| Is there flow, and how fast? | TOF or phase-contrast MRA | inflow; phase | velocity aliasing |
| What is the tissue made of? | MRS | chemical shift | shimming sets linewidth |
| Patient has an implant | screening first, then MR-conditional protocol | static field | not a contrast question at all |
Sequence \(\rightarrow\) property \(\rightarrow\) disease. DWI/ADC \(\rightarrow\) water mobility \(\rightarrow\) acute stroke and tumour cellularity. DTI \(\rightarrow\) fibre orientation \(\rightarrow\) tract mapping. BOLD \(\rightarrow\) neural activity \(\rightarrow\) functional localization. Perfusion \(\rightarrow\) blood supply \(\rightarrow\) penumbra and tumour grade. T1 plus gadolinium \(\rightarrow\) barrier breakdown \(\rightarrow\) tumour, inflammation, scar. \(T_2^*\) \(\rightarrow\) susceptibility \(\rightarrow\) haemorrhage.
4.14 Computational Laboratories
The three laboratories build the signature computational skills of MRI. Each is self-contained and runnable.
Lab 1: k-space, reconstruction, and true undersampling
Goal. Build a numerical head phantom, transform it to k-space, and demonstrate that the image is the inverse Fourier transform of k-space; that the centre of k-space carries contrast while the periphery carries detail; and that undersampling halves the field of view and wraps the object. The last point is done two ways, because the usual shortcut and the physically correct simulation differ in an instructive manner.
N <- 128
ax <- seq(-1, 1, length.out = N)
X <- outer(rep(1, N), ax) # x varies across columns
Y <- outer(ax, rep(1, N)) # y varies across rows
add_ellipse <- function(M, cx, cy, a, b, val) {
M[((X - cx)/a)^2 + ((Y - cy)/b)^2 <= 1] <- val; M
}
phantom <- matrix(0, N, N)
phantom <- add_ellipse(phantom, 0.00, 0.00, 0.82, 0.98, 0.50) # skull
phantom <- add_ellipse(phantom, 0.00, 0.00, 0.72, 0.90, 0.32) # brain
phantom <- add_ellipse(phantom, -0.22, 0.12, 0.12, 0.22, 0.60) # ventricle L
phantom <- add_ellipse(phantom, 0.22, 0.12, 0.12, 0.22, 0.60) # ventricle R
phantom <- add_ellipse(phantom, 0.00, -0.35, 0.10, 0.10, 0.90) # lesion
## Centred FFT shift (its own inverse for even N)
fftshift2 <- function(M) {
n <- nrow(M); m <- ncol(M); h <- n %/% 2; k <- m %/% 2
M[c((h+1):n, 1:h), c((k+1):m, 1:k)]
}
K <- fft(phantom) # k-space, unshifted (DC at [1,1])
Kc <- fftshift2(K) # centred, for masking and display
Kdisp <- log(1 + Mod(Kc))
## Low- versus high-frequency reconstruction
cen <- N/2 + 1; r <- 12
keep <- matrix(FALSE, N, N)
keep[(cen-r):(cen+r), (cen-r):(cen+r)] <- TRUE
img_low <- Mod(fft(fftshift2(Kc * keep), inverse = TRUE))/(N*N)
img_high <- Mod(fft(fftshift2(Kc * !keep), inverse = TRUE))/(N*N)
## (a) The usual shortcut: ZERO every other phase-encode line and keep an N x N matrix.
## Correct wrap geometry, but each replica appears at HALF amplitude, because
## zeroing is multiplication by a comb of mean 1/2.
Kzero <- K; Kzero[seq(2, N, 2), ] <- 0
img_zero <- Mod(fft(Kzero, inverse = TRUE))/(N*N)
## (b) True undersampling: ACQUIRE only every other line and reconstruct into the
## half-height matrix the scanner would actually produce. Delta_k doubles, so
## the FOV halves and the replicas overlap at FULL amplitude.
Ksub <- K[seq(1, N, 2), ]
img_true <- Mod(fft(Ksub, inverse = TRUE))/((N/2)*N)
bpad_grid(
plot_img(phantom, "Phantom (image)"),
plot_img(Kdisp, "k-space magnitude (log)", palette = "viridis"),
plot_img(img_low, "Centre only: contrast, no edges"),
plot_img(img_high, "Periphery only: edges only"),
plot_img(img_zero, "Zero-filled alternate lines"),
plot_img(img_true, "True undersampling (N/2 rows)"),
ncol = 3)
Figure 16: k-space and reconstruction. Top row: the phantom, its k-space magnitude on a logarithmic scale, and the reconstruction from the central 25 by 25 region only, which retains all the contrast but no edges. Bottom row: the reconstruction from the periphery only, which is an edge map; the zero-filled undersampling shortcut, which produces the correct wrap geometry at half amplitude; and true undersampling, in which only every second phase-encode line is acquired and reconstructed into a half-height matrix, which is what a scanner actually produces.
## phantom peak = 0.900## zero-filled shortcut peak = 0.610 (replicas at half amplitude)cat(sprintf("true undersampling peak = %.3f in a %d x %d matrix (full amplitude)\n",
max(img_true), nrow(img_true), ncol(img_true)))## true undersampling peak = 1.220 in a 64 x 128 matrix (full amplitude)cat(sprintf("energy retained in the central %d x %d of k-space: %.1f%%\n",
2*r+1, 2*r+1, 100*sum(Mod(Kc[keep])^2)/sum(Mod(Kc)^2)))## energy retained in the central 25 x 25 of k-space: 96.3%## Note that a tiny fraction of k-space carries almost all the energy, which is## exactly the sparsity that compressed sensing exploits (Section 4.6.3).Discussion. Reconstruction is literally the inverse 2D FFT of k-space, Eq. (21). The central lobe dominates the energy and sets broad contrast; the high-frequency periphery encodes sharp boundaries, and an image built from it alone is an edge map with no anatomy. Both undersampling panels show the same wrap geometry — the Nyquist violation of Eq. (22) made visible — but only the second is what a scanner produces: acquiring half the lines doubles \(\Delta k\), halves the FOV, and folds the object onto itself at full amplitude in a half-height matrix.
Try it yourself. Change r from 12 to 4 and to 40 and watch contrast and edges trade
against each other. Then undersample by 4 instead of 2 and confirm that the object wraps
twice. Finally, zero the columns instead of the rows and observe that the wrap now occurs
along the other axis — a reminder that “phase-encode direction” is a choice the operator
makes, and one of the few artifact controls available at the console.
Lab 2: Truncation ringing and motion ghosting
Goal. Produce two classic MRI artifacts purely by manipulating k-space, reinforcing that artifacts are Fourier-domain phenomena rather than properties of the anatomy.
## Gibbs ringing: keep only a small central k-space window
w <- 14
trunc <- matrix(FALSE, N, N)
trunc[(cen-w):(cen+w), (cen-w):(cen+w)] <- TRUE
img_gibbs <- Mod(fft(fftshift2(Kc * trunc), inverse = TRUE))/(N*N)
## Motion ghosting: periodic phase error varying by phase-encode line (row)
jj <- 1:N
ghost_img <- function(amp) {
ph <- amp*sin(2*pi*6*jj/N)
Mod(fft(K*exp(1i*outer(ph, rep(1, N))), inverse = TRUE))/(N*N)
}
img_motion <- ghost_img(1.2)
## Profile through the lesion showing the ringing explicitly
row_id <- which.min(abs(ax - (-0.35)))
prof <- rbind(
data.frame(x = seq_len(N), v = phantom[row_id, ], trace = "original"),
data.frame(x = seq_len(N), v = img_gibbs[row_id, ], trace = "truncated k-space"))
p_prof <- ggplot(prof, aes(x, v, colour = trace)) +
geom_line(linewidth = 0.8) +
scale_colour_manual(values = c("original" = bpad_pal[7],
"truncated k-space" = bpad_pal[2])) +
labs(x = "Column index", y = "Intensity", title = "Gibbs ringing profile")
## How ghost energy grows with the motion amplitude
amps <- seq(0, 2, by = 0.1)
ghost_frac <- vapply(amps, function(a) {
im <- ghost_img(a)
mask <- phantom > 0
sum(im[!mask])/sum(im) # fraction of energy outside the true object
}, numeric(1))
p_gf <- ggplot(data.frame(a = amps, f = 100*ghost_frac), aes(a, f)) +
geom_line(linewidth = 0.95, colour = bpad_pal[1]) +
geom_point(size = 1.6, colour = bpad_pal[1]) +
labs(x = "Motion phase amplitude (radians)",
y = "Energy outside the object (%)",
title = "Ghosting displaces signal from the object")
bpad_grid(
plot_img(phantom, "Phantom (reference)"),
plot_img(img_gibbs, "Truncated k-space: Gibbs ringing"),
p_prof,
plot_img(img_motion, "Phase error: motion ghosting"),
p_gf,
plot_img(ghost_img(2.0), "Larger motion amplitude"),
ncol = 3)
Figure 17: Artifacts made in k-space. Top: the phantom, the result of truncating k-space to a small central window, and an intensity profile through the lesion showing the characteristic overshoot and decaying oscillation of Gibbs ringing beside a sharp edge. Bottom: a periodic phase error across phase-encode lines produces discrete ghost copies displaced along the phase-encode axis, and increasing the amplitude of that error moves signal out of the true object and into the ghosts.
cat(sprintf("Truncating to a %d x %d window keeps %.2f%% of k-space samples\n",
2*w+1, 2*w+1, 100*sum(trunc)/(N*N)))## Truncating to a 29 x 29 window keeps 5.13% of k-space samplescat(sprintf("and %.1f%% of the energy, yet the image is visibly degraded by ringing.\n",
100*sum(Mod(Kc[trunc])^2)/sum(Mod(Kc)^2)))## and 97.0% of the energy, yet the image is visibly degraded by ringing.cat(sprintf("At a 1.2 rad motion amplitude, %.1f%% of image energy sits outside the object.\n",
100*ghost_frac[which.min(abs(amps - 1.2))]))## At a 1.2 rad motion amplitude, 7.9% of image energy sits outside the object.Discussion. Truncating k-space is multiplication by a box window, which in the image domain is convolution with a sinc (Chapter 1), hence the overshoot and decaying fringes beside every high-contrast edge — visible directly in the profile panel. The periodic phase error spreads each point into evenly spaced ghosts along the phase-encode direction, which is exactly the structure of real respiratory and pulsatile artifacts, and the reason gating and navigators target that axis. Note that the ghost energy is removed from the object: ghosting does not merely add spurious signal, it degrades the true anatomy too.
Lab 3: Diffusion-tensor FA, and why it is rotation-invariant
Goal. Compute fractional anisotropy from diffusion tensors using the Chapter 1 eigenvalue decomposition; verify that FA and MD are invariant to patient orientation; and build a synthetic FA map of a white-matter-like tract.
fa_md <- function(D) {
ev <- eigen(D, symmetric = TRUE)$values
lb <- mean(ev)
c(FA = sqrt(1.5)*sqrt(sum((ev - lb)^2))/sqrt(sum(ev^2)), MD = lb)
}
D_csf <- diag(c(3.00, 3.00, 3.00))*1e-3 # isotropic, free water
D_gm <- diag(c(0.95, 0.90, 0.85))*1e-3 # nearly isotropic
D_wm <- diag(c(1.70, 0.30, 0.30))*1e-3 # strongly anisotropic
## Rotation invariance: rotate the white-matter tensor arbitrarily and recompute
rot <- function(a, b, g) {
Rz <- matrix(c(cos(a),-sin(a),0, sin(a),cos(a),0, 0,0,1), 3, 3, byrow = TRUE)
Ry <- matrix(c(cos(b),0,sin(b), 0,1,0, -sin(b),0,cos(b)), 3, 3, byrow = TRUE)
Rx <- matrix(c(1,0,0, 0,cos(g),-sin(g), 0,sin(g),cos(g)), 3, 3, byrow = TRUE)
Rz %*% Ry %*% Rx
}
set.seed(4); R1 <- rot(0.7, -1.1, 0.35)
D_wm_rot <- R1 %*% D_wm %*% t(R1)
knitr::kable(rbind(CSF = fa_md(D_csf), `Gray matter` = fa_md(D_gm),
`White matter` = fa_md(D_wm),
`White matter, rotated` = fa_md(D_wm_rot)),
digits = 4,
caption = "FA and mean diffusivity (mm^2/s) for representative tensors. The rotated white-matter tensor has completely different matrix entries but identical FA and MD, because both are functions of the eigenvalues alone.")| FA | MD | |
|---|---|---|
| CSF | 0.0000 | 3e-03 |
| Gray matter | 0.0555 | 9e-04 |
| White matter | 0.7990 | 8e-04 |
| White matter, rotated | 0.7990 | 8e-04 |
cat(sprintf("FA before rotation %.6f, after rotation %.6f, difference %.2e\n",
fa_md(D_wm)["FA"], fa_md(D_wm_rot)["FA"],
abs(fa_md(D_wm)["FA"] - fa_md(D_wm_rot)["FA"])))## FA before rotation 0.799022, after rotation 0.799022, difference 3.33e-16## This invariance is why FA is a reportable biomarker: it does not depend on how the## patient happened to be lying in the scanner.## Synthetic FA map with a curved tract
n <- 96
gx <- seq(-1, 1, length.out = n)
band_center <- function(x) 0.35*sin(2.2*x)
D0 <- 0.9e-3
FA_map <- matrix(0, n, n); ang_map <- matrix(0, n, n)
for (i in 1:n) for (jx in 1:n) {
yy <- gx[i]; xx <- gx[jx]
p <- 0.9*exp(-((yy - band_center(xx))/0.12)^2)
l1 <- D0*(1 + 2*p); l2 <- D0*(1 - p); l3 <- D0*(1 - p)
lb <- (l1 + l2 + l3)/3
FA_map[i, jx] <- sqrt(1.5)*sqrt((l1-lb)^2 + (l2-lb)^2 + (l3-lb)^2)/
sqrt(l1^2 + l2^2 + l3^2)
ang_map[i, jx] <- atan(0.35*2.2*cos(2.2*xx)) # local tract tangent
}
p_fa <- plot_img(FA_map, "Synthetic FA map", palette = "viridis")
step <- 6
idx <- seq(3, n, by = step)
vec <- expand.grid(i = idx, j = idx)
vec$FA <- FA_map[cbind(vec$i, vec$j)]
vec$ang <- ang_map[cbind(vec$i, vec$j)]
vec$len <- 2.6*pmax(vec$FA, 0.05)
p_vec <- ggplot(vec, aes(j, -i)) +
geom_segment(aes(xend = j + len*cos(vec$ang), yend = -i + len*sin(vec$ang),
colour = FA), linewidth = 0.6) +
scale_colour_viridis_c() + coord_equal() +
theme_void() +
theme(plot.title = element_text(face = "bold", size = 10),
legend.position = "right", legend.title = element_text(size = 8)) +
labs(title = "Principal eigenvector field", colour = "FA")
bpad_grid(p_fa, p_vec, ncol = 2)
Figure 18: Diffusion-tensor analysis. Left: a synthetic fractional-anisotropy map in which an anisotropic curved tract stands out against a near-isotropic background, mirroring how DTI highlights white-matter pathways. Right: the principal eigenvector field over the same region, whose direction is the seed for tractography; note that the vectors align along the tract precisely where FA is high, and are essentially arbitrary where FA approaches zero, which is why tractography fails in isotropic tissue.
Discussion. The table shows the discriminating power of FA: isotropic CSF gives \(\mathrm{FA} \approx 0\), grey matter is low, and white matter is high — all at broadly similar mean diffusivity, so it is FA and not MD that reveals microstructural organization. The rotated tensor demonstrates the property that makes FA usable as a biomarker at all: its matrix entries change completely under rotation while FA and MD do not, because both are functions of the eigenvalues, which are rotation invariants.
Try it yourself. Construct a tensor with \(\lambda = (1.0, 1.0, 0.1)\times10^{-3}\), a planar rather than linear diffusion profile, and compute its FA. You should find a substantial FA even though there is no single dominant direction — which is the central limitation of DTI in regions of crossing fibres, and the motivation for higher-order models such as constrained spherical deconvolution.
4.15 Conclusions and Discussion
Magnetic resonance imaging is the most versatile instrument in this book. From a single physical phenomenon — nuclear spins precessing in a magnetic field — it builds anatomical images of every organ, movies of the beating heart, maps of water diffusion that detect stroke within minutes, reconstructions of white-matter pathways, maps of brain function, and spectra of tissue chemistry. The unifying insight is that MRI cleanly separates contrast from position: relaxation and pulse timing decide what is bright, while gradients and the Fourier transform decide where it is. Because those two levers are independent, the same hardware and the same protons can be re-tasked by reprogramming the pulse sequence to answer questions that in other modalities would require entirely different machines.
Three themes recur.
MRI is Chapter 1 made physical. The Bloch equations are an ODE system; relaxation is exponential; the detected magnetization is a complex phasor measured in quadrature; k-space and the image are a two-dimensional Fourier pair; field of view, aliasing, and Gibbs ringing are sampling-and-windowing phenomena; magnitude noise is Rician; and the diffusion tensor is analyzed by eigen-decomposition. A student who sees these connections can reason about MRI rather than memorize it: why higher field gives more signal, why undersampling wraps the image, why truncation rings, why DTI needs at least six directions, and why FA is reportable but the raw tensor entries are not.
Every convenient formula has a validity condition. \(TI = T_1\ln2\) nulls a tissue only when \(TR \gg T_1\), and fails for FLAIR by 500 ms at 3 T. “Halving the voxel costs a factor of 8 in SNR” holds only at fixed sample count, and the realistic answer is a factor of 2. The magnitude spectrum of an FID is \(\sqrt3\) wider than the absorption-mode Lorentzian. A bright DWI does not establish restricted diffusion, because of \(T_2\) shine-through. The habit of asking “under what conditions is this true?” is worth more than any individual formula.
Non-ionizing is not the same as hazard-free. SAR scales as \(B_0^2\alpha^2\), gradient switching stimulates nerves and generates over 100 dB of acoustic noise, gadolinium carries renal and retention concerns, and the always-on static field remains the leading cause of serious MRI incidents. These are physics problems with physical limits, and they belong in a physics chapter rather than in a separate safety appendix.
The field continues to move quickly. Higher field strengths push signal and resolution while making SAR and \(B_1\) inhomogeneity the binding constraints. Quantitative MRI replaces qualitative weighted images with reproducible parameter maps. Deep-learning reconstruction now accelerates acquisition beyond classical compressed sensing, with the generalization and hallucination risks treated in Chapter 8. Hybrid PET/MRI fuses molecular and structural contrast with the nuclear methods of Chapter 6. And low-field, portable, point-of-care MRI is beginning to bring the modality to the bedside — trading SNR for access, and reopening every trade-off in this chapter at a different operating point.
Key points
- MRI detects precessing \(^1\)H magnetization; \(\omega_0 = \gamma B_0\), and only about 10 protons per million contribute at 3 T.
- \(M_0 \propto N\gamma^2\hbar^2B_0/4k_BT\): signal grows with field, which is the whole argument for 3 T and 7 T.
- In the rotating frame on resonance, \(B_0\) vanishes and \(\alpha = \gamma B_1\tau_p\).
- The Bloch equations give \(M_z(t) = M_0 + [M_z(0)-M_0]e^{-t/T_1}\) and \(M_{xy}(t) = M_{xy}(0)e^{-t/T_2}\); \(T_2^* \le T_2 \le T_1\).
- \(T_1\) rises substantially with field while \(T_2\) does not, so protocols must be retuned.
- \(S \propto \rho(1-e^{-TR/T_1})e^{-TE/T_2}\); grey–white contrast reverses sign between T1- and T2-weighting.
- Inversion nulling is \(T_1\ln2\) only when \(TR \gg T_1\); FLAIR needs the finite-TR form.
- Gradients encode position; \(\mathrm{FOV} = 1/\Delta k\) and \(\Delta x = 1/(2k_{\max})\); reconstruction is the inverse FFT.
- Spin echo refocuses static dephasing; gradient echo is fast, \(T_2^*\)-weighted, and optimized at the Ernst angle \(\cos\alpha_E = e^{-TR/T_1}\).
- \(\mathrm{SNR}\propto V_{\rm voxel}\sqrt{T_{\rm acq}}\); magnitude noise is Rician with a floor \(\sigma\sqrt{\pi/2}\) that biases high-b ADC.
- Chemical-shift displacement is \(\Delta f/\mathrm{BW}_{\rm px}\) pixels and grows with field.
- \(\mathrm{SAR}\propto B_0^2\alpha^2D\); the static field is a permanent projectile hazard.
- DWI/ADC detect restricted diffusion; DTI/FA come from tensor eigen-decomposition and are rotation invariant; BOLD is a 1–5% \(T_2^*\) effect requiring statistical inference.
Connections to other BPAD chapters
- Chapter 1 (Mathematical and Statistical Foundations). ODE systems (Bloch); exponentials (\(T_1\), \(T_2\), diffusion); complex phasors and quadrature detection; rotating frames; the Boltzmann factor; the 2D Fourier transform (k-space); Nyquist sampling (FOV, aliasing); windowing and convolution (Gibbs); Rician noise; and eigen-decomposition (DTI).
- Chapter 2 (Optical and Thermal Methods). BOLD fMRI measures the same neurovascular haemodynamic response as functional NIRS, read out magnetically rather than optically; MRS recovers a spectrum from a free-induction decay exactly as FTIR does from an interferogram; and the quarter-wave and resonance reasoning of Chapter 2 reappears in RF coil design.
- Chapter 3 (Ultrasound and Photoacoustics). A complementary non-ionizing modality. Both are sampling-limited: pulsed Doppler obeys \(v_{\max}d_{\max}\le c^2/8f_0\) and MRI obeys \(\mathrm{FOV} = 1/\Delta k\), and both trade acquisition time against the extent of the measurement. MRI trades real-time speed for unmatched soft-tissue contrast.
- Chapter 5 (X-ray Imaging and CT). Both reconstruct images from non-pixel measurements — projections in CT, spatial frequencies in MRI — and both are inverse problems. The contrast is dose: CT must manage ionizing radiation, while MRI must manage RF heating and the static field.
- Chapter 6 (Nuclear Medicine). PET/MRI fuses molecular and structural contrast. Note the sensitivity contrast: nuclear imaging counts individual decays, while MRI images a \(10^{-5}\) population excess — which is why PET is picomolar-sensitive and MRI is not.
- Chapter 7 (General Medical Image Processing). Registration, segmentation, bias-field correction, and quantitative parameter mapping all operate on the images produced here.
- Chapter 8 (Data Modeling, AI, and Machine Learning). Quantitative maps are candidate imaging biomarkers; deep-learning reconstruction is an inverse problem with generalization risk; and the difference between a weighted image and a parameter map is precisely the acquisition-shift problem that chapter treats at length.
Glossary
| Term | Meaning |
|---|---|
| ADC | Apparent diffusion coefficient, from \(S/S_0 = e^{-b\,\mathrm{ADC}}\) |
| ASL | Arterial spin labelling; magnetically tagged blood as an endogenous tracer |
| Bloch equations | ODE system combining precession with \(T_1\) and \(T_2\) relaxation |
| BOLD | Blood-oxygen-level-dependent contrast; a \(T_2^*\) effect of 1–5% |
| b-value | Diffusion weighting, \(b = (\gamma G\delta)^2(\Delta - \delta/3)\) |
| Chemical shift | Fat–water resonance offset of 3.5 ppm; displaces fat along the readout axis |
| DTI | Diffusion-tensor imaging; requires \(\ge6\) directions |
| Ernst angle | Flip angle maximizing steady-state GRE signal, \(\cos\alpha_E = e^{-TR/T_1}\) |
| FA | Fractional anisotropy; a rotation-invariant function of tensor eigenvalues |
| FID | Free-induction decay; a damped complex oscillation |
| FLAIR | Fluid-attenuated inversion recovery; nulls CSF |
| Flip angle | \(\alpha = \gamma B_1\tau_p\); tip produced by an RF pulse |
| g-factor | Spatially varying noise amplification in parallel imaging |
| Gibbs ringing | Sinc-convolution artifact from truncated k-space |
| Gyromagnetic ratio | \(\gamma\); 42.577 MHz/T for \(^1\)H |
| k-space | Spatial-frequency domain; the Fourier transform of the image |
| Larmor frequency | \(\omega_0 = \gamma B_0\) |
| MD | Mean diffusivity, the average of the tensor eigenvalues |
| MRS | MR spectroscopy; chemical-shift spectrum from a voxel |
| NSA | Number of signal averages; SNR grows as \(\sqrt{\mathrm{NSA}}\) |
| Rician | Distribution of magnitude data; biased upward at low SNR |
| Rotating frame | Coordinate frame at \(\omega_{\rm RF}\) in which \(B_0\) vanishes on resonance |
| SAR | Specific absorption rate; scales as \(B_0^2\alpha^2D\) |
| Shimming | Correction of \(B_0\) inhomogeneity; lengthens \(T_2^*\) |
| Spin echo | \(90^\circ\)–\(180^\circ\) sequence refocusing static dephasing |
| STIR | Short-tau inversion recovery; nulls fat |
| Susceptibility | Local field distortion near air, bone, or metal |
| \(T_1\) | Longitudinal (spin–lattice) relaxation time |
| \(T_2\), \(T_2^*\) | Transverse relaxation; \(T_2^*\) includes static inhomogeneity |
| TE, TR, TI | Echo, repetition, and inversion times |
| TOF | Time-of-flight angiography; inflow of unsaturated blood |
Review Questions
- Why is \(^1\)H the dominant MRI nucleus, and how does the Larmor frequency depend on \(B_0\)?
- Compute the Boltzmann polarization at 3 T and explain how MRI produces usable signal from so small an excess.
- Explain the rotating frame and derive \(\alpha = \gamma B_1\tau_p\) from it.
- Write the full Bloch equations and identify which terms decouple in the absence of RF.
- Distinguish \(T_1\), \(T_2\), and \(T_2^*\), and state which dephasing a spin echo reverses and why.
- Why does \(T_1\) rise with field while \(T_2\) does not, and what does this imply for transferring a 1.5 T protocol to 3 T?
- Using the spin-echo signal equation, explain how to make an image T1-, T2-, or PD-weighted, and how grey–white contrast can reverse.
- Derive the inversion null time, and explain why FLAIR requires the finite-TR correction while STIR does not.
- What does the centre of k-space control, and what does the periphery control?
- How does the field of view relate to k-space sampling, and what happens if it is too small?
- Compute a slice thickness from a gradient amplitude and an RF bandwidth.
- Derive the Ernst angle from the steady-state gradient-echo signal equation.
- State the SNR relation carefully, and explain why “halving the voxel costs a factor of 8” is usually wrong.
- Why is magnitude MRI noise Rician, and what practical measurement does the noise floor corrupt?
- Why do motion artifacts appear along the phase-encode direction while chemical shift appears along the readout direction?
- Compute the fat–water displacement in pixels at 3 T for a 200 Hz per pixel bandwidth.
- State the four MRI hazard mechanisms and the physical quantity that governs each.
- Why does SAR rise as \(B_0^2\), and what protocol change compensates at 3 T?
- Explain why acute stroke is bright on DWI and dark on ADC, and how \(T_2\) shine-through can mimic it.
- What does fractional anisotropy measure, why does DTI need at least six directions, and why is FA rotation invariant?
- How does the BOLD signal relate neural activity to \(T_2^*\), and why is fMRI inherently a statistical method?
- Why does gadolinium brighten tissue on T1-weighted images, and what does enhancement indicate?
Problems
Problem 1. Larmor frequency
Using \(\bar\gamma = 42.577\) MHz/T for \(^1\)H and \(11.262\) MHz/T for \(^{23}\)Na, find (a) the \(^1\)H resonance frequency at 1.5 T and 3 T, and (b) the \(^{23}\)Na frequency at 3 T.
Problem 2. Polarization
Compute the Boltzmann polarization for \(^1\)H at 3 T and 310 K using Eq. (3). (a) Express it in parts per million. (b) If a voxel contains \(10^{21}\) protons, how many contribute net signal?
Problem 3. Flip angle
A rectangular RF pulse has \(B_1 = 15\) µT. (a) What duration produces a \(90^\circ\) flip? (b) A \(180^\circ\) flip? (c) By what factor does the SAR of the \(180^\circ\) pulse exceed that of the \(90^\circ\) pulse?
Problem 4. \(T_1\) recovery
After a \(90^\circ\) pulse, \(M_z(t) = M_0(1 - e^{-t/T_1})\). (a) What fraction has recovered at \(t = T_1\) and \(t = 2T_1\)? (b) How many \(T_1\) are needed for 95% recovery?
Problem 5. \(T_2\) decay
For \(M_{xy}(t) = M_{xy}(0)e^{-t/T_2}\), find the fraction remaining at \(T_2\), \(2T_2\), \(3T_2\).
Problem 6. Engineering contrast
Tissue A has \(T_1 = 800\) ms, \(T_2 = 80\) ms; tissue B has \(T_1 = 1200\) ms, \(T_2 = 100\) ms. Using \(S \propto (1-e^{-TR/T_1})e^{-TE/T_2}\) with TR = 500 ms and TE = 15 ms, compute the relative signals, state which is brighter, and name the weighting. Then find a TE at which the ranking reverses at TR = 3000 ms.
Problem 7. \(T_2^*\)
A tissue has \(T_2 = 100\) ms and local inhomogeneity contributes \(T_2' = 50\) ms. Find \(T_2^*\). What field inhomogeneity \(\Delta B_0\) does that \(T_2'\) correspond to?
Problem 8. Inversion nulling with finite TR
Using \(T_1(\text{CSF}) = 4300\) ms at 3 T, compute the FLAIR inversion time (a) from \(T_1\ln2\) and (b) from Eq. (16) with TR = 9000 ms. (c) Which matches clinical practice, and what would using the wrong one do to the image?
Problem 9. Slice thickness
A slice-select gradient \(G_z = 10\) mT/m is used with an RF bandwidth of 2 kHz. Find the slice thickness. What gradient would be needed for a 1 mm slice at the same bandwidth?
Problem 10. Field of view and aliasing
A k-space sampling step gives \(\Delta k = 1/(0.20\ \text{m})\). (a) What is the FOV? (b) If the object is 0.26 m wide, what happens and by how much? (c) Name two fixes and their costs.
Problem 11. Resolution and SNR
- At fixed FOV and slice thickness, the in-plane matrix is doubled in both directions. By what factor does SNR change? (b) Repeat for a 3D acquisition in which all three dimensions are halved at fixed FOV. (c) By what factor must NSA increase to double SNR, and what does it cost in time?
Problem 12. Ernst angle
A fast gradient-echo scan uses TR = 20 ms on tissue with \(T_1 = 800\) ms. (a) Find the Ernst angle. (b) Compute the ratio of the signal at the Ernst angle to that at \(90^\circ\).
Problem 13. Chemical shift
At 3 T with a receiver bandwidth of 150 Hz per pixel, compute the fat–water displacement in pixels. (a) What bandwidth would reduce it below one pixel? (b) What does that cost?
Problem 14. SAR scaling
A fast-spin-echo protocol at 1.5 T with \(180^\circ\) refocusing is moved to 3 T unchanged. (a) By what factor does SAR increase? (b) If the refocusing angle is reduced to \(120^\circ\), what is the new factor? (c) Is that sufficient to return to the original SAR?
Problem 15. ADC from DWI
A diffusion scan gives \(S/S_0 = 0.30\) at \(b = 1000\) s mm\(^{-2}\). (a) Compute the ADC. (b) If the true signal at \(b = 3000\) is 9.1 but the Rician noise floor with \(\sigma = 5\) raises the measured magnitude to 10.6, compute the ADC that would be fitted and the percentage error.
Problem 16. Fractional anisotropy
A white-matter voxel has eigenvalues \(\lambda = (1.70, 0.30, 0.30)\times10^{-3}\) mm\(^2\)s\(^{-1}\). (a) Compute MD and FA. (b) Repeat for a planar tensor \((1.00, 1.00, 0.10)\times10^{-3}\) and comment on what FA alone fails to distinguish.
Solutions
| Problem | Computed answer |
|---|---|
| 1(a) | 63.9 MHz at 1.5 T, 127.7 MHz at 3 T |
| 1(b) | 33.8 MHz |
| 2(a) | 9.9 ppm |
| 2(b) | 9.89e+15 protons |
| 3(a) | 0.391 ms |
| 3(b) | 0.783 ms |
| 3(c) | 4x (SAR ~ alpha^2) |
| 4(a) | 63.2% and 86.5% |
| 4(b) | 3.00 T1 |
| 5 | 36.8%, 13.5%, 5.0% |
| 6 | A = 0.385, B = 0.293 -> A brighter (T1-weighted) |
| 7 | T2* = 33.3 ms |
| 8(a) | 2981 ms |
| 8(b) | 2481 ms |
| 9 | 4.70 mm; need Gz = 47.0 mT/m for 1 mm |
| 10(a) | FOV = 0.20 m |
| 11(a) | SNR falls by 2 |
| 11(b) | SNR falls by 2.83 (= 8/sqrt(8)) |
| 11(c) | NSA x4, so 4x the scan time |
| 12(a) | 12.8 degrees |
| 12(b) | 4.53 |
| 13 | 2.98 pixels; need BW > 447 Hz/px |
| 14(a) | 16x |
| 14(b) | 7.1x |
| 14(c) | no: still 7x the original |
| 15(a) | 1.204e-03 mm^2/s |
| 15(b) | 7.481e-04 mm^2/s (-6.4% error) |
| 16(a) | MD = 7.667e-04, FA = 0.799 |
| 16(b) | MD = 7.000e-04, FA = 0.635 |
Solution 2
From Eq. (3), the polarization at 3 T and 310 K is \(\gamma\hbar B_0/(2k_BT) = 9.9\times10^{-6}\), i.e. 9.9 ppm. In a voxel of \(10^{21}\) protons that is about \(10^{16}\) spins in excess. The lesson is that MRI is a low-sensitivity, high-abundance technique: it images a one-in-a-hundred-thousand population imbalance and succeeds only because the population is astronomically large. Contrast this with PET (Chapter 6), which counts individual nuclear decays and is therefore picomolar-sensitive but has far poorer spatial resolution.
Solution 3
- \(\tau_p = (\pi/2)/(\gamma B_1) = 1.5708/(2.675\times10^8\times1.5\times10^{-5}) = 3.91\times10^{-4}\) s \(= 0.39\) ms. (b) Twice that, 0.78 ms. (c) SAR scales as \(\alpha^2\), so the \(180^\circ\) pulse deposits four times the energy — the reason spin-echo trains dominate the SAR budget of a protocol (Section 4.8.1).
Solution 6
Tissue A: \((1-e^{-0.625})e^{-0.1875} = 0.4647\times0.8290 = 0.385\). Tissue B: \((1-e^{-0.4167})e^{-0.15} = 0.3408\times0.8607 = 0.293\). Tissue A is brighter, and short TR with short TE is T1-weighted. At TR = 3000 ms the \(T_1\) factors become 0.976 and 0.918, nearly equal, so the \(T_2\) term dominates; B overtakes A once \(e^{-TE/100}/e^{-TE/80} > 0.918/0.976\), i.e. \(TE > 0.0025^{-1}\ln(1.063) \approx 24\) ms. By TE = 100 ms the ranking has firmly reversed — the same reversal shown in the contrast map of Section 4.3.1.
Solution 7
\(1/T_2^* = 1/100 + 1/50 = 0.03\) ms\(^{-1}\), so \(T_2^* = 33.3\) ms — shorter than both, as Eq. (13) requires. Since \(1/T_2' \approx \gamma\Delta B_0\), \(\Delta B_0 \approx 1/(\gamma T_2') = 1/(2.675\times10^8\times0.050) = 7.5\times10^{-8}\) T, about 0.075 µT — roughly 0.025 ppm at 3 T. Field uniformity of a few parts per hundred million is what shimming must achieve.
Solution 8
- \(T_1\ln2 = 4300\times0.693 = 2981\) ms. (b) With TR = 9000 ms, \(TI = 4300\ln[2/(1+e^{-9000/4300})] = 2481\) ms. (c) The finite-TR value matches clinical practice, where 3 T FLAIR uses TI near 2400–2500 ms. Using 2981 ms would invert past the null: CSF would have recovered to a positive value and would appear as residual bright signal, exactly the artefact FLAIR exists to eliminate, potentially obscuring periventricular lesions.
Solution 11
- Doubling the matrix in both in-plane directions at fixed FOV quarters the voxel area while quadrupling the samples: \(\mathrm{SNR} \propto (1/4)\sqrt{4} = 1/2\), so SNR falls by 2. (b) Halving all three dimensions in a 3D acquisition gives \(V/8\) and \(8\times\) the samples: \((1/8)\sqrt8 = 0.354\), so SNR falls by 2.83. (c) \(\mathrm{SNR}\propto\sqrt{\rm NSA}\), so doubling SNR requires four times the averages and four times the scan time — which is why averaging is the last resort rather than the first.
Solution 14
- SAR \(\propto B_0^2\), so moving from 1.5 T to 3 T multiplies it by \((3/1.5)^2 = 4\); the refocusing angle is unchanged, so the total factor is 4. (b) Reducing refocusing from \(180^\circ\) to \(120^\circ\) multiplies by \((120/180)^2 = 0.444\), giving \(4\times0.444 = 1.78\). (c) No — the protocol still deposits 1.78 times the original SAR. Further reduction requires some combination of a longer TR, fewer slices per TR, a shorter echo train, or a still lower refocusing angle, each of which costs scan time or contrast.
Solution 15
- \(\mathrm{ADC} = -\ln(0.30)/1000 = 1.20\times10^{-3}\) mm\(^2\)s\(^{-1}\), a normal-brain value. (b) The true signal 9.1 would give \(-\ln(0.091)/3000 = 7.99\times10^{-4}\), but the Rician-lifted measurement 10.6 gives \(-\ln(0.106)/3000 = 7.48\times10^{-4}\), an error of \(-6.4\)%. The bias is always downward in ADC, because the noise floor makes tissue look less diffusive than it is — and it worsens rapidly with b-value, as the table in Section 4.6.2 shows.
Solution 16
- \(\mathrm{MD} = 2.30/3 = 0.767\times10^{-3}\); deviations \((+0.933, -0.467, -0.467)\) give \(\sum = 1.306\) and \(\sum\lambda_i^2 = 3.07\), so \(\mathrm{FA} = 1.2247\sqrt{1.306}/\sqrt{3.07} = 0.80\): high, consistent with organized white matter. (b) For \((1.00, 1.00, 0.10)\), \(\mathrm{MD} = 0.70\times10^{-3}\) and \(\mathrm{FA} = 0.64\) — still substantial, yet this tensor is planar, with no single dominant direction. FA alone cannot distinguish linear from planar anisotropy, which is exactly the failure mode in voxels containing crossing fibres, and the reason tractography based on a single principal eigenvector produces spurious tracts there.
References and Further Reading
- Brown, R. W., Cheng, Y.-C. N., Haacke, E. M., Thompson, M. R., and Venkatesan, R. (2014). Magnetic Resonance Imaging: Physical Principles and Sequence Design, 2nd ed. Wiley.
- Nishimura, D. G. (2010). Principles of Magnetic Resonance Imaging. Stanford University.
- McRobbie, D. W., Moore, E. A., Graves, M. J., and Prince, M. R. (2017). MRI from Picture to Proton, 3rd ed. Cambridge University Press.
- Bernstein, M. A., King, K. F., and Zhou, X. J. (2004). Handbook of MRI Pulse Sequences. Elsevier.
- Haacke, E. M. et al. (2004). Susceptibility weighted imaging (SWI). Magnetic Resonance in Medicine, 52(3), 612–618.
- Gudbjartsson, H. and Patz, S. (1995). The Rician distribution of noisy MRI data. Magnetic Resonance in Medicine, 34(6), 910–914.
- Jones, D. K. and Basser, P. J. (2004). Squashing peanuts and smashing pumpkins: how noise distorts diffusion-weighted MR data. Magnetic Resonance in Medicine, 52(5), 979–993.
- Le Bihan, D. et al. (2001). Diffusion tensor imaging: concepts and applications. Journal of Magnetic Resonance Imaging, 13(4), 534–546.
- Ogawa, S. et al. (1990). Brain magnetic resonance imaging with contrast dependent on blood oxygenation. PNAS, 87(24), 9868–9872.
- Lustig, M., Donoho, D., and Pauly, J. M. (2007). Sparse MRI: the application of compressed sensing for rapid MR imaging. Magnetic Resonance in Medicine, 58(6), 1182–1195.
- Pruessmann, K. P. et al. (1999). SENSE: sensitivity encoding for fast MRI. Magnetic Resonance in Medicine, 42(5), 952–962.
- International Electrotechnical Commission (2015). IEC 60601-2-33: Particular requirements for the safety of magnetic resonance equipment for medical diagnosis, 3rd ed.
- Kanal, E. et al. (2013). ACR guidance document on MR safe practices. Journal of Magnetic Resonance Imaging, 37(3), 501–530.
- Stanisz, G. J. et al. (2005). T1, T2 relaxation and magnetization transfer in tissue at 3 T. Magnetic Resonance in Medicine, 54(3), 507–512.
- Prince, J. L. and Links, J. M. (2014). Medical Imaging Signals and Systems, 2nd ed. Pearson.
- Bushberg, J. T., Seibert, J. A., Leidholdt, E. M., and Boone, J. M. (2020). The Essential Physics of Medical Imaging, 4th ed. Wolters Kluwer.
- Dinov, I. D. (2021). Data Science: Time Complexity, Inferential Uncertainty, and Spacekime Analytics. De Gruyter.
- SOCR Probability and Statistics EBook. See also BPAD Chapters 1, 2, 3, 5, 7, and 8.