3D free-breathing cardiac magnetic resonance fingerprinting

Purpose: To develop a novel respiratory motion compensated three-dimensional (3D) cardiac magnetic resonance fingerprinting (cMRF) approach for whole-heart myocardial T 1 and T 2 mapping from a free-breathing scan. Methods: Two-dimensional (2D) cMRF has been recently proposed for simultaneous, co-registered T 1 and T 2 mapping from a breath-hold scan; however, coverage is limited. Here we propose a novel respiratory motion compensated 3D cMRF approach for whole-heart myocardial T 1 and T 2 tissue characterization from a free-breathing scan. Variable inversion recovery and T 2 preparation modules are used for parametric encoding, respiratory bellows driven localized autofocus is proposed for beat-to-beat translation motion correction and a subspace regularized reconstruction is employed to accelerate the scan. The proposed 3D cMRF approach was evaluated in a standardized T 1 / T 2 phantom in comparison with reference spin echo values and in 10 healthy subjects in comparison with standard 2D MOLLI, SASHA and T 2 -GraSE mapping techniques at 1.5 T. Results: 3D cMRF T 1 and T 2 measurements were generally in good agreement with reference spin echo values in the phantom experiments, with relative errors of 2.9% and 3.8% for T 1 and T 2 ( T 2 < 100 ms), respectively. in vivo left ventricle (LV) myocardial T 1 values were 1054 ± 19 ms for MOLLI, 1146 ± 20 ms for SASHA and 1093 ± 24 ms for the proposed 3D cMRF; corresponding T 2 values were 51.8 ± 1.6 ms for T2-GraSE and 44.6 ± 2.0 ms for 3D cMRF. LV coefficients of variation were 7.6 ± 1.6% for MOLLI, 12.1 ± 2.7% for SASHA and 5.8 ± 0.8% for 3D cMRF T 1 , and 10.5 ± 1.4% for T2-GraSE and 11.7 ± 1.6% for 3D cMRF T 2 . Conclusion: The proposed 3D cMRF can provide whole-heart, simultaneous and co-registered T 1 and T 2 maps with accuracy and precision comparable to those of clinical standards in a single free-breathing scan of about 7 min. T 3D following field of ; ; slice echo )/repetition = ms; sinusoidal flip angle in the range 5-10 (cid:3) ; gradient echo readout (gradient spoiling, ie FFE/FISP); variable density spiral sampling, 51 (cid:5) 35 × undersampled in the periphery of k -space ( (cid:5) 18 × undersampled the central 10% radius of k -space); spiral acquisition window 3.6 ms; 30 spirals per shot; cardiac acquisition ms; time-points acquisition


