Freddy Odille
email: freddy.odille@inserm.fr

Anne Menini
email: anne.menini@ge.com

J.-M Escanyé
email: jm.escanye@chu-nancy.fr

P.-Y Marie
email: py-marie@chu-nancy.fr

Jean-Marie Escanyé

Pierre-André Vuissoz
email: pa-vuissoz@chu-nancy.fr

Pierre-Yves Marie

Marine Beaumont
email: m.beaumont@chu-nancy.fr

Member, Jacques Felblinger
email: j-felblinger@chu-nancy.fr

Joint reconstruction of multiple images and motion in MRI: application to free-breathing myocardial T 2 quantification

Keywords: compressed sensing, deformable registration, multi-contrast, relaxometry, T 2 mapping

Exploiting redundancies between multiple images of an MRI examination can be formalized as the joint reconstruction of these images. The anatomy is preserved indeed so that specific constraints can be implemented (e.g. most of the features or spatial gradients should be in the same place in all these images) and only the contrast changes from one image to another need to be encoded. The application of this concept is particularly challenging in cardiovascular and body imaging due to the complex organ deformations, especially with the patient breathing. In this study a joint optimization framework is proposed for reconstructing multiple MR images together with a nonrigid motion model. The motion model takes into account both intra-image and inter-image motion and therefore can correct for most ghosting/blurring artifacts and misregistration between images. The framework was validated with freebreathing myocardial T 2 mapping experiments from nine heart transplant patients at 1.5 T. Results showed improved image quality and excellent image alignment with the multi-image reconstruction compared to the independent reconstruction of each image. Segment-wise myocardial T 2 values were in good agreement with the reference values obtained from multiple breath-holds (62.5 ± 11.1 ms against 62.2 ± 11.2 ms which was not significant with p=0.49).

I. INTRODUCTION

RI examinations consist of acquiring multiple images which often contain a large amount of redundant information. This is because the anatomy is preserved provided the patient does not move significantly between the different imaging sequences. Exploiting these redundancies can be formalized as the joint reconstruction of these images using dedicated regularization schemes. Constraints that can take advantage of the anatomical redundancy can be based, for instance, on the co-occurrence of the spatial gradients (i.e. the gradients should be located in the same place in all images) [START_REF] Haldar | Improved diffusion imaging through SNR-enhancing joint reconstruction[END_REF], [START_REF] Bilgic | Multicontrast reconstruction with Bayesian compressed sensing[END_REF] or on joint statistics [START_REF] Somayajula | PET Image Reconstruction Using Information Theoretic Anatomical Priors[END_REF] (i.e. joint entropy, mutual information etc…). The principle of joint image reconstruction has also been applied to multi-modal imaging such as the joint reconstruction of PET and MRI images [START_REF] Knoll | Simultaneous MR-PET Reconstruction using Multi Sensor Compressed Sensing and Joint Sparsity[END_REF].

The application of this concept is particularly challenging in cardiovascular and body MRI applications. All images should be perfectly aligned in order to benefit best from the redundancy constraints. Organs may be misaligned due to involuntary motion, due to physiological changes (respiratory, cardiac, peristaltic motion etc…) or due to geometric distortions caused by hardware imperfections (e.g. in echo planar imaging sequences). In addition to correcting misalignments from image to image, it is desirable to correct for ghosting/blurring artifacts caused by motion during the kspace sampling (i.e. between the acquisition of different phase encoding steps), which can affect each image individually. These two types of motion will be termed inter-image and intra-image motion in the remainder. In the case of a single image acquisition, only intra-image motion correction needs to be considered. This can be formalized as the joint optimization of the image and a nonrigid motion model [START_REF] Odille | Generalized reconstruction by inversion of coupled systems (GRICS) applied to free-breathing MRI[END_REF]. Such an approach has also been applied to CT reconstruction [START_REF] Hinkle | 4D CT image reconstruction with diffeomorphic motion model[END_REF].

Myocardial T 2 quantification is an example application that might benefit from a joint reconstruction framework. Myocardial T 2 is altered in a number of cardiopathies including heart transplant rejection [START_REF] Marie | Detection and prediction of acute heart transplant rejection with the myocardial T2 determination provided by a black-blood magnetic resonance imaging sequence[END_REF]- [START_REF] Bonnemains | Diagnostic and prognostic value of MRI T2 quantification in heart transplant patients[END_REF]. Patient motion is commonly addressed by using a single-shot acquisition technique consisting of a breath-held T 2 -prepared balanced steady-state free precession (T 2 -prep bSSFP) sequence with three different echo times (TE) [START_REF] Giri | T2 quantification for improved detection of myocardial edema[END_REF]- [START_REF] Van Heeswijk | Self-navigated isotropic three-dimensional cardiac T2 mapping[END_REF]. Such acquisitions are fast but have limited spatial resolution. Furthermore they are sensitive to artifacts because only three data points (three values of TE) are available for the two parameter exponential fit of the form 𝑎 exp (-𝑇𝐸/𝑏). The latter point can be addressed by covering a wider range of TE values. However the acquisition becomes significantly longer and has to be performed during free breathing. Such an approach, combined with nonrigid image registration, has been shown to improve the accuracy of T 2 measurements in [START_REF] Roujol | Impact of motion correction on reproducibility and spatial variability of quantitative myocardial T2 mapping[END_REF]. Finally in order to increase the spatial resolution, single-shot techniques have to be replaced by multi-shot ones. A black-blood fast spin echo sequence can be used in that case as in [START_REF] Marie | Detection and prediction of acute heart transplant rejection with the myocardial T2 determination provided by a black-blood magnetic resonance imaging sequence[END_REF], [START_REF] Bonnemains | Diagnostic and prognostic value of MRI T2 quantification in heart transplant patients[END_REF]. This sequence provides good contrast between blood and myocardium due to the blood signal cancellation.

In this paper a framework is proposed for the joint reconstruction of multiple images and motion, which is designed for a wide range of thoracic and abdominal MRI applications. An application to free-breathing T 2 quantification is presented with the data acquisition consisting of a series of 10 fast spin echo sequences with varying effective times of echo (TE eff ), resulting in 10 images with varying contrast. The objectives of the study are (i) to evaluate the benefit of the proposed framework compared to the independent reconstruction of each image in terms of intra-image motion correction (ghosting/blurring artifact reduction) and interimage motion correction (alignment) and (ii) to validate the method clinically by comparing the quantitative myocardial T 2 values obtained from free-breathing to the conventional breath-held protocol in a population of heart-transplant patients.

II. THEORY

In this section we first review the joint reconstruction of one image and motion and then describe its extension to the joint reconstruction of N images and motion.

A. Joint reconstruction of an image and a motion model

We first consider the reconstruction of a single image which is corrupted by motion occurring during the acquisition. The reconstruction of a motion-corrected image can be formulated as a joint optimization problem where both the image and motion are unknown [START_REF] Odille | Generalized reconstruction by inversion of coupled systems (GRICS) applied to free-breathing MRI[END_REF]:

Here m denotes the acquired k-space data (a vector of 𝑁 𝑚 elements), 𝜌 the motion-corrected image (a vector of 𝑁 𝜌 elements) in a reference motion state and 𝑢 the displacement fields at each k-space sample time (a vector of 𝑁 𝑚 𝑁 𝜌 𝑁 𝑑 elements with 𝑁 𝑑 the number of dimensions of the image, i.e. 𝑁 𝑑 = 2 or 𝑁 𝑑 = 3), E the forward model describing the MRI acquisition process (a matrix of size 𝑁 𝑚 × 𝑁 𝜌 ); 𝑅(𝜌) and 𝑆(𝑢) are regularizers, in this study chosen to be 𝑅(𝜌) = ‖𝜌‖ 2 in order to get the minimal norm solution image and 𝑆(𝑢) = ‖𝛻𝑢‖ 2 in order to impose spatial smoothness on the displacement fields. In this formulation the joint reconstruction is highly underdetermined since the displacement fields need to be described at each pixel in the image and at the acquisition times corresponding to each kspace measurement. It is therefore necessary to parameterize the motion vector 𝑢.

A .

(2)

Assuming the motion sensor signals 𝑠 𝑘 (t) are known for each k-space measurement, Eq. ( 1) is rewritten as:

min (𝜌,𝛼) ‖𝐸(𝛼, 𝑠)𝜌 -𝑚‖ 2 + 𝜆𝑅(𝜌) + 𝜇 𝑆(𝛼), (3)
where 𝛼 = [𝛼 1 … 𝛼 𝑁 𝑠 ] 𝑇 is a vector of 𝑁 𝑠 𝑁 𝜌 𝑁 𝑑 elements. Formulation (3) can be made over-determined by repeating the acquisition a sufficient number a times, so that each k-space measurements is acquired in different motion states (i.e. with different values of 𝑠 𝑘 (t)), rendering the optimization with respect to 𝛼 tractable. In previous work, it was found that 2 to 3 repetitions provide satisfactory results [START_REF] Odille | Generalized reconstruction by inversion of coupled systems (GRICS) applied to free-breathing MRI[END_REF], [START_REF] Odille | Generalized MRI reconstruction including elastic physiological motion and coil sensitivity encoding[END_REF]. In the present work only one repetition was used for each image. Therefore the application of the conventional GRICS is expected to yield suboptimal results. The optimization in (3) is solved using Algorithm 1. This solver involves an alternating optimization (alternating between 𝜌 and 𝛼 optimization); a Gauss-Newton scheme is introduced to linearize the cost function in (3) around the current estimate of 𝛼 (i.e. each Gauss-Newton iteration searches for an optimal refinement 𝛿𝛼); a multi-resolution strategy is used to deal with the initialization of 𝛼 in case of large displacements; the conjugate gradient is used to solve the two linear least squares problems (steps 2 and 4). Algorithm 1 is also termed GRICS (generalized reconstruction by inversion of coupled systems). Detailed implementation of the operators 𝐸 and 𝐽 involved in Algorithm 1 are given in Appendix 1 and Appendix 2 respectively.

B. Joint reconstruction of multiple images and a motion model

We now consider the joint reconstruction of 𝑁 𝑖 images 𝜌 1 … 𝜌 𝑁 𝑖 and a motion model. The framework described in the previous section can be extended by grouping the min (𝜌,𝑢) ‖𝐸(𝑢)𝜌 -𝑚‖ 2 + 𝜆𝑅(𝜌) + 𝜇 𝑆(𝑢).

(1) where all variables have been modified with a subindex 𝑖 referring to the acquisition of the 𝑖 th image. Formulation (4) is equivalent to reconstructing each image and its motion model independently. It needs to be modified in order to take advantage of the redundant information.

Redundancies can be exploited in two ways [START_REF] Menini | Motion correction of multicontrast images applied to T1 and T2 quantification in cardiac MRI[END_REF]: (i) the motion model coefficients 𝛼 𝑖 can be assumed to be constant over the examination and shared by all the images, i.e. 𝛼 𝑖 = 𝛼 for 𝑖 ∈ {1 … 𝑁 𝑖 }; (ii) the anatomical redundancies can be exploited by imposing a joint regularization term 𝑅 that is minimized when images share most information. The sharing of the motion model coefficients means that the relation between the displacement fields and the sensor values, as defined by Eq. ( 2), is assumed to be true during the acquisition of all the images to be reconstructed jointly. This means we are searching for the motion model coefficients that best fit the whole dataset instead of separate coefficients that best fit the acquisition of each individual image. We therefore rewrite the reconstruction as:

For the joint regularization term 𝑅 it was proposed to use a function that is minimized when spatial gradients are perfectly matched in the image series. In this study with 2D images we use a constraint similar to that in Refs [START_REF] Haldar | Improved diffusion imaging through SNR-enhancing joint reconstruction[END_REF], [START_REF] Menini | Motion correction of multicontrast images applied to T1 and T2 quantification in cardiac MRI[END_REF]:

𝑅(𝜌 1 … 𝜌 𝑁 𝑖 ) = ∑‖W∇ρ 𝑖 ‖ 2 𝑁 𝑖 𝑖=1 = ∑ ‖W 𝑥 𝜕ρ 𝑖 𝜕𝑥 ‖ 2 + ‖W 𝑦 𝜕ρ 𝑖 𝜕𝑦 ‖ 2 𝑁 𝑖 𝑖=1 , 𝑤𝑖𝑡ℎ W 𝑞 = 𝑑𝑖𝑎𝑔 (max (∑ 𝛽 𝑖 2 | 𝜕ρ 𝑖 (𝑗-1) 𝜕𝑞 | 2 𝑁 𝑖 𝑖=1 , 𝜀 2 )) - 1 2 𝑓𝑜𝑟 𝑞 ∈ {𝑥, 𝑦} . (6)
The diagonal weighting matrix 𝑊 is computed using image estimates ρ 𝑖 (𝑗-1) from the previous iteration. The 𝛽 𝑖 are scaling factors and are chosen so that all images have approximately the same contribution to the regularization term; here we set 𝛽 𝑖 to be the mean absolute value of the image 𝜌 𝑖 . The threshold 𝜀 is introduced to avoid division by 0 and here it was set to be: 𝜀 = 0.01 × 𝑚𝑎𝑥(|∇ρ|). Formulation (5) has two major advantages over [START_REF] Knoll | Simultaneous MR-PET Reconstruction using Multi Sensor Compressed Sensing and Joint Sparsity[END_REF]. First it inherently deals with the ill-conditioning of the motion model optimization step due to the sharing of the motion model coefficients. In the previous section this was dealt with by repeating the acquisition of each k-space sample at least twice (in different motion states). This is no longer necessary since the optimization with respect to 𝛼 in Eq. ( 5) inherently contains data acquired in different motion states and will be sufficiently well conditioned as soon as 𝑁 𝑖 ≥ 2. The second advantage is that the ill-conditioning of the image reconstruction step is counterbalanced by the use of the shared anatomical information. This ill-conditioning comes from the nonrigid motion operators contained in the 𝐸 𝑖 matrices. This is because the columns of the 𝐸 𝑖 matrices do not form a set of linearly independent equations any longer (as would happen if some parts of the k-space were left unacquired [START_REF] Atkinson | Reconstruction after rotational motion[END_REF]). Solving for ( 5) is a straight-forward adaptation of Algorithm 1 and is a special case of Algorithm 2 described hereafter.

