Magnetic Resonance Imaging (MRI)

Magnetic resonance imaging (MRI, also called MR imaging, and in its first years NMR imaging or “zeugmatography”) produces cross-sectional images and volumes of the body from the radio signal emitted by hydrogen nuclei in a strong magnetic field. It differs from a camera, and from computed tomography, in what the scanner actually records: not intensities at pixels, and not line integrals, but samples of the spatial Fourier transform of the magnetization, a space the MR community calls k-space. The image is reconstructed from those samples, so how k-space is sampled decides resolution, field of view, scan time and many of the artifacts. The brightness of a tissue is not a calibrated physical unit either: it depends on the tissue’s proton density and relaxation times and on the timing parameters the operator chooses, which is how one scanner produces T1-weighted, T2-weighted and many other contrasts of the same anatomy, and why MR intensities have, in the words of West et al. (2013), an “arbitrary grey-scale”.

How It Works

Magnetization and precession

Hydrogen nuclei (protons) carry a magnetic moment. In a strong static field $B_0$ the moments tend to align with the field, giving a net magnetization $\mathbf{M}$ along it, and when tipped away they precess around the field direction at the Larmor frequency

$$f_\ell = \bar{\gamma}\, B_0 ,$$

where $\bar{\gamma}$ is the gyromagnetic ratio expressed in hertz per tesla; Lugauer and Wetzl give 42.576 MHz/T for hydrogen, so protons resonate at about 64 MHz at 1.5 T (and, by the same formula, at about 128 MHz at 3 T). Many texts write $\omega = \gamma B_0$ with $\gamma = 2\pi\bar{\gamma}$ in radians per second per tesla. A radio-frequency (RF) pulse at the Larmor frequency, applied through a transmit coil, tips $\mathbf{M}$ away from $B_0$ (excitation); the precessing transverse component then induces a voltage in receive coils near the body. Water makes up about two thirds of body weight, and the Nobel Assembly’s 2003 press release gives this high water content as the reason MRI became so widely applicable in medicine.

Relaxation: T1, T2, T2* and proton density

After excitation the magnetization returns to equilibrium by two processes with tissue-specific time constants (Lugauer & Wetzl, 2018). After a 90° pulse, the longitudinal component recovers and the transverse component decays as

$$M_z(t) = M_0 \left(1 - e^{-t/T_1}\right), \qquad M_{xy}(t) = M_0\, e^{-t/T_2} .$$

$T_1$ (spin-lattice or longitudinal relaxation) is the time constant of the recovery, and $T_2$ (spin-spin or transverse relaxation) the time after which the transverse signal has fallen to $1/e$ of its initial value. In practice the received signal decays faster than $T_2$ predicts, because inhomogeneity of $B_0$ from magnet imperfections and tissue magnetic susceptibilities dephases the spins; this faster decay has the time constant $T_2^*$, which is always shorter than $T_2$. The equilibrium magnetization $M_0$ is proportional to the density of mobile protons, the proton density (PD).

These three tissue parameters are what MR contrast is built from, and they depend on the main field strength. West et al. (2013) measured them in ten healthy adults on a 1.5 T and a 3.0 T scanner from the same vendor with the same quantitative sequence. At 3.0 T their average white matter had $R_1 = 1/T_1 = 1.38\ \text{s}^{-1}$ ($T_1 \approx 725$ ms), $R_2 = 1/T_2 = 13.76\ \text{s}^{-1}$ ($T_2 \approx 73$ ms) and PD 65.6%; frontal cortex had $T_1 \approx 1.30$ s, $T_2 \approx 94$ ms and PD 86.9%; and cerebrospinal fluid (CSF) had $T_1 \approx 3.1$ s and $T_2 \approx 1.3$ s, with a large spread between subjects ($R_2 = 0.79 \pm 0.59\ \text{s}^{-1}$). White matter $T_1$ was shorter at 1.5 T ($R_1 = 1.62\ \text{s}^{-1}$, about 617 ms). At the much lower field of 0.064 T, Jordanova et al. (2023) measured white matter $T_1$ of about 0.25 s in manually placed regions. Relaxation times are therefore quoted with their field strength and measurement method, and a contrast tuned for one field does not transfer unchanged to another.

Spatial encoding and k-space

Relaxation explains contrast but not location. Location comes from gradient coils that add a field varying linearly with position, $\mathbf{G}\cdot\mathbf{r}$, to $B_0$, so that the precession frequency depends on position (Lugauer & Wetzl, 2018). Lauterbur described this idea in his Nobel lecture: since NMR frequencies depend on the local field, a non-uniform field can locate the signal (Lauterbur, 2003).

