324	P.J. Bones and B. Wu

14.2 Sparsity in MRI Images

In this section, we explore the properties of MR images, which make them amenable to compressed sensing, show examples of how some common transforms can be used to exploit the sparsity, and introduce a novel nonlinear approach to promoting sparsity.

14.2.1 Characteristics of MR Images (Prior Knowledge)

As mentioned in the introduction, MRI measurements are made in the k-space domain. In some cases, the measurements may be at positions constrained by a regular Cartesian grid. Since an inverse Fourier transform is required to generate an image from the sampled k-space data, the regular Cartesian sample positions allow the straightforward use of the efficient FFT algorithm. However to achieve faster scanning or some signal processing advantages non-Cartesian sampling is frequently employed. Radial and spiral sampling [2] are quite common, for example. Some sampling strategies involve a higher density of samples near the center of k-space (i.e., concentrated in the area of lower spatial frequencies). In any case, there is a direct relationship between the extent of the k-space domain within which measurements are distributed and the resolution of the image obtained. Likewise, there is a direct relationship between the field-of-view (FOV) in the spatial domain and the effective spacing of samples in k-space [1].

The main source of noise in MR imaging is due to thermal fluctuations of electrolytes in the region being imaged which are detected by the receiver coil or coils. Electronic noise is inevitably present as well but may usually be of lesser order. Generally, the SNR increases as the square root of the acquisition time and linearly with the voxel size. Thus any moves to increase imaging speed and/or imaging resolution inevitably lead to a loss of SNR. Importantly, the noise statistics of each k-space sample is essentially equal [1]. Since the amplitude of samples near the origin of k-space is much greater than near the periphery, the SNR of these center samples is much better. This consideration often influences the design of a sampling scheme.

In many significant imaging situations, the object is known not to extend throughout the FOV. For example, in making a 2D axial plane image of the brain, the FOV is usually chosen to be a square or rectangular region that entirely contains the outline of the head. There is therefore part of the FOV which lies outside the head and which is therefore known not to contribute a significant signal to the measurements made. This support constraint is explicit. A support constraint may also be implicit: for example if it is known that a transformed version of the image is known to be nonzero only within an unspecified region which spans some known proportion of the transformed space.

14 Sparse Sampling in MRI

325

Fig. 14.2 Compressibility of an MR brain slice. The wavelet and DCT approximations of the original MR brain slice shown in (a), using only the 10% highest amplitude and the 5% highest amplitude coefficients of the transforms, are shown in (b) to (e). There appears to be little loss of information under such high levels of image compression

Constraints may also be temporal. While the relatively slow acquisition of MR data restricts its use in dynamic imaging tasks, the facts that measurements are made at precise timing instants and that objects under observation move relatively smoothly allows temporal constraints to be formulated and exploited.

Biological tissues comprise a large number of types. While at a microscopic level these tissues are generally inhomogeneous, at the resolution observed by MR techniques, each tissue type is relatively homogeneous. It is the difference in MR signal between tissue types which allows such useful anatomical information to be gleaned. The image as a whole therefore exhibits an approximately piecewise homogeneity, which can be exploited.

The forms of prior knowledge about the MR images discussed above can all be seen as evidence for expecting the images to be sparse.

14.2.2 Choice of Transform

The term implicit support was introduced in Sect. 14.2.1. This represents the property that under some transformation an image can be shown to be nonzero only within some unspecified part of the domain. The success of lossy image compression schemes based on the DCT and wavelet transforms indicate that these are useful sparsifying transforms. In Fig. 14.2, we illustrate the degree to which a typical MR image can be compressed by the two transforms. The image in Fig. 14.2a is formed for an axial slice of the brain with the full resolution afforded by the MRI sequence employed (256 ×256). The other parts of the figure show reconstructions with only a fraction of the transform coefficients retained. It is clear that under either of the two transforms a substantial reduction of data volume is possible before

326	P.J. Bones and B. Wu

