A Fast Wavelet-Based Reconstruction Method for Magnetic Resonance Imaging

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 9, SEPTEMBER 2011 1649 A Fast Wavelet-Based Reconstruction Method for Magnetic Resonance Imaging M. Guerquin-Kern*, M. Häberlin, K. P. Pruessmann, and M. Unser Abstract—In this work, we exploit the fact that wavelets can represent magnetic resonance images well, with relatively few coefficients. We use this property to improve magnetic resonance imaging (MRI) reconstructions from undersampled data with arbitrary k-space trajectories. Reconstruction is posed as an optimization problem that could be solved with the iterative shrinkage/thresholding algorithm (ISTA) which, unfortunately, converges slowly. To make the approach more practical, we propose a variant that combines recent improvements in convex optimization and that can be tuned to a given specific k-space trajectory. We present a mathematical analysis that explains the performance of the algorithms. Using simulated and in vivo data, we show that our nonlinear method is fast, as it accelerates ISTA by almost two orders of magnitude. We also show that it remains competitive with TV regularization in terms of image quality. Index Terms—Compressed sensing, fast iterative shrinkage/ thresholding algorithm (FISTA), fast weighted iterative shrinkage/ thresholding algorithm (FWISTA), iterative shrinkage/thresh- olding algorithm (ISTA), magnetic resonance imaging (MRI), non-Cartesian, nonlinear reconstruction, sparsity, thresholded Landweber, total variation, undersampled spiral, wavelets. I. INTRODUCTION M AGNETIC RESONANCE IMAGING scanners providedata that are samples of the spatial Fourier transform (also know as k-space) of the object under investigation. The Shannon–Nyquist sampling theory in both spatial and k-space domains suggests that the sampling density should correspond to the field-of-view (FOV) and that the highest sampled fre- quency is related to the pixel width of the reconstructed images. However, constraints in the implementation of the k-space tra- jectory that controls the sampling pattern (e.g., acquisition du- ration, scheme, smoothness of gradients) may impose locally reduced sampling densities. Insufficient sampling results in re- constructed images with increased noise and artifacts, particu- larly when applying gridding methods. Manuscript received January 05, 2011; revised March 22, 2011; accepted March 26, 2011. Date of publication April 07, 2011; date of current version August 31, 2011. This work was supported by the Swiss National Competence Center in Biomedical Imaging (NCCBI) and the Center for Biomedical Imaging of the Geneva-Lausanne Universities and EPFL. Asterisk indicates corresponding author. *M. Guerquin-Kern is with the Biomedical Imaging Group, École Polytech- nique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland. M. Unser is with the Biomedical Imaging Group, École Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland. M. Häberlin and K. P. Pruessmann are with the Institute for Biomedical En- gineering, ETH Zürich, CH-8092 Zürich. Switzerland. Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TMI.2011.2140121 The common and generic approach to alleviate the recon- struction problem is to treat the task as an inverse problem [1]. In this framework, ill-posedness due to a reduced sampling density is overcome by introducing proper regularization constraints. They assume and exploit additional knowledge about the ob- ject under investigation to robustify the reconstruction. Earlier techniques used a quadratic regularization term, leading to solutions that exhibit a linear dependence upon the measurements. Unfortunately, in the case of severe undersam- pling (i.e., locally low sampling density) and depending on the strength of regularization, the reconstructed images still suffer from noise propagation, blurring, ringing, or aliasing errors. It is well known in signal processing that the blurring of edges can be reduced via the use of nonquadratic regularization. In particular, -wavelet regularization has been found to outper- form classical linear algorithms such as Wiener filtering in the deconvolution task [2]. Indicative of this trend as well is the recent advent of Com- pressed Sensing (CS) techniques in MRI [3], [4]. These let us draw two important conclusions. • The introduction of randomness in the design of trajecto- ries favors the attenuation of residual aliasing artifacts be- cause they are spread incoherently over the entire image. • Nonlinear reconstructions—more precisely, -regulariza- tion—outperform linear ones because they impose con- straints that are better matched to MRI images. Many recent works in MRI have focused on nonlinear re- construction via total variation (TV) regularization, choosing finite differences as a sparsifying transform [3], [5]–[7]. Non- quadratic wavelet regularization has also received some atten- tion [3], [8]–[11], but we are not aware of a study that compares the performance of TV against -wavelet regularization. Various algorithms have been recently proposed for solving general linear inverse problems subject to -regularization. Some of them deal with an approximate reformulation of the regularization term. This approximation facilitates re- construction sacrificing some accuracy and introducing extra degrees of freedom that make the tuning task laborious. Instead, the iterative shrinkage/thresholding algorithm [2], [12], [13] (ISTA) is an elegant and nonparametric method that is mathe- matically proven to converge. A potential difficulty that needs to be overcome is the slow convergence of the method when the forward model is poorly conditioned (e.g., low sampling density in MRI). This has prompted research in large-scale convex optimization on ways to accelerate ISTA. The efforts so far have followed two main directions. • Generic multistep methods that exploit the result of past iterations to speed up convergence, among them: two-step 0278-0062/$26.00 © 2011 IEEE 1650 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 9, SEPTEMBER 2011 iterative shrinkage/thresholding [14] (TwIST), Nesterov schemes [15]–[17], fast ISTA [18] (FISTA), and mono- tonic FISTA [19] (MFISTA). • Methods that optimize wavelet-subband-dependent param- eters with respect to the reconstruction problem: multilevel thresholded Landweber (MLTL) [20], [21] and subband adaptive ISTA (SISTA) [22]. In this work, we exploit the possibility of combining and tai- loring the two generic types of accelerating strategies to come up with a new algorithm that can speedup the convergence of the reconstruction and that can accomodate for every given k-space trajectory. Here, we first consider single-coil reconstructions that do not use sensitivity knowledge. In a second time, we con- firm the results with SENSE reconstructions [1]. We propose a practical reconstruction method that turns out to sensibly outperform linear reconstruction methods in terms of reconstruction quality, without incurring the protracted re- construction times associated with nonlinear methods. This is a crucial step in the practical development of nonlinear algo- rithms for undersampled MRI, as the problem of fixing the reg- ularization parameter is still open. We also provide a mathemat- ical analysis that justifies our algorithm and facilitates the tuning of the underlying parameters. This paper is structured as follows. In Section II, we describe the basic data-formation model for MRI and derive the dis- crete forward model. The representation of the object by wavelet bases is considered in Section II-C2; in particular, we legitimate the use of a wavelet regularization term (which promotes spar- sity) to distinguish the solution from other possible candidates. In Section III, we propose a fast algorithm for solving the non- linear reconstruction problem and present theoretical arguments to explain its superior speed of convergence. Finally, we present in Section IV an experimental protocol to validate and compare our practical method with existing ones. We focus mainly on reconstruction time and signal-to-error ratio (SER) with respect to the reference image. II. MRI AS AN INVERSE PROBLEM In this section, we present the MR acquisition model and the representation of the signal that is used to specify the reconstruc- tion problem. The main acronyms and notations are summarized in Table 1. We motivate the sparsity assumptions in the wavelet domain and the variational approach with regularization that is used to solve the inverse problem of imaging. A. Model of Data Formation 1) Physics: We consider MRI in two dimensions, in which case a 2D plane is excited. The time-varying magnetic gradient fields that are imposed define a trajectory in the (spatial) Fourier domain that is often referred to as k-space. We denote by the coordinates in that domain. The excited spins, which behave as radio-frequency emitters, have their precessing frequency and phase modified depending on their positions. The modulated part of the signal received by a homogeneous coil is given by (1) It corresponds to the Fourier transform of the spin density that we refer to as object. The measurements, concatenated in the vector , correspond to sampled values of this Fourier transform at the frequency locations along the k-space trajectory. 2) Model for the Original Data: a) Spatial discretization of the object: From here on, we consider that the Fourier domain and, in particular, the sampling points , are scaled to make the Nyquist sampling interval unity. This can be done without any loss of generality if the space domain is scaled accordingly. Therefore, we model the object as a linear combination of pixel-domain basis functions that are shifted replicates of some generating function , so that (2) with (3) In MRI, the implicit choice for is often Dirac’s delta. Different discretizations have been proposed, for example by Sutton et al. [23] with as a boxcar function or later by Delattre et al. [24] with B-splines. But it has not been worked out in detail how to get back the image for general that are noninterpolating, which is the case, for instance for B-splines of order greater than 1. The image to be reconstructed [i.e., the sampled version of the object ] is obtained by filtering the coefficients with the discrete filter (4) where denotes the Fourier transform of . b) Wavelet discretization: In the wavelet formalism, some constraints apply on . It must be a scaling function that sat- isfies the properties for a multiresolution [25]. In that case, the wavelets can be defined as linear combinations of the and the object is equivalently characterized by its coefficients in the orthonormal wavelet basis. We refer to Mallat’s book [26] for a full review on wavelets. There exists a discrete wavelet transform (DWT) that bijectively maps the coefficients to the wavelet coefficients that represent the same object in a con- tinuous wavelet basis. In the rest of the paper, we represent this DWT by the synthesis matrix . Note that the matrix multi- plications and have efficient filterbank implementations. B. Matrix Representation of the Model Since a FOV determines a finite number of coefficients , we handle them as a vector , keeping the discrete coordi- nates as implicit indexing. By simulating the imaging of the object (2), and by evaluating (1) for , we find that the noise-free measurements are given by (5) GUERQUIN-KERN et al.: A FAST WAVELET-BASED RECONSTRUCTION METHOD FOR MAGNETIC RESONANCE IMAGING 1651 where , the MRI system matrix, is decomposed as (6) There, is a space-domain vector such that . A more realistic data-formation model is (7) or (8) with and a residual vector representing the effect on measurements of noise and scanner imprecisions. The inverse problem of MRI is then to recover the coefficients (or ) from the corrupted measurements . Its degree of difficulty depends on the magnitude of the noise and the conditioning of the matrix (or ). C. Variational Formulation 1) General Framework: The solution is defined as the minimizer of a cost function that involves two terms: the data fidelity and the regularization that penalizes unde- sirable solutions. This is summarized as (9) where the regularization parameter balances the two constraints. In MRI, is usually assumed to be a realization of a white Gaussian process, which justifies the choice as a proper log-likelihood term. The ill-conditioning in- herent to undersampled trajectories imposes the use of the suit- able regularization term . Standard Tikhonov regularization corresponds to the quadratic term and leads to the closed-form linear solution that is tractable both theoretically and numerically. When the recon- struction problem is sufficiently overdetermined to make noise propagation negligible, regularization is dispensable and the least-squares solution, which corresponds to Tikhonov with , is adequate. The approach is statistically optimal if the object can also be considered as a realization of a Gaussian process. Unfortunately, this assumption is hardly justified for typical MR images. The quality of the solution obtained by this means quickly breaks down when undersampling increases. TV reconstruction is related to the sum of the Euclidean norms of the gradient of the object. In practice, it is defined as , where the operator returns pixelwise the -norm of finite differences. The use of TV regularization is particularly appropriate for piecewise-constant objects such as the Shepp–Logan (SL) phantom. 2) Sparsity-Promoting Regularization: The main idea in this work is to exploit the fact that the object can be well repre- sented by few nonzero coefficients (sparse representation) in an orthonormal basis of functions . Formally, we write that • , (sparse support) and • (small error). It is well documented that typical MRI images admit sparse representation in bases such as wavelets or block DCT [3]. The -norm is a good measure of sparsity with interesting mathematical properties (e.g., convexity). Thus, among the can- didates that are consistent with the measurements, we favor a so- lution whose wavelet coefficients have a small -norm. Specif- ically, the solution is formulated as (10) with (11) This is the general solution for wavelet-regularized inverse problems considered by [13] as well as many other authors. III. WAVELET REGULARIZATION ALGORITHMS In this section, we present reconstruction algorithms that handle constraints expressed in the wavelet domain while solving the minimization problem (10). By introducing weighted norms instead of simple Lipschitz constants, we re- visit the principle of the standard ISTA algorithm and simplify the derivation and analysis of this class of algorithms. We end up with a novel algorithm that combines different acceleration strategies and we provide a convergence analysis. Finally, we propose an adaptation of the fast algorithm to implement the random-shifting technique that is commonly used to improve results in image restoration. A. Weighted Norms Let us first define the weighted norm corresponding to a pos- itive-definite symmetric weight matrix as . The requirement that the eigenvalues of —denoted —must be positive leads to the norm property B. Principle of ISTA Revisited An important observation to understand ISTA is to see that the nonlinear shrinkage operation, sometimes called soft-thresh- olding, solves a minimization problem [27], with (12) By separability of norms, this applies component-wise to vec- tors of . This means that the -regularized denoising problem (i.e., when in (11) is the identity matrix) is precisely solved by a shrinkage operation. The iterative shrinkage/thresholding algorithm (ISTA) [2], [13], also known as thresholded Landweber (TL), generates a sequence of estimates that converges to the minimizer of (11) when it is unique. The idea is to define at each step a new functional whose minimizer will be the next estimate (13) 1652 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 30, NO. 9, SEPTEMBER 2011 Two constraints must be considered for the definition of . 1) It is sufficient for the convergence of the algorithm that is an upper bound of and matches it at ; this guarantees that the sequence is monotonically decreasing. 2) The inner minimization (13) should be performed by a simple shrinkage operation to ensure the rapidity and accuracy of the algorithm. In accordance with Constraint (1), can take the generic quadratically augmented form (14) with the constraint that is positive definite, where the weighting matrix plays the role of a tuning parameter. Then, ISTA corresponds to the trivial choice , with the value of chosen to be greater or equal to the Lipschitz constant of the gradient of , so that . Let us define , , and (15) Then, using standard linear algebra, we can write (16) (17) This shows that Constraint (2) is automatically satisfied. Note that both the intermediate variable in (15) and the threshold values will vary depending on . Algorithm 1: ISTA Repeat ; Beck and Teboulle [18, Thm. 3.1] showed that this algorithm decreases the cost function in direct proportion to the number of iterations , so that . Here, we present a slightly extended version of their original result which is valid for any reference point during iterations. Proposition 1: Let be the sequence generated by Al- gorithm 1 with . Then, for any , (18) Proof: [18, Thm. 3.1] gives the result for . Con- sider a sequence such that . As the iteration does not depend on , we get . The result fol- lows immediately. Selecting as small as possible will clearly favor the speed of convergence. It also raises the importance of a “warm” starting point. Among the variants of ISTA, FISTA, proposed by Beck and Teboulle [18], ensures state-of-the-art convergence properties while preserving a comparable computational cost. Thanks to a controlled over-relaxation at each step, FISTA quadratically decreases the cost function, with . C. Subband Adaptive ISTA (SISTA) SISTA is an extension of the ISTA that was introduced by Bayram and Selesnick [22]. Here, we propose an interpretation of SISTA as a particular case of (14) with a weighting matrix that replaces advantageously the step size in (15). The idea is to use a diagonal weighting matrix —if it is not diagonal, Constraint (2) would not be fulfilled—with coef- ficients that are constant within a wavelet subband. 1) SISTA: In the same fashion as for ISTA, (15) and (17) can be adapted to subband-dependent steps and thresholds.Accord- ingly, SISTA is described in Algorithm 2. Algorithm 2: SISTA Repeat ; 2) Convergence Analysis: By considering the weighted scalar product instead of , we can adapt the convergence proof of ISTA by Beck and Teboulle (see Propo- sition 1). This result is new, to the best of our knowledge. Proposition 2: Let be the sequence generated by Al- gorithm 2 with . Then, for any , (19) Proof: We rewrite the cost function (11) with the change of variable . We then apply ISTA to solve that problem with , , and (Note that is positive-definite if is positive-definite). The iteration can be rewritten, in terms of the original variable, as . The latter is an iteration of SISTA (see Algorithm 2). According to Proposition 1, we have , which translates directly into the proposed result. Therefore, by comparing Propositions 1 and 2, an improved convergence is expected. The main point is that, for a “warm” starting point or after few iterations , the weighted norm in (19) can yield significantly smaller values than the one weighted by in (18). 3) Selection of Weights: Bayram and Selesnick [22] provide a method to select the values of for SISTA. To present this result, let us introduce some notations. We denote by an index that scans all the wavelet