Three uses of gradients produce an image:

• Slice selection. With a gradient along $z$ switched on, an RF pulse containing a narrow band of frequencies excites only the nuclei whose Larmor frequency falls in that band, which is a slab of tissue; the bandwidth sets the slice thickness, a trade-off between resolution and signal-to-noise ratio (SNR) (Lugauer & Wetzl, 2018). Mansfield’s group defined thin slices in this way with a technique they called selective irradiation (Mansfield, 2003).
• Phase and frequency encoding. Within the slice, a gradient applied for a time makes the phase of the spins vary linearly with position, and the coil measures the sum of all spins weighted by that phase pattern. Lugauer and Wetzl show that this sum is a Fourier coefficient of the magnetization: applying different gradients and measuring the net signal performs a Fourier decomposition of the image. Texts distinguish a frequency-encoding (readout) direction, sampled continuously while a gradient is on, from phase-encoding directions, stepped between excitations, but the principle is the same.

Ignoring relaxation during the readout, coil sensitivity and field errors, the demodulated signal from a slice with magnetization $\rho(\mathbf{r})$ is

$$s(t) = \int \rho(\mathbf{r})\, e^{-i 2\pi\, \mathbf{k}(t) \cdot \mathbf{r}}\, d\mathbf{r}, \qquad \mathbf{k}(t) = \bar{\gamma} \int_0^t \mathbf{G}(\tau)\, d\tau ,$$

so the signal at time $t$ is the Fourier transform of $\rho$ at spatial frequency $\mathbf{k}(t)$, in cycles per metre, and the gradient waveforms steer a k-space trajectory. Texts differ in the sign of the exponent and in whether $\mathbf{k}$ includes the factor $2\pi$: Mansfield’s lecture and Lugauer and Wetzl write a positive exponent, Knoll et al. (2020) a negative one with angular $k$. The sign here matches the forward transform of NumPy’s FFT, which the example below uses; with the opposite sign, the image comes out flipped along every axis.

By the sampling theorem with space and frequency exchanged, sampling k-space on a Cartesian grid with spacing $\Delta k$ makes the reconstructed image periodic with period $1/\Delta k$, the field of view (FOV), so an object larger than the FOV wraps around; $N$ samples spanning $\pm k_\text{max}$ give a pixel size of $1/(N\Delta k) \approx 1/(2k_\text{max})$. Sampling to a lower $k_\text{max}$ lowers resolution, and leaving out samples within $\pm k_\text{max}$ causes aliasing (Zbontar et al., 2018).

In a basic Cartesian sequence each excitation acquires one line of k-space along the readout and successive excitations step through the phase-encoding lines; for a typical spin echo sequence the time for one slice is $T_R \cdot N_y \cdot N_\text{ex}$, with repetition time $T_R$, $N_y$ phase-encoding steps and $N_\text{ex}$ excitations (Lugauer & Wetzl, 2018). The readout takes milliseconds while the whole acquisition takes seconds to minutes, which matters for motion. Echo-planar imaging (EPI), introduced by Mansfield in 1977, instead covers the whole 2D k-space plane after a single excitation by oscillating the readout gradient and blipping the phase-encoding gradient, which takes 20–50 ms per image (Mansfield, 2003). Radial, spiral and 3D trajectories exist as well; 3D acquisitions excite a whole volume and phase-encode along two directions.

Sequences and contrast

A pulse sequence is the timed pattern of RF pulses, gradients and readouts that is repeated to fill k-space (Lugauer & Wetzl, 2018). Two families are basic:

• Spin echo (SE). A 90° pulse tips the magnetization into the transverse plane; after $T_E/2$ a 180° pulse flips the dephasing spins so that the faster ones fall behind the slower ones; at the echo time $T_E$ they come back in phase, and the echo is free of the $T_2^*$ loss from static field imperfections. Spin echo images are $T_2$-weighted rather than $T_2^*$-weighted.
• Gradient echo (GRE). A smaller flip angle $\alpha$ leaves $M_0\cos\alpha$ along $z$, so the repetition time can be short, and the echo is formed by reversing the readout gradient instead of a 180° pulse. Without refocusing, GRE is more sensitive to field inhomogeneity and gives $T_2^*$ weighting.

The operator’s main contrast controls are $T_R$, the time between successive excitations, and $T_E$, the time from excitation to the echo. For a spin echo sequence, combining the two relaxation laws above gives the familiar signal model

