Computed Tomography: The Canonical Inverse Problem

Section Roadmap: CT as the Canonical Inverse Problem

Computed tomography is the cleanest example of a Fourier-based imaging inverse problem. The Radon transform has an analytical inverse; the Fourier Slice Theorem connects projections to the 2D spectrum; and filtered back-projection (FBP) provides a non-iterative reconstruction. We develop CT in full, then use it as the baseline against which to measure the additional difficulties of RF imaging: incomplete angular coverage, near-field propagation, and an ill-conditioned forward operator.

Definition:

The Radon Transform

Let f∈L1(R2)f \in L^1(\mathbb{R}^2) be a compactly supported object function. The Radon transform of ff is

Rf(θ,s)=∫−∞∞f(scos⁡θ−tsin⁡θ,  ssin⁡θ+tcos⁡θ) dt,\mathcal{R}f(\theta, s) = \int_{-\infty}^{\infty} f(s\cos\theta - t\sin\theta,\; s\sin\theta + t\cos\theta)\,dt,

where θ∈[0,π)\theta \in [0, \pi) is the projection angle and s∈Rs \in \mathbb{R} is the signed perpendicular distance from the origin to the integration line. The function g(θ,s)≜Rf(θ,s)g(\theta, s) \triangleq \mathcal{R}f(\theta, s) is the sinogram of ff.

Compact operator interpretation: R:L2(R2)→L2([0,π)×R)\mathcal{R}: L^2(\mathbb{R}^2) \to L^2([0,\pi) \times \mathbb{R}) is a bounded linear operator. Its adjoint R∗\mathcal{R}^* is the back-projection operator: smearing each projection back along its line of integration.

The term "sinogram" arises because a point scatterer at (x0,y0)(x_0, y_0) traces the curve s=x0cos⁡θ+y0sin⁡θs = x_0 \cos\theta + y_0 \sin\theta — a sinusoid in (θ,s)(\theta, s) space.

Sinogram

The 2D function g(θ,s)=Rf(θ,s)g(\theta, s) = \mathcal{R}f(\theta, s) collecting all line-integral projections of ff. Each column of the sinogram is one projection at a fixed angle θ\theta.

Related: The Radon Transform

Radon transform

The integral transform R\mathcal{R} mapping a function f(x,y)f(x,y) to its line integrals parameterized by angle θ\theta and offset ss. Introduced by Johann Radon in 1917.

Related: The Radon Transform, The Fourier Slice Theorem (Central Slice Theorem)

Theorem: The Fourier Slice Theorem (Central Slice Theorem)

Let f~(kx,ky)\tilde{f}(k_x, k_y) denote the 2D Fourier transform of ff. The 1D Fourier transform of the Radon projection at angle θ\theta equals a radial slice of f~\tilde{f} through the origin at angle θ\theta:

∫−∞∞Rf(θ,s) e−j2πνs ds=f~(νcos⁡θ,  νsin⁡θ).\int_{-\infty}^{\infty} \mathcal{R}f(\theta, s)\, e^{-j2\pi \nu s}\,ds = \tilde{f}(\nu\cos\theta,\;\nu\sin\theta).

Equivalently: each projection fills one radial line in the 2D Fourier plane of ff.

A projection "integrates out" one spatial dimension. Integration in the spatial domain corresponds to evaluation at zero frequency in the conjugate direction — which picks out a line in 2D frequency space.

Proof
Substitution into the 2D Fourier integral

Write f~(νcos⁡θ,νsin⁡θ)=∬f(x,y) e−j2πν(xcos⁡θ+ysin⁡θ) dx dy\tilde{f}(\nu\cos\theta, \nu\sin\theta) = \iint f(x,y)\,e^{-j2\pi\nu(x\cos\theta + y\sin\theta)}\,dx\,dy. Introduce rotated coordinates s=xcos⁡θ+ysin⁡θs = x\cos\theta + y\sin\theta and t=−xsin⁡θ+ycos⁡θt = -x\sin\theta + y\cos\theta.

Separate the double integral

The Jacobian of the rotation is unity, so

