Medical Image Reconstruction

Medical image reconstruction computes an image from measurements that are not themselves an image. A CT scanner records X-ray attenuation along lines through the body, and an MRI scanner records samples of the Fourier transform of the magnetization; in both cases a computer turns those numbers into the cross-sections that radiologists read and that vision models take as input. How many measurements are taken, how noisy they are and what the algorithm assumes about the image together decide resolution, noise texture and artifacts.

The task matters beyond image quality. Fewer measurements mean a lower radiation dose in CT and a shorter scan in MRI, so much of the research in both modalities is about reconstructing acceptable images from less data. The methods range from closed-form inversion formulas through iterative optimization to trained networks, and the last of these raise a question the earlier ones did not: whether a plausible-looking structure in the output was actually measured.

Problem Definition

Adler and Öktem (2018) formalize an inverse problem as recovering a signal $f_\text{true}$ from data

$$g = T(f_\text{true}) + \delta g,$$

where the forward operator $T$ models how a signal gives rise to data in the absence of noise, and $\delta g$ is one sample of the noise. They name the components of model-driven reconstruction as a forward model derived from the physics, the statistical properties of the data, and prior information about the image.

• Forward model. In CT, after the logarithm, $T$ is the ray (Radon) transform $P$ of the attenuation map. Adler and Öktem also write the physically more correct pre-log model from the Beer–Lambert law, $T(f) = e^{-\mu P(f)}$ with $\mu$ a mass attenuation coefficient, and note that applying $-\log(\cdot)/\mu$ to data and operator linearizes it back to the ray transform. In MRI, Knoll et al. (2020) write the data of coil $j$ as $f_j = F_\Omega C_j u + n_j$, with $u$ the image, $C_j$ a diagonal matrix holding the coil’s sensitivity, $F_\Omega$ the Fourier transform sampled at the acquired locations $\Omega$ and $n_j$ measurement noise; stacking the coils gives $f = Eu + n$ with an encoding operator $E$.
• Noise model. CT counts photons. Kak and Slaney (1988, §5.2.2) model each detector count as a Poisson variable, which makes the variance of the measured log projection about $1/\bar N$ for a mean count $\bar N$, so rays through dense material, which receive few photons, are the noisiest. MRI noise enters k-space additively as $n$ in the model above.
• Prior. Knowledge of what images look like, used to choose among the many images consistent with noisy or incomplete data: smoothness, sparsity in a transform, or a model learned from examples.

Fitting the data alone means minimizing the negative data log-likelihood $L(T(f), g)$. Adler and Öktem point out that for typical forward operators this problem is ill-posed: small changes in the data produce large changes in the reconstruction, and the maximum-likelihood solution over-fits the noise. Variational regularization, which they note is also called model-based iterative reconstruction in medical imaging, adds a regularization functional $S$ that penalizes unlikely images:

$$\min_{f}\ \Big[\, L\big(T(f), g\big) + \lambda\, S(f) \,\Big],$$

with $\lambda \ge 0$ weighing the prior against the fit to the data. In MRI the plain data fit is the least-squares problem $\min_u \tfrac{1}{2}\|Eu - f\|_2^2$, which Knoll et al. (2020) describe as ill-posed because the coil elements do not measure completely independent information, so its inversion amplifies noise; stopping an iterative solver early already acts as a form of regularization. Least squares is the negative log-likelihood for Gaussian noise, and for Poisson counts the corresponding term is the Poisson log-likelihood (our reading of the general formula). Every method below can be read as a choice of these three components, made explicitly or implicitly.

Input and Output

The input is the scanner’s raw data with its acquisition geometry. For CT this is a sinogram, one projection per view angle, either as photon counts (pre-log) or as line integrals (post-log); for MRI it is complex k-space for each receive coil, the sampling mask, and coil sensitivity maps estimated from a calibration scan or the fully sampled centre of k-space (Knoll et al., 2020). The output is an image or volume: attenuation coefficients that are then calibrated to Hounsfield units in CT, and a complex image per coil, usually combined into a magnitude image, in MRI.

Learned methods also need training pairs of measurements and reference images. As Magnetic Resonance Imaging explains, archived clinical images are not a substitute for raw data, because vendor reconstruction has already processed them (Zbontar et al., 2018); simulating measurements from them leaves out the real noise and acquisition physics.