C. Joint reconstruction of multiple images, intra-and interimage motion

The motion model 𝛼 defined in Eq. 2 provides plausible estimates of both intra-and inter-image motion. At this stage intra-image motion estimates can hardly be improved without any additional information (e.g. with more sensors, navigators or rapid images from another imaging modality). However inter-image motion can still be refined by simply applying an image registration technique to correct for the residual misalignment between the reconstructed images 𝜌 𝑖 . Therefore we propose to introduce additional displacement fields 𝑢 ̃𝑖 and their associated spatial transformation matrices 𝑇(𝑢 ̃𝑖) in the acquisition model, such that the images formed by 𝑇(𝑢 ̃𝑖)𝜌 𝑖 , for 𝑖 = {1 … 𝑁 𝑖 }, are registered. By doing so, the sharing of anatomical data through the regularization term 𝑅 will also be improved. Provided the residual motion fields 𝑢 ̃𝑖 are small (see detailed derivation in Appendix 3), ( 5) can be rewritten as:

{ min (𝜌 ̃1…𝜌 ̃𝑁𝑖 ,𝛼,𝑢 ̃1…𝑢 ̃𝑁𝑖 ) ∑‖𝐸 𝑖 (𝛼, 𝑠 𝑖 )𝜌 ̃𝑖 -𝑚 𝑖 ‖ 2 𝑁 𝑖 𝑖=1 +𝜆𝑅( 𝑇(𝑢 ̃1)𝜌 ̃1 … 𝑇(𝑢 ̃𝑁𝑖 )𝜌 ̃𝑁𝑖 ) + 𝜇 𝑆(𝛼) 𝜌 𝑖 = 𝑇(𝑢 ̃𝑖)𝜌 ̃𝑖 (7)
Here we have introduced intermediate variables 𝜌 ̃𝑖 which represent the images corrected by the motion model (𝛼) only; now the 𝜌 𝑖 denote the fully corrected images (including residual inter-image motion correction). Algorithm 2 is proposed for solving [START_REF] Marie | Detection and prediction of acute heart transplant rejection with the myocardial T2 determination provided by a black-blood magnetic resonance imaging sequence[END_REF] and will be termed multi-image GRICS (MI-GRICS). A nonrigid image registration step is introduced prior to the image reconstruction step. A reference image has to be chosen for the registration which corresponds to choosing a reference motion state.

III. METHODS

Nine heart transplant patients were recruited for this study which was approved by the local ethics committee. Written informed consent was obtained from all patients. These patients were scanned with a 1.5T Signa HDxt system

(General Electric Healthcare, Milwaukee, USA). A myocardial T 2 mapping protocol was used which has been validated for detecting and predicting allograft rejection [START_REF] Marie | Detection and prediction of acute heart transplant rejection with the myocardial T2 determination provided by a black-blood magnetic resonance imaging sequence[END_REF], [START_REF] Bonnemains | Diagnostic and prognostic value of MRI T2 quantification in heart transplant patients[END_REF]. The sequence was an ECG-triggered black-blood fast spin echo (FSE) and was repeated 10 times with different effective times of echo (TE eff ): 10, 15, 20, 30, 40, 50 , 55, 60, 70 and 80 ms. One short-axis slice was chosen in the middle of the left ventricle and each of the 10 images was acquired during a breath-hold period (gold standard, i.e. 10 separate breath-hold periods were necessary to obtain a T 2 dataset) and during free breathing (for MI-GRICS reconstructions). The matrix size was 256x160 for breath-held sequences and 256x256 for freebreathing sequences. Images were eventually reconstructed to 512x512. Other parameters were: 10 mm slice thickness, typical field of view of 320 mm, pixel size of 1.25x2 mm² (interpolated to 0.625x0.625 mm²), echo train length=16, number of excitations 𝑁𝑒𝑥 = 1, trigger time set for acquisition in end-diastole (depending on the heart rate), repetition time TR = 2 heart beats. The MR system's body coil was used for transmission and reception in breath-held acquisitions. For the free-breathing scans an 8-element cardiac coil array was used for signal reception. The body coil is used for the clinical diagnostic sequence in our center because the intensity weighting obtained with surface coils depends on the vendor's calibration and coil combination algorithms which are not always sufficiently documented. For the free-breathing sequences, we used our own calculation of coil sensitivity maps from the raw data of the vendor's calibration sequence (i.e. dividing each surface coil image by the body coil image). The coil combination was embedded in our reconstruction framework through the use of these sensitivity maps in the acquisition model. Physiological signals were collected using an in-house realtime system [START_REF] Odille | Noise Cancellation Signal Processing Method and Computer System for Improved Real-Time Electrocardiogram Artifact Correction during MRI Data Acquisition[END_REF]. Respiration signals were provided by two pneumatic belts placed on the patient's chest and abdomen. The system also recorded acquisition windows so that it was possible to relate each acquired k-space line to the corresponding motion signal samples. The two breathing signals as well as their first-order temporal derivatives were used as the 𝑠 𝑘 (𝑡) signals (i.e. 𝑁 𝑠 = 4 signals) needed by GRICS and MI-GRICS reconstructions.

Free-breathing data were reconstructed using the conventional GRICS following Algorithm 1 (independent reconstruction of each image and motion model). Here conventional GRICS is expected to yield suboptimal results since we used 𝑁𝑒𝑥 = 1 (severe ill-conditioning) as explained in the Theory section. They were also reconstructed using multi-image GRICS following Algorithm 2: (i) firstly we turned off the inter-image registration (i.e. 𝑇(𝑢 ̃𝑖) = 𝐼𝑑𝑒𝑛𝑡𝑖𝑡𝑦) and the joint image regularization (i.e. 𝑅 = 𝐼𝑑𝑒𝑛𝑡𝑖𝑡𝑦), i.e. the only difference with conventional GRICS was the use of a shared motion model; (ii) then we only turned off the interimage registration (i.e. using 𝑇(𝑢 ̃𝑖) = 𝐼𝑑𝑒𝑛𝑡𝑖𝑡𝑦); (iii) finally we used the full version of Algorithm 2 providing the complete inter-and intra-motion correction. The latter three reconstructions will be termed MI-GRICS1, MI-GRICS2 and MI-GRICS3 respectively. For the nonrigid registration step required by MI-GRICS3 we used a technique which was extensively validated for T 2 mapping in heart transplant patients [START_REF] Odille | Nonrigid registration improves MRI T2 quantification in heart transplant patient follow-up[END_REF]. It consists of a pre-processing by histogram matching in order to reduce intensity differences in the image series; then a sum-of-squared-differences similarity criterion is optimized using a multi-resolution Gauss-Newton scheme.

A. Implementation of the reconstruction

Reconstructions were implemented in C++ language. The reconstruction code was run on a cluster of workstations comprised of 16 machines on an InfiniBand network. Each machine was equipped with Intel® Xeon® Quad Core CPU (X550, 2.67 GHz) and 24 GB RAM. Our reconstruction library supports both multi-node parallelization (MPI standard) and multicore parallelization (OpenMP standard). The most computationally intensive part lies in the conjugate gradient linear solvers (step 2 and 4 of Algorithm 2).The application of the underlying operators (E H E and J H J) comprises an outer loop on motion states and an inner loop on coil receivers, which were distributed respectively on the 16 nodes (MPI parallelization) and on 8 cores (OpenMP parallelization).

The regularization tuning parameters were set as follows: 𝜆 = 10 -8 , 𝜇 = 0.005, 𝜈 = 0.01 (𝜆, 𝜇, 𝜈 were set relative to the norm of the right hand-side of the linear system in which they are involved). They were optimized empirically on one dataset and then the same values were used for all other datasets.

B. Assessment of the commutation assumption used to derive MI-GRICS

In order to derive Eq. ( 7) and the subsequent Algorithm (2), the assumption has been made that, in our application, two transformation operators can be commuted with little error (Appendix 3). The impact of this assumption was quantified retrospectively. Images and motion fields generated by MI-GRICS3 reconstruction were used for testing the commutation assumption 1 . For each patient dataset the displacement fields with maximal amplitude were extracted from both the motion model and the inter-image motion fields. Those displacement fields, 𝑈 𝑚𝑎𝑥 (1) and 𝑈 𝑚𝑎𝑥 (2) respectively, were then applied to each of the 10 images (corresponding to the 10 echo times) reconstructed by MI-GRICS3. After application of 𝑇 𝑈 𝑚𝑎𝑥

(1) 𝑇 𝑈 𝑚𝑎𝑥 [START_REF] Bilgic | Multicontrast reconstruction with Bayesian compressed sensing[END_REF] or 𝑇 𝑈 𝑚𝑎𝑥 (2) 𝑇 𝑈 𝑚𝑎𝑥 (1) , the normalized root mean squared errors (NRMSE) between the resulting images were computed. Among the 10 images, only the worst case corresponding to maximal NRMSE was recorded. We also tested the commutation assumption when the amplitude of those maximal displacements was divided by 2 and by 4.

C. Assessment of intra-image motion correction (blurring/ghosting artifacts)

The quality of intra-image motion correction was assessed using the gradient entropy as an image quality metric, as expressed by:

Here 𝜌 = [𝜌 1 … 𝜌 𝑁 𝑖 ]
𝑇 is the entire reconstructed image series, ∇ is the spatial gradient and 𝐻 is the Shannon entropy based on the probability 𝑃 of pixel value occurrences. The entropy of an information channel is increased when more uncertainty is associated with the information source and 𝐻(𝑋) is maximal when all intensity values in 𝑋 are equally probable [START_REF] Rafael | Digital Image Processing[END_REF]. Motion artifacts (blurring and ghosting) come from a dispersion of the signal in the image [START_REF] Wood | Mr Image Artifacts from Periodic Motion[END_REF]. This results in a smoothing in the image histogram or in the gradient image histogram (the gradient image highlights the details/features), and thereby an increased entropy. We computed 𝑄 𝑖𝑛𝑡𝑟𝑎 for the different reconstructions (GRICS, MI-GRICS1, MI-GRICS2 and MI-GRICS3) and tested for statistical differences with a Wilcoxon signed rank test, setting the significance level to 5%.

D. Assessment of inter-image motion correction (registration)

Inter-image motion correction was first assessed qualitatively by visualizing a movie of the reconstructed image series. An average image of the spatial gradients (i.e. averaged over the 10 echo times) was also used to represent in a single static image the amount of residual misalignment in the reconstructed series. An objective metric 𝑄 𝑖𝑛𝑡𝑒𝑟 was also computed to quantify the quality of inter-image motion correction, based on the normalized mutual information (NMI) between the first image of the reconstructed series and all other images: 1 We acknowledge a limitation of this experiment since the images and motion fields were generated by MI-GRICS3 which already uses the commutation assumption. However, the MI-GRICS3 images, as shown in the Results section, are very close to the reference breath-hold images. This suggests that the images, the motion model and the inter-image motion fields are plausible solutions and can be used as surrogate data for this experiment. ).