f~(νcos⁡θ,νsin⁡θ)=∫−∞∞e−j2πνs[∫−∞∞f(s,t) dt]⏟Rf(θ,s)ds.\tilde{f}(\nu\cos\theta, \nu\sin\theta) = \int_{-\infty}^{\infty} e^{-j2\pi\nu s} \underbrace{\left[\int_{-\infty}^{\infty} f(s,t)\,dt\right]}_{\mathcal{R}f(\theta,s)} ds.

The inner integral is exactly the Radon projection at angle θ\theta. The outer integral is its 1D Fourier transform. ■\blacksquare

,

Historical Note: Johann Radon and the Transform That Waited 55 Years

1917-1979

Johann Radon published the transform and its inversion formula in 1917, motivated by pure mathematics — the problem of recovering a function from its line integrals. The result was largely forgotten by applied scientists until Allan Cormack (1963) independently rediscovered the inversion for medical imaging, and Godfrey Hounsfield built the first clinical CT scanner in 1971. Cormack and Hounsfield shared the 1979 Nobel Prize in Physiology or Medicine. The 55-year gap between Radon's theorem and its practical realization is a cautionary tale about the lag between mathematical theory and engineering application — and a reminder that the inverse problems we solve today in RF imaging may find their "Hounsfield moment" in the future.

Definition:

Filtered Back-Projection (FBP)

The filtered back-projection formula reconstructs ff from its sinogram g(θ,s)=Rf(θ,s)g(\theta, s) = \mathcal{R}f(\theta, s):

f^FBP(x,y)=∫0π[∫−∞∞∣ν∣ g~(θ,ν) ej2πνs dν]s=xcos⁡θ+ysin⁡θdθ,\hat{f}_{\mathrm{FBP}}(x, y) = \int_0^{\pi} \left[\int_{-\infty}^{\infty} |\nu|\,\tilde{g}(\theta, \nu)\, e^{j2\pi\nu s}\,d\nu\right]_{s = x\cos\theta + y\sin\theta} d\theta,

where g~(θ,ν)\tilde{g}(\theta, \nu) is the 1D Fourier transform of g(θ,⋅)g(\theta, \cdot) and ∣ν∣|\nu| is the ramp filter (Ram-Lak filter).

Interpretation: (1) filter each projection with ∣ν∣|\nu| to compensate for the non-uniform radial density in polar Fourier coordinates; (2) back-project by smearing each filtered projection along its line of integration.

The ramp filter amplifies high frequencies and hence noise. Practical variants — Shepp-Logan, cosine, Hamming — taper ∣ν∣|\nu| at high frequencies, trading resolution for noise suppression.

Definition:

Algebraic Reconstruction Technique (ART)

Discretize the object on an N×NN \times N grid: f∈RN2\mathbf{f} \in \mathbb{R}^{N^2}. Each ray ii gives a linear equation aiTf=gi\mathbf{a}_i^T \mathbf{f} = g_i, forming the system ACTf=g\mathbf{A}_{\mathrm{CT}} \mathbf{f} = \mathbf{g}. The algebraic reconstruction technique (ART) solves this via Kaczmarz iteration:

f(k+1)=f(k)+gi(k)−ai(k)Tf(k)∥ai(k)∥2 ai(k),\mathbf{f}^{(k+1)} = \mathbf{f}^{(k)} + \frac{g_{i(k)} - \mathbf{a}_{i(k)}^T \mathbf{f}^{(k)}} {\|\mathbf{a}_{i(k)}\|^2}\,\mathbf{a}_{i(k)},

where i(k)=k mod Mi(k) = k \bmod M cycles through the MM ray equations. The simultaneous variant SIRT updates using all rays per iteration.

The Key Parallel: CT is Well-Conditioned, RF is Not

The Fourier Slice Theorem guarantees that NvN_v uniformly spaced projections over [0,π)[0, \pi) fill the Fourier plane on a star-shaped grid with known density. The condition number of ACT\mathcal{A}_{\mathrm{CT}} grows modestly with NvN_v — the inverse is stable and FBP provides it analytically.