$$S \propto \text{PD} \left(1 - e^{-T_R/T_1}\right) e^{-T_E/T_2} .$$

Each repetition starts with a 90° pulse that leaves no longitudinal magnetization, which recovers for $T_R$ before the next 90° pulse turns it into transverse signal, and that signal then decays with $T_2$ until the echo. This is our combination of the textbook’s two laws, valid when $T_E \ll T_R$ and no transverse magnetization survives to the next excitation; it neglects the 180° pulse’s own inversion of $M_z$. It reproduces the weighting rules of Lugauer and Wetzl:

	short $T_E$	long $T_E$
short $T_R$ (about $T_1$)	$T_1$ weighting	(mixed, weak signal)
long $T_R$	PD weighting	$T_2$ weighting

A short $T_R$ lets tissues with short $T_1$ recover more and appear brighter; a long $T_E$ lets tissues with short $T_2$ decay more and appear darker. Inversion recovery adds a 180° preparation pulse and waits an inversion time $T_I$ before the imaging sequence; by the recovery law above, now starting from $-M_0$, a tissue whose magnetization crosses zero at $T_I$ gives no signal (our reasoning); DICOM records $T_I$ as the time from the middle of the inverting pulse to the middle of the excitation pulse (NEMA, 2026, C.8.3.1). Many other preparations exist: Lugauer and Wetzl describe $T_2$ preparation for coronary imaging, time-of-flight angiography without contrast agent, and the BOLD effect used in functional MRI.

Data It Produces

The scanner measures complex-valued k-space for each receive coil. Modern scanners use arrays of small coil elements, each most sensitive near itself (Knoll et al., 2020; Zbontar et al., 2018). Fully sampled data are reconstructed by an inverse Fourier transform per coil, and coil images are often combined by the root sum of squares, $\left(\sum_i |m_i|^2\right)^{1/2}$, the ground truth used by the fastMRI dataset (Zbontar et al., 2018). Displayed images are usually such magnitudes; the phase is discarded except in methods that use it. The magnitude of a complex value with Gaussian noise in its real and imaginary parts is not Gaussian: it follows a Rice distribution, which is biased upward where the signal is low. That is a property of the distribution, and it matters wherever image noise in dark regions is modelled; Jordanova et al. (2023), for example, name the Rician noise of magnitude images as a limitation of their relaxation fits.

Raw k-space is rarely what vision code receives. Clinical archives store reconstructed images, and these are not interchangeable with raw data: the fastMRI authors warn that the Fourier transform of a DICOM image does not correspond to the measured k-space, because vendor reconstruction adds post-processing and most clinical images come from accelerated acquisitions reconstructed with parallel imaging (Zbontar et al., 2018). Simulating undersampling by transforming archived magnitude images therefore leaves out the noise, phase and multi-coil structure of real acquisitions.

Geometry depends on the acquisition. Slice-selective 2D acquisitions are typically non-isotropic, with in-plane resolution below 1 mm and slices several millimetres thick, while 3D acquisitions typically have isotropic voxels between the two (Lugauer & Wetzl, 2018). Slices need not touch: West et al. acquired 3 mm slices with a 0.5 mm gap. A volume is therefore a stack of samples with a spacing that differs by axis, and its intensities have no absolute scale, unlike CT’s Hounsfield units.

DICOM

Clinical MR images are stored and exchanged in DICOM. The standard’s MR Image Module (NEMA, 2026, C.8.3.1) defines one sample per pixel, 16 bits allocated and a monochrome photometric interpretation, and records the acquisition: Scanning Sequence (SE spin echo, IR inversion recovery, GR gradient recalled, EP echo planar, RM research mode), Repetition Time and Echo Time in milliseconds, Inversion Time, Echo Train Length, Imaging Frequency in MHz, Imaged Nucleus and the nominal Magnetic Field Strength in tesla. These fields let a pipeline group series by contrast and field strength instead of trusting series descriptions.

Geometry comes from the Image Plane Module (NEMA, 2026, C.7.6.2). Image Position (Patient) gives the x, y, z coordinates in millimetres of the centre of the first transmitted pixel; Image Orientation (Patient) gives the direction cosines of the first row and the first column; Pixel Spacing gives the physical distance between pixel centres, with the spacing between adjacent rows first and the spacing between adjacent columns second. With these, the pixel in column $i$ and row $j$ (both counted from zero) lies at

$$\mathbf{P} = \mathbf{S} + i\,\Delta_c\,\mathbf{X} + j\,\Delta_r\,\mathbf{Y},$$