| INTRODUCTION
Quantitative myocardial tissue characterization has emerged as an important tool for the diagnosis of cardiovascular disease. Tissue parameters such as T 1 , T 2 , T 2 * relaxation times and extracellular volume (ECV) have demonstrated sensitivity to several diseases. 1 These parameters detect underlying mechanisms of disease such as necrosis (T 1 ), edema (T 2 ), iron overload (T 2 *) or fibrosis (ECV), allowing disease progression to be monitored. Conventional myocardial mapping techniques map a single parameter per acquisition using two-dimensional (2D) breath-hold sequences such as modified Look-Locker inversion recovery (MOLLI) 2 and saturation recovery single-shot acquisition (SASHA) 3 for T 1 , or T 2 prepared balanced steady-state free precession (T2prep-bSFFP) 4 and T 2 gradient and spin echo (T 2 -GraSE) 5 for T 2 . Multiple tissue parameters are generally sensitive to a given disease, and recently it has been suggested that simultaneous multi-parametric mapping could improve tissue characterization, as recommended in the updated diagnostic criteria for myocardial inflammation. 6 For example, simultaneous T 1 and T 2 can lead to improved characterization of myocardial infarction, 7,8 help distinguish between acute and chronic myocarditis (leading to improved diagnostic specificity) 6,9 and potentially improve the understanding of the underlying disease. 10 With conventional approaches, multi-parametric mapping requires sequential acquisitions under several breath-holds, often leading to non-registered maps as well as potential patient discomfort. Moreover, conventional mapping techniques provide limited coverage of the heart, and different methods have different confounding factors, including crossover dependencies between T 1 and T 2 if the two parameters are not simultaneously included in the model. 4,[11][12][13] Due to T 1 's and T 2 's complementary diagnostic information and co-dependences in the signal model, joint parametric mapping methods have become of interest. An interleaved 2D T 1 /T 2 encoding strategy based on inversion recovery (IR) and T 2 preparation (T2prep) pulses in freebreathing has been developed. 14 Most cardiac mapping methods acquire data in a small cardiac acquisition window; however, 2D joint T 1 /T 2 mapping has been achieved using IR interrupted b-SSFP readouts in CABIRIA. 15 A combination of T2prep and saturation recovery pulses have also been employed for 2D joint T 1 /T 2 mapping. 16 Similarly, T2preps together with IR pulses have also been developed for this purpose. 17,18 A hybrid T2prep-IR pulse is employed in MR multitasking for motion resolved simultaneous T 1 /T 2 mapping. 19 Extension to three dimensions for wholeheart coverage T 1 /T 2 mapping has been considered using T2prep and IR pulses in three-dimensional (3D-)QALAS in a breath-hold; however, breath-hold duration can impose a limit on spatial resolution and/or coverage. 20 More recently, 3D free-breathing T 1 /T 2 mapping has been achieved using T2prep, saturation recovery pulses and respiratory diaphragmatic navigators; however, this can lead to long and unpredictable scan times. 21 All of the above are steady-state approaches, often sampling discrete points along the relaxation curves and relying on simplified exponential models.
Magnetic resonance fingerprinting (MRF) 22 is an established research technique for simultaneous and multi-parametric quantitative MR where system parameters (such as flip angle and repetition time) vary in time to sensitize the sequence to tissue parameters of interest (such as T 1 and T 2 ). MRF attempts to capture the parametric encoding within the transient state using directly the Bloch equations and can be extended to include not only additional parameters of interest (eg flow, 23 chemical exchange saturation transfer (CEST) 24 ), but also additional model corrections (eg B 1 , 25 slice profile 26 ). Initially proposed for brain acquisitions, MRF has since been extended to 2D cardiac imaging (cMRF) 27 for simultaneous, co-registered T 1 and T 2 myocardial mapping in a single breath-hold. 2D breath-hold cMRF produces simultaneous and co-registered T 1 /T 2 maps comparable to conventional methods, albeit in reduced scan time. Model corrections such as B 1 and slice profile for 2D cMRF have also been evaluated, 28 leading to reduced bias in the estimated parameters. Simultaneous multi-slice imaging has been combined with 2D cMRF 29 to increase the coverage of the heart; however, this is still limited to three slices and is sensitive to through-plane motion.
In this work we develop a novel respiratory motion compensated 3D cMRF approach for whole-heart coverage of T 1 /T 2 myocardial mapping in a single free-breathing scan. Similarly to 2D cMRF, parametric encoding is primarily achieved with IR and T2prep pulses. The extension to three dimensions increases the scan time beyond a common breath-hold duration and free-breathing acquisitions are required. In-plane translational motion due to breathing is corrected with respiratory bellows driven localized autofocus 30,31 , and a patch-based low-rank tensor regularized reconstruction 32 is employed to accelerate the 3D cMRF scan. The proposed 3D cMRF was evaluated in a standardized phantom and in 10 healthy subjects, in comparison with gold standard spin echo measurements (phantom) and conventional MOLLI, SASHA and T 2 -GraSE mapping techniques (in vivo).

