14 Sparse Sampling in MRI

331

14.3.2 Signal Recovery

Assuming that we have chosen a sampling strategy for the k-space data and a transform Φ under which the true image is expected to be sparse, we seek a solution x as close as possible to x, which is constrained in two ways:

1.The solution is consistent with the data dM

2.The solution is sparse under the transformation Φ

Condition 1 can be achieved in principle at least by minimizing the power of the error between the measurements and the values at those measurement points, which are predicted by the imaging model for the current image estimate, that is by minimizing the squared norm ||dM − WM x ||2. Such squared norm minimizations have formed the backbone of image recovery for many years [14].

Condition 2 above implies a minimization of the quantity ||Φ x ||0, that is the number of nonzero elements in Φ x . However, this minimization is computationally intractable [4, 15]. It turns out that a minimization of the quantity ||Φ x ||1, that is the first norm of the transformed image estimate, can achieve Condition 2 remarkably well [4,16]. The l1 norm applied here has the effect of pushing negligible coefficients toward zero while retaining larger components accurately. This is in contrast with a squared norm which tends to penalize large coefficients.

As explained in the previous section, the random sampling patterns which offer advantages in CS do generate noise-like artifacts. Therefore in our experience, it is also useful to apply a further constraint:

3. The solution is piecewise smooth

Minimizing the total variation (TV), that is the sum of the magnitudes of differences between each pixel and its immediate neighbors, has been shown to be effective at meeting Condition 3. We denote the total variation for image vector y, TV(y).

The minimization problem can now be posed: Find an estimate for the required image x by minimizing

||dM − WM x ||2 + λ ||Φ x ||1 + β TV(x )

where λ and β are positive constants used to control the relative importance of the constraints being placed on the solution. A method such as conjugate gradient minimization is suitable to solve the problem. We have employed the SparseMRI package provided by [13] as part of the very comprehensive resource on compressed sensing provided by Rice University (see http://dsp.rice.edu).

14.4 Progress in Sparse Sampling for MRI

In this section, we briefly review the progress made to date in applying the principles of sparse sampling to MRI. We first review the important developments that have appeared in the literature. We believe that the biggest single contribution came

332	P.J. Bones and B. Wu

from Lustig, Donoho, and Pauly [13]. This group has continued to make valuable contributions. Our own contributions, in the form of two new algorithms for applying sparse sampling in MRI, are then presented.

14.4.1 Review of Results from the Literature

Prior to the introduction of compressed sensing, exploiting signal sparseness by utilizing the l1 norm constraint started in mid-1980s when Santosa, Symes, and Raggio [17] utilized an l1 norm to recover a sparse series of spikes from the aliased representation that resulted from sub-Nyquist sampling in the Fourier domain. A similar experiment was implemented by Donoho [18] using soft thresholding, where the individuals in a sparse series of spikes were recovered sequentially in the order of the descending magnitude: the strongest component was first recovered and its aliasing effects were then removed to reduce the overall aliasing artifacts to allow the next strongest component to be recovered, and so on. These simple numerical experiments in fact have the same nature as the application of the modern compressed sensing technique in contrast-enhanced magnetic resonance angiography (CE-MRA). In CE-MRA, the contrast-enhanced regions to be recovered can be usefully approximated as isolated regions residing within a 2D plane, and hence the use of simple l1 norm suffices in recovering the contrast-enhanced angiogram.

Another application of the l1 norm before compressed sensing is in the use of TV filter [19], which imposes a l1 norm in gradient magnitude images (GMI), or the gradient of the image. As discussed previously, l1 norm promotes the strong components while penalizing weak components. In the operations on the GMI, the TV operator suppresses small gradient coefficients, whereas it preserves large gradient coefficients. The former are considered as noise to be removed, whereas the latter are considered to be part of the image features (edges) that need to be retained; hence, TV can serve as an edge-preserving denoising tool. TV itself can be employed as a powerful constraint for recovering images from undersampled data sets. In [20], TV is employed to recover MR images from undersampled radial trajectory measurements; Sidky and Pan [21] investigated the use of TV in recovering computed tomography images from limited number of projection angles.

The formal introduction of compressed sensing into MRI methods was made by Lustig, Donoho, and Pauly in 2007 [13]. Their key contribution is the explicit use of a different transform domain for appropriate application of the l1 norm. Both the sparse set of spikes and the TV filter mentioned previously are special instances of the general transform-based compressed sensing setup. The authors identified the use of DWT and DCT as suitable transform bases for application in MR images, as evidenced by their sparse representation under DWT and DCT. A reconstruction framework was given, which converts the CS formulation into a convex optimization problem and hence allows for computational efficiency. The authors also spelt out that a key requirement in data measurement for successful compressed sensing recovery is to achieve incoherent aliasing. In MRI, such a