(9)

e NMI was chosen because it is one of the most popular similarity criterions used in medical image registration [START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast MR images[END_REF], [START_REF] Studholme | An overlap invariant entropy measure of 3D medical image alignment[END_REF] and is fairly robust to contrast changes that occur in multi-contrast or multi-modal imaging. It is fair to use this metric for validation in this study because our optimization criterions are not based on NMI (our motion model optimizes the fidelity of the acquisition model to the data; our registration optimizes the sum-of-squared-differences of the images after histogram matching). Wilcoxon signed rank tests were performed again to compare the different reconstructions.

E. Quantitative myocardial T 2 analysis

Two of the nine patients were excluded from the quantitative T 2 analysis. The trigger delay for these two patients was improperly set (with echo trains falling in the systolic contraction) leading to unreliable signal intensity.

T 2 quantification was performed in the remaining seven patients. The breath-held image series was considered to be the reference standard and was processed for T 2 quantification in 3 steps: (i) nonrigid registration of the images using the same technique as described earlier, which was previously validated for myocardial T 2 quantification [START_REF] Odille | Nonrigid registration improves MRI T2 quantification in heart transplant patient follow-up[END_REF]; (ii) manual segmentation of endocardium and epicardium (J.-M. E., 30 years of experience in experimental and clinical T 2 quantification) and separation into 6 segments according to American Heart Association (AHA) guidelines [START_REF] Cerqueira | Standardized myocardial segmentation and nomenclature for tomographic imaging of the heart: a statement for healthcare professionals from the Cardiac Imaging Committee of the Council on Clinical Cardiology of the American Heart Association[END_REF]; (iii) segment-wise calculation of T 2 values by an exponential fitting of the averaged segmental signals using a Levenberg-Marquardt algorithm. For the motion-corrected images (GRICS, MI-GRICS1, MI-GRICS2 and MI-GRICS3 reconstructions) the post-processing only consisted of steps (ii) and (iii).

Segment-wise T 2 values obtained from free-breathing data (uncorrected, GRICS, MI-GRICS1, MI-GRICS2 and MI-GRICS3 reconstructions) were compared with breath-held data. We first compared T 2 values from all myocardial segments of our patient population (N=42 segments). Then we also compared T 2 values from myocardial segments where the exponential fit was good, as defined by the coefficient of determination R 2 >0.97. This criterion has been suggested in previous studies [START_REF] Bonnemains | Diagnostic and prognostic value of MRI T2 quantification in heart transplant patients[END_REF] as a quality control index for excluding low quality datasets. The data were tested for statistical significance with Wilcoxon signed rank tests again.

Finally we compared the diagnostic accuracy of the proposed reconstruction (MI-GRICS3) with the reference technique (breath-hold): in our institution, suspicion of allograft rejection is reported if the average T 2 of the 6 myocardial segments is above 62 ms.

𝑄 𝑖𝑛𝑡𝑟𝑎 = 𝐻(|∇𝜌|), 𝑤𝑖𝑡ℎ 𝐻(𝑋) = -∑ 𝑃(𝑋 = 𝑖)log 2 (𝑃(𝑋 = 𝑖)) 𝑖 ). ( 8
)
IV. RESULTS