In RF imaging (Ch 07--15), the forward operator A\mathbf{A} samples k-space on Ewald arcs whose extent depends on bandwidth, aperture, and the near-field geometry. The coverage is fundamentally incomplete: limited angular range, limited bandwidth, and non-uniform density. The condition number of A\mathbf{A} can be orders of magnitude larger than that of ACT\mathcal{A}_{\mathrm{CT}}, and no analytical inverse exists — we must resort to regularized methods (Ch 13--18).

This conditioning gap is the single most important structural difference between CT and RF imaging, and it explains why RF imaging benefits more from learned priors than CT does.

,

Radon Transform and FBP Reconstruction

Left: The Shepp-Logan phantom. Center: Its sinogram (Radon transform over [0,π)[0, \pi)). Right: FBP reconstruction using the selected filter.

Reducing the number of projections introduces streak artifacts that grow worse as the angular sampling becomes sparser. Compare the three filters: Ram-Lak (sharpest, noisiest), Shepp-Logan (moderate), and cosine (smoothest, most blurred).

Parameters
90
0

Definition:

Total Variation Regularized CT

For limited-view CT, total variation (TV) regularization recovers piecewise-constant objects by solving

f^=arg⁡min⁡f≥012∥ACTf−g∥22+λ∥f∥TV,\hat{\mathbf{f}} = \arg\min_{\mathbf{f} \geq 0} \frac{1}{2}\|\mathbf{A}_{\mathrm{CT}}\mathbf{f} - \mathbf{g}\|_2^2 + \lambda \|\mathbf{f}\|_{\mathrm{TV}},

where ∥f∥TV=∑i,j(∇xfi,j)2+(∇yfi,j)2\|\mathbf{f}\|_{\mathrm{TV}} = \sum_{i,j} \sqrt{(\nabla_x f_{i,j})^2 + (\nabla_y f_{i,j})^2} is the isotropic total variation. This is the same regularizer used in RF imaging (Ch 14, [?ch14:def-tv-regularization]), but applied to the Radon forward model instead of A\mathbf{A}.

Theorem: Exactness of FBP Under Full Angular Coverage

If f∈L1(R2)∩L2(R2)f \in L^1(\mathbb{R}^2) \cap L^2(\mathbb{R}^2) is compactly supported and projections Rf(θ,s)\mathcal{R}f(\theta, s) are available for all θ∈[0,π)\theta \in [0, \pi), then the FBP formula recovers ff exactly (in the L2L^2 sense):

f^FBP=f.\hat{f}_{\mathrm{FBP}} = f.

No regularization is needed. This is in stark contrast to the RF imaging problem, where even with "full" multi-view coverage, the bandwidth-limited Ewald arcs leave spectral gaps that prevent exact inversion.

Full angular coverage means every radial line in the Fourier plane is sampled. The ramp filter compensates for the 1/∣ν∣1/|\nu| density of radial sampling, giving uniform coverage. In RF imaging, the Ewald arcs never cover the full Fourier plane — there is always a "cone" of missing frequencies near the origin and beyond the bandwidth cutoff.

Proof
Polar Fourier inversion

Express the 2D inverse Fourier transform in polar coordinates (ν,θ)(\nu, \theta):

f(x,y)=∫0π∫−∞∞∣ν∣ f~(νcos⁡θ,νsin⁡θ) ej2πν(xcos⁡θ+ysin⁡θ) dν dθ.f(x,y) = \int_0^{\pi}\int_{-\infty}^{\infty} |\nu|\,\tilde{f}(\nu\cos\theta, \nu\sin\theta)\, e^{j2\pi\nu(x\cos\theta + y\sin\theta)}\,d\nu\,d\theta.

The Jacobian ∣J∣=∣ν∣|J| = |\nu| provides the ramp filter.

Apply the Fourier Slice Theorem

Substitute f~(νcos⁡θ,νsin⁡θ)=g~(θ,ν)\tilde{f}(\nu\cos\theta, \nu\sin\theta) = \tilde{g}(\theta, \nu) from the Fourier Slice Theorem. The result is precisely the FBP formula. ■\blacksquare

Example: Limited-View CT and the Emergence of Streak Artifacts

