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.
Figure 1. Compressed sensing applied to medical-resonance imaging (MRI). The mid-sagittal view (top left) is sparsely represented in the wavelet domain (top right). The MRI scanner takes samples in the 2D Fourier domain (bottom left). Compressed sensing only measures a subset of this data, such as a small number of radial slices (bottom right). Reconstruction algorithms search for the sparsest representation in the wavelet domain consistent with the observations. Inverse wavelet transformation enables image reconstruction.
Overview of one gradient-pursuit iteration. First, new nonzero coefficients are selected by looking through each element in turn and checking by how much a change in each coefficient reduces the discrepancy between the actual (y
) and estimated (Φx
)observation. We optimize the values of the nonzero elements by taking into account the direction of previous updates, similarly to the conjugate-gradient algorithm.6
The third step uses our knowledge of the observation system to calculate what it would look like with the newly estimated coefficients. x: Vector to represent the data in a domain in which it is sparse.
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
Figure 3. A frame from an MRI movie of a beating mouse heart, synchronized to the heart beat (top left), with magnified detail to the right. (Middle) Reconstruction using our approach based on only one-fifth of the total number of measurements. (Bottom) Reconstruction using a standard linear approach, clearly showing aliasing artifacts and additional blurring in the magnified region.
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.
Thomas Blumensath, Mike Davies
Institute for Digital Communications
The University of Edinburgh