Chapter 8: Iterative Reconstruction

FBP inverted the projection analytically in a single pass. There is another way to look at reconstruction: as solving a system of linear equations. Stack the $N^2$ pixel values into a vector $x$ and the $M$ measurements (projections × detectors) into a vector $p$; forward projection then becomes a linear map by a huge sparse matrix $A$, called the system matrix, whose element $a_{ij}$ is the weight with which ray $i$ crosses pixel $j$, and reconstruction reduces to solving $Ax = p$ for $x$. With few projections the system is underdetermined (more unknowns than independent equations), and with measurement noise no exactly consistent solution exists at all. This is where iterative solvers come in, starting from an initial image and gradually reducing the mismatch with the measurements. That approach is iterative reconstruction.

$$A\mathbf{x} = \mathbf{p}$$

The oldest iterative method is ART (Algebraic Reconstruction Technique), known in mathematics as the Kaczmarz method. It takes the equations one at a time (one per ray) and projects the current image $x$ onto the hyperplane $a_i \cdot x = p_i$ defined by that equation, cycling through all rays. The numerator $p_i - a_i \cdot x$ is the residual, the mismatch between the measured value of that ray and the forward projection of the current image, and it is smeared back along the ray's path to cancel the discrepancy. The relaxation factor $\lambda$ sets how far each update moves toward the hyperplane; smaller values average out noise and stabilize the behavior. For a consistent noise-free system, this sequence of projections provably converges to a solution.

$$\mathbf{x} \;\leftarrow\; \mathbf{x} + \lambda\,\frac{p_i - \mathbf{a}_i\!\cdot\!\mathbf{x}}{\lVert \mathbf{a}_i \rVert^2}\,\mathbf{a}_i$$

Because ART rewrites the image ray by ray, it starts fast but is easily jerked around by a single noisy ray. SIRT (Simultaneous Iterative Reconstruction Technique) instead computes the residuals of all rays at once and applies a single combined correction, backprojected with $A^{\mathsf T}$ and normalized by the row sums $R$ (per ray) and column sums $C$ (the per-pixel sensitivity). Each iteration advances gently, but averaging over many rays makes it robust to noise. SART sits in between: it updates with the rays of one projection angle (view) at a time using the same formula as SIRT, balancing speed and stability.

$$\mathbf{x} \;\leftarrow\; \mathbf{x} + \lambda\, C^{-1} A^{\mathsf T} R^{-1}\,(\mathbf{p} - A\mathbf{x}), \qquad R = \operatorname{diag}\Bigl(\sum_j a_{ij}\Bigr),\quad C = \operatorname{diag}\Bigl(\sum_i a_{ij}\Bigr)$$

The methods so far treat reconstruction geometrically, as equation solving. A statistical method instead starts from the measurement distribution. In emission tomography, measured photon counts are modeled as Poisson variables. MLEM is a multiplicative maximum-likelihood update that backprojects the ratio $p/(Ax)$, normalizes by sensitivity, and preserves non-negativity. OSEM divides the projections into subsets to accelerate the updates and is standard in PET and SPECT.

$$x_j \;\leftarrow\; \frac{x_j}{\sum_i a_{ij}} \sum_i a_{ij}\, \frac{p_i}{(A\mathbf{x})_i}$$

So when should each be used? With plenty of projections and little noise, single-pass FBP is fast and excellent, and there is little reason to iterate. Under harsh conditions, however, with few projections, low dose (heavy noise), or incomplete data, the $|\omega|$ filter of FBP amplifies undersampling and noise directly, and the image quality collapses. Iterative methods re-enforce consistency with the measurements step by step and remain stable with far fewer artifacts under the same conditions. The computational cost is orders of magnitude higher, but the dose reductions of modern clinical CT rest on the practical deployment of iterative reconstruction and its successors.

Simulation: driving the iterations

Pick a method, then advance the iteration one step at a time with "Step" or automatically with "Run", and watch the reconstruction rise out of the zero image. ART shows an image within a few iterations, SIRT converges slowly and smoothly, and OSEM converges roughly a subset-count times faster than MLEM. Changing the method, $\lambda$, the number of subsets, the number of projections, or the noise restarts the iteration from scratch (mandatory when switching methods: MLEM/OSEM require a non-negative image). The display window is fixed to the true phantom, and RMSE is evaluated inside the central 90% circle of the FOV.

Simulation: watching convergence

The left plot traces the relative residual $\|Ax - p\|/\|p\|$ of the method running in the simulator above, iteration by iteration (log scale). The residual measures how much the forward projection of the current image still disagrees with the measurements. The comparison mode on the right runs all five methods for 10 iterations each under the identical conditions (projections, noise, $\lambda$, subsets) and overlays their residual curves. Compare the fast early convergence of ART, the gentleness of SIRT, and the subset acceleration of OSEM. With noise present, the residual stops decreasing at the noise floor; pushing it lower would mean overfitting the noise.

Simulation: when to prefer iterative over FBP

Press one of the harsh-condition presets to reconstruct the same sinogram with FBP (Ram-Lak) and with 10 iterations of the method selected above, side by side. With sparse projections (24) FBP shows strong radial streaks while the iterative result stays much cleaner. With heavy noise ($I_0 = 3\times10^3$) Ram-Lak FBP amplifies the noise while SIRT or OSEM deliver smoother images. Confirm that the RMSE also favors the iterative method (changing the method or $\lambda$ clears the result).

Choosing a method in practice

Reconstruction can be posed as solving the linear system $Ax = p$. ART and SIRT/SART reduce residuals through additive updates, while the simplified MLEM/OSEM example demonstrates multiplicative updates and subsets. Iterative methods can remain stable under sparse-view and low-dose conditions where FBP struggles, at a much higher computational cost. Image quality depends on the statistical model, regularization, iteration count, relaxation factor, and subset count.

References

• Andersen AH, Kak AC. Simultaneous algebraic reconstruction technique (SART). Ultrasonic Imaging 6, 81–94 (1984).
• Shepp LA, Vardi Y. Maximum likelihood reconstruction for emission tomography. IEEE Transactions on Medical Imaging 1, 113–122 (1982).
• Hudson HM, Larkin RS. Accelerated image reconstruction using ordered subsets of projection data. IEEE Transactions on Medical Imaging 13, 601–609 (1994).
• Thibault JB et al. A three-dimensional statistical approach to improved image quality for multislice helical CT. Medical Physics 34, 4526–4544 (2007).