where $\mathbf{S}$ is the image position, $\mathbf{X}$ and $\mathbf{Y}$ are the row and column direction cosines, and $\Delta_c$ and $\Delta_r$ the column and row spacings. The patient-based coordinate system is right-handed, with x increasing toward the patient’s left, y toward the posterior, and z toward the head. Slice Thickness and Spacing Between Slices are separate attributes, which is how a gap is recorded.

NIfTI

Research pipelines, especially in neuroimaging, often convert DICOM series to NIfTI-1, a format adapted from ANALYZE 7.5 that stores a volume as a single .nii file or a .hdr/.img pair (NIfTI documentation). Its 348-byte header holds the array dimensions (dim), the grid spacings (pixdim) and their units, a data scaling (scl_slope, scl_inter), slice-timing fields for functional MRI, and two ways of relating voxel indices $(i, j, k)$ to spatial coordinates $(x, y, z)$: a rotation stored as a quaternion with offsets (the “qform”, flagged by qform_code) and a general affine matrix stored row by row in srow_x, srow_y, srow_z (the “sform”, flagged by sform_code). Vision code that drops the header and keeps only the voxel array loses the spacing and orientation, so left and right, voxel size and alignment with other scans all become guesses. Comparing the affine after conversion with the DICOM geometry above is a simple check that the conversion kept orientation.

Characteristics

Property	What sets it	Notes
Field strength	Magnet	1.5 T and 3 T are the conventional systems; Jordanova et al. (2023) describe renewed interest in systems at 0.55 T and below for point-of-care use. Relaxation times and contrast change with field.
Contrast	$T_R$, $T_E$, $T_I$, flip angle, preparation	One anatomy, many contrasts; no absolute intensity scale.
Field of view and resolution	k-space spacing and extent	FOV $= 1/\Delta k$; pixel size about $1/(2k_\text{max})$.
Slice geometry	2D slice-selective or 3D	Anisotropic 2D stacks or near-isotropic 3D volumes.
Scan time	$T_R$, number of phase encodes, excitations	Seconds to minutes per series.
Acceleration	Coils, sampling pattern, reconstruction	Clinical parallel imaging factors are typically about two to three (Zbontar et al., 2018).

Calibration

MRI has no equivalent of a camera’s geometric calibration target that the user runs, but several calibrations happen in every exam. The field is made uniform by shimming: Lauterbur (2003) describes the correction or “shim” fields on NMR machines that cancel field non-uniformities term by term. Coil sensitivities are estimated for each scan, either from a separate low-resolution reference scan or from a fully sampled region at the centre of k-space acquired with the undersampled data; one simple estimate divides each coil’s low-resolution image by the root sum of squares of all of them (Lugauer & Wetzl, 2018; Knoll et al., 2020). Quantitative MRI, which maps $T_1$, $T_2$ or PD in physical units, is validated against reference samples: Jordanova et al. (2023) compared their in vivo protocol with measurements of test samples on a separate NMR system at the same field, and found $T_1$ within 10% for nine of ten samples and long $T_2$ values underestimated.

Accelerated Acquisition

Because time per image grows with the number of phase-encoding lines, much of MR engineering is about acquiring fewer samples. Undersampling violates the sampling condition, so a plain inverse FFT aliases, as the first script below shows; acceleration methods remove the aliasing with extra information. Parallel imaging uses the different sensitivity maps of the coils in an array: SENSE (Pruessmann et al., 1999) unfolds the aliased pixels in image space with known sensitivities, and GRAPPA (Griswold et al., 2002) estimates each missing k-space line from acquired neighbouring lines of all coils, with kernels fitted to a fully sampled autocalibration region (Knoll et al., 2020). Knoll et al. write that essentially all clinical scanners include parallel imaging and that it is the default for many protocols. Compressed sensing (Lustig et al., 2007) draws phase encodes pseudo-randomly with a density that falls with distance from the k-space centre, so that the aliasing is incoherent, and reconstructs with a nonlinear method that enforces sparsity in a known transform and consistency with the data. Learned reconstruction replaces hand-designed regularizers with trained networks, developed on shared raw multi-coil data such as the fastMRI dataset (Zbontar et al., 2018). Medical image reconstruction describes these methods, their benchmarks and their evaluation, including the risk of hallucinated detail.

Practical Example

The three scripts below simulate k-space acquisition with NumPy 2.3.5 and scikit-image 0.25.2. They use the Shepp–Logan phantom bundled with scikit-image (skimage.data.shepp_logan_phantom, a 400 × 400 image read from the package with no download) and a synthetic tissue map.

