L3 · Chapter 6

Undersampled reconstruction: why random skipping can be repaired and regular skipping cannot

Sandbox Physics · 2026-10-29

Both ways of recording fewer lines in Chapter 5 have a price: skipping lines folds the image and keeping only the centre blurs it. This chapter describes two ways to repair the damage. The first uses several receiver coils: each is sensitive mainly to tissue near it, and the extra information is enough to unfold the image. The second skips lines at random and then fills the missing lines using the prior knowledge that the image has only a few large coefficients in some transform. The two approaches rest on different assumptions and fail in different ways.

Several coils unfold the image

The MIT handout considers two coils, one on each side of the head. Coil 1 is sensitive only to the left half, so its image shows only the left half; if only the even columns are measured, the image from coil 1 is added to a shifted copy of itself. Because coil 1 sees only the left half, the copy lands in the empty right half and the two do not overlap. Taking the left half of the folded coil 1 image and the right half of the folded coil 2 image rebuilds the full image, and since the two coils record at the same time, the scan takes half as long.

The handout then changes the sensitivities. If one coil covers the whole field of view and the other only half, subtracting one image from the other returns to the first case. The last pair, $c_5$ and $c_6$ , fails: $c_6$ gives exactly the same weight to a pair of pixels that fold onto each other, so its image holds no information that could separate them, and the same holds for $c_5$ .

In general, skipping every other line stacks each pixel $y_1$ on $y_2=y_1+\mathrm{FOV}/2$ . The folded image from coil $j$ at that point is the sum of the two, weighted by the coil sensitivity (the factor $\tfrac12$ of Chapter 5 is absorbed into the sensitivities). Two coils give two equations in two unknowns, which can be solved when the sensitivity matrix is invertible. In the $c_5$ , $c_6$ case the two columns of the matrix are equal and the determinant is zero.

\begin{pmatrix}g_1\\ g_2\end{pmatrix}=\begin{pmatrix}c_1(y_1)&c_1(y_2)\\ c_2(y_1)&c_2(y_2)\end{pmatrix}\begin{pmatrix}f(y_1)\\ f(y_2)\end{pmatrix}

One

2\times2

system per pair of folded pixels. If the matrix is invertible, the fold can be undone.

Source: MIT 6.3000 Signal Processing · M. Lustig

Worked example: unfolding one pixel pair and the noise it costs

Problem: an example made up for this course. Coil 1 has sensitivities 0.9 and 0.3 at pixels $y_1$ and $y_2$ , coil 2 has 0.2 and 0.8, and the true values are $f(y_1)=0.8$ and $f(y_2)=0.3$ . Find the folded values of both coils, solve back for the true values, and find the noise in the recovered $f(y_1)$ when each coil carries independent noise of standard deviation $\sigma$ .

The folded values are $g_1=0.9\times0.8+0.3\times0.3=0.81$ and $g_2=0.2\times0.8+0.8\times0.3=0.40$ . The determinant is $0.9\times0.8-0.3\times0.2=0.66$ . The first row of the inverse is $(0.8,\,-0.3)/0.66$ , so $f(y_1)=(0.8\times0.81-0.3\times0.40)/0.66=0.8$ , the true value; in the same way $f(y_2)=0.3$ .

Noise propagates through the rows of the inverse: the noise standard deviation of $f(y_1)$ is $\sigma\sqrt{0.8^2+0.3^2}/0.66=1.29\,\sigma$ and that of $f(y_2)$ is $1.40\,\sigma$ . The more alike the two coil sensitivities, the smaller the determinant and the larger the amplification. When the determinant is zero, as for $c_5$ and $c_6$ , the amplification is infinite, which is what the handout means by no useful information.

\sigma_{f(y_1)}=\frac{\sigma\sqrt{c_2(y_2)^2+c_1(y_2)^2}}{\left|c_1(y_1)c_2(y_2)-c_1(y_2)c_2(y_1)\right|}

From the first row of the inverse. The determinant sits in the denominator.

Source: MIT 6.3000 Signal Processing · M. Lustig