Variants

• Sparse-view CT. Fewer projections than the sampling rule asks for, the setting FBPConvNet was trained for (Jin et al., 2017).
• Low-dose CT. A full set of views with fewer photons per ray, hence noisier data, as in Adler and Öktem’s (2018) simulations and the AAPM Low Dose CT Grand Challenge described below.
• Accelerated MRI. Fewer k-space lines than full sampling needs, with one coil or many. The fastMRI dataset defines a multi-coil task and a single-coil task with simulated single-coil data, the latter to lower the barrier for researchers new to MRI and because single-coil machines are still in use (Zbontar et al., 2018).
• Image-domain or data-domain. A method can work on an initial reconstruction, as image denoising does, or on the measurements themselves; the AAPM challenge accepted both, either denoising filtered back-projection images or reconstructing from projection data with a noise map (AAPM, Low Dose CT Grand Challenge).

Challenges

• Ill-posedness. Inverting the forward model amplifies noise, and with undersampled data many images fit the measurements equally well (Adler & Öktem, 2018; Knoll et al., 2020).
• Model mismatch. The linear models above leave out beam hardening, scatter and motion in CT, and field inhomogeneity and motion in MRI; the sensor articles describe the resulting artifacts.
• No true ground truth. The reference image is itself a reconstruction: a full-dose or full-view FBP image in CT, a fully sampled reconstruction in MRI. Some acquisitions, such as imaging of the moving heart or diffusion imaging of the brain, cannot be acquired fully sampled because of constraints on spatio-temporal resolution, so supervised training has no fully sampled labels for them (Knoll et al., 2020).
• Scarce raw data. Public raw k-space data were, in Knoll et al.’s survey, on the order of hundreds of datasets, limited to specific exams.
• Generalization. A learned reconstruction is trained for one scanner, protocol, sampling pattern and acceleration, and may not transfer to others, as the open problems below describe.

Approaches

Analytic reconstruction: filtered back-projection

The analytic route starts from the central slice theorem, which Computed Tomography states: the one-dimensional Fourier transform $S_\theta(w)$ of a parallel projection at angle $\theta$ equals the object’s two-dimensional Fourier transform $F(u, v)$ along a radial line at the same angle, $S_\theta(w) = F(w\cos\theta,\; w\sin\theta)$ (Kak & Slaney, 1988, §3.2). The proof is a change of variables: integrating $f$ along lines of constant $t$ and then transforming in $t$ is the same as evaluating the 2D transform on the rotated frequency axis (Kak & Slaney, 1988, ch. 3, Eqs. 14–18). It follows that projections at enough angles fill the 2D frequency plane along radial lines, and an inverse 2D Fourier transform recovers the object. Direct Fourier inversion works, but it needs interpolation from the polar grid of radial samples to a Cartesian one, and the radial samples thin out at high frequencies, so high frequencies are estimated less accurately (Kak & Slaney, 1988, §3.2).

Rewriting the inverse 2D Fourier transform in polar coordinates gives filtered back-projection, which Kak and Slaney describe as used in almost all applications of straight-ray tomography in 1988. The Jacobian of the polar change of variables, $|w|$, ends up multiplying each projection’s spectrum (Kak & Slaney, 1988, ch. 3, Eqs. 24–35):

$$f(x, y) = \int_0^{\pi} Q_\theta(x\cos\theta + y\sin\theta)\, d\theta, \qquad Q_\theta(t) = \int_{-\infty}^{\infty} S_\theta(w)\, |w|\, e^{j2\pi w t}\, dw.$$

Each projection is filtered with the ramp $|w|$, and the filtered projection $Q_\theta$ is then back-projected: its value at $t$ is added to every pixel on the line $x\cos\theta + y\sin\theta = t$, smeared back across the image. Summing over angles completes the reconstruction. Kak and Slaney (1988, §3.3.1) give the intuition: each projection supplies one radial line of the object’s spectrum but has to stand for a pie-shaped wedge of it, and the wedge’s width at frequency $w$ is proportional to $|w|$. It follows that back-projection without the filter over-weights the low frequencies, where all the lines meet, and returns a blurred image.

