MODEL-BASED ITERATIVE RECONSTRUCTION FOR ... - IEEE Xplore

MODEL-BASED ITERATIVE RECONSTRUCTION FOR MAGNETIC RESONANCE FINGERPRINTING Bo Zhao Department of Electrical and Computer Engineering, University of Illinois at Urbana-Champaign

Index Terms— Model-based image formation, parameter estimation, low-rank modeling, Bloch equation, sparse sampling, magnetic resonance imaging.

tion exactly produces the conventional MRF reconstruction, while additional iterations push the solution toward a statistical optimal solution. In this work, we present a new model-based MRF reconstruction method within the statistical estimation framework, incorporating a low-rank image model with the Bloch-equation based MR physical model. It enables direct parameter estimation from highly undersampled, noisy k-space data. To solve the resulting optimization problem, we describe an efficient algorithm based on variable splitting (VS) [4–6], the alternating direction method of multipliers (ADMM) [4, 5, 7], and the variable projection method (VARPRO) [4, 8, 9]. The performance of the proposed method has been evaluated using numerical experiments. Both qualitative and quantitative improvements of the proposed method over the conventional MRF reconstruction and recently proposed ML reconstruction will be demonstrated. The rest of the paper is organized as follow. In Section 2, we describe the proposed formulation and solution algorithm in detail. Section 3 shows representative results obtained with the proposed method, followed by the discussion and conclusion in Section 4.

1. INTRODUCTION

2. PROPOSED METHOD

ABSTRACT Magnetic resonance fingerprinting (MRF) is an emerging quantitative magnetic resonance (MR) imaging technique that simultaneously acquires multiple tissue parameters (e.g., spin density, T1 , and T2 ) in an efficient imaging experiment. A statistical estimation framework has recently been proposed for MRF reconstruction. Here we present a new model-based reconstruction method within this framework to enable improved parameter estimation from highly undersampled, noisy k-space data. It features a novel mathematical formulation that integrates a low-rank image model with the Bloch equation based MR physical model. The proposed formulation results in a nonconvex optimization problem, for which we develop an efficient iterative algorithm based on variable splitting, the alternating direction method of multipliers, and the variable projection method. Representative results from numerical experiments are shown to illustrate the performance of the proposed method.

Magnetic resonance fingerprinting (MRF) [1] is a recent breakthrough in quantitative magnetic resonance (MR) imaging, which utilizes a specialized fast pulse sequence to simultaneously obtain multiple tissue parameters, e.g., spin density, spin-lattice relaxation time T1 , and spin-spin relaxation time T2 . The conventional MRF reconstruction [1] employs a pattern-matching based non-iterative procedure: for each voxel, it selects a signal evolution (from a pre-specified dictionary) that yields the maximum correlation with the one obtained from the gridding reconstruction, and then assigns the underlying parameters generating the selected signal evolution as its reconstruction. Although this procedure is computationally efficient, it is suboptimal from a statistical estimation standpoint. To achieve accurate parameter estimation, a relatively large number of time points is often needed, leading to prolonged acquisition for applications requiring high-resolution and/or broad volumetric coverage. A number of advanced MRF reconstruction methods have recently been proposed. For example, Ref. [2] utilized sparse representations to reduce acquisition time and improve SNR, whereas Ref. [3] formulated the reconstruction problem as a low-dimensional manifold recovery problem, which was then solved by a greedypursuit algorithm. In our early work in [4], a statistical estimation framework was introduced to MRF reconstruction, within which a maximum likelihood (ML) reconstruction formulism was proposed. In particular, it has been shown in [4] that with the gridding reconstruction as an initialization, the first iteration of the ML reconstruc-

978-1-4799-8339-1/15/$31.00 ©2015 IEEE

3392

2.1. Formulation 2.1.1. Data Model In MRF, the contrast-weighted images {Im (x)}M m=1 can be parameterized as [3, 4]: Im (x) = ρ (x) ϕm (T1 (x) , T2 (x) , f0 (x)) ,

(1)

for m = 1, · · · , M , where ρ (x) denotes the spin density map, and ϕm denotes the contrast-weighting function that depends on the tissue parameters (e.g., T1 (x) and T2 (x)), imaging system parameters (e.g., the off-resonance frequency f0 (x)), and applied acquisition parameters (e.g., flip angles and repetition times). Spatial discretization of (1) yields Im = Φm (T1 , T2 , f0 ) ρ,

(2) N ×N

where ρ, Im ∈ C , T1 , T2 , f0 ∈ R , and Φm ∈ C diagonal matrix with [Φm ]n,n = ϕm (T1 [n], T2 [n], f0 [n]). Incorporating (2) into the imaging equation, i.e., N

N

is a