Consider a Shepp-Logan phantom (256×256256 \times 256) imaged with NvN_v projections uniformly distributed over [0,π)[0, \pi). For each Nv∈{180,36,18,9}N_v \in \{180, 36, 18, 9\}, compute the FBP reconstruction and its PSNR relative to the ground truth.

Task: Explain the origin of the streak artifacts and relate them to the Fourier Slice Theorem.

Solution
Fourier coverage analysis

With NvN_v projections, the Fourier plane is sampled on NvN_v radial lines. The angular gap between adjacent lines is Δθ=π/Nv\Delta\theta = \pi / N_v. At spatial frequency ∣ν∣|\nu|, the arc-length gap between samples is π∣ν∣/Nv\pi |\nu| / N_v. When this gap exceeds the Nyquist spacing, aliasing occurs.

Streak artifact mechanism

Missing angular data leaves wedge-shaped gaps in the Fourier plane. FBP, which assumes uniform angular density, misinterprets the missing data as zero, producing oscillatory artifacts (streaks) oriented along the missing projection angles.

Quantitative comparison
NvN_v Angular gap PSNR (dB) Visual quality
180 1.0 deg ∼\sim38 Near-perfect
36 5.0 deg ∼\sim28 Mild streaks
18 10.0 deg ∼\sim22 Severe streaks
9 20.0 deg ∼\sim17 Heavily degraded

Each halving of NvN_v removes half the Fourier data. The PSNR degrades roughly logarithmically, and streak severity grows with the angular gap.

Example: TV-Regularized CT vs FBP for Sparse Views

For the same Shepp-Logan phantom with Nv=18N_v = 18 projections, compare (a) FBP reconstruction and (b) TV-regularized reconstruction (solved via ADMM with 200 iterations).

Task: Explain why TV regularization suppresses streak artifacts while preserving edges.

Solution
FBP result

FBP with Nv=18N_v = 18 yields PSNR ≈22\approx 22 dB with prominent streak artifacts. The ramp filter amplifies the angular aliasing because it has no mechanism to distinguish signal from artifacts.

TV-regularized result

TV regularization penalizes total variation, favouring piecewise-constant solutions. The Shepp-Logan phantom is piecewise constant, so TV is a near-perfect prior for this object. ADMM converges to PSNR ≈34\approx 34 dB — a 12 dB improvement over FBP.

Why TV works here

The Shepp-Logan phantom has a sparse gradient: most pixels have zero gradient, with non-zero gradient only at the ellipse boundaries. TV regularization fills the missing Fourier data by enforcing gradient sparsity, effectively "interpolating" the missing angular information.

Parallel to RF imaging: The same TV prior improves RF reconstructions when the scene is piecewise constant (e.g., walls, metallic objects). But natural RF scenes are rarely piecewise constant, motivating learned priors (Ch 18, Ch 22).

Definition:

FBPConvNet — Learned Post-Processing for CT

FBPConvNet (Jin et al., 2017) applies a convolutional neural network (CNN) as a post-processor to the FBP output:

f^CNN=Gθ(f^FBP),\hat{f}_{\mathrm{CNN}} = \mathcal{G}_\theta(\hat{f}_{\mathrm{FBP}}),

where Gθ\mathcal{G}_\theta is a U-Net trained to map noisy/artifact- corrupted FBP images to clean ground truth. The key insight: FBP provides a good initial estimate even when incomplete, and the CNN only needs to remove structured artifacts rather than reconstruct from scratch.

This is the CT analogue of the matched-filter + U-Net pipeline for RF imaging (Ch 26, Chapter 26): replace FBP with c^BP\hat{\mathbf{c}}^{\text{BP}}, and the architecture transfers directly.

Quick Check

According to the Fourier Slice Theorem, a single projection at angle θ=0\theta = 0 provides information about which region of the 2D Fourier plane?

The horizontal line ky=0k_y = 0

The vertical line kx=0k_x = 0

A circle of radius ν\nu centered at the origin

The entire Fourier plane

Correction:
The horizontal line ky=0k_y = 0

At θ=0\theta = 0, the projection integrates along yy. The Fourier Slice Theorem states that the 1D FT of this projection equals f~(ν,0)\tilde{f}(\nu, 0) — the horizontal axis of the 2D Fourier plane.

