MMANUSCRIPT SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE 1
Maxwell Parallel Imaging
Matteo Alessandro Francavilla, Stamatios Lefkimmiatis, Jorge F. Villena, and Athanasios G. Polimeridis
Abstract — Purpose : To develop a general framework for Parallel Imaging(PI) with the use of Maxwell regularization for the estimation ofthe sensitivity maps (SMs) and constrained optimization for theparameter-free image reconstruction.
Theory and Methods : Certain characteristics of both the SMs andthe images are routinely used to regularize the otherwise ill-posedoptimization-based joint reconstruction from highly acceleratedPI data. In this paper we rely on a fundamental property ofSMs–they are solutions of Maxwell equations– we construct thesubspace of all possible SM distributions supported in a givenfield-of-view, and we promote solutions of SMs that belong inthis subspace. In addition, we propose a constrained optimizationscheme for the image reconstruction, as a second step, once anaccurate estimation of the SMs is available. The resulting method,dubbed Maxwell Parallel Imaging (MPI), works seamlessly forarbitrary sequences (both 2D and 3D) with any trajectory andminimal calibration signals.
Results : The effectiveness of MPI is illustrated for a wide range ofdatasets with various undersampling schemes, including radial,variable-density Poisson-disc, and Cartesian, and is comparedagainst the state-of-the-art PI methods. Finally, we include somenumerical experiments that demonstrate the memory footprintreduction of the constructed Maxwell basis with the help of tensordecomposition, thus allowing the use of MPI for full 3D imagereconstructions.
Conclusions : The MPI framework provides a physics-inspiredoptimization method for the accurate and efficient image recon-struction from arbitrary accelerated scans.
Index Terms —constrained optimization, electromagnetic basis,Maxwell regularization, parallel imaging, tensor decomposition.
I. I
NTRODUCTION
Parallel Imaging (PI) is admittedly one of the most dis-ruptive technologies in modern magnetic resonance imaging(MRI) and probably the best example of a successful tran-sition from academic research to widespread usage in clinic.Essentially, PI exploits the multi-physic nature of MRI and theubiquitous use of sophisticated spatially-distributed receivingcoils in order to significantly reduce the scan time. Indeed,the interplay of electrodynamics and spin-dynamics in thespatiotemporal encoding, as evinced by the bilinear form ofthe MR signal equation, suggests that the spatial selectivity ofthe receivers could be harnessed in order to reduce the time-consuming gradient encoding.There is a plethora of PI reconstruction methods that couldbe roughly categorized into two main approaches: the image-space (or spatial-domain) and the k-space (or spectral-domain).As main representatives of the former approach, which callsfor the a-priori knowledge of the associated sensitivity maps(SMs), one can mention the pioneering works of SMASH [1]and SENSE [2]. The k-space methods followed a few years
All authors are with Q Bio Inc., San Carlos, CA 94070, USA.Manuscript submitted to Magnetic Resonance in Medicine. later aiming exactly at breaking the dependence of separatepre-calibration scans, which increase the overall acquisitiontime and are more susceptible to motion artifacts. The be-ginning of those so-called auto-calibrating methods can beidentified with the emergence of GRAPPA [3], which makesuse of some extra auto-calibration signals (ACS) in order tofit the kernels for approximating the missing k-space lines. Amore detailed description of all the methods developed in theearly days of PI can be found in the review paper [4].The first PI techniques, both image-space and k-space,were geared to fast reconstruction times, allowing certainsimplifications at the expense of extra pre-calibration scans orACSs in order to transform the inherently non-linear probleminto a linear one. Naturally, more sophisticated PI methodsfollowed that consider the original bilinear form of the inverseproblem at hand, incorporating the estimation of the coil SMs.The common point of the most notable among them (JSENSE[5] and NLINV [6]) is the use of appropriate regularization,necessary for the otherwise ill-posed inverse problem. Morespecifically, in both methods the authors exploit the smooth-ness of the SMs by making use of a polynomial expansionfor constraining the subspace of the possible solutions of theSMs in the former while applying a smoothness-enforcingregularization term in the latter. Recently, NLINV was furthergeneralised to a method dubbed ENLIVE with the addition ofextra bilinear forms in order to account for the violation ofthe standard model in case of limited field-of-view (FOV) [7].Another aspect of the smoothness and the spatial selectivityof the SMs is that they also favor purely algebraic techniquesbased on modern numerical linear algebra algorithms thatpromote low-rank and subspace-specific solutions [8]–[11].As the above-mentioned iterative PI reconstruction ap-proaches started gaining more traction, the interest shiftedtowards the use of more expressive regularizers, rangingfrom ones readily available in the mathematical optimizationliterature [12]–[14] to more modern data-driven variationalmodels [15]. Again, it became clear that although the jointreconstruction of both the SMs and the images was offeringcertain advantages, there was a strong argument for consider-ing the SMs estimation first and then using those SMs in thesolution of the linear image reconstruction. This justifies thefurther proliferation of numerical methods that are tailored tothe accurate estimation of the SMs [16]–[20]; among themESPIRiT [19] deserves a special mention as it appears tobe a true workhorse and the method of choice for most ofthe recent studies, including the benchmark challenge forthe deep-learning PI reconstruction techniques [21]. Morespecifically, ESPIRiT is based on an eigenvalue decompositionof an image-domain operator, and essentially exploits thesmoothness of the SMs and the rank-deficient properties ofthe calibration matrix. a r X i v : . [ ee ss . I V ] A ug ANUSCRIPT SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE 2
Evidently, the modern PI reconstruction techniques havegone a long way from the first days of accelerated MR scansand today it is quite common to use more sophisticated meth-ods of linear and bilinear numerical optimization as well asdeep-learning for reconstructing both the SMs and the images.Nevertheless, even the most effective PI methods availabletoday are based on regularizers that are oversimplified and/orrequire case-dependent fine tuning of the penalty parameters.In this work, we develop a general PI framework that relies onphysics-inspired regularization for the estimation of the SMsand parameter-free constrained optimization for the imagereconstruction. More specifically, we note that smoothness isonly one of the characteristics of SMs that depends, amongother, on the scanner’s main field strength. Foremost, SMsare solutions of Maxwell equations; they correspond to themagnetic fields collected by the receiving coils in the presenceof the patient. Hence we choose to generate the subspaceof the associated SMs (i.e. a complete numerical basis ofmagnetic fields in the FOV) in a patient-agnostic fashion andwe proceed to the solution of the regularized bilinear opti-mization problem, where the SMs are expressed as arbitrarylinear combinations of the elements of the Maxwell basis. Inaddition, we make use of a tensor compression scheme forreducing the memory footprint of the Maxwell basis in thecase of 3D reconstruction. Finally, we appreciate the needfor more expressive regularizers and we propose a parameter-free, constrained optimization scheme for improving the imagequality when SMs are available. The effectiveness of theproposed general PI framework, dubbed Maxwell ParallelImaging (MPI), is demonstrated for a wide range of typicalsequences (both 2D and 3D) with various reduction factors(R) and ACSs. II. T
HEORY AND M ETHODS
A. Problem Formulation
We consider the discretized form of the PI problem, whichcan be described by the forward model y = FSp + n , (1)where p is an N -dimensional vector that contains in rasterizedform the samples of the unknown density to be reconstructed(e.g., a 2-D MRI slice, a 3-D MRI volume, or a 4-D MRImulti-contrast tensor), and y , n ∈ C KC are column vectorscorresponding to the k -space samples obtained from the C re-ceiver coils and i.i.d Gaussian noise, respectively. Furthermore, S ∈ C NC × N is a matrix composed as S = (cid:2) S H . . . S H C (cid:3) H ,where S k ∈ C N × N is a diagonal matrix constructed by theSM s k ∈ C N of the k th coil, k = 1 , . . . , C , ( · ) H denotes theHermittian transpose, and F ∈ C KC × NC is a block diagonalmatrix obtained as I C ⊗ F s , where I C ∈ R C × C is the identitymatrix, ⊗ denotes the Kronecker product, and F s ∈ C K × N is the undersampled operator that provides a mapping fromimage space to k -space, with K ≤ N . The nominal R isdefined as N/K and corresponds to the undersampling rate ofthe k -space.The recovery of the underlying density p from the acquired k -space data y belongs to the category of inverse problems. Due to the presence of noise n , whose exact realization isunknown, and since the operator F is singular, it is an ill-posed problem [22]. This implies that in order to obtain astatistically or physically meaningful solution, we need toexploit any prior knowledge we might have about the solution.Another complicating factor that makes the recovery of p even more challenging, is that the SMs embedded in S aretypically unknown and need to be also recovered. This resultsin an observation model that is not anymore linear w.r.t. theunknown quantities, but instead has the following bilinearform: y = G ( p , s , . . . , s C ) + n , (2)where G (cid:16) x ≡ (cid:2) p H , s H , . . . , s H C (cid:3) H (cid:17) = F s ( s (cid:12) p ) ... F s ( s C (cid:12) p ) (3)and (cid:12) indicates element-wise multiplication of vectors. B. Regularized Nonlinear Inversion
One popular way to tackle the joint recovery problem of p and { s k } Ck =1 is to employ the Iteratively Regularized Gauss-Newton (IRGN) method that was introduced in [23] and waslater used for PI reconstruction in [6]. The underlying idea ofthis approach consists of (a) considering the linearized approx-imation of the nonlinear operator G ( x ) around some currentestimate of the solution, G ( x n + ∆ x ) ≈ G ( x n )+ J G ( x n ) ∆ x ,where J G ( x n ) is the Jacobian of G evaluated at x n , (b)minimizing an objective function of the form: ∆ x ∗ = arg min ∆ x (cid:107) ( y − G ( x n )) − J G ( x n ) ∆ x (cid:107) + R ( x n + ∆ x ) , (4)where R ( · ) is a regularization functional and (c) updatingthe current estimate as x n +1 = x n + γ ∆ x ∗ , where γ can becomputed using a line-search strategy.Initially, in [6] the authors considered using the regularizer R ( x n ) = α n (cid:107) p n (cid:107) + β n C (cid:80) k =1 (cid:13)(cid:13)(cid:13) W ˜ Fs nk (cid:13)(cid:13)(cid:13) , where ˜ F ∈ C N × N is the DFT matrix, W ∈ R N × N is a weighting diagonal matrixand α n , β n ≥ . Note that, given that the SMs are expected tobe smooth, the second term of this regularizer penalizes theirhigh-frequency content. Later, in [14] the authors replacedthe Tikhonov regularizer on the density, (cid:107) p n (cid:107) , with non-quadratic regularizers that can better model certain propertiesof the underlying density, at the cost of a more involvedminimization strategy; the minimizer of Eq. (4) cannot bederived anymore as the solution of a system of linear equationsand more advanced convex optimization techniques must beemployed.In this work, we also rely on the IRGN method as an initialstep that provides an estimate of the unknown SMs. Then ina second step, as we describe later, we use the estimated SMsin order to recover a high-quality estimate of the underlyingdensity by solving the linear inverse problem of Eq. (1).We note that such a two-step strategy has been regularly ANUSCRIPT SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE 3 followed in other image processing applications, such as blinddeconvolution [24], where apart from the underlying imagethe degradation operator is also unknown. Unlike Refs. [6],[14], we consider a modified version of the objective functionin (4), where instead of regularizing directly the k th SM s k ,we penalize its expansion coefficients α k on a predefinedsubset of basis vectors U ∈ C N × q . Details on the constructionof an appropriate physics-inspired basis are provided in thefollowing Sections. In particular, we express each SM as s k ≈ U α k and thus, our observation model takes the form: y = G ( p , α , . . . , α C ) + n , (5)where G (cid:16) ˆ x ≡ (cid:2) p H , α H , . . . , α H C (cid:3) H (cid:17) = (cid:34) F s ( U α (cid:12) p ) ... F s ( U α C (cid:12) p ) (cid:35) . Then, we seek for the solution of the following minimizationproblem: ∆ˆ x ∗ = arg min ∆ p , ∆ α ,..., ∆ α C (cid:107) ( y − G (ˆ x n )) − J G (ˆ x n ) ∆ˆ x (cid:107) + α n (cid:107) p n + ∆ p (cid:107) + β n C (cid:88) k =1 (cid:107) α nk + ∆ α k (cid:107) , (6)where we impose an (cid:96) -squared penalty both on the densityand the expansion coefficients of the SMs. It is worth noticingthat by estimating the expansion coefficients α instead ofthe SMs themselves and since the coils are represented ina reduced order model , i.e. N >> q , the solution of Eq. (6)corresponds to that of an over-determined problem. Due to thequadratic form of the objective function to be minimized inEq. (6), the solution can be derived by solving the relevantnormal equations using the conjugate gradient method [25].This requires the ability to compute the matrix-vector productsof the Jacobian of G and its adjoint, with a vector. Theseproducts are computed as follows: J G (ˆ x ) ∆ p ∆ α ... ∆ α C = G (cid:16)(cid:2) ∆ p H , α H , . . . , α H C (cid:3) H (cid:17) + G (cid:16)(cid:2) p H , ∆ α H , . . . , ∆ α H C (cid:3) H (cid:17) = (cid:34) F s ( U α (cid:12) ∆ p + p (cid:12) U ∆ α ) ... F s ( U α C (cid:12) ∆ p + p (cid:12) U ∆ α C ) (cid:35) (7)and J H G (ˆ x ) y ... y C = C (cid:80) k =1 ( U α k ) H (cid:12) F H s y k U H (cid:0) p H (cid:12) F H s y (cid:1) ... U H (cid:0) p H (cid:12) F H s y C (cid:1) . Finally, it’s worth noting that the same procedure can befollowed in the case of dynamic or multi-contrast PI data,expecting that the associated artifacts will be captured by thedensity while the estimation of the (constrained) SMs willremain unaffected, especially when using some average of the k -space data. C. Regularized Density Reconstruction via Constrained Opti-mization
While the regularization applied on p in Eq. (6) is ratherplain and thus, not capable of modeling complex propertiesof the underlying density, it allows us to perform a jointreconstruction of the sensitivities and the density withouthaving to rely upon a computationally heavy and time con-suming minimization scheme. Furthermore, due to the implicitregularization of the SMs, by expressing them in terms of aproper basis expansion, and the explicit Tikhonov regulariza-tion of the corresponding expansion coefficients, we expectthat most of the reconstruction errors will be accumulated inthe recovered density, while the unknown SMs will be moreaccurately restored.Having this in mind, we use the estimated SMs, discardthe estimated density and solve the linear inverse problem de-scribed in Eq. (1). Hence, we obtain the final density estimateas the minimizer of the following constrained optimizationproblem: p ∗ = arg min p ∈ C N (cid:107) y k − F s S k p (cid:107) ≤ ε k , ∀ k R ( p ) , (8)where ε k is a scalar that is proportional to the standarddeviation of the complex Gaussian noise realization thatdegrades the k -space measurements acquired from the k thcoil. While for the experiments that we report in this work,we have considered Total Variation [26] as the regularizationfunctional R ( p ) of choice, the minimization strategy that weoutline next can be also used without modifications whendifferent and more expressive regularizers are considered,such as the Structure Tensor Total Variation (STV) [27] andit’s non-local extension [28] or the Hessian-Schatten normregularizers of [29]. We also note that one particular advantageof the above constrained problem formulation, compared to theunconstrained minimization approach that is typically pursuedin PI reconstruction, is that in this case there is no need of fine-tuning any regularization penalty parameter, which in practiceis not a straightforward task and requires a certain level ofexperience from the user. The only parameters involved in theabove formulation, are the scalars ε k which can be directlyestimated from the k -space measurements.Now, let us first note that the constrained formulation ofEq. (8) can be equivalently expressed in the unconstrainedform p ∗ = arg min p ∈ C N R ( p ) + C (cid:88) k =1 ι C ( y k ,ε k ) ( F s S k p ) , (9)where ι C ( y k ,ε k ) ( z ) = (cid:40) , if (cid:107) y k − z (cid:107) ≤ ε k ∞ , otherwiseis an indicator function which ensures that the imposed con-straints are satisfied by the solution. Next, since the trans-formed problem is still hard to solve directly, we rely on theAlternating Direction Method of Multipliers (ADMM) [30],[31]. The strategy of ADMM is to split the original problemin smaller and easier ones to solve, by decoupling the different ANUSCRIPT SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE 4 terms of the objective function. Based on this idea andfollowing a similar splitting strategy as in [32], we insteadconsider the problem min Ap + Bz = R ( z ) + C (cid:88) k =1 ι C ( y k ,ε k ) ( z k ) , (10)where A = (cid:104) I N , ( F s S ) H , . . . , ( F s S C ) H (cid:105) H , B = − I ( N + KC ) and z = (cid:2) z H . . . z H C (cid:3) H ∈ C N + KC . Then, using the scaled form of ADMM [31] we obtain the solution to our original problemof Eq. (9) in an iterative way, where each iteration involvesthe following update steps: z n +10 = prox /ρ ·R ( z n − ( p n + u n )) , z n +1 k = Π C ( y k ,ε k ) ( F s S k p n + u nk ) , ∀ k = 1 , . . . , C, p n +1 = (cid:32) I N + C (cid:88) k =1 S H k F H s F s S k (cid:33) − (cid:32) z n +10 − u n + C (cid:88) k =1 S H k F H s (cid:0) z n +1 k − u nk (cid:1)(cid:33) , u n +1 = u n + Ap − z . (11)In Eq. (11) we have that Π C ( y ,ε ) ( z ) = y + ε ( z − y )max ( (cid:107) z − y (cid:107) ,ε ) ,prox /ρ R ( z ) = arg min x ρ/ (cid:107) x − z (cid:107) + R ( x ) is theproximal operator [33] of the regularizer R ( · ) /ρ , u = (cid:2) u H . . . u H C (cid:3) H ∈ C N + KC are the dual variables , and ρ isthe ADMM penalty parameter. In order to avoid fine-tuningthe ADMM penalty parameter ρ , whose value can affect theconvergence rate of the minimization algorithm, we adaptivelychoose its value in each iteration so as to balance the primaland dual residuals (see [31] for their definitions), as proposedin [34].The linear reconstruction algorithm that we proposed aboveis general enough to accommodate for different MRI acquisi-tion modalities. In particular, the steps described in Eq. (11)can also be applied when multicontrast or dynamic MRI areconsidered. The main difference is that for multicontrast MRI,the regularizer R ( · ) instead of being applied only on thespatial dimensions of the underlying density, it should also acton the different contrast channels so that it accounts for the de-pendencies that exist among them. A possible regularizer thatcan be used for this task is the Vectorial Total Variation [35].As far as it concerns the dynamic MRI case, the solution canbe expressed as the minimization of a constrained problemvery similar to the one in Eq. (8), p ∗ = arg min p ∈ C N × T (cid:107) y k,t − F ts S k p t (cid:107) ≤ ε k,t , ∀ k,t R ( p ) , (12)where p = [ p . . . p T ] , y k,t , ε k,t correspond to the k -spacemeasurements and a scalar proportional to the standard devi-ation of the noise from the k th coil and the t th time instance,respectively, F tk is the undersampled mapping operator usedat time instance t = 1 , . . . , T and R ( p ) is a spatiotemporalregularizer. Then, one can follow the strategy described aboveto obtain a slightly modified version of the algorithmic stepsprovided in Eq. (11). D. Maxwell Regularization
A key ingredient of the proposed non linear inversionscheme is the physics inspired regularization of the coil model.Because SMs s k are solutions to Maxwell equations, wepropose to constrain the solution space of the imaging problemto a subspace where SMs are indeed solutions to Maxwellequations, and conjecture that it is possible to express each SMas s k ≈ U h α k , where α k ∈ C q is a column vector collectingthe expansion coefficients of the k th coil, and U h ∈ C N × q isa proper change of basis matrix, referred to as Maxwell basisin the following. The dimension q of the basis will play animportant role in controlling the accuracy of the representationand the regularization properties.Because the basis U h collects solutions of Maxwell equa-tions, a scheme for solving Maxwell equations within theFOV is a prerequisite. One approach is based on Love’sform of the field equivalence theorem [36], [37]: fields insidea source-free volume are fully determined if the tangentialelectromagnetic (EM) fields on the boundary of the volume areknown. Following this idea, the problem of finding volumetricfields inside the source-free FOV is conveniently addressedas a two-step procedure: first, solve for equivalent electric ( j )and magnetic ( m ) currents on the boundary. Subsequently, EMfields e and h inside the FOV are expanded as (cid:34) eh (cid:35) = (cid:34) K ej K em K hj K hm (cid:35) (cid:34) jm (cid:35) , (13)where K αβ is a linear integro-differential operator defined as K αβ f = ˆ R G αβ ( r , r (cid:48) ) · f ( r (cid:48) ) d r (cid:48) , (14)and G αβ is the dyadic Green function, mapping β -kind cur-rents to α -kind fields.Note that j and m are only proxies for computing h : we areinterested in finding a basis to represent a basis for all possiblerealizations of h . This reflects into the need of spanning therange of the integral operators K hj and K hm , and not in aparticular solution of Eq. (13).One way to obtain an orthonormal basis U h is to computethe left singular vectors of K = (cid:2) K hj K hm (cid:3) , (15)where K hj and K hm are the discrete representations ofoperators K hj and K hm , respectively. Then, if K = UΣV H is a Singular Value Decomposition (SVD) of K , U is anorthonormal basis for the range of K . A significant advantageof obtaining the basis via SVD is that it also provides theoptimal low-rank approximation of K . If K ∈ C m × n , amongall matrices B ∈ C m × n with rank k , the one obtained bytruncated SVD is the one with minimum error w.r.t. thespectral norm, (cid:13)(cid:13) K − U k Σ k V H k (cid:13)(cid:13) = σ k +1 (cid:107) K (cid:107) , (16)where only the k column vectors of U and V correspondingto the k largest singular values are kept, and σ k +1 is the ( k +1) -th singular value of K . In other terms, K k ≡ U k Σ k V H k approximates K with error σ k +1 . In turn, by defining U h ≡ ANUSCRIPT SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE 5 q=1 q=2 q=10 q=20 q=100 q=200
Fig. 1: Elements of the basis: pictorial representation of a selection of elements of the magnetic field basis U h adoptedto represent coil SMs. The Maxwell basis U h is computed for a 2D FOV and matrix of size 220x220mm and 320x320,respectively. The number of random excitations to sample the range of the integral operator via randomized SVD is 500. Thebasis is computed over a support covering the full FOV (top two rows), and with a circular support with a diameter of 220mm(bottom two rows). For the two scenarios, the magnitude (rows 1 and 3, arbitrary units) and phase (rows 2 and 4, representedon a cyclical HSV color scale to suppress π phase jumps) are plotted for the basis vectors u q , with q ∈ [1 , , , , , .Because of the SVD-based scheme, U h is a spectral basis: vectors are ordered according to their spatial frequency content,higher index vectors modelling faster variations (from left to right). As a consequence, increasing the dimension of the basisincreases the high-frequency components of the SM that can be captured by U h . U k we have an orthonormal basis to approximate the rangeof K with error σ k +1 .Unfortunately, evaluation of the dyadic Green functions inEq. (13) requires knowledge of the object to be imaged: thisimplies that U h is acquisition dependent, which would clearlybe a major limitation. However, in view of the investigationsdocumented in [38] and references therein, at MRI frequenciesthe magnetic field is only slightly perturbed by the biologicaltissue, due to its weakly magnetic properties and the relativelysmall (in terms of electric length) FOV. We then conjecturethat, in the absence of fast spatial variations in the mag-netic field, the problem can be simplified by a homogeneousmedium problem, and one can rewrite the field equation for h in Eq. (13) in terms of the free-space scalar Green function g ( r , r (cid:48) ) = e − jk | r − r (cid:48) | π | r − r (cid:48) | , with k = ω √ (cid:15) µ the wavenumber invacuum: G hj ( r , r (cid:48) ) = ∇ g ( r , r (cid:48) ) × I , (17) G hm ( r , r (cid:48) ) = 1 jωµ (cid:0) ∇∇ + k (cid:1) g ( r , r (cid:48) ) . (18) This is crucial for the practical applicability of the method:because the basis is computed in the absence of the biologicaltissue, it is universally applicable to all imaging problemssharing the same FOV. In other words, the basis is pre-computed offline for a few FOVs of interest, given only thedimensions of the FOV and the target resolution. In practice,this is achieved via a numerical discretization of Eq. (13).More specifically, in this work we obtain the discretized linearoperator in matrix form with the help of the open-sourcepackage MARIE [39], based on the methods presented in [40],[41].Finally, we observe that computing the SVD of K is notfeasible for practical problems, due to the extremely large sizeof K . As a matter of fact, K is only known via its sparsefactorization. A remedy to this is to resort to the so-calledrandomized matrix decompositions [42], [43], numerical tech-niques that have attracted growing interest recently thanks totheir effectiveness in computing low-rank approximations ofvery large matrices. Because the range of a linear operator canbe sampled with arbitrary precision if the images of indepen- ANUSCRIPT SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE 6 dent and random source distributions are known, by excitingdipoles located on the boundary with random amplitudes andphases it is possible to sample the left subspace of K withoutactually building it. Finally, because the detected MRI signalis a circular polarization of the magnetic field h , the subspaceis further restricted to span only circularly polarized fields.Figure 1 exemplifies the elements of a typical Maxwellbasis over square and circular supports. The randomized SVDbased approach guarantees that the basis vectors possess aspatial frequency content growing with the index of the basisvector: increasing the dimension of the basis increases thehigh-frequency components of the SM that can be capturedby U h . Consequently, the low-pass filtering properties behaveas a regularizer for the inverse problem. The representationproperties of the basis are demonstrated in Figure 2, where thecapability to expand known synthetic 2D SMs via the basisis analyzed. The convergence of the error of the expansion ofa known SM is shown as a function of the basis dimension q , proving that by increasing the dimension of the basis it ispossible to control the accuracy of the representation. −6 −5 −4 −3 −2 −1 ℓ - n o r m e rr o r Fig. 2: A set of 8 synthetic SMs, each of size 256x256pixels, is generated via an open-source Python package [44].The magnetic field basis is then used to represent the sameset of SMs, with dimension of the basis (number of basisvectors) ranging from 20 to 500. If U h is the matrix storingthe basis, and s i is the SM of the i th coil unrolled as acolumn vector, we denote the projection onto the basis as (cid:101) s i ≈ U h (cid:0) U Hh U h (cid:1) − U Hh s i . For each coil i and each size q ofthe basis, the projection error is computed as e q,i = (cid:107) (cid:101) c i − c i (cid:107) (cid:107) c i (cid:107) ,and the largest error among all coils max i ( e q,i ) is displayedas a function of q . Inset: the SM magnitude of one of the coilsis shown. E. A Compression Scheme for Maxwell Basis
The proposed method is valid for fully 3D problems, i.e. for3D acquisitions over volumetric FOVs, or can be restricted to2D problems. In the latter case, the range of K , and thusthe support of the basis, is restricted to a single slice. On theother hand, when the problem is fully 3D, storage requirementsfor the basis itself can be a limitation. Because the adopteddiscretization is a finite-element basis, each entry of one basisvector is proportional to the field intensity sampled at thecentroid of a voxel: each column of U h can be reshaped asa three-dimensional tensor representing a three-dimensionalfield distribution U i = r ( u i ) , with r : C N −→ C n × n × n , N = n n n being a reshape operator reordering entries of a column vectoronto a Cartesian grid.Here we follow the idea pioneered by Tucker in [45], whichintroduces a high-order singular value decomposition knownas Tucker decomposition. More specifically, Tucker decompo-sition is used to decompose a tensor T ∈ C n × n ... × n M to acore tensor G ∈ C n × n ... × n M multiplied by a unitary matrix U k ∈ C n k × n k along each mode k . In three dimensions: T = G × U × U × U . (19) A × k B denotes the k -mode product between a tensor A anda matrix B obtained as a convolution along the k th axis. Forinstance, the 1-mode product is defined as: C = A × B , C ijk = n (cid:88) p =1 A pjk B ip For an intuition of the decomposition, we find it useful topictorially visualize the 3D version as in Figure 3. If we acceptan approximation of T , the size of G (the multilinear ranks ofthe decomposition) can be much smaller than the size of T ,hence the compression. Similarly to the truncated SVD, we cantruncate the expansion in Eq. (19) with a reduced core tensor (cid:101) G ∈ C r × r × r and reduced unitary matrices (cid:101) U k ∈ C r k × n k ,with r k ≤ n k : T ≈ (cid:101)
G × (cid:101) U × (cid:101) U × (cid:101) U . (20)A key feature of the expansion in Eq. (20) is that it canbe obtained with controlled accuracy, i.e. given ε > it ispossible to find a Tucker expansion of T such that (cid:13)(cid:13)(cid:13) T − (cid:101)
G × (cid:101) U × (cid:101) U × (cid:101) U (cid:13)(cid:13)(cid:13) < ε (cid:107)T (cid:107) (21)For an overview of the algorithms to obtain a compressedTucker representation, the interested reader is referred to [45]–[47] and references therein.III. R ESULTS
In all the examples, densities and SMs obtained via Reg-ularized Nonlinear Inversion will be labelled as ”MPI-BL”(MPI BiLinear), while densities obtained via Regularized Den-sity Reconstruction via Constrained Optimization as ”MPI-L”(MPI Linear).
ANUSCRIPT SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE 7 n n n = r r r U n U n U n Fig. 3: Pictorial representation of the Tucker decomposition ofa 3D tensor T = G × U × U × U . For a better under-standing of the k − mode product × k , it is useful to realize thatit amounts to multiplying each mode- k fiber of the core tensor G by the matrix U k , i.e. it is a convolution along the k th axisof G . If the tensor T is rank deficient, we have that r k < n k and the decomposition results in a compression. When T isfull rank we can introduce an approximated decomposition of T up to an arbitrary accuracy ε , by truncating the ranks r k such that (cid:13)(cid:13)(cid:13) T − (cid:101)
G × (cid:101) U × (cid:101) U × (cid:101) U (cid:13)(cid:13)(cid:13) < ε (cid:107)T (cid:107) . A. 2D Cartesian sequences
Figures 4 and 5 depict the extracted SMs and density of anMPI reconstructed axial slice from a Cartesian acquisition ofa human brain obtained from the fastMRI database [21], [48].The data is a fully sampled Flash acquisition (TR/TE=250/3.4ms, FA=70 ◦ , matrix size: 320x320, slice thickness: 5mm)with a FOV of 220x220 mm , acquired at 3T using a 16-channel head coil. The data is retrospectively downsampled,according to different Cartesian undersampling patterns andACS regions. The stability of the recovered SMs for differentcombination of R and ACS regions proves the effectivenessof the physics-based regularization scheme. Aliasing artifactsare visible in the final reconstructed image for accelerationfactors R ≥ , and substantial stability of the image isobserved for ACS lines ≥ . Computation time on an IntelXeon CPU E5-2650 with NVIDIA Tesla K80 GPU is 346sand 36s for SMs and image reconstruction (R=2, ACS=16),respectively. Figure 6 investigates the performances of MPIfor simultaneous Cartesian accelerations along phase and slicedirections. The dataset is a fully sampled BRAVO acquisition(TR/TE=9.972/3.92 ms, FA=10 ◦ , matrix size: 192x192x170)with a FOV of 240x240x204 mm , acquired with a 1.5TGE using a 12-channel head coil. The frequency encoding isresolved and one single axial slice is reconstructed: differentCartesian downsampling schemes are retrospectively applied,with fixed number of ACS lines (16) and Maxwell basis with q = 50 . The reconstructed image is free from artifacts forcombined R ≤ . B. 2D synthetic radial sequences
Figure 7 shows the capability of MPI to address non-Cartesian acquisitions. Provided that the operator F of Eq. (1)is available, the described formulation is directly applicable.A synthetic 8 channels acquisition with golden angle radialsampling is generated from a Shepp-Logan model of size R = ACS=32 ACS=16 ACS=8 ACS=4 ACS=2 R = R = R = Fig. 4: A 3T fully sampled Cartesian acquisition of a humanhead is retrospectively downsampled, with different Cartesianundersampling factors R and ACS regions. MPI-BL SMsextraction is carried out with a fixed basis dimension q = 100 :the magnitude of the extracted SM of one coil is shown forincreasing acceleration (top to bottom) and decreasing ACSs(left to right). R = ACS=32 ACS=16 ACS=8 ACS=4 ACS=2 R = R = R = Fig. 5: The same dataset and undersampling strategy as inFig. 4 and the corresponding extracted SMs are used toreconstruct the image with MPI-L for increasing Cartesiandownsampling R (top to bottom) and decreasing number ofACSs (left to right).
ANUSCRIPT SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE 8 R p = Rs=1 Rs=2 Rs=3 Rs=4 R p = R p = R p = Fig. 6: A 1.5T fully sampled Cartesian acquisition of a humanhead is retrospectively downsampled, with different Cartesianundersampling factors along the phase encoding (Rp) and sliceencoding (Rs) dimensions, and a fixed number of 16 ACSs.The frequency encoding is resolved and one single axial sliceis reconstructed as a 2D problem. The reconstructed MPI-Ldensity is shown, with SMs extracted via MPI-BL and basisdimension q = 50 .256x256, in the presence of additive white Gaussian noiseindependent for each coil. C. Comparison with previous studies
Figure 8 explores variable-density Poisson-disc undersam-pled reconstructions of a knee, comparing MPI with ENLIVEand SAKE. All methods provide artifact-free reconstructionsup to acceleration R=3, with the denoising step performedby MPI-L providing a generally cleaner image. For higheraccelerations (R=5) SAKE misses signal from the center ofthe image, ENLIVE and MPI-BL both provide a rather noisyimage, while the MPI-L reconstructed image has significantbetter quality. Figure 9 shows Cartesian reconstructions withCAIPIRINHA patterns with different acceleration factors and24 ACS lines, with comparisons to ESPIRiT and ENLIVE.All images appear free from artifacts even at R=16.
D. 3D sequence
Figure 10 shows the results of a full 3D reconstruction ofthe same dataset of Figure 6, undersampled with a combinedacceleration factor 4 and 16 ACS lines, and reconstructedwith MPI-BL with a variable dimension of the Maxwellbasis. The results show the flexibility of the formulation inseamlessly addressing 3D k -spaces with the same formulation. R = SNR=∞ SNR=25dB SNR=20dB SNR=15dB R = R = R = Fig. 7: A 256x256 pixels Shepp-Logan phantom and a set of8 synthetic SMs with the same size are used to simulate asynthetic golden angle acquisition with readout length 256.Independent white Gaussian noise is then added to the sim-ulated k -space signal of each coil, yielding decreasing SNRvalues [ ∞ , dB , dB , dB ] (left to right). The solution ofMPI-L is shown for different numbers of acquired spokes N ∈ [200 , , , , corresponding to acceleration factors R ∈ [2 , , , (top to bottom). The basis dimension is set to q = 100 for all cases.At q = 50 aliasing artifacts are visible, as highlighted by theyellow arrow. When the basis is enlarged to capture theseartifacts, MPI-BL yields artifact-free images. Additionally,Tucker compression reduces the memory footprint of theMaxwell basis from 19.1GB to 31MB when q = 200 , inturn enabling accelerated computations on GPU (see alsoSupporting Information Table S1). Computation time on anIntel Xeon E5-2686 CPU with NVIDIA Tesla V100 GPU is73 minutes for 9 iterations of MPI-BL and q = 200 .IV. D ISCUSSION
A. Forward Model Extensions
There are cases where the bilinear form of the MR signalfails to capture accurately the underlying physics. More specif-ically, it is well documented that the image-domain methods,with the exception of ENLIVE [7], produce erroneous resultswhen the chosen FOV does not include entirely the objectunder study. As mentioned above, the SMs are essentially thecircularly polarized magnetic fields received by the coils, anddue to the nature of Maxwell equations their values dependstrongly on the EM properties of the entire object, not onlythe portion inside the FOV. Hence the estimation of the actualSMs for small FOVs is an ill-posed problem. Fortunately, MPIallows the extension of the original signal equation, much likeENLIVE, with the addition of extra bilinear terms, resultingin a fairly accurate approximation of the governing physics,
ANUSCRIPT SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE 9 R = . SAKE ENLIVE MPI-BL MPI-L R = . R = . Fig. 8: Variable-density Poisson-disc undersampled data ofa human knee with varying undersampling factors recon-structed with SAKE [9], ENLIVE [7], MPI-BL, and MPI-L. All methods generate images free from artifacts up toacceleration R=3, with MPI-L providing better SNRs. Forhighly accelerated sequences (R=5), ENLIVE and MPI-BLprovide similar images with high noise, while SAKE missessignal in the center of the image. MPI-L can provide a cleanerimage even at R=5. R = ESPIRiT ENLIVE MPI-BL MPI-L R = R = Fig. 9: Comparison of MPI (MPI-BL and MPI-L) with ES-PIRiT [19] reconstruction using 2 maps and ENLIVE [7]reconstruction using 2 maps of a single slice of a humanhead, undersampled with Cartesian CAIPIRINHA patternswith differing undersampling factors R and fixed size of theACS region (24). The MPI-BL joint density (third column)and SMs reconstruction is obtained with a basis of dimension q = 200 ; the reconstructed SMs are subsequently used as aninput for the MPI-L density reconstruction (last column). q=50 q=100 q=200 q=300 Fig. 10: Full 3D reconstruction with MPI-BL of the samedataset of Figure 6, with cartesian undersampling factor 2along the phase encoding and 2 along the slice encodingdimensions, corresponding to a combined acceleration 4, and16 ACSs, for different dimensions q of the basis along thecolumns. The full 3D density is reconstructed at once usingMPI-BL with Tucker compression of the Maxwell basis: slicesof the solution along the axial (top row), coronal (middle row)and sagittal (bottom row) are shown.though in this case the estimated SMs do not correspondanymore to the true magnetic field distributions and shouldbe considered as merely dummy variables. Nevertheless, theimage reconstruction is devoid of artifacts, as evinced by theSupporting Information Figure S2, where MPI with 2 sets ofmaps is applied on a dataset from Ref. [7]. B. Maxwell-Constrained Deep CNNs
In recent years we have witnessed some dramatic devel-opments in the field of machine and deep learning, whereDeep Convolutional Neural Networks (DCNNs) have shownsuperior performance over more traditional methods in variousimage reconstruction tasks, such as denoising [49], [50],demosaicking [51], super-resolution [52], etc. Consequently,this has also lead to an increased interest in the development ofdeep learning methods that could efficiently tackle the problemof MRI reconstruction [15], [53]. While in this work we haveintentionally focused on an optimization-based reconstructionapproach, we are convinced that a very promising futureresearch direction, which could lead to further improvementsin the reconstruction quality and offer additional robustness,is the design of physics-constrained deep reconstruction net-works. The main idea here is that by constraining the solutionspace of a neural network, we can gain more control onthe reconstruction outcome and reduce the risk of introduc-ing erroneous reconstruction artefacts, which are completelyundesirable in medical applications. In this direction, andfollowing the discussion on the construction of the Maxwellbasis, one possible way to enforce such kind of physics-based
ANUSCRIPT SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE 10 constraints is to combine the implicit Maxwell regularizationapproach with a variational-inspired deep network such asthose introduced in [51], [54]. This way, we can learn moremeaningful and accurate representations for the SMs, whichin turn can lead to better and more robust reconstructionresults. At the same time, it is expected that by providingmore information to the network about the space of solutions,we can avoid its overfitting during training and further requireless training data. V. C
ONCLUSIONS
In this work, we described a general framework for the jointreconstruction of PI data. The proposed framework introducesan expressive, physics-based regularizer for the estimation ofthe SMs and a constrained optimization scheme for the sub-sequent parameter-free density reconstruction, for improvedimage quality. In addition, the use of a Maxwell basis for theexpansion of the SMs reduces dramatically the overall numberof the unknowns in the inverse problem and acceleratesthe convergence of the iterative reconstruction. Finally, weutilized some relatively modern tensor decomposition methodsin order to reduce the memory footprint of the Maxwell basis,which can become prohibitively large for high-resolution 3Dscans. We expect this framework to allow MRI scientists andpractitioners to obtain images of higher quality from datasetswith even more aggressive acceleration, while its extensions incombination with deep learning-based reconstructions to offera paradigm shift in next-generation data-driven PI approaches.A
CKNOWLEDGEMENTS
We thank Daniel Sodickson, Riccardo Lattanzi, and ThomasWitzel for useful discussions.R
EFERENCES[1] D. K. Sodickson and W. J. Manning, “Simultaneous acquisition of spatialharmonics (SMASH): Fast imaging with radiofrequency coil arrays,”
Magn. Reson. Med , vol. 38, no. 4, pp. 591–603, 1997.[2] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, and P. Boesiger,“SENSE: Sensitivity encoding for fast MRI,”
Magn. Reson. Med , vol. 42,no. 5, pp. 952–962, 1999.[3] M. A. Griswold, P. M. Jakob, R. M. Heidemann, M. Nittka, V. Jellus,J. Wang, B. Kiefer, and A. Haase, “Generalized autocalibrating partiallyparallel acquisitions (GRAPPA),”
Magn. Reson. Med , vol. 47, no. 6, pp.1202–1210, 2002.[4] D. J. Larkman and R. G. Nunes, “Parallel magnetic resonance imaging,”
Phys. Med. Biol. , vol. 52, no. 7, pp. 15–55, Mar 2007.[5] L. Ying and J. Sheng, “Joint image reconstruction and sensitivityestimation in SENSE (JSENSE),”
Magn. Reson. Med , vol. 57, no. 6,pp. 1196–1202, 2007.[6] M. Uecker, T. Hohage, K. T. Block, and J. Frahm, “Image reconstructionby regularized nonlinear inversion–Joint estimation of coil sensitivitiesand image content,”
Magn. Reson. Med , vol. 60, no. 3, pp. 674–682,2008.[7] H. C. M. Holme, S. Rosenzweig, F. Ong, R. N. Wilke, M. Lustig, andM. Uecker, “ENLIVE: An efficient nonlinear method for calibrationlessand robust parallel imaging,”
Sci Rep , vol. 9, p. 3034, 2019.[8] J. D. Trzasko and A. Manduca, “A Calibrationless parallel MRI usingCLEAR,”
In Conf. Rec. Asilomar Conf. Signals Syst. Comput. , no. 45,pp. 75–79, 2011.[9] P. J. Shin, P. E. Z. Larson, M. A. Ohliger, M. Elad, J. M. Pauly, D. B. Vi-gneron, and M. Lustig, “Calibrationless parallel imaging reconstructionbased on structured low-rank matrix completion,”
Magn. Reson. Med ,vol. 72, no. 4, pp. 959–970, 2014. [10] J. P. Haldar, “Low-rank modeling of local k -space neighborhoods(LORAKS) for constrained MRI,” IEEE Trans. Med. Imag. , vol. 33,no. 3, pp. 668–681, Mar. 2014.[11] J. P. Haldar and J. Zhuo, “P-LORAKS: Low-rank modeling of local k -space neighborhoods with parallel imaging data,” Magn. Reson. Med ,vol. 75, no. 4, pp. 1499–1514, 2016.[12] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application ofcompressed sensing for rapid MR imaging,”
Magn. Reson. Med , vol. 58,no. 6, pp. 1182–1195, 2007.[13] F. Huang, Y. Chen, W. Yin, W. Lin, X. Ye, W. Guo, and A. Reykowski,“A rapid and robust numerical algorithm for sensitivity encoding withsparsity constraints: Self-feeding sparse SENSE,”
Magn. Reson. Med ,vol. 64, no. 4, pp. 1078–1088, 2010.[14] F. Knoll, C. Clason, K. Bredies, M. Uecker, and R. Stollberger, “Par-allel imaging with nonlinear reconstruction using variational penalties,”
Magn. Reson. Med , vol. 67, no. 1, pp. 34–41, 2012.[15] F. Knoll, K. Hammernik, C. Zhang, S. Moeller, T. Pock, D. K. Sodick-son, and M. Akcakaya, “Deep-learning methods for parallel magneticresonance imaging reconstruction: A survey of the current approaches,trends, and issues,”
IEEE Signal Processing Magazine , vol. 37, no. 1,pp. 128–140, 2020.[16] R. L. Morrison, M. Jacob, and M. N. Do, “Multichannel estimationof coil sensitivities in parallel MRI,”
In Proceedings of the 4th IEEEInternational Symposium on Biomedical Imaging: From Nano to Macro ,pp. 117–120, April 2007.[17] J. Jin, F. Liu, E. Weber, Y. Li, and S. Crozier, “An electromagneticreverse method of coil sensitivity mapping for parallel MRI – Theoreticalframework,”
Journal of Magnetic Resonance , vol. 207, no. 1, pp. 59 –68, 2010.[18] M. J. Allison, S. Ramani, and J. A. Fessler, “Regularized MR coil sensi-tivity estimation using augmented Lagrangian methods,”
In Proceedingsof the 9th IEEE International Symposium on Biomedical Imaging,Barcelona, Spain , pp. 394–397, May 2012.[19] M. Uecker, P. Lai, M. J. Murphy, P. Virtue, M. Elad, J. M. Pauly,S. S. Vasanawala, and M. Lustig, “ESPIRiT–an eigenvalue approach toautocalibrating parallel MRI: Where SENSE meets GRAPPA,”
Magn.Reson. Med , vol. 71, no. 3, pp. 990–1001, 2014.[20] Y.-J. Ma, W. Liu, X. Tang, and J.-H. Gao, “Improved SENSE imagingusing accurate coil sensitivity maps generated by a global magnitude-phase fitting method,”
Magn. Reson. Med , vol. 74, no. 1, pp. 217–224,2015.[21] J. Zbontar, F. Knoll, A. Sriram, and et al., “fastMRI: An open dataset andbenchmarks for accelerated MRI,” arXiv:1811.08839 preprint. , 2018.[22] P. C. Hansen, J. G. Nagy, and D. P. O’Leary,
Deblurring Images:Matrices, Spectra, and Filtering . SIAM, 2006.[23] A. B. Bakushinsky and M. Y. Kokurin,
Iterative methods for approximatesolution of inverse problems . Springer, 2005, vol. 577.[24] A. Levin, Y. Weiss, F. Durand, and W. T. Freeman, “Understanding blinddeconvolution algorithms,”
IEEE transactions on pattern analysis andmachine intelligence ∼ jrs/jrspapers.html[26] L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noiseremoval algorithms,” Physica D , vol. 60, pp. 259–268, 1992.[27] S. Lefkimmiatis, A. Roussos, P. Maragos, and M. Unser, “Structuretensor total variation,”
SIAM Journal on Imaging Sciences , vol. 8, no. 2,pp. 1090–1122, 2015.[28] S. Lefkimmiatis and S. Osher, “Nonlocal structure tensor functionals forimage regularization,”
IEEE Transactions on Computational Imaging ,vol. 1, no. 1, pp. 16–29, March 2015.[29] S. Lefkimmiatis, J. Ward, and M. Unser, “Hessian Schatten-normregularization for linear inverse problems,”
IEEE Trans. Image Process. ,vol. 22, no. 5, pp. 1873–1888, 2013.[30] E. Esser, “Applications of Lagrangian-based alternating direction meth-ods and connections to split Bregman,”
CAM report , vol. 9, 2009.[31] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein,
DistributedOptimization and Statistical Learning via the Alternating DirectionMethod of Multipliers . Now Publishers, 2011.[32] S. Lefkimmiatis and M. Unser, “Poisson image reconstruction withHessian Schatten-norm regularization,”
IEEE Trans. Image Process. ,vol. 22, pp. 4314–4327, 2013.[33] P. L. Combettes and V. R. Wajs, “Signal recovery by proximal forward-backward splitting,”
Multiscale Model. Simul. , vol. 4, no. 4, pp. 1168–1200, 2005.
ANUSCRIPT SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE 11 [34] B. He, H. Yang, and S. Wang, “Alternating direction method withself-adaptive penalty parameters for monotone variational inequalities,”
Journal of Optimization Theory and applications , vol. 106, no. 2, pp.337–356, 2000.[35] P. Blomgren and T. F. Chan, “Color tv: total variation methods forrestoration of vector-valued images,”
IEEE transactions on image pro-cessing , vol. 7, no. 3, pp. 304–309, 1998.[36] A. E. H. Love, “I. the integration of the equations of propagation ofelectric waves,”
Philosophical Transactions of the Royal Society ofLondon. Series A, Containing Papers of a Mathematical or PhysicalCharacter , vol. 197, no. 287–299, pp. 1–45, 1901.[37] A. Ishimaru,
Wave propagation and scattering in random media . Aca-demic Press, 1978.[38] M. V. Vaidya, C. M. Collins, D. K. Sodickson, R. Brown, G. C. Wiggins,and R. Lattanzi, “Dependence of b1+ and b1- field patterns of surfacecoils on the electrical properties of the sample and the mr operatingfrequency. concepts in magnetic resonance.”
Concepts Magn Reson PartB Magn Reson Eng , vol. 46, no. 1, pp. 25–40, 2016.[39] A. G. Polimeridis and J. F. Villena, “MARIE: MAgnetic ResonanceIntegral Equation suite.” [Online]. Available: https://github.com/thanospol/MARIE[40] A. Polimeridis, J. Villena, L. Daniel, and J. White, “Stable FFT-JVIEsolvers for fast analysis of highly inhomogeneous dielectric objects,”
J.Comput. Phys. , vol. 269, pp. 280 – 296, 2014.[41] J. F. Villena, A. G. Polimeridis, Y. Eryaman, E. Adalsteinsson, L. L.Wald, J. K. White, and L. Daniel, “Fast electromagnetic analysis ofMRI transmit RF coils based on accelerated integral equation methods,”
IEEE Trans. Biomed. Eng. , vol. 63, no. 11, pp. 2250–2261, Nov. 2016.[42] E. Liberty, F. Woolfe, P.-G. Martinsson, V. Rokhlin, and M. Tygert,“Randomized algorithms for the low-rank approximation of matrices,”
Proceedings of the National Academy of Sciences , vol. 104, no. 51, pp.20 167–20 172, 2007.[43] N. Halko, P. G. Martinsson, and J. A. Tropp, “Finding structurewith randomness: probabilistic algorithms for constructing approximatematrix decompositions,”
SIAM Rev. , vol. 53, no. 2, pp. 217–288, 2011.[44] Sigpy. [Online]. Available: https://sigpy.readthedocs.io/en/latest/mri.html[45] L. R. Tucker, “Some mathematical notes on three-mode factor analysis,”
Psychometrika , vol. 31, pp. 279–311, 1966.[46] I. I. Giannakopoulos, M. S. Litsarev, and A. G. Polimeridis, “Memoryfootprint reduction for the FFT-based volume integral equation methodvia tensor decompositions,”
IEEE Trans. Antennas Propag. , vol. 67,no. 12, pp. 7476–7486, 2019.[47] S. Rabanser, O. Shchur, and S. Gnnemann, “Introduction to tensordecompositions and their applications in machine learning,” 2017.[48] fastmri. [Online]. Available: https://fastmri.med.nyu.edu/[49] U. Schmidt and S. Roth, “Shrinkage fields for effective image restora-tion,”
In Proc. IEEE Int. Conf. Computer Vision and Pattern Recognition(CVPR) , pp. 2774–2781, 2014.[50] S. Lefkimmiatis, “Universal denoising networks : A novel cnn architec-ture for image denoising,”
In Proc. IEEE Int. Conf. Computer Visionand Pattern Recognition (CVPR) , June 2018.[51] F. Kokkinos and S. Lefkimmiatis, “Iterative joint image demosaickingand denoising using a residual denoising network,”
IEEE Transactionson Image Processing , vol. 28, no. 8, pp. 4177–4188, Aug 2019.[52] M. Haris, G. Shakhnarovich, and N. Ukita, “Deep back-projectionnetworks for super-resolution,”
In Proc. IEEE Int. Conference ComputerVision and Pattern Recognition (CVPR) , pp. 1664–1673, 2018.[53] D. Lee, J. Yoo, and J. C. Ye, “Deep artifact learning for compressedsensing and parallel mri,” arXiv preprint arXiv:1703.01120 , 2017.[54] F. Kokkinos and S. Lefkimmiatis, “Iterative residual cnns for burstphotography applications,”
In Proc. IEEE Int. Conf. Computer Visionand Pattern Recognition (CVPR) , June 2019.
ANUSCRIPT SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE 12 S UPPORTING MATERIAL M P I - B L ACS=6 ACS=12 ACS=24 M P I - L Fig. 11: MPI reconstruction of the same dataset of Figure9, undersampled with Cartesian CAIPIRINHA pattern withundersampling factor R=9 and and variable number of ACSlines. The results of MPI-BL (with basis dimension q = 200 )and with MPI-L demonstrate the flexibility of MPI, whichdoes not explicitly depend on the sampling pattern, allowingto reduce the size of the low-frequency portion of k -spacewithout abruptly breaking up. Fig. 12: MPI reconstruction of the dataset from [19], corre-sponding to a truncated FOV acquisition of a factor 2 under-sampled 2D spin-echo dataset acquired at 1.5T. 19 iterationsof the MPI-BL solver with basis dimension q = 50 are usedto generate solutions with 1 and 2 maps: the solution allowing2 maps (right) is free from artifacts, which are clearly visiblein the center of the single map solution (left). ANUSCRIPT SUBMITTED TO MAGNETIC RESONANCE IN MEDICINE 13
Dense Tucker q RAM CPU MVP GPU MVP RAM CPU MVP GPU MVP[MB] [sec] [sec](*) [MB] [sec] [sec](*)75 7,172 0.18 0.1 11 18.9 0.6200 19,125 0.35 – 31 50.4 1.1500 47,812 0.8 – 129 127.1 2.9
TABLE I: Time and memory requirements for the basis stored as a dense matrix and in compressed form with accuracy (cid:15) t = 10 − . The FOV has size 192x192x170 voxels, with corresponding matrix size N × q, N = 6266880= 6266880