All datasets were successfully acquired and reconstructed with GRICS, MI-GRICS1, MI-GRICS2 and MI-GRICS3. The total effective scan time (i.e. time between the beginning of the first sequence and the end of the last sequence) was 460 s on average for the breath-hold protocol, and 520 s for the freebreathing protocol. The reconstruction times for a T 2 dataset were: 507 s with GRICS, 439 s with MI-GRICS1, 465 s with MI-GRICS2 and 524 s with MI-GRICS3 (including 34 s for inter-image registration steps). A complete dataset from patient 1 is shown in Fig. 1 as an illustration, including the 10 reconstructed images and the post-processed T 2 map. A video file for this patient is provided as supplementary multimedia material in order to better assess the image alignment with changing TE in the different reconstructions (FB, GRICS, MI-GRICS1, MI-GRICS2, MI-GRICS3) and in the registered breath-hold images (BH).

A. Commutation assumption

Worst case errors obtained in the images when commuting transformation operators are summarized in Table I. When maximal displacement fields were applied (maximal amplitude of 27.9 mm for the motion model and 11.5 mm for the inter-image registration on average), the commutation resulted in an NRMSE of 11 % on average.

When the amplitude of these displacement fields was half (respectively one quarter) of the maximum, the NRMSE dropped to 4% (respectively 1%). NRMSE( U (1) max , U (2) max ) a 0.11 0.07 0.17 0.15 0.07 0.09 0.08 0.10 0.14 0.11 NRMSE( U (1) max/2, U (2) max /2 ) 0.04 0.03 0.08 0.05 0.03 0.04 0.03 0.03 0.05 0.04 NRMSE( U (1) max/4, U (2) max /4 ) 0.01 0.01 0.03 0.02 0.01 0.01 0.01 0.01 0.02 0.01 a NRMSE(U (1) ,U (2) ): normalized root mean squared error introduced in the image when applying 𝑇 𝑈 (1) 𝑇 𝑈 [START_REF] Bilgic | Multicontrast reconstruction with Bayesian compressed sensing[END_REF] or 𝑇 𝑈 (2) 𝑇 𝑈 [START_REF] Haldar | Improved diffusion imaging through SNR-enhancing joint reconstruction[END_REF].

B. Intra-image motion correction

The visual assessment of the reconstructed images shows a drastic reduction in image artifacts with all tested reconstruction (GRICS, MI-GRICS1, MI-GRICS2 and MI-GRICS3) compared to the uncorrected free-breathing images as can be seen in Fig. 1. Differences between the different reconstructions can be observed in a typical example illustrated in Fig. 2. It can be seen that the conventional GRICS reconstruction shows residual artifacts (annotated by the arrows) which are reduced with MI-GRICS1 and further reduced with MC-GRICS2 and MC-GRICS3.

The image quality metric 𝑄 𝑖𝑛𝑡𝑟𝑎 obtained with the three reconstructions are summarized in Fig. 3 and confirm the visual assessment. The improvement with this metric from MI-GRICS1 to MI-GRICS2 was statistically significant (p=0.02).

C. Inter-image motion correction

The effect of the inter-image motion correction (misalignments) is illustrated in Fig. 4. The average gradient maps over the echo times show that the spatial gradients are more spread out in the GRICS images than in the MI-GRICS1 images; there is also an improvement from MI-GRICS1 to MI-GRICS2 and from MI-GRICS2 to MI-GRICS3. The quality of the registration seemed to be similar between MI-GRICS3 and the reference (registered) breath-held images.

The inter-image motion correction quality objective scores, as assessed by 𝑄 𝑖𝑛𝑡𝑒𝑟 , are given in Fig. 5 and are in good agreement with the visual observations. The improvement from to MI-GRICS1 was statistically significant (p=0.004); the improvement from MI-GRICS2 to MI-GRICS3 was significant as well (p=0.004). The spatial gradient images, averaged over the 10 echo times, are shown for the uncorrected (g), GRICS (h), MI-GRICS1 (i), MI-GRICS2 (j), MI-GRICS3 (k) and registered breath-held reference (l) images. The benefit of the joint reconstruction framework for inter-image motion correction is observed visually by the decreased blurring in the averaged gradient images near the endocardial and epicardial borders and in the papillary muscles (numerical results are shown in Fig. 5). MI-GRICS3 seemed to perform as well as the registered breath-held images.

FB GRICS MI-GRICS1 MI-GRICS2 BH

(a) (b) (c) (d) (e) (g) (h) (i) (k) (k) (f) (l)

MI-GRICS3

D. Quantitative T 2 analysis T 2 quantification results are summarized in Table II. The segment-wise analysis showed that T 2 values obtained with all reconstructions are close to the reference values. For all reconstructions, differences in T 2 values were not statistically significant. Considering the quality of the T 2 fit, the number of myocardial segments that were able to be processed with good confidence (exponential fit with R 2 >0.97) was gradually increased from uncorrected [START_REF] Akçakaya | Improved quantitative myocardial T2 mapping: Impact of the fitting model[END_REF], to GRICS [START_REF] Odille | Noise Cancellation Signal Processing Method and Computer System for Improved Real-Time Electrocardiogram Artifact Correction during MRI Data Acquisition[END_REF], MI-GRICS1 [START_REF] Studholme | An overlap invariant entropy measure of 3D medical image alignment[END_REF], MI-GRICS2 [START_REF] Cerqueira | Standardized myocardial segmentation and nomenclature for tomographic imaging of the heart: a statement for healthcare professionals from the Cardiac Imaging Committee of the Council on Clinical Cardiology of the American Heart Association[END_REF] and MI-GRICS3 [START_REF] Filipovic | Motion compensated generalized reconstruction for free-breathing dynamic contrast-enhanced MRI[END_REF] and approached the reference breath-hold technique [START_REF] Zhao | Low rank matrix recovery for real-time cardiac MRI[END_REF]. T 2 values in those segments only were also close to the reference values and no statistical difference was observed.

The diagnostic accuracy of T 2 quantification with MI-GRICS3 is presented in Table III. Among the seven patients, three were suspected of allograft rejection due to abnormally elevated T 2 . The proposed technique (MI-GRICS3) agreed with the reference technique (breath-hold) in all seven patients regarding the final diagnosis (suspicion of allograft rejection or not).

V. DISCUSSION

A framework has been proposed for the joint reconstruction of multiple images and motion and has been applied to freebreathing myocardial T 2 quantification in heart transplant patients. Our results show that the joint reconstruction of multiple images and motion (MI-GRICS1, MI-GRICS2 and MI-GRICS3) has a benefit over the joint reconstruction of a single image and motion (GRICS). Specifically, the analysis of the intra-and inter-image motion correction metrics revealed the following significant improvements: (i) the sharing of the motion model improved inter-image correction over conventional GRICS (improved 𝑄 𝑖𝑛𝑡𝑒𝑟 ) which might be explained by a more stable fit of the motion model due to dimensionality reduction; (ii) the joint image regularization improved intra-image correction (improved 𝑄 𝑖𝑛𝑡𝑟𝑎 ) which might be explained by the more stable image reconstruction step; (iii) the added inter-image registration improved interimage correction further (improved 𝑄 𝑖𝑛𝑡𝑒𝑟 ). The latter improvement was implemented with a limited extra computational cost of 13% compared to MI-GRICS2, including 9% for the registration steps alone. This was possible due to the commutation assumption described in Appendix 3 and we have checked that this assumption did not introduce large errors in our acquisition model, even in the worst case motion amplitudes. It should be noted that MI-GRICS3 is different from applying MI-GRICS2 followed by registration. This is because the inter-image motion fields ensure that the images are optimally aligned when applying the multi-image regularizer. In our results, although this did not translate into a statistically significant improvement of the gradient entropy metric, a reduction of the artifacts was visible in some cases as shown in Fig. 2.