The ramp grows with frequency, and so does its amplification of noise. Kak and Slaney note that multiplying the ramp by a smooth window, such as a Hamming window, usually gives better results in the presence of noise, because the high frequencies it suppresses are mostly noise (Kak & Slaney, 1988, ch. 3, Eq. 43). The choice of window is the CT version of the trade-off between resolution and noise that Hounsfield (1979) described: filtering out high frequencies reduces noise at the cost of resolution. Implementation details also matter: filtering with the FFT produces a slight dishing and a DC shift, the first from the circular convolution, which zero-padding each projection removes, and the second from zeroing the whole lowest frequency cell rather than the single frequency $w = 0$ (Kak & Slaney, 1988, §3.3.3). In their bibliographic notes, Kak and Slaney credit the idea of filtered back-projection to Bracewell and Riddle (1967) and, independently, to Ramachandran and Lakshminarayanan (1971).

How many views are needed follows from the slice theorem. For parallel projections spread uniformly over 180°, with point sources and detectors, each projection of $N_\text{ray}$ samples contributes one radial line of the frequency disc. Requiring the spacing between neighbouring radial lines at the edge of the disc to match the spacing of samples along each line gives

$$\frac{M_\text{proj}}{N_\text{ray}} \approx \frac{\pi}{2}$$

for $M_\text{proj}$ projections, which Kak and Slaney summarize as needing roughly as many projections as rays per projection (Kak & Slaney, 1988, ch. 5, Eqs. 10–15); in real systems the finite detector aperture and focal spot smooth the projections and soften the aliasing (§5.1.2). Fewer views leave gaps between the radial lines at high frequencies, and the reconstruction shows thin streaks: the filtered back-projection of a point from a few views is a star-shaped pattern (Kak & Slaney, 1988, §5.1.1). This is aliasing in the angular direction.

In terms of the problem definition, filtered back-projection inverts the linear forward model exactly for complete, noiseless data and has no explicit noise model or prior; the window is its only control over noise, applied equally to every ray (our reading). Kak and Slaney’s analysis of noise in FBP images shows the cost: with independent Poisson counts, each ray adds to a pixel’s variance in proportion to the squared filter and inversely to the ray’s mean count (1988, §5.2.2, Eq. 61). The MRI counterpart of filtered back-projection is the inverse FFT of fully sampled Cartesian k-space, which Magnetic Resonance Imaging uses; applied to undersampled data with the missing samples set to zero, it aliases.

Example: reconstruction filters under photon noise

The example simulates a parallel-beam scanner with scikit-image 0.25.2, NumPy 2.3.5 and SciPy 1.16.3, using the Shepp–Logan phantom bundled with scikit-image (skimage.data.shepp_logan_phantom, loaded from a file inside the package; Computed Tomography identifies it as the higher-contrast variant). The phantom’s values have no physical unit. To simulate photon counting, the code adopts an illustrative scale of its own: pixel values are read as attenuation coefficients in cm⁻¹ over a 25.6 cm field, so that the value 0.2 plays the role of water. It then forms line integrals with radon, draws Poisson counts from $10^5$ incident photons per ray, takes the logarithm, and reconstructs with iradon and each of its filters (ramp, the default, shepp-logan, hann, and None for unfiltered back-projection).

import numpy as np
from skimage.data import shepp_logan_phantom
from skimage.transform import radon, iradon

phantom = shepp_logan_phantom()            # 400 x 400, values in [0, 1]

# Illustrative physical scale (our assumption, not part of the phantom):
# read pixel values as attenuation in 1/cm and give the 400-pixel field 25.6 cm.
dx = 25.6 / 400                            # cm per pixel
mu = phantom                               # 1/cm; 0.2 plays the role of water

theta = np.linspace(0.0, 180.0, 360, endpoint=False)
p = radon(mu, theta=theta) * dx            # line integrals of mu (dimensionless)
print("sinogram (detector bins, views):", p.shape,
      " largest line integral: %.2f" % p.max())

# Photon counting: N = N0 exp(-p), Poisson noise, then the log transform.
rng = np.random.default_rng(0)
N0 = 1e5                                   # incident photons per ray
counts = rng.poisson(N0 * np.exp(-p))
p_noisy = -np.log(np.maximum(counts, 1) / N0)

