Low-rank reconstruction

The idea

Some acquisitions repeat the same volume many times, with the image changing in a structured way between repetitions: an MR elastography wave moving through the tissue, for example. Reconstructing each repetition on its own throws away what the repetitions share. If the k-space trajectory also changes between repetitions, as in OSCILLATE acquisitions where the spiral is rotated, each repetition samples k-space differently and the repetitions can fill in for each other.

oscillateRecon uses a partial separability model. The image series, arranged as a matrix with one column per repetition, is assumed to have low rank:

\[X = U\,V_L^H\]

V_L holds L temporal basis vectors (NRep × L), fixed in advance. U holds L spatial images (voxels × L), and these are what the solver estimates. With L much smaller than the number of repetitions, there are far fewer unknowns than in separate reconstructions, and the low-rank constraint removes noise that doesn't follow the temporal pattern.

How FORGE solves it

The forward operator (LRobj) forms the full series U V_L^H, then applies each repetition's own pcSenseTimeSeg operator to its image. Each repetition can therefore have its own trajectory, shot phase maps and time segmentation. The adjoint runs each repetition's adjoint and projects the result back onto the basis.

The quadratic roughness penalty (R_lowRank) is applied to each repetition's image in the full series, not to the spatial basis images, so -B means the same as in the single-repetition programs.

The solver is the same conjugate gradient as elsewhere. Because LRobj works on Armadillo vectors only, the solver's own vector arithmetic runs on the CPU, while each repetition's operator still uses the GPU in Metal and CUDA builds.

Where the basis comes from

FORGE reads the basis from the v array in the input file and uses its first --Rank columns. It is usually computed during data preparation from a fully sampled low-resolution version of the series: the right singular vectors of its Casorati matrix (voxels × repetitions). The singular values show how many components carry signal above the noise, which is a starting point for the rank.

Trade-offs

• A rank that is too low forces real temporal changes into the wrong shape, and biases the time course.
• A rank that is too high lets noise back in and gives the solver more unknowns.
• Image-quality scores that measure noise between repetitions can look unrealistically good after a low-rank reconstruction, because the model itself removes that noise. Compare against a per-repetition reconstruction instead.

Further reading

• Z.-P. Liang, "Spatiotemporal imaging with partially separable functions," Proceedings of the IEEE International Symposium on Biomedical Imaging, 2007.