Conventions matter in FFT code. The scripts use NumPy’s fft2 and ifft2 with norm="ortho", which makes the transform unitary, so image and k-space have the same energy and zero-filled reconstructions keep their scale. They also apply ifftshift before and fftshift after each transform, so that $k = 0$ is at index $N/2$ of the k-space array and the centre of the object is at index $N/2$ of the image. Axis 0 is the phase-encoding direction and axis 1 the readout.

One simplification applies throughout: the phantom’s k-space is the DFT of a pixel grid, whereas a real object is continuous, so simulating data and reconstructing on the same grid can look better than real acquisitions.

The first script compares full sampling with uniform undersampling:

import numpy as np
from skimage.data import shepp_logan_phantom

# Centred, orthonormal 2D DFT: k = 0 sits at index N//2, the image centre at N//2.
def to_kspace(img):
    return np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(img), norm="ortho"))

def to_image(ksp):
    return np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(ksp), norm="ortho"))

def nrmse(a, ref):
    return np.linalg.norm(a - ref) / np.linalg.norm(ref)

img = shepp_logan_phantom()          # 400 x 400, values 0..1, bundled with scikit-image
N = img.shape[0]
K = to_kspace(img)                   # axis 0 = phase encoding (ky), axis 1 = readout (kx)

print("full sampling, max error:", f"{np.abs(to_image(K) - img).max():.1e}")
centre = np.abs(K[N//2-8:N//2+8, N//2-8:N//2+8])**2
print("energy in central 16x16 of k-space:", f"{centre.sum() / (np.abs(K)**2).sum():.3f}")

# Uniform undersampling: keep every R-th phase-encoding line, fill the rest with zeros.
for R in (2, 4):
    mask = np.zeros(N, bool); mask[N//2 % R::R] = True
    zf = to_image(K * mask[:, None])
    copies = sum(np.roll(img, s * N // R, axis=0) for s in range(R)) / R
    print(f"R={R} uniform: {mask.sum()} lines, NRMSE {nrmse(np.abs(zf), img):.3f}, "
          f"max |zero-filled - mean of {R} shifted copies| {np.abs(zf - copies).max():.1e}")

full sampling, max error: 8.9e-16
energy in central 16x16 of k-space: 0.602
R=2 uniform: 200 lines, NRMSE 0.586, max |zero-filled - mean of 2 shifted copies| 5.6e-16
R=4 uniform: 100 lines, NRMSE 0.712, max |zero-filled - mean of 4 shifted copies| 4.5e-16

Full sampling inverts exactly, and 60% of the phantom’s energy sits in the central 16 × 16 samples, 0.16% of k-space, which is why every practical mask samples the centre densely. Keeping every $R$-th line makes the zero-filled image equal, to rounding error, to the mean of $R$ copies of the object shifted by $N/R$ pixels along the phase-encoding axis (exactly so because the kept lines include $k = 0$ and $R$ divides $N$; otherwise the copies carry phase factors): the field of view has shrunk to FOV$/R$ and the object wraps into itself. This is the folding that SENSE undoes pixel group by pixel group. Medical image reconstruction continues the script with variable-density random masks of the same number of lines, whose interference is spread into a noise-like background instead of a coherent ghost.

The second script truncates k-space to a smaller acquired matrix, as every acquisition does, with and without a tapering window:

import numpy as np
from skimage.data import shepp_logan_phantom

def to_kspace(img):
    return np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(img), norm="ortho"))

def to_image(ksp):
    return np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(ksp), norm="ortho"))

img = shepp_logan_phantom()          # non-negative, piecewise constant, 400 x 400
N = img.shape[0]
K = to_kspace(img)

# Keep only the central n x n samples (n odd, symmetric about k = 0), with or
# without a Hann taper. The object is real and non-negative, so any negative
# value in the real part of the result is ringing, not blur.
print(" kept   | taper | most negative value | edge 10-90% width (px)")
for n in (201, 101, 51):
    lo = N//2 - n//2
    for taper in ("none", "Hann"):
        w = np.ones(n) if taper == "none" else np.hanning(n + 2)[1:-1]
        win = np.zeros(N); win[lo:lo + n] = w
        rec = to_image(K * np.outer(win, win))
        assert np.abs(rec.imag).max() < 1e-12
        rec = rec.real
        # 10-90% rise of the left skull edge on the central row (0 outside, 1 in the skull)
        row = rec[N//2, :N//2]
        i90 = np.argmax(row > 0.9 * row.max())
        i10 = np.flatnonzero(row[:i90] < 0.1 * row.max())[-1]
        print(f"{n:>3}x{n:<3} | {taper:<5} | {rec.min():+.3f}              | {i90 - i10}")

 kept   | taper | most negative value | edge 10-90% width (px)
201x201 | none  | -0.148              | 3
201x201 | Hann  | -0.006              | 5
101x101 | none  | -0.135              | 4
101x101 | Hann  | -0.008              | 7
 51x51  | none  | -0.141              | 6
 51x51  | Hann  | -0.012              | 12

Truncation multiplies k-space by a box, which convolves the image with an oscillating kernel, and the result is Gibbs ringing: ripples parallel to sharp edges. Since the true object is non-negative, the negative values measure the ringing directly. Their depth hardly changes with the matrix size; fewer samples make the ripples wider rather than weaker. A Hann taper suppresses them almost entirely but roughly doubles every edge width. This is the ringing-versus-blur trade-off of Fourier filtering; in MRI it is settled by the acquisition matrix and the reconstruction filter, before any vision code sees the image.

The third script builds a synthetic slice with white matter, a cortical grey matter ribbon, CSF around it and two CSF-filled ventricles, and renders it with the spin echo model and the 1.5 T and 3.0 T tissue parameters of West et al. (2013):

import numpy as np

# Tissue parameters measured with QRAPMASTER in 10 healthy adults (West et al., 2013,
# Table 1): R1 = 1/T1 and R2 = 1/T2 in 1/s, PD in percent. "WM" is their average
# white matter, "GM" frontal cortex, "CSF" cerebrospinal fluid.
TISSUE = {  #        R1    R2     PD    (3.0 T)        R1    R2     PD    (1.5 T)
    "WM":  {"3.0": (1.38, 13.76, 65.60), "1.5": (1.62, 12.84, 64.51)},
    "GM":  {"3.0": (0.77, 10.69, 86.88), "1.5": (0.81, 10.29, 88.06)},
    "CSF": {"3.0": (0.32,  0.79, 94.40), "1.5": (0.24,  0.87, 103.76)},
}

# A synthetic axial "slice": CSF around a cortical GM ribbon, WM inside, two CSF ventricles.
n = 256
y, x = np.mgrid[:n, :n] - n / 2
r = np.hypot(x / 1.0, y / 1.25)
labels = np.full((n, n), "", dtype=object)
labels[r < 110] = "CSF"
labels[r < 104] = "GM"
labels[r < 92] = "WM"
for cx in (-14, 14):
    labels[np.hypot((x - cx) / 7, (y + 5) / 22) < 1] = "CSF"

def spin_echo(TR, TE, field):
    """S = PD (1 - exp(-TR/T1)) exp(-TE/T2), times in ms (90-180 spin echo, TE << TR)."""
    img = np.zeros((n, n))
    for t, p in TISSUE.items():
        R1, R2, PD = p[field]
        img[labels == t] = PD * (1 - np.exp(-TR / 1000 * R1)) * np.exp(-TE / 1000 * R2)
    return img

protocols = [("T1-weighted", 500, 15), ("PD-weighted", 4000, 15), ("T2-weighted", 4000, 100)]
print("3.0 T          TR/TE (ms)   WM     GM     CSF   (signal / brightest tissue)")
for name, TR, TE in protocols:
    img = spin_echo(TR, TE, "3.0")
    s = {t: img[labels == t].mean() for t in TISSUE}
    top = max(s.values())
    print(f"{name:<13} {TR:>5}/{TE:<4}    " + "  ".join(f"{s[t] / top:.2f}" for t in TISSUE))

for f in ("1.5", "3.0"):
    img = spin_echo(500, 15, f)
    print(f"{f} T, T1-weighted 500/15: WM/GM signal ratio "
          f"{img[labels == 'WM'].mean() / img[labels == 'GM'].mean():.2f}")

3.0 T          TR/TE (ms)   WM     GM     CSF   (signal / brightest tissue)
T1-weighted     500/15      1.00  0.89  0.52
PD-weighted    4000/15      0.75  1.00  0.95
T2-weighted    4000/100     0.26  0.45  1.00
1.5 T, T1-weighted 500/15: WM/GM signal ratio 1.18
3.0 T, T1-weighted 500/15: WM/GM signal ratio 1.12

The same label map gives three different images. With a short $T_R$, white matter, whose $T_1$ is shortest, is brightest and CSF darkest; with a long $T_E$ the order reverses, because CSF’s long $T_2$ keeps its signal while the tissues decay. The “PD-weighted” protocol still carries some $T_1$ weighting: a 4 s $T_R$ is only about 1.3 times CSF’s $T_1$, so CSF has not fully recovered and grey matter comes out slightly brighter despite its lower proton density. The ordering of intensities, not only their scale, depends on the protocol, so a threshold or a model fitted to one contrast does not apply to another. The last two lines show the same for field strength: an identical protocol gives less white–grey contrast at 3.0 T than at 1.5 T; in our reading, because both $T_1$ values are longer at 3.0 T and their ratio is smaller (about 1.8 against 2.0 in West et al.’s measurements). The CSF values vary widely between subjects and the model leaves out much of real sequences, so the numbers illustrate the mechanism rather than predict a scanner’s output.

Artifacts

• Aliasing (wrap-around). An object larger than the FOV along a phase-encoding direction folds back into the image, as in the uniform undersampling above. Along the readout, scanners avoid it by oversampling, which is standard in clinical exams; fastMRI crops its images to the central region for this reason (Zbontar et al., 2018).
• Gibbs ringing. Ripples at sharp edges from the finite acquired matrix, as in the second script.
• Motion. Patient motion during the acquisition is inevitable, and breath-holds and short $T_R$ are used to limit respiratory and cardiac motion (Lugauer & Wetzl, 2018). Because lines are acquired over seconds to minutes, motion makes successive lines inconsistent, and since a line is a set of frequencies along the readout, the inconsistency spreads along the phase-encoding direction. That is our reading of the sampling model: periodic motion such as breathing puts ghost copies along the phase-encoding axis, at offsets set by the motion period in units of $T_R$.
• Field inhomogeneity and susceptibility. Imperfect $B_0$ and susceptibility differences dephase spins and shorten $T_2^*$, which gradient echoes show more than spin echoes (Lugauer & Wetzl, 2018). Off-resonance also shifts position. A spin precessing $\Delta f$ away from the frequency the encoding assumes picks up extra phase that grows with time, which a Cartesian reconstruction reads as a displacement of $\Delta f$ divided by the bandwidth per pixel, along the readout; in EPI, whose phase-encoding direction has a very low bandwidth per pixel, the displacement lies along that direction instead. This is our derivation from the signal equation. West et al. (2013) name EPI readout distortion and a larger water–fat shift (fat signal displaced because fat protons precess at a slightly different frequency) among the differences between their 3.0 T and 1.5 T data.
• Intensity non-uniformity. Receive coils are most sensitive near themselves (Zbontar et al., 2018), and transmit-field (B1) inhomogeneity is larger at 3.0 T than at 1.5 T (West et al., 2013). The result is a smooth shading across the image that has nothing to do with tissue, and that segmentation by intensity has to correct; West et al. note that normalization strategies or ad hoc filters may be needed.
• Noise. Because magnitude noise is Rician, image denoising methods designed for Gaussian noise are biased in low-signal regions unless they account for it (our inference).

Processing the Data

For vision work, the sensor’s properties translate into a few habits:

• Read the geometry, not just the pixels. Use the DICOM image position, orientation and spacing, or the NIfTI affine, to place voxels in millimetres. Resample anisotropic stacks deliberately, and measure distances and volumes in physical units rather than voxels.
• Normalize intensities. Since intensities depend on protocol, field strength, coil and position in the image, a model trained on raw intensities from one scanner sees different numbers from another; pipelines normalize per image and correct the shading before learning across scans. Quantitative maps of $T_1$, $T_2$ or PD are an alternative where they are acquired, though they still depend on field strength (West et al., 2013).
• Treat contrasts as different inputs. T1-, T2- and PD-weighted series have different intensity orderings; a pipeline can register them and stack them as channels, and image registration between contrasts uses similarity measures, such as mutual information, that do not assume equal intensities.
• Segment in 3D with spacing in mind. Tissue, organ and lesion delineation is a central vision task on MR volumes, often with U-Net-style networks for semantic segmentation; medical image segmentation adapts them to anisotropic voxels, class imbalance and small structures.
• Reconstruct from raw data when the task is reconstruction. Since archived images do not contain the measured data, learned reconstruction, super-resolution and artifact correction are best trained and evaluated on raw multi-coil k-space, with care for hallucinated detail (Zbontar et al., 2018; Knoll et al., 2020).

Strengths

• No ionizing radiation, in contrast to radiography and CT (Nobel press release, 2003; Lugauer & Wetzl, 2018).
• Soft-tissue contrast that can be tuned: the same anatomy can be imaged with $T_1$, $T_2$, PD, diffusion, flow and functional contrasts (Lugauer & Wetzl, 2018).
• Oblique slice orientations and true 3D acquisition, set by combining the three gradient coils (Lugauer & Wetzl, 2018).
• Access to the raw complex measurements, which lets reconstruction be posed as an inverse problem and improved by parallel imaging, compressed sensing and learned priors.

Limitations

• Slow acquisition. The physics of acquisition makes MRI inherently slower than X-ray, CT or ultrasound (Knoll et al., 2020), which brings motion artifacts and limits throughput (Lugauer & Wetzl, 2018).
• No calibrated intensity scale, which complicates thresholds, transfer between sites and comparisons over time.
• Scarce raw data. Public raw k-space datasets are few and specific to certain exams (Knoll et al., 2020).

Typical Applications

The 2003 Nobel press release describes MRI as used to examine almost all organs, and as especially valuable for detailed imaging of the brain and spinal cord. For vision and machine learning, the main uses are:

• Segmentation and quantification of tissues, organs and lesions, including brain tissue classes, which West et al. (2013) segmented from quantitative $R_1$, $R_2$ and PD maps.
• Registration across time, contrasts and modalities, including MR to CT, and to anatomical atlases (medical image registration).
• Reconstruction from undersampled data, the task fastMRI was built for (Zbontar et al., 2018).
• Image enhancement (denoising, super-resolution, artifact correction), and classification and detection on MR images and volumes.

History

Nuclear magnetic resonance in protons was demonstrated in 1946 by Felix Bloch and Edward Purcell, who received the 1952 Nobel Prize in Physics; in the following decades it was used mainly to study chemical structure (Nobel press release, 2003). In his Nobel lecture, Paul Lauterbur recounts that in 1971 Raymond Damadian observed that some malignant tissue from implanted rat tumours had longer NMR relaxation times than many normal tissues, and that in September 1971, after watching such tissue measurements, he conceived of locating NMR signals by applying linear field gradients in many directions and reconstructing the object from the resulting projections. He tested the idea on capillaries of ordinary water inside tubes of heavy water, coined the word “zeugmatography” for the new kind of image, and after an initial rejection published it in Nature in 1973 (Lauterbur, 1973; Lauterbur, 2003). The press release credits Lauterbur with the discovery that gradients make two-dimensional images possible.

At Nottingham, Peter Mansfield and Peter Grannell began in 1972 to use gradients to reveal structure in solids; their 1973 papers, the journal one titled “NMR ‘diffraction’ in solids?”, emphasized a Fourier transform approach (Mansfield & Grannell, 1973; Mansfield, 2003). The group then defined thin slices by selective irradiation. Early line-scanning methods took 10–20 minutes for a 64 × 64 image; echo-planar imaging in 1977 reduced acquisition to tens of milliseconds, and actively screened gradient coils later made such fast switching practical inside superconducting magnets (Mansfield, 2003). The first MRI equipment reached health care at the beginning of the 1980s, and by 2002 about 22,000 scanners were in use worldwide, performing more than 60 million examinations a year (Nobel press release, 2003). Lauterbur and Mansfield shared the 2003 Nobel Prize in Physiology or Medicine “for their discoveries concerning magnetic resonance imaging”.

Since then, speed has driven much of the engineering: parallel imaging (SENSE, Pruessmann et al., 1999; GRAPPA, Griswold et al., 2002), compressed sensing (Lustig et al., 2007) and learned reconstruction on shared raw data such as fastMRI (Zbontar et al., 2018).

Comparison with Other Sensors

• Computed tomography measures X-ray attenuation along lines; its voxels are calibrated in Hounsfield units, so thresholds transfer between scans far better than MR intensities do, though only approximately across scanners and tube voltages, since CT numbers depend on the effective beam energy and scanner calibration. CT is faster (Knoll et al., 2020) but uses ionizing radiation; MRI offers many soft-tissue contrasts per exam. Both relate to the Fourier domain: CT through the central slice theorem, MRI by sampling k-space directly.
• Cameras record intensities at pixels, so their raw data are already an image; MRI’s raw data are Fourier coefficients, and many of its artifacts (wrap-around, ringing, ghosts) have no direct counterpart in photographs.
• Ultrasound, PET and radiography are other medical modalities with their own physics; ultrasound and radiography acquire faster (Knoll et al., 2020), while MRI gives cross-sectional soft-tissue detail without radiation.