dm = Fm Im + nm , where dm ∈ C denotes the measured data, Fm ∈ CPm ×N undersampled discrete Fourier transform matrix, nm ∈ CPm the measurement noise, we can form the following data model for MRF: Pm

dm = Fm Φm (T1 , T2 , f0 ) ρ + nm .

ICIP 2015

With further knowledge of the noise statistics (e.g., nm is often assumed to be complex white Gaussian noise), the unknown parameters {T1 , T2 , f0 , ρ} can be determined using the ML estimation [4].

LA (ρ, θ, gm , Q, Λgm , ΛQ ) =

2.1.2. Image Prior Model

arg

ˆ 1, T ˆ 2, ρ ˆ , ˆf0 T min

}

M ∑

T1 ,T2 ,f0 ,ρ

M ∑

∥dm − Fm gm ∥22 +

m=1

We can incorporate an image prior model with the above data model to potentially improve the performance of parameter estimation. Similar to conventional MR parameter mapping (e.g., [10, 11]), there are a number of models that can be used for the parameter maps {ρ, T1 , T2 , f0 } and/or the contrast-weighted images I = [I1 , · · · , IM ] ∈ CN ×M . Here we impose a low-rank image model for I, motivated by the following observations: 1) there are often very few types of tissues involved in an MRF experiment, and 2) magnetization evolutions associated with the same type of tissue can be highly correlated. Integrating the data model in (3) and a low-rank image model, we can formulate the following penalized maximum likelihood (PML) reconstruction problem: {

We can form the augmented Lagrangian functional LA (·) for (4) as follows:

= ∥dm − Fm Φm (T1 , T2 , f0 ) ρ∥22 + λψr (I) ,

m=1

(3)

M ∑

λψr (Q) +

Re (< Λgm , gm − Φm (θ) ρ >) +

m=1 M ∑ µg µQ ∥gm − Φm (θ) ρ∥22 + ∥G − Q∥2F + 2 2 m=1

Re (< ΛQ , G − Q >) ,

(5)

where Λgm ∈ CN and ΛQ ∈ CN ×M are the Lagrangian multipliers, µg > 0, µQ > 0 are the penalty parameters, and Re(·) denotes the operation of taking the real part of a complex number. Second, we apply the ADMM algorithm [5,7] to iteratively minimize (5). Specifically, at the (k + 1)th iteration, it consists of minimizing LA (·) with respect to ρ, θ, gm , and Q in an alternating manner, and updating the Lagrangian multipliers as follows: ( ) (k+1) (k) ΛQ = ΛQ + µQ G(k+1) − Q(k+1) , [ ( ) ] (k+1) Λ(k+1) = Λ(k) − Φ θ(k+1) ρ(k+1) . gm gm + µgm gm

where ψr (·) denotes a low-rank promoting regularization functional, and λ > 0 is the regularization parameter. Note that thanks to the parametric model in (2), the proposed formulation enables direct parameter estimation from sparsely-sampled k-space data, although the rank constraint is imposed on the contrast-weighted image sequence. Regarding ψr (·), we can either impose the low-rank constraint explicitly through matrix factorization (e.g., [12]), or implicitly with the use of surrogate functions (e.g., [13, 14]). Here, we use a nonconvex surrogate function based on the Schatten p-norm (p < 1) for (∑ p )1/p ψr (·). Mathematically, it can be written as ψr (I) = , i σi where σi denotes the ith singular value of I.

In the following, we describe the solution to each sub-problem in the alternating minimization.

2.2. Solution Algorithm

2.2.2. Updating gm

In this subsection, we describe an iterative algorithm to address the PML reconstruction problem in (3). Note that (3) results in a nonconvex and nonlinear optimization problem. Besides the difficulty arising from the data fidelity term [4], the penalty term ψr (I) is not nondifferentiable, making the problem even more challenging. Here we extend our previously proposed algorithm in [4], based on variable splitting, the alternating direction method of multipliers, and the variable projection method, to efficiently solve (3). For notational simplicity, let θ = {T1 , T2 , f0 } denote the unknown parameters in the contrast weighting function. In the proposed algorithm, we first apply an VS strategy that splits the parametric physical model from the data fidelity term, and the image sequence from the penalty term, i.e., gm = Φm (θ) ρ and G = Q, where G = [g1 , · · · , gM ] ∈ CN ×M . As a consequence, (3) can be converted into the following equivalent constrained optimization problem:

The sub-optimization problem with respect to {gm }M m=1 can be written as {

2 } ∑M µg (k) 2 min ∥F g − d ∥ + − y +

g m m m 2 m m m=1 2

min

∑M

{gm },ρ,θ,Q

m=1

The sub-optimization problem with respect to {ρ, θ} can be written as

( ) 2 M (k)

∑ Λgm

(k) min

Φm (θ) ρ − gm +

. ρ,θ

µgm m=1

2

Note that this is a separable nonlinear least-squares problem, which also appears in the ML reconstruction [4]. Thus, we can use the VARPRO [4, 8, 9] algorithm to update {ρ, θ}.

{gm }M m=1

µQ 2

2

2

G − X(k) ,

( ) (k) where ym = Φm θ(k+1) ρ(k+1) −

(6)

F

(k)

Λgm µgm

, and X(k) = Q(k) −

(k)

ΛQ

. Eq. (6) is a quadratic optimization problem that is separable with respect to gm . Here we can apply the preconditioned conjugate gradient method to solve it. µQ

2.2.3. Updating Q The sub-optimization problem with respect to Q can be written as

∥dm − Fm gm ∥22 + λψr (Q) ,

gm = Φm (θ) ρ and G = Q.

2.2.1. Updating ρ and θ

(4)

3393

2 λ 1

min Q − Tk + ψr (Q) , Q 2 µQ F

(7)

Fig. 1. Ground truth spin density map (a), reconstructed spin density maps (b), and corresponding relative error maps (c), obtained by the conventional MRF, ML-MRF, and MBIR-MRF at the two acquisition lengths (i.e., M = 200 and 400). Note that the associated NRMSE is labeled at the lower right corner of each error map. (k)

Λ

where T(k) = G(k) + µQQ . Its solution can be obtained by the following singular value thresholding [14]: Q

(k+1)

) ∑( λ p−1 = σi − σ ui viH , µQ i + i

(8)

where σi is the ith singular value of Tk , ui and vi are respectively the ith left and right singular vectors of Tk , and (a)+ = max {a, 0} denotes the thresholding operation. 2.2.4. Algorithm initialization and termination Since (3) is a nonconvex nonlinear optimization problem, there is generally no theoretical guarantee for the global convergence of the proposed algorithm. In practice, we observed that the algorithm performance depends on the use of initialization. Similar to [4], we use the gridding reconstruction as an initialization of the algorithm, which consistently yields good empirical performance. Regarding the stopping criteria, we terminate the algorithm when a prespecified maximum number of iterations is reached, or the relative change between the consecutive iterations is smaller than some prespecified tolerance. 3. RESULTS In this section, we show representative results from computer simulations to illustrate the performance of the proposed method. The simulations were performed using a brain phantom from the Brainweb [15] database that contains the spin density, T1 and T2 maps. For simplicity, f0 is assumed to be zero for all the voxels. The Bloch simulation was performed to generate a sequence of contrastweighted images using the balanced steady-state free precession sequence with the same set of flip angles and repetition times as in [1]. For data acquisition, k-space data were acquired using one simulated spiral interleave at each time instant (the fully sampled data

3394

contain 48 interleaves), and different interleaves were used for different time instants. Moreover, complex white Gaussian noise was added to measured data such that SNR = 35 dB. We performed the conventional MRF reconstruction [1], ML reconstruction (called ML-MRF) [4], and the proposed model-based iterative reconstruction (called MBIR-MRF) at two different acquisition lengths, i.e., M = 200 and M = 400. We initialized both MLMRF and MBIR-MRF with the gridding reconstruction. Moreover, we empirically tuned the regularization parameter λ and penalty parameters µQ and µgm to optimize the performance of the proposed method. Fig. 1, 2, and 3 show the reconstructed spin density maps, T1 maps, and T2 maps along with the relative error maps and the overall normalized root mean square error (NRMSE)1 . As expected, the performance of all three methods improves with the increased number of measurements. At each acquisition length, MBIR-MRF outperforms the conventional MRF and ML-MRF both qualitatively and quantitatively. In particular, when the acquisition length is short (i.e., M = 200), the improvements are substantial for the spin density map and T2 map (e.g., significantly reduced image artifacts and NRMSE), illustrating the benefits of introducing low-rank image model into the proposed method. 4. DISCUSSION AND CONCLUSION In this paper, we proposed a novel model-based MRF reconstruction method, incorporating the low-rank image model with the Bloch equation based MR physical model. It performs parameter estimation directly from highly undersampled, noisy k-space data with a PML formulism. An efficient algorithm based on the VS, ADMM, and VARPRO method was developed to solves the resulting optimization problem. The performance of the proposed method was evaluated using computer simulations, and the results have demonstrated that the proposed method leads to improved performance, both qualita1 Note that the NRMSE was calculated with respect to all brain tissues except the skull, scalp, and cerebrospinal fluid, which are typically not regionof-interest.

Fig. 2. Ground truth T1 map (a), reconstructed T1 maps (b), and corresponding relative error maps (c), obtained by the conventional MRF, ML-MRF, and MBIR-MRF at the two acquisition lengths (i.e., M = 200 and 400). Note that the associated NRMSE is labeled at the lower right corner of each error map.

Fig. 3. Ground truth T2 map (a), reconstructed T2 maps (b), and corresponding relative error maps (c), obtained by the conventional MRF, ML-MRF, and MBIR-MRF at the two acquisition lengths (i.e., M = 200 and 400). Note that the associated NRMSE is labeled at the lower right corner of each error map.

tively and quantitatively, over the conventional MRF reconstruction and recently proposed ML reconstruction. Within the PML reconstruction framework, a number of extensions can be performed. For example, we can incorporate other image prior models (e.g., joint sparse structure of contrast-weighted image sequences [11, 16]) into the current formulation, which may lead to better performance. This extension is currently under investigation, and will be reported our future work. Furthermore, it is worthwhile to integrate the proposed method with parallel MR imaging [17], which has the potential of significantly enhancing the practical utility of the proposed method.

3395

5. REFERENCES [1] D. Ma, V. Gulani, N. Seiberlich, K. Liu, J. L. Sunshine, J. L. Duerk, and M. A. Griswold, “Magnetic resonance fingerprinting,” Nature, vol. 495, pp. 187–192, 2013. [2] Z. Wang, Q. Zhang, J. Yuan, and X. Wang, “MRF denoising with compressed sensing and adaptive filtering,” in Proc. IEEE Int. Symp. Biomed. Imaging, 2014, pp. 870–873. [3] M. Davies, G. Puy, P. Vandergheynst, and Y. Wiaux, “Compressed quantitative MRI: Bloch response recovery through iterated projection,” in Proc. IEEE Int. Conf. Acoust. Speech Signal Process., 2014, pp. 6899–6903.

[4] B. Zhao, F. Lam, B. Bilgic, H. Ye, and K. Setsompop, “Maximum likelihood reconstruction for magnetic resonance fingerprinting,” in Proc. IEEE Int. Symp. Biomed. Imaging, 2015, pp. 905–909. [5] S. Ramani and J. A. Fessler, “Parallel MR image reconstruction using augmented Lagrangian methods,” IEEE Trans. Med Imaging, vol. 30, pp. 694–706, 2011. [6] S. G. Lingala, E. DiBella, and M. Jacob, “Deformation corrected compressed sensing (DC-CS): A novel framework for accelerated dynamic MRI,” IEEE Trans. Med Imaging, vol. 34, pp. 72–85, 2015. [7] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Found. Trends Mach. Learn., vol. 3, pp. 1–122, 2011. [8] G. Golub and V. Pereyra, “Separable nonlinear least squares: the variable projection method and its applications,” Inverse Probl., vol. 19, pp. R1–R26, 2003. [9] J. P. Haldar, J. Anderson, and S. W. Sun, “Maximum likelihood estimation of T1 relaxation parameters using VARPRO,” in Proc. Int. Symp. Magn. Reson. Med., 2007, p. 41. [10] B. Zhao, F. Lam, and Z.-P. Liang, “Model-based MR parameter mapping with sparsity constraints: Parameter estimation and performance bounds,” IEEE Trans. Med Imaging, vol. 33, pp. 1832–1844, 2014. [11] B. Zhao, W. Lu, T. K. Hitchens, F. Lam, C. Ho, and Z.-P. Liang, “Accelerated MR parameter mapping with low-rank and sparsity constraints,” Magn. Reson. Med., 2014, In press. [12] B. Zhao, J. P. Haldar, C. Brinegar, and Z.-P. Liang, “Low rank matrix recovery for real-time cardiac MRI,” in Proc. IEEE Int. Symp. Biomed. Imaging, 2010, pp. 996–999. [13] B. Recht, M. Fazel, and P. Parrilo, “Guaranteed minimumrank solutions of linear matrix equations via nuclear norm minimization,” SIAM Rev., vol. 52, pp. 471–501, 2010. [14] S. G. Lingala, Yue Hu, E. DiBella, and M. Jacob, “Accelerated dynamic MRI exploiting sparsity and low-rank structure: kt SLR,” IEEE Trans. Med Imaging, vol. 30, pp. 1042–1054, 2011. [15] D. L. Collins, A. P. Zijdenbos, V. Kollokian, J. G. Sled, N. J. Kabani, C. J. Holmes, and A. C. Evans, “Design and construction of a realistic digital brain phantom,” IEEE Trans. Med Imaging, vol. 17, pp. 463–468, 1998. [16] A. Majumdar and R. K. Ward, “Joint reconstruction of multiecho MR images using correlated sparsity,” Magn Reson Imaging, vol. 29, pp. 899 – 906, 2011. [17] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, and P. Boesiger, “SENSE: Sensitivity encoding for fast MRI,” Magn. Reson. Med., vol. 42, pp. 952–962, 1999.

3396