We propose a motion and contrast enhancement separation

model in dynamic magnetic resonance imaging (MRI). Furthermore, the

reconstruction is done from partial measurements to achieve faster dynamic

MR imaging. The algorithm minimizes a linear combination of

three terms, a data fitting functional and two regularization functionals

corresponding to the nuclear and ℓ1 norm. The proposed method

is tested on simulated and real dynamic datasets. This paper suggests

an image reconstruction model that directly induces clinically-relevant

informations from partial measurements

Arridge, SR

Atkinson, D

Tremoulheac, B

UCL Discovery

Motion and Contrast Enhancement SeparationModel Reconstruction from PartialMeasurements in Dynamic MRIBenjamin Tre´moulhe´ac, David Atkinson, and Simon R. ArridgeCentre for Medical Image Computing, University College London, U.K.Abstract. We propose a motion and contrast enhancement separationmodel in dynamic magnetic resonance imaging (MRI). Furthermore, thereconstruction is done from partial measurements to achieve faster dy-namic MR imaging. The algorithm minimizes a linear combination ofthree terms, a data fitting functional and two regularization function-als corresponding to the nuclear and ℓ1 norm. The proposed methodis tested on simulated and real dynamic datasets. This paper suggestsan image reconstruction model that directly induces clinically-relevantinformations from partial measurements.Keywords: low-rank plus sparse model, dynamic MRI, motion, under-sampling, split Bregman-based method1 IntroductionDynamic magnetic resonance imaging (MRI) consists of acquiring multiple im-ages (frame by frame) of a dynamic object in time. Since this is a sequential andslow acquisition process, motion artefacts are often present due to unintentionalmovements of the patient. Faster dynamic MR imaging, by reducing the numberof measurements (under-sampling), can minimize the chance of motion affectinga frame, giving a better temporal resolution for dynamic processes, or tradingthe time saving for higher spatial resolution.A few years ago, a more efficient signal acquisition than the traditional Shan-non sampling theory was introduced for sparse or compressible signals. This the-ory, called Compressed Sensing (CS), states that one can recover images fromonly a few random measurements. This is possible providing that (1) the imageis sparse in some domain or in its own, (2) sensing and sparsifying matricessatisfies the restricted isometry property and are incoherent, and (3) non-linearreconstruction techniques is used to recover the signal. Since these assumptionscan generally be met in medical imaging, CS has been applied in MRI as a rapidimaging technique [1].Recently, the recovery of a matrix from incomplete information has attractedinterest. This research topic is known as low-rank matrix completion [2] and hasstrong parallels with CS framework. Indeed, one can see the notion of matrices2 Benjamin Tre´moulhe´ac, David Atkinson, and Simon R. Arridgewith a low-rank as an equivalent of signals with only a few significant coeffi-cients. Similarly, analogous to ℓ1 norm minimization promoting sparse signalsis the notion of nuclear norm minimization promoting low-rank matrices. Com-bining low-rank and sparsity constraints for matrix problems has recently beeninvestigated in mathematics and computational science because it has a largespectrum of potential applications. The idea is to recover low-rank and sparsecomponents from a matrix which is supposed to have an underlying low-rankstructure [3–5].In medical imaging, these type of priors have been recently studied. In par-ticular, it has been used by Majumdar and Ward [6] for accelerating MRI, byLingala et al. for accelerating dynamic MRI [7] and by Gao et al. [8] to reduceradiation dose in computed tomography (CT). While these works mostly focuson achieving better image quality from under-sampled data than state of theart methods, we take a different approach in this paper. The novelty is thatregularization is viewed as providing a useful decomposition of dynamic MRimages that can be used technically or is shown to be clinically relevant. Thealgorithm used, a split Bregman based method proposed by Gao et al. in CT[8], has strong links with alternating direction method for solving the low-rankand sparse decomposition problem [4, 5].The rest of the article is organized as the following. Section 2 provides abrief review of the low-rank and sparse decomposition problem. In section 3,we discuss under-sampling strategies for the proposed method. Reconstructionmodel and algorithm are explained in section 4. Numerical simulations are shownin section 5 and section 6 concludes this paper.2 Low-rank and Sparse DecompositionIn this section, we introduce the decomposition from a mathematical and com-putational perspective. We assume we have a matrix X ∈ Rn×m which is thesuperposition of a low-rank and a sparse component. The goal is to decom-pose the given matrix into its low-rank L and sparse S components such thatX = L+ S. This can be cast naturally as the following optimization problem,minL,Srank(L) + λ‖S‖0 such that X = L+ S , (1)but reveals to be intractable as it encloses rank and ℓ0 minimization problems,which both are known to be NP-hard and consequently difficult to solve. From astatistical point of view, this method can be seen as a robust variant of principalcomponent analysis (PCA), where one seeks a low-rank matrix from grosslycorrupted measurements.To make the problem easier to solve, one can reformulate it as a convex one[2–4],minL,S‖L‖∗ + λ‖S‖1 such that X = L+ S . (2)Here ‖L‖∗ and ‖S‖1 denote respectively the nuclear norm of L (the sum of its sin-gular values), ‖L‖∗ =∑i σi(L) with σi(L) the ith singular value of L, and the ℓ1Motion and contrast enhancement separation model 3norm of S defined as the sum of the absolute values of all entries, ‖S‖1 =∑i |si|with si the ith entries of S viewed as a vector. λ can be viewed as a trade-offparameter between low-rank and sparse components. Theoretical analysis in [5]leads to the trade-off parameter value λ = max(n,m)−1/2 for recovering anylow-rank, incoherent matrix with high probability. The convex relaxation ap-proach (2) is usually solved using alternating direction method, an improvementof the augmented Lagrangian method, which exploits the favourable structureappearing in both the objective function and the constraint [4, 5].3 Dynamic MR Data and Sampling SchemeIn this study, the low-rank and sparse decomposition model operates as a spatio-temporal one. The matrixX to decompose is formed as the concatenation of eachimage frames {xi}Ji=1 ∈ RN=n×n casted as column vectors of X , i.e.,X = [x1, x2, ..., xJ ] ∈ RN×J . (3)Hence, columns of X correspond to spatial pixels and rows correspond to tem-poral intensity variations of the pixels in the image domain. Since each framesxi exhibits significant correlations, the whole sequence lies in fact in a low-complexity model, an underlying low-rank structure, and thus motivates thisdecomposition. We highlight that the method relies on an appropriate numberof frames since the matrix X can be said low-rank only if rank(X)≪ min(N, J).To accelerate MRI, the idea is to reduce acquisition time by sampling onlypartially the k-space and reconstruct images from these incomplete data. Creat-ing better random k-space sampling trajectories for under-sampled MRI acquisi-tion has been subject to a lot of research in the CS MRI literature, as in practiceit is conditioned by hardware and physiological considerations. Here, the pro-posed method is more closely related to low-rank matrix completion than CS,since we do not employ any sparsifying transform. Under-sampling is possiblebecause of the spatio-temporal matrix formulation and the use of nuclear normand ℓ1 minimization in the model.Since we do not want to use non uniform Fourier transform in order to avoidlonger computational cost, we limit our trajectory choices to Cartesian sampling.Experiments showed that the reconstruction quality is sensitive to the samplingstrategy. In this paper, we choose to densely sample the low spatial frequenciesand randomly select phase encoding lines elsewhere in each time-frame. Thisresults in a sampling scheme random in time, except for the center of the k-space. We make this choice not only because it is easily done in practice, butalso because it is common practice in CS MRI and because we could achievemodest acceleration rate with accurate reconstruction. Here, acceleration ratesare computed by dividing the number of k-space points acquired over the totalnumber of points in the k-space. For example, a 2-fold acceleration implies thatonly 50% of the k-space is used.4 Benjamin Tre´moulhe´ac, David Atkinson, and Simon R. Arridge4 Reconstruction Model and AlgorithmIn this paper, we consider the following unconstrained optimization problem forthe reconstruction,minL,S12‖A(L+ S)− y‖22 + µ(‖L‖∗ + λ‖S‖1) (4)where A : RN×J → Cp is the Fourier operator providing a vector y of p < N×Jmeasurements. The first term in the objective function is a data fitting criterion,and the second term encloses the original decomposition problem.A slightly different problem has been proposed by Gao et al. in dynamic CT[8] and solved using a split Bregman based method. We rely on this efficientapproach to solve (4). The iterative algorithm, described by equations in (5),consists of three major steps : a quadratic problem solved using conjugate gra-dient, an ℓ1 problem solved by shrinkage formula and a nuclear norm problemsolved by singular value thresholding.(Ln+1, Sn+1) = argminL,S‖A(L+ S)− y + fn‖22 + α‖L− dnL + vnL‖22+ β‖S − dnS + vnS‖22dn+1S = argmindS12‖Sn+1 + vnS − dS‖22 + λµβ‖dS‖1vn+1S = vnS + Sn+1 − dn+1Sdn+1L = argmindL12‖Ln+1 + vnL − dL‖22 +µα‖dL‖∗vn+1L = vnL + Ln+1 − dn+1Lfn+1 = fn +A(Ln+1 + Sn+1)− y(5)These three major steps have closed form solutions, hence the iterative schemecan be rewritten as described in algorithm 1, where S and D are respectivelythe shrinkage and singular value thresholding operators, defined asSτ (X) = sign(X)max(|X | − τ, 0) , (6)Dτ (X) = USτ (Σ)V∗ . (7)Here X = UΣV ∗ represents any singular value decomposition. To ensure goodconvergence, stopping criterion is met when a maximum number of iteration(typically 500 or 1000 iterations) is achieved or when ‖A(L+S)−y‖/‖y‖ ≤ 10−5.α and β are set equal to µ.We emphasize the main goal of this work is not to compare image recon-struction quality with state of the art methods in dynamic MR reconstructionfrom under-sampled data. In fact, it is to provide the direct reconstruction ofthe low-rank and sparse components, yielding the exploitation of the proposeddecomposition. Accordingly, we suggest that the parameter λ should be tunedMotion and contrast enhancement separation model 5Algorithm 1 Low-rank and sparse decomposition from partial measurementsInput: Measurements y, parameters µ (regularization) and λ (decomposition)Output: L, SL0 = S0 = v0L = v0S = d0L = d0S ← 0 ; α = β = µ ; f0 ← 0 ; n← 0while stopping criterion is not met doSn+1 ← (ATA+ µI)−1(AT (y − fn) + β(dnS − vnS)− ATA(Ln))Ln+1 ← (ATA+ µI)−1(AT (y − fn) + α(dnL − vnL)− ATA(Sn))dn+1S ← S λµβ(Sn+1 + vnS)vn+1S ← vnS + Sn+1 − dn+1Sdn+1L ← D µα(Ln+1 + vnL)vn+1L ← vnL + Ln+1 − dn+1Lfn+1 ← fn + A(Ln+1 + Sn+1)− yn← n+ 1end whileby visual inspection, to obtain the best possible motion and contrast enhance-ment separation. More precisely, a correct scaling for this parameter should beof order of λ = max(N, J)−1/2 given by Cande`s et al. [5].5 Numerical SimulationsTo evaluate the method, we tested the proposed model on simulated and invivo datasets by targeting cardiac imaging and dynamic contrast enhanced(DCE) MRI applications. Numerical simulations were implemented on MAT-LAB R2012a environment and run on a 2.2GHz Intel 8-core processor laptop.Simulated dataset were created in image domain and real datasets were acquiredfully-sampled on MRI scanner. Both were then retrospectively under-sampledwith a corresponding acceleration rate. The decomposition is computed usingalgorithm 1 with µ found by empirical study and λ tuned by visual inspection.5.1 Dynamic phantomThe method is first studied on simulated data. A 2D plus time numerical phan-tom of size 128× 128× 40 is created to mimic three time-varying components indynamic MRI. In particular, the phantom simulates periodic respiratory motion,periodic local motion and abrupt local changes of intensity.To our knowledge there is no commonly used method that provides this jointseparation and reconstruction procedure. Hence, we test our approach againstthe zero filling Fourier transform and dynamic MRI acceleration method k-tFOCUSS [9] for the reconstruction stage, followed up by a principal componentanalysis decomposition for the separation stage. More precisely, the PCA stage isdone by taking only the first three principal components to model a componentwith a low-rank. The difference between reconstructed data and this low-rankcomponent is done to obtain an equivalent of the sparse component of the pro-posed method. Results in figure 1 show that the suggested approach is able to6 Benjamin Tre´moulhe´ac, David Atkinson, and Simon R. Arridgereconstruct with very good accuracy the original sequence, and provide a morerobust decomposition than the classical PCA for time-varying components suchas motion and intensity changes. We provide quantitative results in table 1.Zero filling Fourier transform and PCAk-t FOCUSS and PCAProposed methodPCA DifferenceTemporal evolutionExtracted frameOriginalOriginal dataL STemporal evolutionExtracted frameL+STemporal evolutionExtracted framek-t FOCUSSPCA DifferenceTemporal evolutionExtracted frameZero filling FTPCA DifferenceFig. 1. Dynamic phantom dataset. Original data and different reconstruction methodsfrom partial measurements with a 2-fold acceleration along with PCA decomposition.The sampling scheme used is the one described in section 3. For the proposed method,parameters are set to µ = 10 and λ = 2max(N, J)−1/2.Zero filling FT k-t FOCUSS Proposed methodReconstruction error 0.1072 0.0449 0.0020‖AX − y‖22 1.5010e-28 0.5869 3.3621e-05Nuclear norm ‖.‖∗ 291.2141 299.7161 289.1930Sparsity ‖.‖1 1.1675e+04 5.7848e+03 3.5804e+03Table 1. Quantitative results for phantom reconstruction. Reconstruction error is com-puted as relative reconstruction error, i.e., ‖X˜−X‖2/‖X‖2, where X˜ are reconstructeddata.5.2 In vivo dataCardiac imaging dataset consists of a cine MRI sequence of the heart with asize of 224×155×50. The original dataset and its reconstruction/decompositionMotion and contrast enhancement separation model 7from partial k-space is presented in figure 2. We tune λ to obtain a rank-1 com-ponent, acting as an average of the whole sequence, while completely isolatingmovement of the heart in the sparse frames. Presentation of these images fromthis perspective may be useful from a clinical point of view by isolating heartmovements that may help detect motion abnormalities.Full FOV L+S L SFully-sampled Under-sampled Close upFig. 2. Cardiac cine MRI dataset. Left figures are fully-sampled data with full field ofview, close up and temporal evolution along a dynamic profile (corresponding dashedline). Right figures are reconstructed images from partial measurements. Top figuresshow specific frame with the separation into low-rank and sparse components. Bottomfigures present temporal evolution along corresponding dynamic profiles. Accelerationrate is 2, parameters are set to µ = 10 and λ = 12max(N, J)−1/2. Relative reconstruc-tion error is 0.0792.In DCEMRI, acquisition of multiple MR images are taken before, during, andafter the administration of an MR contrast agent. As in dynamic MRI, artefactsdue to patient movements (such as breathing, heartbeat, etc.) can appear duringthe acquisition. A common approach is to register images to reduce importantmisalignments. Here, dataset is of size 256× 256× 10 with respiratory motion.The decomposition from partial k-space, shown in figure 3, is able to isolatecontrast enhancement in the heart through the sparse component. Since thecontrast uptake and washout is clearly separated in this decomposition, it hasbeen shown that it could help in the registration of the time frames of the low-rank component to correct for respiration-induced motion [10].6 ConclusionWe have proposed reconstruction of low-rank and sparse components of dynamicMR data from partial k-space measurements. While providing faster dynamicMR imaging, this method leads to a flexible separation of motion and contrastenhancement that can be exploited in different ways.More generally, this approach suggests a joint reconstruction model with aclinically-relevant separation process. This is an interesting direction for fur-8 Benjamin Tre´moulhe´ac, David Atkinson, and Simon R. ArridgeL+S L SFully-sampled Under-sampledFig. 3. DCE MRI dataset. Left figure is fully-sampled extracted frame from the se-quence. Right figures present extracted frames of under-sampled reconstructed se-quence along with low-rank and sparse decomposition. 2-fold acceleration, parametersare set to µ = 10 and λ = 2max(N, J)−1/2. Relative reconstruction error is 0.0587.ther research, where reconstruction models would extract informations relateddirectly to clinically-relevant measures. This might be possible through incorpo-ration of appropriate priors and by targeting specific applications in medicine.Acknowledgments. The authors thank Dr. Muhammad Usman and ValentinHamy for useful discussions and providing respectively the cardiac and DCEMRI datasets. This work was funded by EPSRC (grant EP/H046410/1).References1. Lustig, M., Donoho, D.L., Pauly, J.M., Sparse MRI: The Application of CompressedSensing for Rapid MR Imaging. Magnetic Resonance in Medicine 58, pp. 1182–1195(2007)2. Cande`s, E.J., Recht, B., Exact Matrix Completion via Convex Optimization. Foun-dations of Computational Mathematics 9, pp. 717–772 (2009)3. Chandrasekaran, V., Sanghavi, S., Parrilo, P.A., Willsky, A.S., Rank-sparsity Inco-herence for Matrix Decomposition. SIAM Journal on Optimization 21 (2011)4. Yuan X. M., Yang J.F., Sparse and Low-rank Matrix Decomposition via AlternatingDirection Methods, Pacific Journal of Optimization, Preprint (2009)5. Cande`s, E.J., Li, X., Ma, Y., Wright, J., Robust Principal Component Analysis?Journal of the ACM 58 (2011)6. Majumdar A., Ward R.K., Exploiting Rank Deficiency and Transform Domain Spar-sity for MR Image Reconstruction. Magnetic Resonance Imaging 30, pp. 9–18 (2012)7. Lingala, S.G., Hu, Y., DiBella, E., Jacob, M., Accelerated Dynamic MRI ExploitingSparsity and Low-rank Structure: k-t SLR, IEEE Transactions on Medical Imaging30, pp. 1042–1054 (2011)8. Gao H., Cai, J.-F., Shen Z., Zhao H., Robust Principal Component Analysis-basedFour-dimensional Computed Tomography. Physics in Medicine and Biology 56(2011)9. Jung H., Sung K., Nayak K.S., Kim E.Y., Ye J.C., k-t FOCUSS: A General Com-pressed Sensing Framework for High Resolution Dynamic MRI. Magnetic Resonancein Medicine 61, pp. 103–116 (2009)10. Hamy, V., Melbourne A., Tre´moulhe´ac B., Punwani S., Atkinson D., Registrationof DCE-MRI using Robust Data Decomposition. Proc. ISMRM, p. 749 (2012)

Motion and contrast enhancement separation model reconstruction from partial measurements in dynamic MRI

http://discovery.ucl.ac.uk/1371089/1/Atkinson_Final.pdf

Motion and contrast enhancement separation model reconstruction from partial measurements in dynamic MRI

Abstract

Similar works

Full text

Available Versions

UCL Discovery