Common Mistake: Do Not Assume the CT Forward Model Applies to RF

Mistake:

A common error is to treat the RF imaging forward operator as a Radon transform and apply FBP directly. This works only in the far-field, high-frequency limit where the Born approximation reduces to straight-line propagation.

At RF frequencies (λ∼\lambda \sim cm), the wavelength is comparable to the object size, diffraction is significant, and the forward operator is the scattering integral (Ch 05-06), not the Radon transform. FBP applied to RF data produces blurred, artifact-ridden images because it ignores the curvature of the Ewald arcs.

Correction:

Use the RF-specific forward model y=Ac+w\mathbf{y} = \mathbf{A}\mathbf{c} + \mathbf{w} from Ch 07 and the diffraction-tomography-aware FBP from Ch 15. For best results, use iterative or learned methods from Ch 13-18 that account for the exact structure of A\mathbf{A}.

Common Mistake: The Ramp Filter Amplifies Noise in FBP

Mistake:

The ramp filter ∣ν∣|\nu| in FBP grows linearly with frequency. If the sinogram data contain noise (inevitable in practice), FBP amplifies the high-frequency noise components. Students sometimes assume FBP is always superior because it is "exact" — but exactness holds only for noiseless, fully sampled data.

Correction:

In practice, use a windowed ramp filter (Shepp-Logan, cosine, Hamming) that tapers the high-frequency response. Alternatively, use iterative methods (ART, TV-CT) that can incorporate noise models.

⚠️Engineering Note

Radiation Dose and the Few-View CT Problem

Clinical CT typically uses 500-2000 projections per rotation. Reducing NvN_v lowers the X-ray dose (important for patient safety) but degrades image quality. The "few-view CT" problem — reconstructing from Nv≪200N_v \ll 200 projections — is the medical imaging analogue of the limited-aperture problem in RF imaging.

Regulatory context: Dose reduction is mandated by the ALARA (As Low As Reasonably Achievable) principle. The FDA evaluates AI-based CT reconstruction methods under the same framework as conventional algorithms, requiring demonstration of diagnostic equivalence at reduced dose.

Practical Constraints
  • •

    Projection count reduction is limited by SNR: below ∼\sim20 projections, even TV-regularized methods fail

  • •

    Learned methods must be validated against clinical ground truth, not just phantom experiments

Why This Matters: CT vs RF Forward Operator Conditioning

The conditioning gap between CT and RF imaging is the central lesson of this section. In CT, the Fourier Slice Theorem guarantees dense, uniform k-space coverage from projections over [0,π)[0, \pi) — the condition number of ACT\mathcal{A}_{\mathrm{CT}} is modest and FBP provides a stable analytical inverse.

In RF imaging, the forward operator A\mathbf{A} samples k-space on bandwidth-limited Ewald arcs with angular gaps, near-field curvature, and element-dependent weighting. The condition number can be 10-100x larger. This structural difference explains why (1) FBP suffices for CT but not for RF, (2) regularization/learning is essential for RF, and (3) architectures borrowed from medical imaging must be modified to handle the worse conditioning of A\mathbf{A}.

See full treatment in Chapter 13

Key Takeaway

CT provides an analytically invertible forward model (Radon transform

  • Fourier Slice Theorem + FBP) that serves as the gold standard for tomographic reconstruction. The RF imaging forward operator is structurally harder — incomplete k-space coverage and ill-conditioning prevent analytical inversion and demand regularized or learned methods. Every architecture we borrow from CT (FBPConvNet, TV-CT, ART) must be adapted to the RF forward model y=Ac+w\mathbf{y} = \mathbf{A}\mathbf{c} + \mathbf{w}.

The Radon Transform: Projections Build a Sinogram

A projection line rotates around a Shepp-Logan phantom, integrating the attenuation along each ray. The resulting 1D projection is stacked into a sinogram — the raw data of computed tomography. The Fourier Slice Theorem connects each projection angle to a radial line in k-space, directly paralleling how each Tx-Rx pair in RF imaging samples a point in the spatial-frequency domain.
Prerequisites & NotationMagnetic Resonance Imaging and Compressed Sensing