Fast compressed-sensing reconstruction for magnetic-resonance imaging
To digitally process continuous physical phenomena, signals must be described using a finite set of numbers. This sampling is traditionally done by taking equally spaced signal readings, where the required number of measurements is directly related to the highest frequency present.1,2 However, for some applications this would require an impractically large number of measurements. For example, in magnetic-resonance imaging (MRI)3 the goal is to acquire detailed images of a patient's anatomy. To reduce scanning time and the patient's exposure to electromagnetic radiation it is, however, desirable to take as few measurements as possible.
Compressed sensing opens up a new paradigm that allows signal acquisition using far fewer samples than predicted on the basis of traditional theory whenever the signal is sparse, i.e., whenever most of the contributing elements are approximately zero.4,5 The technique is applicable to a wide range of signal-acquisition problems. For example, images are often sparsely represented in the wavelet domain, and compressed-sensing approaches can be applied to a wide range of imaging modalities such as hyperspectral imaging, projection tomography, and synthetic-aperture radar imaging.
In Figure 1, a patient's mid-sagittal view (top left) is sparsely represented in the wavelet domain (top right). An MRI scanner samples in the 2D Fourier domain (bottom left), but, contrary to traditional sampling, in compressed sensing only a subset of the frequency information is measured (bottom right). To approximate the missing data, we exploit the sparsity in the wavelet domain. Importantly, if the sampling procedure satisfies certain conditions and if the wavelet representation is sparse enough, any other wavelet representation leading to the same observations will have far more nonzero elements. In order to estimate the wavelet coefficients, we need to search for the sparsest representation that is consistent with the measurements.
More generally, let the vector x represent the data in a domain in which it is sparse, with the observations collected in a vector y. If a linear transformation exists from the sparse domain to the observations, we can model the measurement process as
where e accounts for observational or quantization noise. In the MRI example, the matrix Φ represents the sequential inverse wavelet and Fourier transforms, followed by subsampling. Since in compressed sensing only a few measurements are taken, y contains far fewer elements than x. Therefore, for each y there are infinitely many matching signals x. From among all possible solutions, we search for x with the fewest nonzero elements such that e is small. Since this is computationally expensive, much effort has gone into the development of suboptimal but computationally efficient algorithms.7–12
We have developed one such iterative approach, ‘gradient pursuit’ (see Figure 2).7,9 This algorithm is characterized by very low computational complexity and performs well for many applications. The main idea is to search for a few nonzero wavelet coefficients and estimate their coefficient values on the basis of a least-squares procedure. Their contribution to the observation is then subtracted from the actual measurements, and the algorithm continues with the modified reference data set until a sparse representation is found that accounts for most of the observation. While it is not guaranteed that our approach will generally find the sparsest matching vector x, we can show that—under certain conditions—discovery of good estimates is certain.8
The advantages of the compressed-sensing approach can be demonstrated using an MRI example. Figure 3 (top left) shows one frame taken from an eight-frame sequence of a dynamic-MRI mouse scan. Each original data frame was acquired by traditional sampling synchronized to the mouse's heartbeat, taking a small number of measurements each time. We then artificially downsampled the data by a factor of five and used the gradient-pursuit algorithm for reconstruction (middle panels). For comparison, standard linear-reconstruction algorithms show clear aliasing artifacts and blurring (bottom panels).
The gradient-pursuit algorithm was developed to efficiently calculate reconstructions on the basis of sparse data representations, in particular for large data sets—as found in many imaging problems. While sparsity is a powerful constraint, Figure 1 shows that the wavelet representation exhibits additional structure. For example, large wavelet coefficients on one scale typically correspond to large coefficients on the next scale too. We are currently developing a more general class of signal models and new algorithms that can optimally exploit this tree structure. This may enable us to further reduce the number of measurements and improve reconstruction accuracy.
This research was supported by grants D000246/2 and EP/F039697/1 from the UK's Engineering and Physical Sciences Research Council. Mike Davies acknowledges support for his position from the Scottish Funding Council and their support of the Joint Research Institute with Heriot-Watt University as a component part of the Edinburgh Research Partnership. We thank Ian Marshall and Yuehui Tao for providing access to the MRI data.