The results also show the feasibility of using the technique for myocardial T 2 quantification in heart transplant patients. Motion-compensation strategies have been proposed previously for such applications but have been restricted to the registration of relatively low-resolution single-shot images [START_REF] Giri | Myocardial T2 mapping with respiratory navigator and automatic nonrigid motion correction[END_REF], [START_REF] Roujol | Impact of motion correction on reproducibility and spatial variability of quantitative myocardial T2 mapping[END_REF]. The proposed method should be particularly useful in patients who have difficulties holding their breath repeatedly and when lengthy high-resolution (multi-shot) images are desired. Another benefit is that the images generated by the reconstruction are already registered and are readily available for segmentation. Importantly segment-wise myocardial T 2 values were found to be in good agreement with the reference breath-held images. The bias (62.5 ms for MI-GRICS3 against 62.2 ms for the reference) was small compared to the inter-observer variability of T 2 measurement with this technique which was evaluated to be 2.7 ms in Ref [START_REF] Odille | Nonrigid registration improves MRI T2 quantification in heart transplant patient follow-up[END_REF]. The diagnosis with MI-GRICS3 was in agreement with that of the reference protocol in all seven patients. The acquisition time was nearly the same for the free-breathing and for the breath-held protocols although the k-space size was increased in the free-breathing experiments (from 160 to 256 in the phase encoding dimension). This was possible because the breath-held protocol comprises many idle periods so that the patient can rest between successive apneas, unlike the free-breathing protocol. Finally the reconstruction time for a T 2 dataset was 9 min which is close to the scan time, making the technique useable clinically.

In this study we used a generally applicable constraint to take advantage of anatomical redundancies in the images based on the spatial gradients (i.e. gradient co-occurrence constraint). Other redundancy constraints would be of interest, in particular those based on joint statistics of the images such as the joint entropy or the mutual information [START_REF] Somayajula | PET Image Reconstruction Using Information Theoretic Anatomical Priors[END_REF]. As an alternative in specific applications it might be possible to impose an analytical MR signal model [START_REF] Sumpf | Model-based nonlinear inverse reconstruction for T2 mapping using highly undersampled spin-echo MRI[END_REF], [START_REF] Dikaios | Direct parametric reconstruction from undersampled (k, t)-space data in dynamic contrast enhanced MRI[END_REF] (exponential models for T 2 or T 1 mapping, pharmacokinetic models for perfusion imaging etc…), mathematical models such as Bsplines [START_REF] Filipovic | Motion compensated generalized reconstruction for free-breathing dynamic contrast-enhanced MRI[END_REF], low-rank signal evolution models [START_REF] Liang | Dynamic imaging by model estimation[END_REF]- [START_REF] Otazo | Low-rank plus sparse matrix decomposition for accelerated dynamic MRI with separation of background and dynamic components[END_REF], or sparsity constraints in other domains such as the x-f space [START_REF] Jung | Improved k-t BLAST and k-t SENSE using FOCUSS[END_REF]- [START_REF] Usman | Motion corrected compressed sensing for freebreathing dynamic cardiac MRI[END_REF] or x-p space with p the varying acquisition parameter space [START_REF] Doneva | Compressed sensing reconstruction for magnetic resonance parameter mapping[END_REF]. Other models based on the singular value decomposition have also been proposed recently for accelerated parameter mappings [START_REF] Velikina | Accelerating MR parameter mapping using sparsitypromoting regularization in parametric dimension[END_REF]. Here the multi-image reconstruction framework was applied to multi-contrast MRI datasets. It might also be adapted to dynamic imaging (with or without significant contrast changes) such as perfusion imaging or cardiac cine imaging and is therefore conceptually close to the recent work in [START_REF] Lingala | Deformation Corrected Compressed Sensing (DC-CS): A Novel Framework for Accelerated Dynamic MRI[END_REF]. Another possible use of the framework would be in the context of multi-modal imaging. Finally other motion sensors might be used as surrogates of (or in combination with) the external sensors used here. In particular, partial MRI data (navigators or self-gated approaches) [START_REF] Pang | ECG and navigator-free fourdimensional whole-heart coronary MRA for simultaneous visualization of cardiac anatomy and function[END_REF], [START_REF] Feng | XD-GRASP: Golden-angle radial MRI with reconstruction of extra motion-state dimensions using compressed sensing[END_REF] might give signals better correlated with the motion of organs and improve the results further.

The joint reconstruction of image and motion, as defined in Eqs. [START_REF] Haldar | Improved diffusion imaging through SNR-enhancing joint reconstruction[END_REF][START_REF] Knoll | Simultaneous MR-PET Reconstruction using Multi Sensor Compressed Sensing and Joint Sparsity[END_REF][START_REF] Marie | Detection and prediction of acute heart transplant rejection with the myocardial T2 determination provided by a black-blood magnetic resonance imaging sequence[END_REF], is a non-convex optimization procedure. The proposed approach relies on a multi-scale implementation and a linearization with respect to the motion parameters (Appendix 2). The latter step involved the optic flow equation which is known to be ill-posed (one equation with two unknowns, known as the aperture problem). Occlusions are also known to render motion estimation difficult. Here occlusions might result from large motion as structures move in and out of the imaging plane. This might cause the algorithms to be trapped in a local minimum. When applied in 2D, the method may fail when the slice thickness is small compared to the through-plane motion. In practice the clinical slice thicknesses of 8-10 mm used in cardiac MRI generally give good results. Otherwise the method might be applied in 3D to resolve through-plane motion.

One limitation of the method is the choice of the registration technique which might not be suitable for all applications. The use of a mutual information based similarity criterion might improve the registration step in Algorithm 2 slightly. More sophisticated techniques might still be necessary in case of complex contrast changes from one image to another, as can occur for instance in first-pass perfusion imaging. Another limitation of this study is the small patient population.

In conclusion we have presented a joint optimization technique that allows reconstructing N images together with both intra-image and inter-image motion models. The method allows taking advantage of anatomical redundancies in the image dataset despite the complex organ deformations during the acquisition. The feasibility of the technique was demonstrated with free-breathing myocardial T 2 quantification in heart transplant patients providing equivalent diagnostic outcome compared to the reference technique.

Algorithm 1 : 4 .

 14 Joint reconstruction of one image and motion (GRICS) min 𝜌 ‖𝐸(𝛼, 𝑠)𝜌 -𝑚‖ 2 + 𝜆𝑅(𝜌) min 𝛿𝛼 ‖𝐸(𝛼 + 𝛿𝛼, 𝑠)𝜌 -𝑚‖ 2 + 𝜇 𝑆(𝛼 + 𝛿𝛼) = min 𝛿𝛼 ‖𝐽(𝛼, 𝑠, 𝜌)𝛿𝛼 -𝜀‖ 2 + 𝜇 𝑆(𝛼 + 𝛿𝛼) 𝛼 = 𝛼 + 𝛿𝛼 Input: the k-space data 𝑚; the motion sensor signals 𝑠 Output: one image 𝜌 and one motion model 𝛼 Initialize the motion model : 𝛼 = 0 for resolution level = coarse to fine for Gauss-Newton iteration = 1 to N 1. Image reconstruction 2. Compute residual error 𝜀 = 𝑚 -𝐸(𝛼, 𝑠)𝜌 3. Motion model optimization (Gauss-Newton step with 𝐽 the corresponding Jacobian): Motion model update end for end for reconstruction of the 𝑁 𝑖 images and their corresponding motion models 𝛼 1 … 𝛼 𝑁 𝑖 into a single optimization problem:

min(

  𝜌 1 …𝜌 𝑁 𝑖 ,𝛼 1 …𝛼 𝑁 𝑖 ) ∑‖𝐸 𝑖 (𝛼 𝑖 , 𝑠 𝑖 )𝜌 𝑖 -𝑚 𝑖 ‖

Algorithm 2 :

 2 Joint reconstruction of multiple images and motion (𝛼, 𝑠 𝑖 )𝜌 ̃𝑖 -𝑚 𝑖 ‖ 2 𝑁 𝑖 𝑖=1 + 𝜆𝑅( 𝑇(𝑢 ̃1)𝜌 ̃1 … 𝑇(𝑢 ̃𝑁𝑖 )𝜌 ̃𝑁𝑖 ) min 𝛿𝛼 ∑‖𝐸 𝑖 (𝛼 + 𝛿𝛼, 𝑠 𝑖 )𝜌 ̃𝑖 -𝑚 𝑖 ‖ 2 𝑁 𝑖 𝑖=1 + 𝜇 𝑆(𝛼 + 𝛿𝛼) = min 𝛿𝛼 ‖𝐽(𝛼, 𝑠, 𝜌 ̃)𝛿𝛼 -𝜀‖ 2 + 𝜇 𝑆(𝛼 + 𝛿𝛼) 𝛼 = 𝛼 + 𝛿𝛼 Input: the k-space data 𝑚 1 … 𝑚 𝑁 𝑖 ; the motion sensor signals 𝑠 1 … 𝑠 𝑁 𝑖 Output: images 𝜌 1 … 𝜌 𝑁 𝑖 , a motion model 𝛼 and residual inter image motion fields 𝑢 ̃1 … 𝑢 ̃𝑁𝑖 Initialize the motion model : 𝛼 = 0 Initialize the images 𝜌 ̃𝑖 with a naïve reconstruction and choose one as the reference 𝜌 ̃𝑟𝑒𝑓 for resolution level = coarse to fine for Gauss-Newton iteration = 1 to N 1. Inter-image registration 2. Multi-image reconstruction 3. Compute residual error for each image 𝜀 𝑖 = 𝑚 𝑖 -𝐸 𝑖 (𝛼, 𝑠 𝑖 )𝜌 ̃𝑖 4. Motion model optimization (Gauss-Newton step with 𝐽 the corresponding Jacobian): 5. Motion model update end for end for Finally apply the last inter-image warping fields for each image: 𝜌 𝑖 = 𝑇(𝑢 ̃𝑖)𝜌 ̃𝑖

  𝑋, 𝑌) = -∑ 𝑃(𝑋 = 𝑖, 𝑌 = 𝑗)log 2 (𝑃(𝑋 = 𝑖, 𝑌 = 𝑗)) 𝑖,𝑗

Fig. 1 .

 1 Fig. 1. Example myocardial T2 quantification image series (effective echo times TE ranging from 10 to 80 ms) and resulting T2 map obtained from freebreathing experiments: uncorrected (a); with conventional GRICS reconstruction, i.e. independent reconstruction of each image (b); with MI-GRICS1 reconstruction, i.e. GRICS with a shared motion model (c); with MI-GRICS2 reconstruction, i.e. joint reconstruction of all images and a motion model (d); with MI-GRICS3, i.e. joint reconstruction of all images, intra-and inter-image motion (e). The breath-held reference images and T2 maps are shown in (f) (the noise region including blood has been masked out).

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. Comparison of the image quality in patient 2 (image at TEeff=70 ms) obtained from free-breathing data, including uncorrected (a), GRICS (b), MI-GRICS1 (c),MI-GRICS2 (d) and MI-GRICS3 (e) correction, with the breath-held reference data (f). Noticeable ghosting artifacts in the uncorrected image are annotated by the arrows. A gradual improvement can be observed in these areas from (b) to (e) illustrating the benefit of the joint reconstruction framework in terms of reducing ghosting/blurring artifacts.

Fig. 5 .

 5 Fig.5. Quality of the inter-image motion correction, as assessed by a metric based on normalized mutual information, obtained with the proposed reconstructions. The improvement from GRICS to MI-GRICS1, as well as the improvement from MI-GRICS2 to MI-GRICS3, were statistically significant.GRICS MI-GRICS1 MI-GRICS2

APPENDIX 1 FORWARD

 1 MODEL OF MRI ACQUISITION IN THE PRESENCE OF PATIENT MOTIONThe forward model of the MRI acquisition process in the presence of motion is described by the following block matrix:

  𝑢(𝑥, 𝑦, 𝑧, 𝑡) ≃ ∑ 𝛼 𝑘 (𝑥, 𝑦, 𝑧)𝑠 𝑘(𝑡)

	𝑁 𝑠
	𝑘=1

separable motion model has been proposed in order to express 𝑢 as a linear combination of a small number 𝑁 𝑠 of signals obtained from motion sensors 𝑠 𝑘 (t) (e.g. respiratory belts, navigators…):

