The solver and regularization¶
What is being solved¶
Every iterative reconstruction in FORGE solves a penalized weighted least-squares problem:
The first term asks the image to explain the data through the encoding operator A (see Operators). W weights the data; the command-line programs use uniform weights. The second term, R, penalizes images that are implausible, here ones that are rough.
Conjugate gradient¶
solve_pwls_pcg() minimizes this with nonlinear conjugate gradient, using the Polak-Ribière-Polyak update for the search direction and a step size from a quadratic surrogate of the penalty. Each iteration costs one forward and one adjoint application of A, which dominate the run time.
Conjugate gradient builds up the image from the most strongly encoded components first. Early iterations recover the coarse, low-frequency structure, and later ones add fine detail and, eventually, noise. Stopping after a fixed number of iterations (-n) is therefore itself a form of regularization: too few iterations give a blurred image, too many give a noisy one. Your first reconstruction shows the effect.
The solver stops early if the normalized gradient falls below 10−10, and between iterations if the program receives Ctrl-C. It starts from a zero image.
When A and the penalty support forgeCol and the build has a GPU backend, the solver keeps its vectors on the GPU for the whole solve and runs its vector arithmetic there. Otherwise it uses Armadillo on the CPU.
Regularization¶
Penalties derive from Robject, which applies a finite-difference operator C between neighbouring voxels and a potential function ψ to each difference:
| Class | ψ(d) | Effect |
|---|---|---|
QuadPenalty |
½ d2 | Smooths everywhere, edges included. Keeps the problem quadratic, so the solver's line search is exact in one step. |
TVPenalty |
δ2(√(1+(d/δ)2) − 1) | Smoothed total variation: quadratic for differences smaller than δ, close to linear above it, so it smooths noise but preserves edges better |
The command-line programs use QuadPenalty. -B sets β and -D chooses whether differences are taken in 2 (x, y) or 3 (x, y, z) dimensions. Use -D 2 when slices are far apart compared with the in-plane voxel size, or reconstructed independently.
To write a new penalty, subclass Robject and override wpot, dpot and pot, which return ψ′(d)/d, ψ′(d) and ψ(d).
Choosing β¶
β balances data against smoothness, and its useful range depends on the scale of the data. If the data values are large, a β that is meaningful for one scanner or preparation can be negligible for another. Look at the penalty and error_norm values the solver reports each iteration (shown by forgeview, and in the JSON lines): if the penalty is many orders of magnitude below the error norm, β has no effect. A sweep over a few orders of magnitude on one representative dataset is usually the quickest way to settle it.
Further reading¶
- J. A. Fessler, "Penalized weighted least-squares image reconstruction for positron emission tomography," IEEE Transactions on Medical Imaging, 13(2), 1994. doi:10.1109/42.363108