| Acquisition
The proposed 3D cMRF uses a free-breathing, ECG-triggered acquisition with a variable density stack of spirals, where a different preparation pulse is applied in each heartbeat ( Figure 1A). The preparation scheme was chosen heuristically, inspired by recent works in 2D cMRF. 27,28 The acquisition is divided into blocks of 18 heartbeats with the following preparation scheme: IR15, NP, NP, NP, NP, IR150, NP, NP, NP, NP, IR300, NP, T2prep20, T2prep30, T2prep40, T2prep50, T2prep60, T2prep80, where NP denotes "no preparation pulse" and the numerical values after IR or T2prep denote the corresponding IR time delay (T I ) or T2prep echo time (T ET2p ), respectively. One block was used per slice encoding (with linear F I G U R E 1 A, Acquisition sequence for the proposed 3D cMRF framework. IR pulses are used to encode T 1 , T2prep pulses are used to encode T 2 and SPIR modules provide fat suppression; T I denotes the inversion time delay and T ET2p denotes the T2prep echo time. The diagram depicts the encoding for one slice (repeated for all slices). The same sinusoidal flip angle pattern varying from 5 to 10 degrees was employed at each heart beat (as indicated in heartbeats 2, 3, 4 and 5 of the diagram). B, The acquired MRF k-space and respiratory bellows are used to drive an autofocus algorithm in order to correct for beat-to-beat translational motion of the heart (five potential translationally corrected images x α shown for the mid-slice; the autofocus solution is highlighted with the red box). Subsequently, the data is reconstructed with 3D LRI-HDPROST (six singular images shown for the mid-slice) and standard template matching to produce 3D whole-heart T 1 and T 2 maps (basal, mid-and apical slices shown for one representative subject) ordering) and a minimum of a 4 s recovery period was employed between blocks to recover longitudinal magnetization (after 4 s the sequence resumed when the subsequent R wave was detected). A spectral presaturation with inversion recovery (SPIR) 33 pulse is used before data acquisition to suppress aliasing artefacts originating from the fat signal. The effective flip angle of the SPIR (which determines the zero crossing of the fat) was varied to account for the preparation pulse employed in each heartbeat. For IR delay < 200 ms, no SPIR pulse was employed, as the longitudinal fat magnetization is still negative. Otherwise, the SPIR flip angle θ at each heartbeat was given by where T 1f and T 2f are the T 1 and T 2 relaxation times of fat, T 2prep denotes the echo time of the T 2 preparation employed, T If is the inversion time delay of the SPIR and HR is the nominal heartrate for a given subject. Adiabatic pulses were used for all 180 pulses in IR and T2prep modules to reduce sensitivity to potential B 1 errors 34,35 and subsequent errors in the parametric maps. 28 In addition to ECG, respiratory bellows were used to record respiratory motion. The one-dimensional (1D) respiratory signal obtained from the bellows was later used for respiratory motion correction.