n = phantom.shape[0]
yy, xx = np.mgrid[:n, :n] - n // 2
inside = yy**2 + xx**2 <= (n // 2) ** 2
flat = (slice(276, 316), slice(238, 278))  # 40 x 40 patch of value 0.2

def rmse(rec):
    # RMSE after the least-squares gain, so that unfiltered backprojection,
    # which has no meaningful scale, can be compared on shape alone
    a = (rec[inside] @ mu[inside]) / (rec[inside] @ rec[inside])
    return np.sqrt(np.mean((a * rec - mu)[inside] ** 2)), a

print("filter       RMSE noiseless  RMSE noisy  noise SD (flat patch)")
for name in ["ramp", "shepp-logan", "hann", None]:
    clean = iradon(p, theta=theta, filter_name=name) / dx
    noisy = iradon(p_noisy, theta=theta, filter_name=name) / dx
    (e0, _), (e1, a) = rmse(clean), rmse(noisy)
    print("%-12s %14.4f %11.4f %22.4f" % (name, e0, e1, (a * noisy)[flat].std()))

sinogram (detector bins, views): (400, 360)  largest line integral: 6.80
filter       RMSE noiseless  RMSE noisy  noise SD (flat patch)
ramp                 0.0345      0.0368                 0.0156
shepp-logan          0.0378      0.0393                 0.0124
hann                 0.0511      0.0514                 0.0054
None                 0.2141      0.2141                 0.0031

The longest ray loses all but about $e^{-6.8} \approx 0.1\%$ of its photons. Going from the ramp to the Hann window, the noise in a uniform patch falls by a factor of about three, while the error on noiseless data rises by half: the window trades resolution, which matters at the phantom’s sharp edges, for lower noise. Most of the RMSE here comes from edges rather than noise, which is why the noisy and noiseless columns are close. Unfiltered back-projection has the lowest noise and by far the largest error even after the best scaling, because it is the blurred image that the ramp exists to correct. Computed Tomography continues the simulation with fewer views and a dense insert.

Algebraic and statistical iterative reconstruction

A different approach discretizes the slice into $N$ pixels $f_j$ and writes each measured ray sum $p_i$ as a weighted sum of the pixels it crosses, $\sum_j w_{ij} f_j = p_i$, where $w_{ij}$ is the fraction of pixel $j$ intercepted by ray $i$ (Kak & Slaney, 1988, ch. 7, Eq. 1). For a 256 × 256 image the system has about 65,000 unknowns and as many equations, far too many for direct inversion, but it is sparse. Iterative methods based on Kaczmarz’s method of projections solve it by projecting the current estimate onto the hyperplane of one equation after another (Kak & Slaney, 1988, §7.1). The algebraic reconstruction technique (ART) updates after every ray; SIRT averages the corrections over all rays before updating, giving smoother images at slower convergence; and SART, introduced by Andersen and Kak in 1984, updates per projection with bilinear basis functions and gives good images in one iteration (Kak & Slaney, 1988, §7.2–7.4). Kak and Slaney credit the earliest expositions of algebraic reconstruction to Gordon, Herman and colleagues from 1970, and note that Hounsfield’s 1972 patent already ordered projections so that successive ones were well separated, to speed convergence (§7.5).

In their 1988 assessment, algebraic methods lacked the accuracy and speed of filtered back-projection for medical use with complete data, and earned their place where projections are few, unevenly distributed, or follow curved paths (Kak & Slaney, 1988, ch. 7). Statistical, model-based reconstruction, in the framework of Adler and Öktem (2018), adds the other two components of the problem definition to the forward model: a likelihood that matches the noise, such as the Poisson model of photon counts, and a regularizer, as in the variational formulation above. Jin et al. (2017) call regularized iterative algorithms the standard approach to ill-posed inverse problems, and total-variation regularization is the iterative baseline against which both learned CT methods below are compared.

Parallel imaging

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; the methods below remove the aliasing with extra information.

Each coil of an array sees the object multiplied by its own sensitivity map $c_j(\mathbf{r})$, so with $n_c$ coils the measurements are $f_j = F_\Omega C_j u + n_j$ for $j = 1, \dots, n_c$, where $F_\Omega$ samples the Fourier transform at the acquired locations $\Omega$ (Knoll et al., 2020). Parallel imaging uses this redundancy to fill in skipped lines. Knoll et al. sort the methods into two families. SENSE (Pruessmann et al., 1999) works in image space: with every $R$-th line acquired, each pixel of a coil’s aliased image is a sensitivity-weighted sum of $R$ pixels of the true image, and with known sensitivities a small linear system per pixel group unfolds them (Lugauer & Wetzl, 2018). GRAPPA (Griswold et al., 2002) works in k-space: it estimates each missing line of each coil as a linear combination of acquired neighbouring lines from all coils, with convolution kernels fitted to a fully sampled autocalibration region at the centre of k-space (Knoll et al., 2020). Inverting the coil system amplifies noise where coil sensitivities overlap, which the SENSE paper describes by the g-factor, and this limits practical acceleration (Knoll et al., 2020; Zbontar et al., 2018). Knoll et al. write that essentially all clinical scanners include parallel imaging and that it is the default for many protocols.

Compressed sensing

Lustig, Donoho and Pauly (2007) adapted compressed sensing to MRI with three requirements: the image must have a sparse representation in a known transform (such as wavelets or finite differences), the aliasing from undersampling must be incoherent, that is noise-like, in that transform domain, and a nonlinear reconstruction must enforce both sparsity and consistency with the data. They reconstructed by solving

$$\min_{m}\ \|\Psi m\|_1 \quad \text{subject to} \quad \|F_u m - y\|_2 < \epsilon ,$$

with $\Psi$ the sparsifying transform, $F_u$ the undersampled Fourier operator, $y$ the measured k-space data and $\epsilon$ set below the expected noise level, often adding a total-variation term. Regular undersampling produces coherent copies, so they drew phase encodes pseudo-randomly with a density that falls with distance from the k-space centre, where most image energy lies, and measured incoherence by the point spread function’s largest sidelobe-to-peak ratio. In their simulations, reconstructions degraded with acceleration by losing low-contrast features rather than resolution. Knoll et al. (2020) describe the general form, a data-fidelity term plus weighted regularizers, as the basis of almost all modern MRI reconstruction, and combined parallel imaging and compressed sensing as common. In the terms of the problem definition, the $\ell_1$ norm of the transform coefficients is the prior and the constraint is the data fit.

Example: incoherent undersampling

The second example follows Magnetic Resonance Imaging’s simulation, with NumPy 2.3.5 and scikit-image 0.25.2: the k-space of the bundled phantom by a centred, orthonormal FFT (norm="ortho"), axis 0 the phase-encoding direction. It compares the uniform fourfold mask, whose zero-filled image is folded copies of the object, with variable-density random masks of the same number of lines and a fully sampled centre, scoring the zero-filled reconstruction by its normalized error and each mask by the largest sidelobe of its point spread function. As in the MRI article, simulating and reconstructing on the same pixel grid can look better than real acquisitions.

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)