serious degradation of the image occurs. Note that the compression here is “lossy,” in that there is always some degradation generated by setting small coefficients to zero, but not so much degradation that the image usefulness is seriously impaired.

The discrete cosine transform (DCT) was the transform of choice for many years for image compression [7]. It was the basis of the original JPEG standard. The properties of the DCT which make it a useful choice for image compression include:

1.It lends itself to representing small rectangular blocks of pixels

2.It can be shown to have a faster fall off in coefficient amplitude as frequency increases in comparison with the DFT

3.It is relatively efficient to compute via the FFT

The DWT has taken over from the DCT in certain regards [8]. The properties of the DWT which make it a useful choice for image compression include:

1.It naturally distributes information over a number of scales and localizes that information in space

2.It offers a wide range of basis function (wavelet families) from which to choose

3.It is inherently efficient to compute

The decision between the transforms is unlikely to be critical. The nature of the DWT, however, does render it better at representing boundaries in the image between two tissue types where the image function exhibits a step change. With the DCT, such a boundary necessarily injects some energy in high frequency components and adversely affects the sparsity in the transformed representation. The DWT with an appropriate choice of wavelet may perform better in this regard.

14.2.3 Use of Data Ordering

A quite distinctly different approach for increasing sparsity has recently been proposed. In 2008, Adluru and DiBella [9] and Wu et al. [10] independently proposed performing a sorting operation on the signal or image as part of the reconstruction process. The principle is presented in Fig. 14.3 for a 2D axial brain image. In Fig. 14.3a, the situation is shown whereby the image is transformed by the 2D DCT and then a compression occurs by setting all coefficients less than a given threshold to zero. The resulting reconstruction is similar to the original, but noticeably smoother due to the loss of some small amplitude high frequency components. In Fig. 14.3b, the image pixels are sorted from largest amplitude in the lower right to highest amplitude in the upper left to make the resulting function monotonic and the mapping required to do this is retained (denoted “R”). The same transformation and recovery operation after thresholding as in (a) is performed and a re-sorting (denoted “R−1”) is performed. Because the compression retains much of the shape of the image after sorting, the result has much higher fidelity than in

14 Sparse Sampling in MRI

327

Fig. 14.3 Illustration of how a data ordering can achieve a higher sparsity for a 2D image. In (a), the signal is compressed by retaining only those DCT coefficients with amplitudes higher than a threshold. In (b), the image pixels are sorted to generate a monotonic function and then the same recovery operation is performed before a final resorting. Because the sorted data in (b) is more sparse, the recovery is of higher quality

Fig. 14.3a. We argue that many fewer coefficients need to be retained in the DCT of the sorted image than in the original, hence the more successful reconstruction.

Clearly, the process depicted in Fig. 14.3 requires knowledge of the original signal to derive R. The practical utility of what has been demonstrated is likely therefore to be questioned. However, we show in Sect. 14.4.2 that several methods to derive an approximate R are possible and that they lead to useful and practical algorithms for MR image recovery.

The advantage claimed for data ordering depends on the sorted data function being more sparse than the original. It is difficult if not impossible to prove this,

328	P.J. Bones and B. Wu

Fig. 14.4 Examples which show the higher sparsity that data ordering can achieve for 2D images: (a) a set of 5 original images; (b) the number of DCT coefficients needed to achieve a reconstruction of each image with a relative mean square error (RMSE) ≤ 1%; (c) the image after a data ordering operation; and (d) the number of DCT coefficients needed to achieve a reconstruction of the ordered image with RMSE≤ 1%

but experiments indicate that in virtually all cases it is true. In Fig. 14.4, we show a set of 2D example images; three are typical brain axial slices (128 × 128 pixels) while the others are popular test images (200 ×200 pixels). A data sorting operation was performed on each image such that the highest intensity pixel was positioned in the top left corner and the lowest intensity pixel in the bottom right corner. A simple algorithm arranged the others in a type of raster which achieved the “wedge” images shown in column (c). Clearly, other algorithms for arranging the sorted pixels are possible. To the right of each image in columns (b) and (d) is the number of DCT