| Motion correction and reconstruction
The proposed 3D cMRF respiratory motion correction and reconstruction is summarized in Figure 1B. Respiratory motion of the heart is commonly monitored with diaphragmatic navigators 36 or self-navigation 37 ; however, these methods have generally been validated in steady-state imaging applications. Diaphragmatic navigator gating can lead to long and unpredictable scan time, since data is acquired only within a giving gating window. Moreover, the parametric encoding can deviate from the prescribed non-gated sequence (particularly for IR based sequences where T 1 is encoded across multiple shots), as data outside the gating window is rejected and re-acquired. Here, respiratory bellows are employed to be less sensitive to the variable contrasts employed in 3D cMRF, which could otherwise affect motion estimation. Respiratory bellows provide only a relative 1D signal of motion in arbitrary units, r(t). In order to correct for translational respiratory motion this signal is used to drive a localized autofocus algorithm, 30 similar to previous work in abdominal 31 and cardiac 38 imaging. This approach reconstructs a set of 3D images x α from all the acquired data and evaluates an image quality metric on these images to determine the optimal motion correction. To simplify the autofocus search space, we assume a linear relationship between the bellows signal and the respiratory motion of the heart, similar to previous works. 31,38 Motion estimation is performed for each spatial dimension sequentially. Initially a set of αr(t) scaled respiratory bellows signals are considered. A corresponding bank of translationally corrected images x α is obtained by reconstructing k-space translationally corrected with αr(t) (using the inverse non-uniform Fourier transform). Following, the optimal scalingα (which determines the estimated respiratory translational motion in milli- where h α is the normalized spatial gradient and x α (i) is the ith pixel intensity. 31 The local gradient entropy is computed around a small region of interest, manually selected around the heart for each subject. The resulting optimal scaling for each dimension is used for global translation motion correction (ie corresponding phase shifts in k-space) prior to MRF reconstruction.
Although MRF initially used zero-filled reconstructions, 22 it has been shown that some aliasing artefacts can propagate into the parametric maps for highly accelerated acquisitions such as those required in cMRF. More recently it has been shown that advanced reconstructions of the MRF time-series can improve the resulting T 1 /T 2 estimates. 32,[39][40][41][42][43] In addition to parallel imaging, 44,45 subspace modelled reconstructions 41,42 and compressed sensing regularization 32,43,46 have been applied to MRF reconstruction, improving accuracy and precision. For 3D cMRF, we use a low-rank inversion 41 (LRI) regularized with a high-dimensionality undersampled patch based reconstruction 32 (HDPROST). The LRI-HDPROST reconstruction solves the following problem: where A, U R , F and C are sampling, temporal compression (obtained from a truncated singular value decomposition of the MRF dictionary), nonuniform fast Fourier transform and coil sensitivity operators; P b constructs a 3D local tensor T b around voxel b by concatenating local (within a patch), non-local (between similar patches) and contrast (in the singular value domain) voxels along each dimension; x are the singular images reconstructed and k 0 is the translationally corrected k-space. The problem above was solved via the alternating direction method of multipliers (ADMM) 47 ; additional details can be found in Reference 32.

| Simulations
The sensitivity of the proposed sequence to variable heart rates was investigated in simulations. Extended phase graphs (EPGs) 48,49 were employed to simulate the response of the magnetization to a sequence of 18 heartbeats, corresponding to a slice encoding (as depicted in Figure 1A). Four RR intervals were considered for the simulations: RR = [1500, 1000, 750, 600]ms (corresponding to [40,60,80, 100] beats per minute (bpm)). In each case, the RR intervals were further perturbed by random RR variations drawn for a Gaussian distribution of zero mean and standard deviation of [50,300] ms. Finally, each simulated fingerprint was corrupted by Gaussian noise of zero mean and standard deviation of

| Phantom acquisition
The proposed 3D cMRF was evaluated in a standardized phantom and in 10 healthy subjects (five males, age 30 ± 2 years) on a 1.5 T Ingenia MR system (Philips, Best, The Netherlands) using a 28-channel cardiac coil. The study was approved by the institutional review board and written informed consent was obtained from all subjects according to institutional guidelines.

| In vivo acquisition
Ten healthy subjects were scanned with the same 3D cMRF protocol as used in the phantom, in short axis orientation and with subject specific mid-diastolic triggering (determined by preceding CINE imaging). ECG and respiratory bellows signals were recorded for dictionary simulation and translational motion correction, respectively. Conventional 2D MOLLI, SASHA and T 2 -GraSE acquisitions were performed in basal, mid-and apical slices (for a total of nine breath-holds). All conventional mapping methods used FOV = 300 × 300 mm 2 , resolution = 2 × 2 mm 2 and slice thick-

| Motion correction and reconstruction
Motion estimation was performed in vivo via respiratory bellows driven autofocus, as described above. The bellows signal (r(t)) was normalized to unity and a set of trial motion signals αr(t) was obtained using α = [0 : 0.1 : 1]β; consequently, the target motion state for correction was end-expiration. The value of β represents the maximum expected physiological motion along a given dimension and was set to 8 mm based on a preliminary evaluation (not shown) of the data acquired in this study. Localized gradient entropy was evaluated on the data subset with T2prep for improved contrast. Preliminary results (Supporting Material Figure S6) estimated negligible motion along the slice dimension (presumably due to the 8 mm acquired resolution), and thus only in-plane motion correction was considered in this study.
All acquired data was reconstructed with LRI-HDPROST and solved with ADMM, 47

| Image analysis
Left ventricular myocardial T 1 and T 2 values were evaluated in regions of interest defined by the American Heart Association's (AHA's) 16-segment model. 53 3D cMRF parametric maps were interpolated to 1 × 1 × 4 mm 3 , for a total of 38 slices (2D reference methods were correspondingly interpolated to 1 × 1 mm 2 ) using cubic interpolation. 3D cMRF slices corresponding to the basal, mid-and apical slices in the conventional 2D acquisitions were manually selected for analysis and comparison. The mean and coefficient of variation (CoV) over all subjects were evaluated for each of these slice positions, for T 1 and T 2 . Segments where the myocardium was not well depicted (eg due to remaining cardiac motion, respiratory motion or field inhomogeneities) were excluded from analysis. These included three segments for MOLLI and 20 segments for SASHA (out of 160 segments), primarily in apical and/or inferior regions. No segments were excluded for the proposed 3D cMRF.