TABLE II MYOCARDIAL

 II T2 QUANTIFICATION RESULTS FROM 7 HEART TRANSPLANT PATIENTS (42 MYOCARDIAL SEGMENTS IN TOTAL) Statics performed on myocardial segments where a good fit was found with either the breath-hold or the tested reconstruction, i.e. segments such that R 2 (breathhold)>0.97 or R 2 (reconstruction)>0.97.

		Uncorrected	GRICS	MI-GRICS1	MI-GRICS2	MI-GRICS3	Breath-hold (reference for T2)
	T2 in ms (all myocardial segments)	64.9 ± 13.2 (p=0.50)	62.6 ± 11.1 (p=0.66)	62.8 ± 12.1 (p=0.73)	62.4 ± 11.8 (p=0.46)	62.5 ± 11.1 (p=0.49)	62.2 ± 11.2
	Number of
	segments with good	12	18	23	24	27	29
	T2 fit (R 2 >0.97)
	T2 in ms (segments with good fit) †	64.1 ± 11.7 (p=0.97)	62.5 ± 10.8 (p=0.44)	62.9 ± 11.8 (p=0.43)	62.6 ± 11.4 (p=0.24)	63.0 ± 11.2 (p=0.32)	64.2 ± 9.9
	†

  𝑇(𝑢 𝑡 1 ) ⋮ 𝜉 𝑡 𝑁 𝐹𝐶 1 𝑇(𝑢 𝑡 𝑁 ) ⋮ 𝜉 𝑡 1 𝐹𝐶 𝑁 𝑐 𝑇(𝑢 𝑡 1 ) ⋮ 𝜉 𝑡 𝑁 𝐹𝐶 𝑁 𝑐 𝑇(𝑢 𝑡 𝑁 ) ] , (A1) where the acquisition has been split into 𝑁 motion states 𝑡 1 … 𝑡 N and 𝑁 𝑐 parallel receiver coils. Each motion state 𝑡 n is described by a dense displacement vector field 𝑢 𝑡 𝑁 associated with a spatial transformation matrix 𝑇(𝑢 𝑡 𝑛 ) (sparse matrix of size 𝑁 𝜌 × 𝑁 𝜌 ) warping an image from the reference motion state to the 𝑛 th motion state. 𝑇(𝑢 𝑡 𝑛 ) are interpolation operators; if a linear interpolation kernel is chosen, it has 2 𝑁 𝑑 elements per row (i.e. 4 in 2D and 8 in 3D). The 𝜉 𝑡 𝑛 matrices are sparse and describe the sampling of k-space data during the 𝑛 th motion state. 𝐹 is the Fourier transform operator (of size 𝑁 𝜌 × 𝑁 𝜌 ) and 𝐶 1 … 𝐶 𝑁 𝑐 are the coil sensitivity 𝐽(𝛼, 𝑠, 𝜌)𝛿𝛼 = 𝐸(𝛼 + 𝛿𝛼)𝜌 -𝐸(𝛼)𝜌 = [ 𝜉 𝑡 1 𝐹𝐶 1 (𝑇(𝛼 + 𝛿𝛼, 𝑠 𝑡 1 )𝜌 -𝑇(𝛼, 𝑠 𝑡 1 )𝜌) ⋮ 𝜉 𝑡 𝑁 𝐹𝐶 1 (𝑇(𝛼 + 𝛿𝛼, 𝑠 𝑡 𝑁 )𝜌 -𝑇(𝛼, 𝑠 𝑡 𝑁 )𝜌) ⋮ 𝜉 𝑡 1 𝐹𝐶 𝑁 𝑐 (𝑇(𝛼 + 𝛿𝛼, 𝑠 𝑡 1 )𝜌 -𝑇(𝛼, 𝑠 𝑡 1 )𝜌) ⋮ 𝜉 𝑡 𝑁 𝐹𝐶 𝑁 𝑐 (𝑇(𝛼 + 𝛿𝛼, 𝑠 𝑡 𝑁 )𝜌 -𝑇(𝛼, 𝑠 𝑡 𝑁 )𝜌) ] 𝐹𝐶 1 s 𝑘 (t 1 )∇(𝑇(𝛼, 𝑠 𝑡 1 )ρ) 𝑇 ⋮ 𝜉 𝑡 𝑁 𝐹𝐶 1 s 𝑘 (t 𝑁 )∇(𝑇(𝛼, 𝑠 𝑡 𝑁 )ρ) 𝜉 𝑡 1 𝐹𝐶 𝑁 𝑐 s 𝑘 (t 1 )∇(𝑇(𝛼, 𝑠 𝑡 1 )ρ) 𝑇 ⋮ 𝜉 𝑡 𝑁 𝐹𝐶 𝑁 𝑐 s 𝑘 (t 𝑁 )∇(𝑇(𝛼, 𝑠 𝑡 𝑁 )ρ)

					(A2)
	= ∑ 𝑁 𝑠 𝑘=1	[	𝜉 𝑡 1 𝑇	]	𝛿𝛼 𝑘 .

𝑇 ⋮

Copyright (c) 2015 IEEE. Personal use of this material is permitted. However, permission to use this material for any other purposes must be obtained from the IEEE by sending a request to pubs-permissions@ieee.org.

weightings (diagonal matrices of size 𝑁 𝜌 × 𝑁 𝜌 ) of the 𝑁 𝑐 receiver coils.

In order to reduce the complexity of the problem the number of motion states is reduced based on a clustering of the motion sensor signal values so that k-space data acquired in similar motion states are grouped into a single motion state. Further implementation details can be found in Ref [START_REF] Odille | Generalized reconstruction by inversion of coupled systems (GRICS) applied to free-breathing MRI[END_REF].

APPENDIX II JACOBIAN MATRIX DESCRIBING THE PROPAGATION OF SMALL ERRORS IN THE MOTION MODEL

The Jacobian of the data fidelity term in Eq. (3) needs to be computed with respect to 𝛼 in order to solve this nonlinear least squares optimization problem. The Jacobian describes the effect of a small perturbation 𝛿𝛼 of the motion model. Recalling the optic flow equation (given an image 𝐼, 𝐼(𝑢 + 𝛿𝑢) -𝐼(𝑢) = 𝛻𝐼 𝑇 • 𝛿𝑢), the Jacobian is written and expanded below:

APPENDIX III INTEGRATION OF RESIDUAL INTER-IMAGE WARPING IN THE JOINT RECONSTRUCTION FRAMEWORK

Residual inter-image motion fields are estimated by registration of the reconstructed images 𝜌 ̃𝑖 so that the registered images are given by 𝜌 𝑖 = 𝑇(𝑢 ̃𝑖)𝜌 ̃𝑖. In order to integrate this residual motion in the forward acquisition model, 𝑇(𝑢 ̃𝑖) -1 should be applied after the estimate (𝛼, 𝑠 𝑖 ) is applied to the unknown solution image. The new forward acquisition model 𝐸 ̃𝑖 is accurately written as:

⋮ 𝜉 𝑖,𝑡 𝑁 𝐹𝐶 1 𝑇(𝑢 ̃𝑖) -1 𝑇(𝛼, 𝑠 𝑡 𝑁 ) ⋮ 𝜉 𝑖,𝑡 1 𝐹𝐶 𝑁 𝑐 𝑇(𝑢 ̃𝑖) -1 𝑇(𝛼, 𝑠 𝑡 1 ) ⋮ 𝜉 𝑖,𝑡 𝑁 𝐹𝐶 𝑁 𝑐 𝑇(𝑢 ̃𝑖) -1 𝑇(𝛼, 𝑠 𝑡 𝑁 ) ]

𝜌 𝑖 (A3)

Here we propose to simplify (A3). Under the assumption that 𝑢 ̃𝑖 is small and 𝛼 is sufficiently smooth, the two successive transformation operators can be commuted. Indeed, if 𝑣 and 𝑤 are two displacement fields and 𝜌 an image, for 𝑣 small, using the optic flow equation we have:

Furthermore if the Jacobian of 𝑤 is close to 1 (which is true for w sufficiently smooth), ∇(𝑇(𝑤)𝜌) ≃ 𝑇(𝑤)∇𝜌 and therefore 𝑇(𝑣)𝑇(𝑤) ≃ 𝑇(𝑤)𝑇(𝑣).

Finally, under the assumption that the transformations can be commuted (this assumption is termed commutation assumption in the remainder of the paper), the residual interimage motion operators can be split away from the acquisition matrix:

Using the intermediate variables 𝜌 ̃𝑖 for the optimization (instead of 𝜌 𝑖 ), the inter-image motion correction can be incorporated with a simple change of the regularization term as described in Eq. ( 7). This implementation was chosen because: (i) it was found to be more stable than introducing an additional motion operator in the forward model, which was found to lead to faster convergence and (ii) it saves computation time since the additional operators only need to be applied in the regularization step.