# Variable-density random undersampling at about R=4, with a fully sampled centre
# (8% of the lines, as in the fastMRI masks); density falls with |ky|.
def random_mask(seed, R=4, centre_frac=0.08):
    rng = np.random.default_rng(seed)
    m = np.zeros(N, bool)
    c = int(round(centre_frac * N))
    m[N//2 - c//2 : N//2 + c//2] = True
    p = (1 - np.abs(np.arange(N) - N//2) / (N/2))**2
    p[m] = 0
    m[rng.choice(N, N//R - m.sum(), replace=False, p=p / p.sum())] = True
    return m

# Incoherence: the point spread function of a mask along ky and its largest
# sidelobe-to-peak ratio (the measure used by Lustig et al., 2007), counted at
# offsets of 5 or more pixels (nearer offsets belong to the main lobe).
def sidelobe(m):
    psf = np.abs(np.fft.ifft(np.fft.ifftshift(m.astype(float))))
    return psf[5:N-4].max() / psf[0]

uniform4 = np.arange(N) % 4 == N//2 % 4
zf = to_image(K * uniform4[:, None])
print(f"uniform R=4: {uniform4.sum()} lines, NRMSE {nrmse(np.abs(zf), img):.3f}, "
      f"max sidelobe-to-peak {sidelobe(uniform4):.2f}")
mask = random_mask(seed=0)
zf = to_image(K * mask[:, None])
print(f"random R~4 (seed 0): {mask.sum()} lines, NRMSE {nrmse(np.abs(zf), img):.3f}, "
      f"max sidelobe-to-peak {sidelobe(mask):.2f}")
res = np.array([(nrmse(np.abs(to_image(K * m[:, None])), img), sidelobe(m))
                for m in map(random_mask, range(20))])
print(f"random R~4, seeds 0-19: NRMSE {res[:,0].min():.3f}-{res[:,0].max():.3f}, "
      f"max sidelobe-to-peak {res[:,1].min():.2f}-{res[:,1].max():.2f}")

uniform R=4: 100 lines, NRMSE 0.712, max sidelobe-to-peak 1.00
random R~4 (seed 0): 100 lines, NRMSE 0.279, max sidelobe-to-peak 0.19
random R~4, seeds 0-19: NRMSE 0.244-0.287, max sidelobe-to-peak 0.15-0.21

The random mask uses the same number of lines, but its point spread function has no sidelobe higher than about a fifth of the peak over all 20 seeds, against sidelobes as tall as the peak for the uniform mask. Its interference is spread into a noise-like background instead of a coherent ghost, the incoherence compressed sensing needs. Part of its lower error comes simply from sampling the centre fully. Zero filling removes none of the interference; that takes the nonlinear reconstruction above.

Learned reconstruction

Adler and Öktem (2018) sort learned tomographic reconstruction into three directions. The first, learned post-processing or learned denoisers, applies a trained network to an approximate inverse such as filtered back-projection; it is easy to implement, since the initial reconstruction can be computed before training. The second learns a regularizer and uses it in a classical variational scheme. The third, learned iterative schemes, resembles classical optimization but lets a network choose each update from the last iterate and the results of applying the forward operator and its adjoint.

Jin et al. observed that unrolled iterative methods resemble a convolutional neural network when the forward model’s normal operator is a convolution, a class they note includes X-ray CT, and proposed FBPConvNet: filtered back-projection followed by a U-Net with a residual connection, trained to remove the artifacts of sparse-view reconstruction (Jin et al., 2017). In the arXiv version, they used synthetic ellipse images, parallel-beam sinograms simulated from 500 Mayo Clinic images of the Low-dose CT Grand Challenge, and real rat-brain sinograms from a synchrotron beamline. Taking the FBP of the full 1,000-view (rat: 721-view) sinogram as ground truth and subsampling the views, they reported higher SNR than total-variation-regularized iterative reconstruction on the Mayo images at 143 and 50 views and on the rat data at 145 and 51 views (Table III; the text gives 52), while TV did better on the ellipses.

Applied to filtered back-projection images of low-dose scans, learned post-processing is low-dose CT denoising: a network trained on pairs of noisy FBP images and low-noise reference images learns a mapping between them (Adler & Öktem, 2018), the image-domain route that the AAPM Low Dose CT Grand Challenge accepted alongside reconstruction from projection data. Generic image denoising methods apply, with a caveat: CT noise concentrates at high frequencies (Hounsfield, 1979) and, by Kak and Slaney’s analysis above, its level varies with the photon counts of the rays through each pixel.

Adler and Öktem instead unrolled a primal–dual optimization algorithm for a fixed number of iterations and replaced its proximal operators with small CNNs, so the network applies the forward projector and its adjoint inside its layers and works from the raw sinogram, without an initial FBP (Adler & Öktem, 2018). On simulated low-dose fan-beam data from the same Mayo Clinic scans (nine patients for training, one held out for evaluation, the full-dose image as reference), with Poisson noise for $10^4$ incident photons per pixel, the arXiv version reports a gain of 6.6 dB in PSNR over total-variation reconstruction and 2.2 dB over learned post-processing of FBP, with higher SSIM. Both papers measure agreement with reference images by image-quality metrics, on simulated or retrospectively subsampled data; neither reports a reader study.

In MRI, deep networks now take the place of hand-designed regularizers. Knoll et al. (2020) review image-space methods that unroll an iterative reconstruction into a network of a fixed number of stages, each combining a learned regularizer with a data-consistency step that uses the coil sensitivities and sampling mask, and k-space methods that learn the interpolation of missing lines, some trained only on a scan’s own calibration data. The variational network of Hammernik et al. (2018) follows the image-space design. In the arXiv version, it embeds a compressed-sensing reconstruction posed as a variational model in ten steps of unrolled gradient descent, and learns the filter kernels and activation functions of the regularizer and the data-term weight of every step. The authors scanned 20 patients for each of five sequences of a clinical knee protocol on a 3 T scanner, without acceleration, and undersampled the data retrospectively. They trained one network per contrast, sampling pattern and acceleration on 200 slices from 10 patients, and tested on 10 other patients whose scans included pathology absent from the training set, with the fully sampled, coil-combined reconstruction as the reference. Against CG-SENSE and a combined parallel imaging and compressed sensing reconstruction with total generalized variation, each tuned by grid search to minimize error on the reference, they reported lower mean squared error and higher SSIM in all cases at accelerations of 3 and 4, for five sequences with regular sampling and two with variable-density random sampling (Table 1). Random sampling helped at acceleration 4 and not at 3.

Datasets and Benchmarks

• fastMRI. The fastMRI dataset of the NYU School of Medicine and Facebook AI Research released raw multi-coil k-space with retrospective undersampling masks: for test data, the lowest-frequency lines are kept (8% of the lines at fourfold acceleration, 4% at eightfold) and the rest are chosen at random for knee and equidistantly for brain (Zbontar et al., 2018). As of October 2026, access requires accepting the fastMRI Data Sharing Agreement and submitting an application, for internal research or educational use (fastMRI website). Its organizers ran public challenges on it; in the 2020 challenge, participants received data from 7,299 clinical brain scans, with the fully sampled data of 894 held back for evaluation, and a new Transfer track tested models on data from vendors outside the main fastMRI dataset (Muckley et al., 2021).
• AAPM Low Dose CT Grand Challenge (2016). Mayo Clinic provided contrast-enhanced abdominal CT scans with Poisson noise inserted into the projection data to simulate quarter-dose acquisitions: ten training cases with full-dose and quarter-dose data, and twenty test cases with quarter-dose data only. Participants returned denoised or reconstructed images, which radiologists read for the detection of liver lesions (AAPM, Low Dose CT Grand Challenge). The Mayo images that Jin et al. and Adler and Öktem used came from this challenge.
• LDCT-and-Projection-data. A larger library of CT projection and image data from Mayo Clinic, curated by The Cancer Imaging Archive: 99 head, 100 chest and 100 abdomen scans from two vendors, each with projection data at the routine dose and a second set at a simulated lower dose (25% of the routine dose for head and abdomen, 10% for chest). As of October 2026 the collection page lists the chest, abdomen (“liver”) and phantom images under CC BY 4.0 and the head scans under the NIH Controlled Data Access Policy; its download table lists 99 chest subjects against the 100 chest scans of its summary (McCollough et al., 2020).
• Simulation. Phantoms and synthetic images with known ground truth remain common, as in Jin et al.’s and Adler and Öktem’s ellipse experiments and the examples above. Knoll et al. (2020) name mridata.org alongside fastMRI as public raw MRI data.

Evaluation Metrics

Reconstruction is evaluated against a reference standard, and in practice the reference is another reconstruction: the full-view FBP for Jin et al. (2017), the full-dose image for Adler and Öktem (2018), and the fully sampled, coil-combined image for Hammernik et al. (2018) and fastMRI, which uses the root sum of squares of the coil images for its multi-coil data. The fastMRI authors note that their fully sampled ground truth may in many cases be of higher quality than the clinical standard, for which full sampling is not routine (Zbontar et al., 2018). The reference’s own flaws enter the scores: in the 2020 fastMRI challenge, the noise in the background of root-sum-of-squares images was treated as ground truth and so influenced the quantitative ranking, which the organizers counterbalanced with radiologist evaluation (Muckley et al., 2021). In each of these studies the evaluation is retrospective: undersampling or dose reduction is simulated from complete acquisitions.

• Pixel-wise errors. fastMRI defines the normalized mean squared error, $\|\hat v - v\|_2^2 / \|v\|_2^2$, computed over whole volumes, and recommends it as the primary measure for its data, while warning that it tends to favour smoothness over sharpness and should be reported with other metrics (Zbontar et al., 2018). PSNR is a log-scaled form of the same squared error.
• Structural similarity. SSIM compares local means, variances and covariances. On fastMRI’s total-variation baseline, NMSE and PSNR were highly correlated and favoured stronger regularization than SSIM did, and stronger regularization gave smoother images that lacked the fine texture of the ground truth (Zbontar et al., 2018). In Adler and Öktem’s human-phantom table, U-Net post-processing of FBP beat total variation by 4.4 dB in PSNR, yet its SSIM was no better than TV’s. The metrics are not interchangeable, and the choice of training loss affects the properties of the trained network (Knoll et al., 2020).
• Reader evaluation. The AAPM challenge scored submissions by radiologists’ detection of liver lesions against a reference diagnosis, with points for correctly marked lesions and cases and penalties for false positives and negatives (AAPM, Low Dose CT Grand Challenge). The 2020 fastMRI challenge ranked submissions by SSIM, then sent the top three of each track to six radiologists, blinded to the teams, who ranked them on how well they depicted pathology, using the fully sampled image as reference, and rated artifacts, sharpness and contrast-to-noise ratio. Averaged over radiologists, the top-ranked submission was also the one with the highest SSIM, while second and third places in the fourfold and eightfold tracks were swapped relative to SSIM (Muckley et al., 2021).

Hallucinated structure can escape image-quality metrics. Knoll et al. (2020) warn that adversarially trained networks have the potential to produce images visually indistinguishable from fully sampled ones, with features that the acquired data do not support. The 2020 fastMRI challenge documented such features in submissions from all three tracks: a false vessel, a linear bright signal mimicking a cleft of cerebrospinal fluid, and a false sulcus or prominent vessel. The organizers describe methods that created artifact-mimics, and models that morphed an abnormality into a more normal-looking structure. These images had high SSIM scores, and the organizers wrote that, as of the challenge’s end, they knew of no automated, quantitative method for detecting such hallucinations (Muckley et al., 2021). Their earlier challenge had found that, while many submissions were of high quality at fourfold acceleration, all of them began to miss pathology at eightfold. It follows that an image can match its reference closely on average while adding or removing exactly the small structure a reader looks for; the 2020 challenge asked its radiologists to rank the depiction of pathology rather than general image quality.

Applications

Reconstruction sits upstream of every vision task on CT and MRI. Its main uses are lowering the dose of CT, by acquiring fewer views (Jin et al., 2017) or fewer photons per view, and shortening MRI acquisitions, where parallel imaging is already routine (Knoll et al., 2020). Its choices, from the CT filter to the MRI sampling pattern and the learned model, change the noise, texture and resolution of the images that medical image segmentation and other analyses receive, so it follows that a model trained on images from one reconstruction may see different inputs from another. Learned super-resolution and artifact correction share the same evaluation concerns when they produce detail the measurements did not contain.

Open Problems

• Failure modes and generalization. Knoll et al. (2020) list open issues: learned models are non-convex and their failure modes are less understood than those of convex methods; one study found that they generalized across image contrast but were sensitive to systematic changes in SNR; and it is open how specific trained models must be to anatomy, sequence, trajectory, acceleration, vendor, field strength and coil. In the 2020 fastMRI challenge, models trained only on fastMRI data from one vendor were tested in the Transfer track on data from two others, where SSIM was lowest and varied most; many participants struggled to adapt their models to one vendor’s data, which lacked the frequency oversampling of the training data, a difference the organizers observed can lower SSIM by 0.1–0.4 if no other adjustments are made (Muckley et al., 2021).
• Hallucination detection. The 2020 fastMRI organizers called automated stability and hallucination analysis a topic of great interest for future challenges (Muckley et al., 2021).
• Higher accelerations. Muckley et al. take the results of the 2019 challenge, and a later validation at fourfold that they cite, to suggest that the upper limit of 2D machine-learning reconstruction lay between fourfold and eightfold acceleration.
• Training without full sampling. Acquisitions that cannot be fully sampled call for scan-specific or unsupervised training, which Knoll et al. (2020) describe as ongoing research.
• Shared raw data. Since the training data become part of a model’s performance, methods trained on unshared data are almost impossible to compare (Knoll et al., 2020).