| Simulations
No significant variations of the (absolute) error in T 1 and T 2 were found when comparing RR intervals between 600 and 1500 ms (

| Phantom
T 1 and T 2 values for the proposed 3D cMRF in the phantom were found to be in agreement with reference spin echo values. A slight negative bias in short T 2 (<100 ms) and a considerable negative bias for long T 2 (>100) was observed ( Figure 2 Representative T 1 and T 2 maps obtained with 3D cMRF, 2D MOLLI, 2D SASHA and 2D T 2 -GraSE are depicted in Supporting Material Figure S4, revealing comparable quality between 3D cMRF and corresponding conventional methods.

| In vivo
3D cMRF was successfully acquired and reconstructed in all healthy subjects; the cohort's average heartrate was 61 ± 10 bpm; the cohort's acquisition time was 6.9 ± 1.1 min. The proposed 3D cMRF produced maps comparable to those of conventional 2D MOLLI, SASHA and T 2 -GraSE, as seen in representative subject A (Figure 3). Myocardium depiction is similar across different maps, except in apical slices, where MOLLI and especially SASHA present more artefacts (primarily at the inferolateral segment). Whole-heart T 1 and T 2 maps were obtained with 3D cMRF as seen in representative subject B (Figures 4 and 5, respectively). Myocardium and papillary muscles are well defined for all slices in the co-registered T 1 and T 2 maps obtained with 3D cMRF. T 1 measurements with the proposed 3D cMRF were generally higher than those with MOLLI and slightly lower than those with SASHA; 3D cMRF T 2 measurements were generally lower than those with T 2 -GraSE. The mean ± standard deviation of the mean over the subject cohort is depicted in the AHA 16-segment model shown in Figure 6A (as a surrogate for accuracy  Figure 6B (as a surrogate for precision). Corresponding measured CoVS for MOLLI, SASHA and 3D cMRF T 1 values were 7.6 ± 1.6%, 12.1 ± 2.7% and 5.8 ± 0.8%, respectively. For T 2 -GraSE and 3D cMRF the measured T 2 values were 10.5 ± 1.4% and 11.7 ± 1.6%, respectively. Septal T 1 and T 2 values over basal, mid-and apical slices revealed slightly higher slice variations for 3D cMRF relative to MOLLI and  Figure 9 compares the proposed 3D cMRF with and without respiratory motion correction for subject C ($72 bpm, 8.3 mm motion amplitude), where residual blurring of the myocardium and papillary muscles is observed in the absence of translation motion correction. 3D translation motion correction was investigated in a preliminary study; however, no substantial differences were observed (Supplementary Material Figure S6). Key statistical results are summarized in Table 1.

| DISCUSSION
A framework for 3D free-breathing cMRF was developed and evaluated, showing comparable performance with 2D MOLLI, SASHA and T 2 -GraSE mapping techniques. In contrast to conventional 2D approaches, 3D cMRF does not require breath-holds, achieving whole-heart coverage and providing simultaneous and co-registered T 1 and T 2 parameter maps from a single scan. This approach helps streamline myocardial tissue F I G U R E 3 3D cMRF T 1 , 2D MOLLI, 2D SASHA, 3D cMRF T 2 and 2D T 2 -GraSE in three slices for representative subject A (RR interval = 930 ± 49 ms) characterization, potentially improving the experience of the radiographer (fewer scans to plan, higher throughput), the patient (fewer breathholds, greater comfort) and the cardiologist (co-registered whole-heart maps, facilitated and comprehensive analysis).
3D cMRF T 1 values were generally higher than those of MOLLI (mean bias of +38 ms) and lower than those of SASHA (mean bias of −55 ms), which has also been observed in previous 2D cMRF studies. 27,54 There are several confounding factors that may account for these differences, F I G U R E 5 A, Whole-heart T 2 mapping (20 short-axis slices) produced with 3D cMRF for representative subject B (RR interval = 1069 ± 62 ms). B, Corresponding basal, mid-and apical slices for T 2 -GraSE F I G U R E 4 A, Whole-heart T 1 mapping (20 short-axis slices) produced with 3D cMRF for representative subject B (RR interval = 1069 ± 62 ms). B, C, Corresponding basal, mid-and apical slices for MOLLI (B) and SASHA (C) such as T 2 dependences, partial volume and B 1 . 12 In particular, lack of magnetization transfer 55 in the current MRF model may account for the underestimation bias relative to SASHA. Similarly, 3D cMRF T 2 values were generally lower than T 2 -GraSE (mean bias −7.3 ms), also in agreement with previous studies in two dimensions. 56 2D cMRF studies comparing against T2prep-bSSFP and T2prep-FLASH (fast low angle shot) have also reported an underestimation bias. [27][28][29] T2prep methods and GraSE based methods have shown a negative and positive bias relative to multi-echo spin echo, 57 respectively. Underestimation was observed for high T 2 values (not relevant for myocardial characterization), likely due to the absence of long T2prep modules (to encode high T 2 values) in the proposed sequence. T 2 underestimation has been previously reported in brain, 58 abdominal 59 and cardiac 60 MRF, and confounding factors that may account for these differences should be further investigated in the future.
Coefficients of variation were slightly lower for 3D cMRF T 1 than for MOLLI (and more so than for SASHA); for 3D cMRF T 2 these values were slightly higher (within one standard deviation) than for T 2 -GraSE. In all cases CoVs were measured on cubic interpolated data: the smoothness effect of the interpolation may lead to small reductions in CoV; however, the effect is consistent for all methods. A general decrease of parameter values from the septal to the lateral wall was observed with 2D MOLLI, T 2 -GraSE and more so with 3D cMRF for both T 1 and T 2 . Some slice variabilities were also observed with 3D cMRF; these were slightly higher than with MOLLI and T 2 -GRASE. These variations are likely due to B 0 and B 1 field inhomogeneities not accounted for in the respective models.
The proposed 3D cMRF achieves T 1 and T 2 encoding with IR pulses and T 2 preparation pulses, similar to the original 2D cMRF. 27 Here, sequence design was made heuristically, guided by results from previous studies (eg low flip angles to reduce B 1 sensitivity) in 2D cMRF. 28 F I G U R E 6 16-segment AHA bull'seye plot for the mean and the CoV for 2D MOLLI, 2D SASHA, 2D T 2 -GraSE and 3D cMRF T 1 and T 2 Preliminary simulations in this work suggest that the flip angle pattern could affect precision at low flip angles and/or low SNR; however, further studies comparing variable and fixed flip angle patterns should be investigated. The proposed sequence differs from previous cMRF 61 by separating IR and T2prep pulses (instead of interleaving them). This approach is expected to improve the dynamic range of the T 1 encoding at the expense of some SNR in T 2 encoding; however, it remains a heuristic solution. Comprehensive analysis on optimal sequence design for conventional (nontriggered) MRF has been developed previously. [62][63][64] A similar analysis for (triggered) cMRF, possibly considering B 0 /B 1 sensitivities, would be of interest and will be investigated in future studies. Optimal sequences would be expected to lead to reduced scan times and/or improved parametric maps. Fat signal can create strong undersampling artefacts and has been shown to bias T 1 estimation, 65 in addition to introducing blurring artefacts due to spiral readouts (no spiral blurring was considered in the proposed framework). Moreover, translational motion correction can introduce artefacts originating from (incorrectly motion corrected) static fat tissue. 66 A SPIR pulse is used in this framework for fat suppression, but alternative approaches such as water selective excitations 67 could be considered for similar effect. Corrections using more complex respiratory motion models can also help in this regard. 66 This study has several limitations. 3D cMRF has a predictable scan time of $7 min, and, while adequate for clinical deployment, future work should consider additional strategies to reduce scan time. Higher resolution will also be desirable in the future and alternative 3D non-Cartesian trajectories could be considered, as recently proposed for T 1 /T 2 mapping. 68 Although 3D cMRF was compared against three conventional approaches, it was not compared against T2prep-bSSFP at this time, due to unavailability on our system. As a proof of concept, no patients were included in this study, but should be considered in future studies. No in vivo repeatability was investigated in this study, primarily due to the long scan times associated with the experiment. Some residual blurring was observed with 3D cMRF, which could be primarily due to spiral (off-resonant) blurring and residual respiratory or cardiac (due to acquisition window or mistriggering) motion. Residual cardiac motion within the acquisition window of 204 ms could potentially lead to artificially enlarged wall thickness or reduced T 2 values 69 ; residual respiratory motion could occur due to errors in the bellows measurements (eg drift). 70 Minor blurring can also arise as an artefact of the undersampled reconstructions performed. SPIR modules will likely decrease the method's sensitivity to detect fat infiltration; however, 3D cMRF could potentially be combined with Dixon based readouts for simultaneous T 1 , T 2 and fat fraction estimation. 56 The proposed approach relies on a respiratory bellows driven autofocus algorithm to estimate the beat-to-beat translational motion of the heart. However, residual motion artefacts could also occur due to F I G U R E 7 A, T 1 values in the septum in each healthy subject (labelled from S1 to S10) for MOLLI, SASHA and 3D cMRF T 1 . B, Corresponding T 2 values for T 2 -GraSE and 3D cMRF T 2 . Reported values are averaged over basal, mid-and apical slices; error bars correspond to inter-slice standard deviations the simplified (translational) motion model. Moreover, as the slice thickness is reduced, through-plane motion correction will increasingly be more important to avoid errors in T 1 and T 2 . 71,72 Future work will consider generalized elastic motion compensated reconstructions to correct the remaining non-rigid components of motion. 73 A major limitation of the current approach, however, is the reconstruction time of approximately 10 h. 3D cMRF requires subject specific dictionaries and long iterative reconstructions. While patient tailored dictionaries are expected to grant higher heart rate insensitivity, they are also associated with long computation times ($3 h in this study) which can impose limitations on dictionary resolution. A GPU implementation of the method will considerably cut computational times. 74 In addition, recent developments in deep learning indicate that dictionary computation, 75 image reconstruction 76 and template matching 77 computational demands could be reduced by orders of magnitude, all of which could be combined with 3D cMRF to considerably reduce computational times. With the current approach, 3D cMRF is limited to simultaneous T 1 /T 2 and has some residual biases in these measurements. cMRF can potentially be extended to map any combination of parameters and consider any combination of model corrections; however, dictionaries increase exponentially (in addition to increased acquisition and reconstruction times). 78 Deep learning solutions could also play a role in enabling acceptable computational times for truly multi-parametric fingerprinting (eg T 1 , T 1ρ , T 2 , T 2 *, B 0 , B 1 , MT and diffusion) from a single scan and should be investigated as future work.

| CONCLUSION
A novel approach for 3D free-breathing respiratory motion compensated cMRF was proposed and validated in healthy subjects. The proposed method produces co-registered T 1 and T 2 maps with whole-heart coverage in a predictable scan time of $7 min. Free-breathing 3D cMRF showed similar accuracy and precision to conventional 2D breath-hold methods; however, a slight underestimation of T 1 and T 2 values with cMRF was observed when compared with 2D SASHA and T2-GraSE, respectively. 2D MOLLI, SASHA and T 2 -GRASE for representative subject C (RR interval = 870 ± 55 ms). Blurring and ghosting motion artefacts in T 1 and T 2 maps are reduced after translation correction (arrows)