CUBO A Mathematical Journal Vol.10, N o ¯ 02, (135–144). July 2008 Iterative Regularization Methods for a Discrete Inverse Problem in MRI A. Leitão Universidade Federal de Santa Catarina, Departamento de Matemática, P.O. Box 476, 88040-900, Florianópolis – Brasil email: aleitao@mtm.ufsc.br and J.P. Zubelli IMPA - Instituto Nacional de Matemática Pura e Aplicada, Estrada Dona Castorina 110, 22460-320, Rio de Janeiro – Brasil email: zubelli@impa.br ABSTRACT We propose and investigate efficient numerical methods for inverse problems related to Magnetic Resonance Imaging (MRI). Our goal is to extend the recent convergence results for the Landweber-Kaczmarz method obtained in [7], in order to derive a con- vergent iterative regularization method for an inverse problem in MRI. RESUMEN Nosotros investigamos y proponemos métodos numéricos eficientes para problemas in- versos relacionados con resonancia Magnética de Imagen (MRI). Nuestro objetivo es extender resultados recientes de convergencia para el método de Landweber-Kaczmarz 136 A. Leitão and J.P. Zubelli CUBO 10, 2 (2008) obtenido en [7], a fin de obtener un método de regularización iterativo convergente para un problema inverso en MRI. Key words and phrases: Magnetic Resonance Imaging, MRI, Tomography, Medical Imaging, Inverse Problems. Math. Subj. Class.: 94A08, 47A52, 65D18 1 Introduction Magnetic Resonance Imaging, also known as MR–Imaging or simply MRI, is a noninvasive tech- nique used in medical imaging to visualize body structures and functions, providing detailed images in arbitrary planes. Unlike X-Ray Tomography it does not use ionizing radiation, but uses a pow- erful magnetic field to align the magnetization of hydrogen atoms in the body. Radio waves are used to systematically alter the alignment of this magnetization, causing the hydrogen atoms to produce a rotating magnetic field detectable by the scanner. More specifically, when a subject is in the scanner, the hydrogen nuclei (i.e., protons, found in abundance in the human body as water) align with the strong magnetic field. A radio wave at the correct frequency for the protons to absorb energy pushes some of the protons out of alignment. The protons then snap back to alignment, producing a detectable rotating magnetic field as they do so. Since protons in different areas of the body (e.g., fat and muscle) realign at different speeds, the different structures of the body can be revealed. The image to be identified in MRI corresponds to a complex valued function P : [0, 1]×[0, 1]→C and the image acquisition process is performed by so-called receivers. Due to the physical nature of the acquisition process, the information gained by the receivers does not correspond to the unknown image, but instead, to P multiplied by receiver dependent sensitivity kernels. In real life applications, the sensitivity kernels are not precisely known and have to be identified together with P. This corresponds to a version of the blind deconvolution problem that has been investigated by many authors. See for example [2, 12, 14] Our main goal in this article is to investigate efficient iterative methods of Kaczmarz type for the identification problem related to MRI. We extend the convergence results for the loping Landweber-Kaczmarz method in [7] and derive a convergent iterative regularization method for this inverse problem. This article is outlined as follows. In Section 2 the description of a discrete mathematical model for Magnetic Resonance Imaging is presented. In Section 3 we derive the corresponding inverse problem for MRI. In Section 4 we investigate efficient iterative regularization methods for this inverse problem. Using a particular hypothesis on the sensitivity kernels, we are able to derive convergence and stability results for the proposed iterative methods. CUBO 10, 2 (2008) Iterative Regularization Methods ... 137 2 The direct problem In what follows we present a discrete model for MRI. In our approach, we follow the notation introduced in [1]. The image to be identified is considered to be a discrete function P : {1, . . . , phor} × {1, . . . , pver} =: B → C , where phor, pver ∈ N0 are known. Therefore, the number of degrees of freedom related to this parameter is pnum := phor × pver (typical values are phor = pver = 256; pnum = 65536). As mentioned above, the image acquisition process is performed by several receivers, denoted here by Rj , j = 0, . . . , rnum − 1, where rnum ∈ N0 is given (typically one faces the situation where rnum << pnum). Due to the physical nature of the acquisition process, the information gained by the receivers does not correspond to the unknown image, but instead, to P multiplied by receiver dependent sensitivity kernels Sj = S(Rj ) : B → C , j = 0, . . . , rnum − 1 . In real life applications, the sensitivity kernels Sj are not precisely known and have to be identified together with P. This fact justifies the following ansatz: (A1) The sensitivity kernels Sj can be written as linear combination of the given basis functions Bn : B → C, for n = 1, . . . , bnum, and bnum ∈ N0. In other words, we assume the existence of coefficients bj,n ∈ C such that Sj (m) = bnum∑ n=1 bj,n Bn(m) , m ∈ B , j = 0, . . . , rnum − 1 . (2.1) In the sequel we make use the abbreviated notations bj := (bj,n) bnum n=1 and (bj ) := (bj ) rnum−1 j=0 . Notice that the coefficient vectors bj = b(Rj ) are receiver dependent. The measured data for the inverse problem is given in a subset of the Fourier space of the image P, i.e. there exists a known subset M ⊂ B (consisting of pproj elements) such that the receiver dependent measurement Mj = M(Rj ) satisfies Mj := P[F(P × Sj )] , j = 0, . . . , rnum − 1 . where F is the Discrete Fourier Transform (DFT) operator defined by F : {f | f : B → C} → {f̂ | f̂ : B → C} f 7→ (F(f ))(m) := pnum−1∑ n=0 f (n) exp ( − 2πi pnum nm ) , 138 A. Leitão and J.P. Zubelli CUBO 10, 2 (2008) and P is the operator defined by P : { f | f : B → C } → { g | g : M → C } =: Y f 7→ (P[f ])(m) := f (m) , m ∈ M . Notice that, due to ansatz (A1) and the linearity of F and P, the measured data Mj ∈ Y can be written in the form Mj = bnum∑ n=1 bj,n P[F(P × Bn)] , j = 0, . . . , rnum − 1 . (2.2) Remark 2.1. The numerical evaluation of the DFT requires naively O(p2num) arithmetical opera- tions. However, in practice the DFT must be replaced by the Fast Fourier Transform (FFT), which can be computed by the Cooley-Tukey algorithm1 and requires only O(pnum log(pnum)) operations. 3 The inverse problem Next we use the discrete model discussed in the previous section as a starting point to formulate an inverse problem for MRI. The unknown parameters to be identified are the discrete image function P and the sensitivity kernels Sj . Due to the ansatz (A1), the parameter space X consists of pairs of the form (P, (bj )), i.e. X := { (P, (bj )) ; P ∈ C pnum , (bj ) ∈ (C bnum ) rnum } . It is immediate to observe that X can be identified with C(pnum+bnum×rnum), while Y can be identified with C pproj . The parameter to output operators Fi : X → Y are defined by Fi : (P, (bj )) 7→ bnum∑ n=1 bi,n P[F(P × Bn)] , i = 0, . . . , rnum − 1 . (3.1) Due to the experimental nature of the data acquisition process, we shall assume that the data Mi in (2.2) is not exactly known. Instead, we have only approximate measured data M δ i ∈ Y satisfying ‖M δ i − Mi‖ ≤ δi , (3.2) with δi > 0 (noise level). Therefore, the inverse problem for MRI can be written in the form of the following system of nonlinear equations Fi(P, (bj )) = M δ i , i = 0, . . . , rnum − 1 . (3.3) It is worth noticing that the nonlinear operators Fi’s are continuously Fréchet differentiable, and the derivatives are locally Lipschitz continuous. 1The FFT algorithm was published independently by J.W. Cooley and J.W. Tukey in 1965. However, this algorithm was already known to C.F. Gauss around 1805. CUBO 10, 2 (2008) Iterative Regularization Methods ... 139 4 Iterative regularization In this section we analyze efficient iterative methods for obtaining stable solutions of the inverse problem in (3.3). 4.1 An image identification problem Our first goal is to consider a simplified version of problem (3.3). We assume that the sensitivity kernels Sj are known, and we have to deal with the problem of determining only the image P. This assumption can be justified by the fact that, in practice, one has very good approximations for the sensitivity kernels, while the image P is completely unknown. In this particular case, the inverse problem reduces to F̃i(P) = M δ i , i = 0, . . . , rnum − 1 , (4.1) where F̃i(P) = Fi(P, (bj )), the coefficients (bj ) being known. This is a much simpler problem, since F̃i : X̃ → Y are linear and bounded operators, defined at X̃ := {f | f : B → C}. We follow the approaches in [7, 5] and derive two iterative regularization methods of Kaczmarz type for problem (4.1). Both iterations can be written in the form P δ k+1 = P δ k − ωkαksk , (4.2) where sk := F̃[k](P δ k ) ∗ (F̃[k](P δ k ) − M δ i ) , (4.3) ωk := { 1 ‖F̃[k](P δ k ) − Mδ i ‖ > τ δ[k] 0 otherwise . (4.4) Here τ > 2 is an appropriately chosen constant, [k] := (k mod rnum) ∈ {0, . . . , rnum − 1} (a group of rnum subsequent steps, starting at some multiple k of rnum, is called a cycle), P δ 0 = P0 ∈ X̃ is an initial guess, possibly incorporating some a priori knowledge about the exact image, and αk ≥ 0 are relaxation parameters. Distinct choices for the relaxation parameters αk lead to the definition of the two iterative methods: 1) If αk is defined by αk := { ‖sk‖ 2/‖F̃[k](P δ k )sk‖ 2 ωk = 1 1 ωk = 0 , (4.5) the iteration (4.2) corresponds to the loping Steepest-Descent Kaczmarz method (lSDK) [5]. 2) If αk ≡ 1, the iteration (4.2) corresponds to the loping Landweber-Kaczmarz method (lLK) [7]. 140 A. Leitão and J.P. Zubelli CUBO 10, 2 (2008) The iterations should be terminated when, for the first time, all Pk are equal within a cycle. That is, we stop the iteration at the index kδ ∗ , which is the smallest multiple of rnum such that P k δ ∗ = P k δ ∗ +1 = · · · = Pkδ ∗ +rnum−1 . (4.6) Convergence analysis of the lSDK method From (3.1) follows that the operators F̃i are linear and bounded. Therefore, there exist M > 0 such that ‖F̃i‖ ≤ M , i = 0, . . . , rnum − 1 . (4.7) Since the operators F̃i are linear, the local tangential cone condition is trivially satisfied (see (4.16) below). Thus, the constant τ in (4.4) can be chosen such that τ > 2. Moreover, we assume the existence of P ∗ ∈ Bρ/2(P0) such that F̃i(P ∗ ) = Mi , i = 0, . . . , rnum − 1 , (4.8) where ρ > 0 and (Mi) rnum−1 i=0 ∈ Y rnum corresponds to exact data satisfying (3.2). In the sequel we summarize several properties of the lSDK iteration. For a complete proof of the results we refer the reader to [5, Section 2]. Lemma 4.1. Let the coefficients αk be defined as in (4.5), the assumption (4.8) be satisfied for some P∗ ∈ X̃, and the stopping index kδ ∗ be defined as in (4.6). Then, the following assertions hold: 1) There exists a constant α > 0 such that αk > α, for k = 0, . . . , k δ ∗ . 2) Let δi > 0 be defined as in (3.2). Then the stopping index k δ ∗ defined in (4.6) is finite. 3) Pδ k ∈ Bρ/2(P0) for all k ≤ k δ ∗ . 4) The following monotony property is satisfied: ‖P δ k+1 − P ∗ ‖ 2 ≤ ‖P δ k − P ∗ ‖ 2 , k = 0, 1, . . . , kδ ∗ , (4.9) ‖P δ k+1 − P ∗ ‖ 2 = ‖P δ k − P ∗ ‖ 2 , k > kδ ∗ . (4.10) Moreover, in the case of noisy data (i.e. δi > 0) we have ‖F̃i(P δ k δ ∗ ) − M δ i ‖ ≤ τ δi , i = 0, . . . , rnum − 1 . (4.11) Next we prove that the lSDK method is a convergent regularization method in the sense of [3]. Theorem 4.2 (Convergence). Let αk be defined as in (4.5), the assumption (4.8) be satisfied for some P∗ ∈ X̃, and the data be exact, i.e. Mδ i = Mi in (3.2). Then, the sequence P δ k defined in (4.2) converges to a solution of (4.1) as k → ∞. CUBO 10, 2 (2008) Iterative Regularization Methods ... 141 Proof. Notice that, since the data is exact, we have ωk = 1 for all k > 0. The proof follows from [5, Theorem 3.5]. Theorem 4.3 (Stability). Let the coefficients αk be defined as in (4.5), and the assumption (4.8) be satisfied for some P∗ ∈ X̃. Moreover, let the sequence {(δ1,m, . . . , δrnum,m)}m∈N (or simply {δm}m∈N) be such that limm→∞(maxi δi,m) = 0, and let M δm i be a corresponding sequence of noisy data satisfying (3.2) (i.e. ‖Mδm i − Mi‖ ≤ δi,m, i = 0, . . . , rnum − 1, m ∈ N). For each m ∈ N, let km ∗ be the stopping index defined in (4.6). Then, the lSDK iterates Pδm k m ∗ converge to a solution of (4.1) as m → ∞. Proof. The proof follows from [5, Theorem 3.6]. Convergence analysis of the lLK method The convergence analysis results for the lLK iteration are analog to the ones presented in Theo- rems 4.2 and 4.3 for the lSDK method. In the sequel we summarize the main results that we could extend to the lLK iteration. Theorem 4.4 (Convergence Analysis). Let αk ≡ 1, the assumption (4.8) be satisfied for some P∗ ∈ X̃, the operators F̃i satisfy (4.7) with M = 1, and the stopping index k δ ∗ be defined as in (4.6). Then, the following assertions hold: 1) Let δi > 0 in (3.2). Then the stopping index k δ ∗ defined in (4.6) is finite. 2) Pδ k ∈ Bρ/2(P0) for all k ≤ k δ ∗ . 3) The monotony property in (4.9), (4.10) is satisfied. Moreover, in the case of noisy data, (4.11) holds true. 4) For exact data, i.e. δi = 0 in (3.2), the sequence P δ k defined in (4.2) converges to a solution of (4.1) as k → ∞. 5) Let the sequence {δm}m∈N, the corresponding sequence of noisy data M δm i , and the stopping indexes km ∗ be defined as in Theorem 4.3. Then, the lLK iterates Pδm k m ∗ converge to a solution of (4.1) as m → ∞. Proof. The proof follows from corresponding results for the lLK iteration for systems of nonlinear equations in [7]. Notice that the assumption M = 1 in Theorem 4.4 is nonrestrictive. Indeed, since the op- erators F̃i are linear, it is enough to scale the equations in (4.1) with appropriate multiplicative constants. 142 A. Leitão and J.P. Zubelli CUBO 10, 2 (2008) 4.2 Identification of image and sensitivity Our next goal is to consider the problem of determining both the image P as well as the sensitivity kernels Sj in (3.3). The lLK and lSDK iterations can be extended to the nonlinear system in a straightforward way (P δ k+1, (bj ) δ k+1) = (P δ k , (bj ) δ k ) − ωkαksk , (4.12) where sk := F ′ [k](P δ k , (bj ) δ k ) ∗ (F[k](P δ k , (bj ) δ k ) − M δ i ) , (4.13) ωk := { 1 ‖F[k](P δ k , (bj ) δ k ) − Mδ i ‖ > τ δ[k] 0 otherwise . (4.14) In the lLK iteration we choose αk ≡ 1, and in the lSDK iteration we choose αk := { ‖sk‖ 2/‖F ′[k](P δ k , (bj ) δ k )sk‖ 2 ωk = 1 1 ωk = 0 . (4.15) In order to extend the convergence results in [7, 5] for these iterations, we basically have to prove two facts: 1) Assumption (14) in [7]. 2) The local tangential cone condition [7, Eq. (15)], i.e. the existence of (P0, (bj )0) ∈ X and η < 1/2 such that ‖Fi(P, (bj )) − Fi(P̄, (b̄j )) − F ′ i (P, (bj ))[(P, (bj )) − (P̄, (b̄j ))]‖Y ≤ η‖Fi(P, (bj )) − Fi(P̄, (b̄j ))‖Y , (4.16) for all (P, (bj )), (P̄, (b̄j )) ∈ Bρ(P0, (bj )0), and all i = 1, . . . , rnum. The first one represents no problem. Indeed, the Fréchet derivatives of the operators Fi are locally Lipschitz continuous. Thus, for any (P0, (bj )0) ∈ X and any ρ > 0 we have ‖F ′ i (P, (bj ))‖ ≤ M = Mρ,P0,(bj )0 for all (P, (bj )) in the ball Bρ(P0, (bj )0) ⊂ X. The local tangential cone condition however, does not hold. Indeed, the operators Fi are second order polynomials of the variables bj,n and P. Therefore, it is enough to verify whether the real function f (x, y) = xy satisfies |f (x, y) − f (x̄, ȳ) − f ′(x, y)((x − x̄, y − ȳ))| ≤ η|f (x, y) − f (x̄, ȳ)| , in some vicinity of a point (x0, y0) ∈ R 2 containing a local minimizer of f . This, however, is not the case. CUBO 10, 2 (2008) Iterative Regularization Methods ... 143 Therefore, the techniques used to prove convergence of the lLK and lSDK iterations in [7, 5] cannot be extended to the nonlinear system (3.3). It is worth noticing that the local tangential cone condition is a standard assumption in the convergence analysis of adjoint type methods (Landweber, steepest descent, Levenberg-Marquardt, asymptotical regularization) for nonlinear inverse problems [3, 4, 6, 9, 10, 11, 13]. Thus, none of the classical convergence proofs for these iterative methods can be extended to system (3.3) in a straightforward way. Motivated by the promising numerical results and efficient performance of the lLK and lSDK iterations for problems known not to satisfy the local tangential cone condition (see [7, 8, 5]), we intend to use iteration (4.12) for computing approximate solutions of system (3.3). This numerical investigation will be performed in a forthcoming article. 5 Conclusions We presented the description of a discrete mathematical model for Magnetic Resonance Imaging and derived the corresponding inverse problem for MRI. We investigate efficient iterative regularization methods for this inverse problem. An iterative method of Kaczmarz type for obtaining approximate solutions for the inverse problem is proposed. Using a particular assumption on the sensitivity kernels, we are able to prove convergence and stability results for the proposed iterative methods. The convergence analysis presented in this article extends the results for the loping Landweber- Kaczmarz method in [7]. Moreover, we prove that our method is a convergent iterative regulariza- tion method in the sense of [3]. Acknowledgments The work of A.L. was supported by the Brazilian National Research Council CNPq, grants 306020/2006- 8 and 474593/2007-0. J.P.Z. was supported by CNPq under grants 302161/2003-1 and 474085/2003- 1. This work was developed during the permanence of the authors in the Special Semester on Quantitative Biology Analyzed by Mathematical Methods, October 1st, 2007 -January 27th, 2008, organized by RICAM, Austrian Academy of Sciences. Received: March 2008. Revised: May 2008. 144 A. Leitão and J.P. Zubelli CUBO 10, 2 (2008) References [1] F. Bauer and S. Kannengiesser, An alternative approach to the image reconstruction for parallel data acquisition in MRI, Mathematical Methods in the Applied Sciences, 30 (2007), 1437–1451 [2] M. Bertero and P. Boccacci, Introduction to inverse problems in imaging, Bristol: IoP, Institute of Physics Publishing, 1998. [3] H. Engl, M. Hanke and A. Neubauer, Regularization of Inverse Problems, Kluwer, Dordrecht, 1996. [4] H. Engl and A. Leitão, A Mann iterative regularization method for elliptic Cauchy problems, Numer. Funct. Anal. Optim., 22 (2001), 861–864. [5] A. De Cezaro, M. Haltmeier, A. Leitão and O. Scherzer, On Steepest-Descent Kaczmarz methods for regularizing nonlinear systems of ill-posed equations, Applied Math- ematics and Computation, to appear, 2008. [6] P. Deuflhard, H.W. Engl and O. Scherzer, A convergence analysis of iterative methods for the solution of nonlinear ill–posed problems under affinely invariant conditions, Inverse Problems, 14 (1998), 1081–1106. [7] M. Haltmeier, A. Leitão and O. Scherzer, Kaczmarz methods for regularizing nonlinear ill-posed equations. I. Convergence analysis, Inverse Probl. Imaging, 1 (2007), 289–298. [8] M. Haltmeier, R. Kowar, A. Leitão and O. Scherzer, Kaczmarz methods for reg- ularizing nonlinear ill-posed equations. II. Applications, Inverse Probl. Imaging, 1 (2007), 507–523. [9] M. Hanke, A. Neubauer and O. Scherzer, A convergence analysis of Landweber iteration for nonlinear ill-posed problems, Numerische Mathematik, 72 (1995), 21–37. [10] B. Kaltenbacher, A. Neubauer and O. Scherzer, Iterative Regularization Methods for Nonlinear Ill–Posed Problems, Springer Verlag, 2008, to appear. [11] R. Kowar and O. Scherzer, Convergence analysis of a Landweber-Kaczmarz method for solving nonlinear ill-posed problems, Ill-posed and inverse problems, 253–270, VSP, Zeist, 2002. [12] F. Natterer and F. Wuebbeling, Mathematical methods in image reconstruction, SIAM Monographs on Mathematical Modeling and Computation, 5, Philadelphia, PA, 2007. [13] U. Tautenhahn, On the asymptotical regularization of nonlinear ill-posed problems, In- verse Problems, 10 (1994), 1405–1418. [14] J. P. Zubelli, R. Marabini, C. Sorzano, and G. Herman, Three-dimensional re- construction by Chanine’s method from electron microscopic projections corrupted by in- strumental aberrations, Inverse Probl. 19, 4 (2003), 933–949. N9