Random skipping: interference that looks like noise

With a single coil, regular skipping cannot be repaired: as Lustig and co-workers put it, the original signal and its shifted replicas are equally likely and cannot be told apart. Random skipping is a very different situation. The artifacts of the zero-filled reconstruction look like random noise, but they are not noise: undersampling leaks energy from each nonzero coefficient into others, including coefficients that were zero.

Lustig and co-workers measure this leakage with the point spread function: the values the undersampled reconstruction of a single point produces elsewhere. When $N$ of $D$ grid points are chosen at random, the off-centre values behave like zero-mean white Gaussian noise, with the standard deviation of Eq. [1], derived from energy conservation. The course benchmark draws 400 sets with $D=256$ and $N=32$ : the measured standard deviation is 0.165683, the exact finite-grid value, 0.2% above the 0.165359 of Eq. [1]. The largest sidelobe averages only 0.38, whereas regular sampling produces replicas of height exactly 1.

This is why random skipping can be repaired: strong coefficients stand above the interference and can be picked out first. Figure 2 of Lustig and co-workers shows the process on an 8-fold undersampled sparse signal: threshold to find the strongest components, compute the interference they cause and subtract it, then look for the next strongest once the interference is lower. The benchmark repeats this stagewise procedure and recovers the signal completely after 5 stages, with an error of $3\times10^{-16}$ .

\sigma_{\mathrm{SPR}}=\sqrt{\frac{p-1}{D}},\qquad p=\frac{D}{N}

Eq. [1] of Lustig and co-workers;

p

is the sampling reduction factor.

Source: M. Lustig · MIT 6.3000 Signal Processing

Sparse reconstruction: the simplest image that fits the data

Lustig and co-workers pose the reconstruction as a constrained optimization, Eq. [3]: among all images that differ from the measured data by less than $\varepsilon$ , find the one with the smallest L1 norm under the transform $\Psi$ . $\varepsilon$ is usually set a little below the expected noise level. The equivalent form in their appendix, Eq. [A1], replaces the constraint with a weight $\lambda$ chosen so that the residual is about $\varepsilon$ . They solve it with nonlinear conjugate gradients, and their examples need 80 to 200 iterations.

The k-space experiment solves the same kind of objective with three changes. $\Psi$ is a three-level Haar wavelet transform and only the detail coefficients are shrunk. The solver is an accelerated soft-threshold iteration with 50 steps. The weight is not set from the noise level; instead one eighth of the recorded lines outside the centre are held out, the fit is made without them, and the result is scored on how well it predicts them. With every line recorded the objective has a closed-form solution, and the iteration matches it to $2\times10^{-15}$ .

In the “Random lines, sparse fit fills gaps” preset (2-fold acceleration), the held-out lines pick a weight of 0.008, while the best weight, which only a simulation can know, is 0.004; their errors on the fit without the held-out lines are 17.2% and 17.0%. Rebuilt from every recorded line, the image error falls from 23.9% with zero filling to 13.3%, and the dot contrast recovers from 71% to 87.3%.

\min_m\ \|\Psi m\|_1\quad \mathrm{s.t.}\quad \|F_u m-y\|_2<\varepsilon

Eq. [3] of Lustig and co-workers.

F_u

is the undersampled Fourier operator and

y

the measured k-space data.

\min_m\ \tfrac12\|F_u m-y\|_2^2+\lambda\|\Psi m\|_1

The Lagrangian form the experiment uses. Eq. [A1] has no

\tfrac12

, which only rescales

\lambda

Source: M. Lustig · MIT 6.3000 Signal Processing

The price: low-contrast detail goes first

Lustig and co-workers write that, as the acceleration grows, the most distinct artifacts of compressed sensing are not the usual loss of resolution or an increase in aliasing interference but the loss of low-contrast features: they are weaker than the interference, and once submerged they cannot be recovered. The method suits applications with high-resolution, high-contrast features where imaging must be fast; their own examples include angiography.

They also point out that single-slice 2D imaging can only undersample the phase encodes, so the interference spreads along one direction only and, by Eq. [1], is $D^{1/4}$ times larger than for pure 2D random sampling, 16 times for a $256\times256$ image. Only modest accelerations can be expected. Their brain example uses 80 of 192 phase encodes, a 2.4-fold acceleration. The head in the experiment is also a single-slice 2D scan, accelerated 2-fold, with one coil; multi-coil unfolding appears only in the worked example above.

Source: M. Lustig · MIT 6.3000 Signal Processing

Try it in the Lab

Let the sparse fit restore the missing lines

Choose “Random lines, gaps left empty”. Nothing folds, but the image is covered in grainy artifacts. Note the 23.9% error and the 71% dot contrast.
Switch to “Random lines, sparse fit fills gaps”. Predict that most of the artifacts vanish and the contrast recovers. Check against the 13.3% and 87.3% readouts.
Open “Fit weight, centre rows and seed” and drag the sparsity weight to 0.064. Predict that a larger weight flattens the low-contrast small dots first; the readouts should show the contrast down to 64.7% and the error up to 17.6%. Then find the lowest held-out error on the weight curve in the data row.
Finally choose “Skip every other line: the head folds” and switch to the sparse wavelet fit. The fold stays, because the replicas of regular skipping are just as sparse as the original, and the prior cannot tell which one is real.

Exercises with explained answers

Two coils have sensitivities $(0.6,\,0.6)$ and $(0.5,\,0.5)$ at a pair of folded pixels. Can the pair be unfolded?
Show the answer
No. The determinant $0.6\times0.5-0.6\times0.5=0$ , and both coils measure only the sum of the two pixels. This is the $c_5$ , $c_6$ case of the MIT handout.
Sixty-four of 256 grid points are chosen at random. What is the interference standard deviation from Eq. [1], and is it larger or smaller than with 32 points?
Show the answer
$p=4$ and $\sigma=\sqrt{3/256}=0.108$ , smaller than the 0.165 for $N=32$ . More samples, less interference.
A real scan has no full scan to compare with. Why can held-out lines be used to choose the weight, and what are their limits?
Show the answer
Held-out lines are real measurements that took no part in the fit. A weight that predicts them well usually also fills the unrecorded lines well. The limits: the score is taken only at recorded positions, and leaving lines out of the fit means fitting with less data. In the experiment the two choices differ by only 0.2 percentage points, but that is one phantom and one seed.

Sources and what we recomputed

MIT 6.3000 Signal Processing, Magnetic Resonance Imaging lecture handout (6 May 2025), data and examples from E. Adalsteinsson
Read in full. Page 9: the proton gyromagnetic ratio rounded to 42.58 MHz/T. Pages 11–13: gradient fields, the phase integral and the definition of k as the time integral of the gradient. Pages 15–18: k-space data, the complex image and its magnitude. Page 19: 256 lines at 2 s each, about 8.5 minutes. Pages 20–21: keeping only even columns of k-space. Pages 24–36: two-coil unfolding, including the coil pair whose sensitivities cannot separate the folded pixels. The handout uses the sign convention of the Fourier kernel that this course follows.
M. Lustig, D. Donoho and J. M. Pauly, Sparse MRI: the application of compressed sensing for rapid MR imaging, Magnetic Resonance in Medicine 58, 1182–1195 (2007), doi:10.1002/mrm.21391
Read in full. Fig. 1 and the three requirements of compressed sensing; Fig. 2 (a 256-sample sparse signal, 8-fold undersampled); Eq. 1 (spread of the interference) and its derivation in Appendix B, which the text calls Appendix II; the remark that single-slice 2DFT interference is D^{1/4} larger; Eq. 2 (transform point spread); Eqs. 3–4 (constrained reconstruction); Appendix A, Eqs. A1–A3 (Lagrangian form and nonlinear conjugate gradients, 80 to 200 iterations); Results and Discussion (2.4-fold brain, 5-fold angiography from 22 s to 4.4 s, loss of low-contrast features). The in-vivo images are described, not reproduced.