Free‐running 3D whole heart myocardial T1 mapping with isotropic spatial resolution

Purpose To develop a free‐running (free‐breathing, retrospective cardiac gating) 3D myocardial T1 mapping with isotropic spatial resolution. Methods The free‐running sequence is inversion recovery (IR)‐prepared followed by continuous 3D golden angle radial data acquisition. 1D respiratory motion signal is extracted from the k‐space center of all spokes and used to bin the k‐space data into different respiratory states, enabling estimation and correction of 3D translational respiratory motion, whereas cardiac motion is recorded using electrocardiography and synchronized with data acquisition. 3D translational respiratory motion compensated T1 maps at diastole and systole were generated with 1.5 mm isotropic spatial resolution with low‐rank inversion and high‐dimensionality patch‐based undersampled reconstruction. The technique was validated against conventional methods in phantom and 9 healthy subjects. Results Phantom results demonstrated good agreement (R2 = 0.99) of T1 estimation with reference method. Homogeneous systolic and diastolic 3D T1 maps were reconstructed from the proposed technique. Diastolic septal T1 estimated with the proposed method (1140 ± 36 ms) was comparable to the saturation recovery single‐shot acquisition (SASHA) sequence (1153 ± 49 ms), but was higher than the modified Look‐Locker inversion recovery (MOLLI) sequence (1037 ± 33 ms). Precision of the proposed method (42 ± 8 ms) was comparable to MOLLI (41 ± 7 ms) and improved with respect to SASHA (87 ± 19 ms). Conclusions The proposed free‐running whole heart T1 mapping method allows for reconstruction of isotropic resolution 3D T1 maps at different cardiac phases, serving as a promising tool for whole heart myocardial tissue characterization.


| INTRODUCTION
Cardiovascular magnetic resonance has been increasingly used to diagnose and monitor different cardiac diseases. 1,2 Quantitative mapping of tissue parameters has demonstrated the potential to characterize subtle pathological changes of the myocardium and is especially valuable in detecting both focal and diffuse fibrosis. 3,4 T 1 is sensitive to the changes of myocardial tissue properties, such as fibrosis, fat, and water content, and therefore is considered a promising novel biomarker to characterize myocardium and assess various cardiomyopathies. [5][6][7][8] In addition, extracellular volume fraction, which is important to detect diffuse myocardial fibrosis, can be estimated from native and post-contrast myocardial T 1 mapping. 9 Myocardial T 1 mapping has been recognized as one of the most valuable quantitative mapping techniques to support diagnostic, therapeutic, and prognostic decision making in ischemic and non-ischemic cardiomyopathies. 4 The modified Look-Locker inversion recovery (MOLLI) 10 and saturation recovery single-shot acquisition (SASHA) 11 are commonly used cardiac T 1 mapping techniques. They acquire multiple images with different T 1 contrasts and estimate pixel-wise T 1 by fitting the acquired signal to an exponential recovery model. Respiratory and cardiac motion may result in blurring artefacts in each T 1 -weighted image and can cause misalignment between different images, influencing the T 1 mapping accuracy. 12 Therefore, breath-holding or respiratory gating and electrocardiography (ECG) triggering are usually required during cardiac imaging. However, besides the requirement of patient compliance, a single breath-hold limits acquisition time and results in low spatial resolution and limited coverage as typically only 1 image is acquired per cardiac cycle. 10,11 Multiple breath-holds may allow longer acquisition time, but may result in different breath-hold positions, which can lead to mis-registration artefacts and inaccuracy of T 1 estimation. 13 Respiratory gating has been used to limit respiratory motion during free-breathing imaging, but has low imaging efficiency and may lead to unpredictable long acquisition time. 14 ECG triggering is usually used to gate the acquisition only at diastole in each cardiac cycle, which has low acquisition efficiency and can be challenging in arrhythmic patients who have shortened diastolic phase.
Recently, the magnetic resonance multitasking technique using 2D radial acquisition has been proposed for motionresolved quantitative myocardial T 1 and T 2 imaging without breath-holding or ECG triggering. 15,16 2D myocardial mapping is typically performed with thick slices to avoid the influence of through-plane motion 13 and low SNR. 3D acquisitions can overcome these limitations and provide whole heart coverage for comprehensive characterization of diffuse diseases of myocardium. 17 However, 3D T 1 mapping is technically challenging, because respiratory gating and ECG-triggering together with the requirement of whole heart coverage, may lead to very long scan time. To cover the whole heart in an acceptable scan duration, current 3D T 1 mapping techniques 14,18 acquire only a few slices in the shortaxis view with a slice thickness of >8 mm, which require complex scan planning and may suffer from partial volume effect along the slice direction.
To address the current technical challenges of myocardial T 1 mapping, we propose a free-running (free-breathing, retrospective cardiac gating) 3D whole heart T 1 mapping technique with high isotropic spatial resolution (1.5 mm 3 ), which is able to provide images in any desired orientation. The proposed sequence consists of an inversion recovery (IR)-prepulse, followed by an efficient 3D golden angle radial acquisition 19 to sample the T 1 relaxation curve. The k-space center is repeatedly sampled within each spoke and used to estimate 1D respiratory motion. ECG signal is recorded and synchronized with data acquisition. Respiratory binning is performed based on the estimated 1D respiratory signal, enabling estimation and correction of 3D translational respiratory motion. For the reconstruction, a recently proposed algorithm that combines dictionary-based low-rank inversion 20 and highdimensionality 3D patch-based undersampled reconstruction (HD-PROST) 21 is used. The proposed approach allows for respiratory motion compensated 3D myocardial T 1 mapping at different cardiac phases. The feasibility, accuracy and precision of the proposed technique was validated against conventional T 1 mapping methods in a standardized T 1 phantom and 9 healthy subjects.

| Free-running 3D T 1 mapping sequence
The proposed IR-prepared 3D radial sequence is shown in Figure 1A. IR-preparation is used to generate T 1 contrast, after which a spoiled gradient-echo (SPGR) readout is performed with small flip angle (θ). Spatial-spectral pulse for selective water excitation is adopted for fat suppression. To achieve pseudo-uniform distribution of radial spokes for retrospective cardiac and respiratory binning, the continuously acquired radial spokes conform to the 3D golden angle distribution. 19 Specifically, the azimuthal angle and polar angle increments (Δυ, Δτ) between adjacent spokes are defined as: where ϕ 1 = 0.4656, ϕ 2 = 0.6823 are the 2D golden angle means. 19 (1) Δ = acos 1 , Δ = 2 ⋅ 2 , | 1333 QI et al.

| Motion extraction and retrospective sorting
The flowchart of retrospective processing, including motion extraction, respiratory motion correction and data sorting for T 1 mapping of a given cardiac phase is shown in Figure 1B. Independent component analysis (ICA) was performed on the k-space center of all spokes from all coils to extract the signal component that is most relevant to the 1D superior-inferior respiratory motion. This approach has been previously validated for cardiac cine and coronary artery MRI using either radial or spiral trajectory. 22,23 To extract the self-navigated respiratory signal from the k-space center, coil compression is first performed using principle components analysis. ICA is applied to the k-space center amplitudes from the compressed coils to obtain 5 independent components (IC). Example of the extracted ICs and the corresponding spectral power are shown in Figure 2A and B. The harmonic corresponding to the contrast change because of IR preparation is removed from the ICs based on the IR repetition time (IRTR in Figure 1A). Then, the IC with the highest spectral power in the respiratory frequency range (0.1-0.5 Hz) is selected (IC 4 in Figure 2A) and band-pass filtered to get the 1D respiratory signal, based on which the k-space data is binned into 5 equally populated respiratory phases from end-inspiration to end-expiration. An example of the estimated respiratory signal validated against the corresponding respiratory bellow signal is shown in Figure 2C. Cardiac motion synchronization is achieved by logging the ECG time stamps, from which the temporal occurrence of each spoke relative to the most recent cardiac R-wave preceding the spoke is calculated and termed as cardiac delay.
After retrospective respiratory binning and cardiac gating, 3D translational respiratory motion is estimated via registration from the intermediate reconstruction of low-resolution respiratory bin images at diastolic cardiac phase by using central part of selected k-space data. Representative diastolic bin images of 5 respiratory motion states and the difference images with the end-expiratory bin are shown in Figure 2D. A rectangular volume of interest for motion estimation is selected around the heart at the end-expiratory bin and propagated to other respiratory bins. Respiratory motion correction is performed by correcting the phase of the k-space data using the estimated 3D translational motion parameters. The respiratory motion corrected k-space data can then be binned according to the cardiac delay and inversion recovery time (TI) for T 1 mapping at a given cardiac phase. preparation, spoiled gradient echo (SPGR) readout with low flip angle (θ) was performed using 3D golden angle radial trajectory. T gap is the time between the IR pulse and the first excitation and T ex is the time interval between the last excitation in the SPGR readout and the next IR pulse. IRTR is the IR repetition time. (B) Data sorting process for reconstruction of multiple T 1 contrasts for a given cardiac phase, which includes 3 steps: 1D respiratory motion estimation from k-space center of all radial spokes and cardiac motion extraction from ECG log; respiratory motion correction of k-space using motion parameters estimated by 3D translational image registration of respiratory bin images at diastole; binning the respiratory motion corrected k-space into different T 1 contrasts for a given cardiac phase

| HD-PROST algorithm
For 3D T 1 mapping at a given cardiac phase, spokes for the corresponding cardiac delay and acquisition window are selected, which are then sorted into a time series of T 1 contrasts according to TI and a given temporal window for each T 1 frame. After data selection, the highly undersampled 3D radial data is reconstructed over multiple T 1 contrasts by combining a dictionary-based low-rank inversion 20 to efficiently reduce the number of T 1 contrasts to be reconstructed with a recently proposed high-dimensionality 3D patch-based undersampled reconstruction (HD-PROST), exploiting local (within a patch), non-local (between similar patches), and contrast redundancies. 21 The reconstruction for the proposed technique can be formulated as the following unconstrained Lagrangian 21 where ||·|| F and ||·|| * denote the Frobenius norm and nuclear norm respectively; I is the compressed T 1 image series to be reconstructed; E = √ DAU r FS is the encoding operator, with S being sensitivity maps, F being Fourier transform, U r being the low-rank operator obtained by truncating the singular value decomposition (SVD) of a dictionary generated by Bloch simulation, 20 A being the convolutional gridding operator, transforming Cartesian data back to 3D radial, and D being the non-Cartesian density compensation function; K is the undersampled data; P p (·) is the patch selection operator at pixel p of a 3D multi-contrast image set. This operator selects patches on local (patch for a given pixel location and contrast), non-local (similar patches within a neighborhood for a given contrast) and contrast (patches from all the contrasts) scales and  p is a 3D tensor built by the selected patches centered at pixel p; T represents the denoised multi-contrast images constructed by folding and aggregating the tensors  p for each pixel p (see the step 2 below); Y is the Augmented Lagrangian multiplier; λ is the sparsity-promoting regularization parameter and µ is the penalty parameter. Equation 2 can be efficiently solved by operator-splitting via alternating direction method of multipliers (ADMM), which divides the optimization process into the following 3 steps: Step 1: joint reconstruction update The first step is a joint reconstruction of the compressed TI image series (i.e., singular images) I by incorporating the denoised multi-contrast images T obtained at the end of step 2 as prior information The above equation can be efficiently solved using the conjugate gradient (CG) algorithm.
Step 2: high-order singular value decomposition based denoising With obtained I and Y from step 1 and step 3 separately, the second step is to minimize the following equation regarding to the high order tensor  p The details to solve the above equation can be found in Bustin et al. 21 Generally speaking,  p is obtained by thresholding the singular values obtained via high-order singular value decomposition of the 3D tensor built by the patches selected from the multi-contrast images. The thresholding parameter is defined by . The denoised  p is then rearranged to form the denoised image patches. This process is repeated for all the pixels in the singular images, and aggregation is then performed to generate the final denoised singular images T.
Step 3: Lagrangian multiplier update Finally, with the optimized I and T from step 1 and step 2, the Lagrangian multiplier is updated by Y = Y + I − T. The above 3 steps are iteratively interleaved in the reconstruction process to improve the reconstructed image quality.
For the retrospective data selection for a given cardiac phase, the cardiac acquisition window was set to 186 ms, similar to the acquisition window of conventional 2D MOLLI and SASHA imaging. 10,11 The temporal window per T 1 frame was also set to 186 ms, resulting in 10 T 1 frames to be reconstructed. A dictionary-based low-rank compression 20 was performed along the T 1 contrast dimension to further reduce the number of T 1 frames to be reconstructed (Equation 2). The dictionary was generated using Bloch simulation for T 1 in the range of 100-3000 ms, with an increment of 1% with respect to the previous T 1 . The contrast of a specific TI frame is given by the average of all the radial spokes binned into the corresponding frame. The signal of each radial spoke is calculated using Bloch equation, which has been previously described for 3D radial-based carotid quantitative mapping 24,25 and is also provided in the Supporting Information file. The singular values by SVD of the dictionary are shown in Figure 3A. The low-rank operator U r in Equation 2 was obtained by keeping the 3 largest singular value vectors, resulting in 3 singular images to be reconstructed ( Figure 3B). The first singular image, with the highest SNR, was used for patch selection in the step 2 of HD-PROST (Equation 4).
In this study, the maximum number of CG iterations in the joint MR reconstruction (step 1) was set to 15, and the regularization parameter µ, balancing the contribution of the prior term, was empirically set to 0.01. Coil sensitivity maps were estimated by using spokes of TI larger than 1000 ms and reconstructed in low resolution by using only the central half of k-space to reduce noise. The regularization parameter λ of HD-PROST was optimized by visual inspection of the reconstruction quality on 1 subject and used for all other data set. For the patch-based denoising step in HD-PROST, 21 the following parameters were used: patch size = 5 × 5 × 5; search window = 10 × 10 × 10; patch offset = 3; number of selected similar patches = 15. The number of ADMM iterations was set to 5. Example singular images after HD-PROST reconstruction are shown in Figure 3B. T 1 maps were generated by a dot product matching between the reconstructed singular images and the dictionary.
For the proposed 3D radial T 1 mapping sequence, the ECG signal used for the retrospective binning was generated by physiological simulation with heart rate varying between 50 bpm and 70 bpm during the acquisition. The acquired 3D radial data was retrospectively binned into 10 T 1 contrasts for a simulated diastolic phase with cardiac delay of 650 ms. The phantom T 1 map acquired with the 3D imaging sequence was generated using the HD-PROST reconstruction as explained above.

| In vivo experiments
The study was approved by the local institutional review board and 9 healthy subjects (5 males, 31.8 ± 3.5 y) were imaged using a 28-channel cardiac coil. All volunteers provided written informed consent before inclusion into this study. First, standard breath-hold 2D cine images were acquired in the mid ventricular short-axis view with retrospective gating to 16 cardiac phases, based on which the cardiac delays for systole and diastole were determined. The T 1 mapping reference consisted of the clinically adopted breath-hold 2D MOLLI (3-3-5) 27 and 2D SASHA 11 sequences that were performed in diastole. The imaging parameters of MOLLI were TR/TE = 2.6/1.3 ms, flip angle = 35°, FOV = 288 × 288 mm 2 , in-plane resolution = 2 × 2 mm 2 , and slice thickness = 8.0 mm. The SASHA sequence was performed with TR/TE = 2.6/1.3 ms, flip angle = 70°, and the same FOV and spatial resolution as used for MOLLI. Parallel imaging with SENSE acceleration factor of 2 in the phase-encoding direction was used in the 2D MOLLI and SASHA mapping acquisitions, resulting in a mid-diastolic acquisition window of ~187 ms. MOLLI and SASHA images were acquired in the same shortaxis location as the cine scan and reconstructed to an in-plane resolution of 1.5 × 1.5 mm 2 . The 2D T 1 maps for MOLLI and SASHA were obtained by fitting the corresponding standard 3-parameter models. 10,11 The proposed free-running 3D T 1 mapping sequence was also acquired in short-axis orientation covering the entire heart and centered around the slice location of the 2D MOLLI and SASHA scans. Besides diastole, T 1 mapping was also reconstructed for systole using the proposed framework, considering that systolic T 1 maps may find applications in patients with thin myocardium or arrhythmias. 28,29

| Image analysis
For phantom data, circular region-of-interests (ROIs) were drawn for each vial and the mean and SD of T 1 values were calculated for each vial. The accuracy of the proposed 3D T 1 mapping technique was evaluated by comparing the estimated T 1 values to the 2D IR-SE method using Pearson's linear correlation and Bland-Altman analyses.
For the analysis of in vivo images, the diastolic T 1 map in the mid short-axis view that was most similar to the MOLLI and SASHA geometry was selected from the 3D T 1 maps of the proposed approach for each subject. Septal ROIs with the same size were drawn for the 3 T 1 mapping methods. The mean and SD of T 1 values in the septum were measured and compared between the proposed 3D method and 2D MOLLI and SASHA separately using the Wilcoxon rank-sum test with Bonferroni correction to evaluate the in vivo accuracy and precision.
For the analysis of 3D in vivo T 1 maps, 3 slices from base, mid to apex were selected from the diastolic and systolic T 1 maps. T 1 values of the myocardium were measured according to the myocardial segments model defined by the American Heart Association (AHA), 30 with 6 segments in the base and middle slices and 4 segments in the apical slice. The mean T 1 values were calculated for all segments to evaluate the spatial T 1 distribution of the proposed method. The T 1 value in each segment was compared between diastole and systole using the paired t-test. Then the diastolic and systolic T 1 values for each segment were determined by averaging across all the volunteers and visualized with bull'seye plots.
All image reconstruction and analysis were performed using MATLAB (The MathWorks, Natick, MA) on a server with a dual 16-core CPU and 256 GB RAM. Statistical analysis was carried out in GraphPad Prism 7 (La Jolla, CA). A P-value of <0.05 was considered statistically significant.

| Phantom study
Phantom T 1 maps with the proposed approach and 2D IR-SE are shown in Figure 4A and B, and the mean and SD of T 1 values of all vials are shown in Figure 4C. The linear correlation indicates good agreement of T 1 estimations between the proposed 3D technique and the 2D IR-SE reference (R = 0.99). Bland-Altman plot demonstrates no correlation between the difference and average of T 1 measurements of the 2 methods ( Figure 4D). Compared with 2D IR-SE, the percentage error of the T 1 measurements with the proposed approach was 1.2 ± 0.7%, ranging from 0.3-2.1%.

| In vivo study
All volunteers completed the scans, and the diastolic and systolic 3D T 1 maps were successfully reconstructed for each subject with recorded heart rates of 62 ± 8 bpm, ranging from 45-73 bpm. For reconstruction of the 3D T 1 map, each ADMM iteration in HD-PROST takes ~34.9 min consisting of 26.3 min for the joint reconstruction step and 8.6 min for the patch-based denoising step, resulting in a total reconstruction time of 174.5 min with 5 ADMM iterations. T 1 mapping results at diastole from a representative healthy subject are shown in Figure 5, consisting of 8 representative short-axis slices and a reformatted long-axis slice. As can be seen, the spatial distribution of myocardium T 1 was uniform over the whole left ventricle.

F I G U R E 4 (A and B) Phantom T 1 maps of 2D inversion recovery spin echo (2D IR-SE) and the proposed 3D T 1 mapping approach.
(C) Linear correlation of phantom T 1 estimation with the proposed 3D method in comparison with reference 2D IR-SE sequence. (D) Bland-Altman plot showing the difference between the 2 methods and their average. The grey solid line indicates the mean difference, and the grey dotted lines indicate the 95% confidence intervals of limits of agreement F I G U R E 5 Proposed 3D T 1 mapping technique at diastolic cardiac phase for a representative healthy subject. Eight short-axis slices, from apex to base of the left ventricle and the reformatted longaxis view are shown. Uniform T 1 distribution across the myocardium can be observed Diastolic T 1 mapping results from another 2 subjects are shown in Figure 6, including the mid short-axis and longaxis slices from the 3D T 1 mapping approach and a middle short-axis slice from 2D MOLLI and SASHA techniques. Diastolic T 1 maps reconstructed by the proposed technique were visually comparable to the reference 2D methods with improved depiction of the right ventricle. The mean and SD of septal T 1 derived from 2D MOLLI, 2D SASHA, and the proposed 3D technique are compared in Figure 7 for all subjects. The mean septal T 1 measured with the proposed approach is 1140 ± 36 ms and is comparable to that of SASHA (1153 ± 49 ms) (P = 0.68), but is higher than that of MOLLI (1037 ± 33 ms) ( P < 0.01). The SD, indicating the precision of septal T 1 measurements, was 42 ± 8 ms with the proposed approach, which is much lower than the SD of SASHA (87 ± 19 ms) (P < 0.01) and similar to that of MOLLI (41 ± 7 ms) (P = 0.79).
Representative systolic and diastolic T 1 maps from 2 subjects are shown in Figure 8, where 3 short-axis slices (base, mid and apex) are included. The systolic T 1 maps also demonstrated a uniform T 1 distribution with thicker myocardium compared with diastole. The diastolic and systolic T 1 values of each AHA segment are shown in the bull's-eye and box plots in Figure 9. The comparison of T 1 values in diastole and systole for all myocardium segments are summarized in Table 1. No significant differences were found between the diastolic and systolic T 1 values estimated with the proposed technique in most of the AHA segments, except for the inferior segments in the apical and mid ventricular slices, and the inferoseptal segment in the basal slice, with significantly longer T 1 values in diastole than in systole.

| DISCUSSION
In this study, free-running 3D whole heart myocardial T 1 mapping with 1.5 mm 3 isotropic spatial resolution has been successfully achieved with predictable 9.5 min scan time.
The proposed sequence features free-breathing, retrospective cardiac gating and continuous data acquisition, so T 1 maps can be compensated for 3D translational respiratory motion and reconstructed at different cardiac phases by retrospective data binning. Using combined dictionary-based low-rank inversion 20 and high-dimensionality 3D patch-based undersampled reconstruction algorithm (HD-PROST) 21 to reconstruct the highly undersampled 3D radial data, feasibility of myocardial T 1 mapping was demonstrated in vivo. Myocardial T 1 measurements with the proposed 3D technique were comparable to conventional breath-hold, cardiac-gated 2D MOLLI and SASHA. In the phantom study, T 1 values estimated with the proposed technique showed good accuracy for T 1 ranging from 250-1600 ms as compared with the reference 2D IR-SE. Septal T 1 values at diastole with the proposed technique were comparable to SASHA and were significantly larger than those for MOLLI, which has been reported to underestimate T 1 . 12 In this study, the dictionary was simulated with a linear T 1 increment of 1%, resulting in 10-12 ms increments in the myocardial T 1 range of 1000-1200 ms, which are small compared with the measured septal T 1 SD (42 ± 8 ms). The precision of septal T 1 was similar to MOLLI and significantly higher than SASHA. The results indicate that the proposed 3D T 1 mapping technique improves accuracy compared to MOLLI and precision compared to SASHA. Therefore, the proposed technique may be a promising tool for myocardial T 1 mapping with the advantages of both high accuracy and precision. In addition, with the proposed 3D high resolution T 1 mapping technique, visualization of the right ventricle was improved, which may facilitate fibrosis quantification in the thin wall of the right ventricle. Nevertheless, T 1 values of the thin right ventricle should be interpreted carefully, considering that there may be partial volume effects resulting from residual fat caused by the imperfection of the water selective excitation pulse.
In this study, respiratory motion was addressed by correcting the phase of k-space data using motion parameters estimated from 3D translational image registration of images acquired at different respiratory states. More sophisticated registration methods, such as affine and non-linear transforms could be considered in future work to account for more complex respiratory-induced cardiac movement. To extract reliable cardiac trigger signal to account for cardiac motion, ECG signal was logged and synchronized with data acquisition for retrospective data selection. Previous studies have proposed several approaches for cardiac self-navigation using radial or spiral trajectories. 22,23 Therefore, cardiac selfnavigation, besides respiratory self-navigation, could be also investigated in future studies to achieve free-breathing, ECGfree 3D myocardial T 1 mapping.
T 1 mapping at systole has special diagnostic value when the diastolic T 1 map cannot be reliably obtained, such as in patients with thin myocardium, where partial volume effects may corrupt myocardial T 1 estimation, and in patients suffering from frequent arrhythmias, where the diastolic onset and duration are changing. 28,29 In this study, systolic T 1 maps were reconstructed along with diastolic T 1 maps for each subject. The AHA segment analysis showed good spatial uniformity of the myocardial T 1 values measured across the left ventricle in both diastole and systole. In general, the systolic myocardium T 1 was comparable to the diastolic myocardium T 1 . The shorter T 1 in the inferoseptal segment of the basal slice and inferior segments of the mid and apical slices may be explained by increased myocardial thickness and therefore reduced partial volume effects from blood during systole compared to diastole. 29 The acquisition window for T 1 mapping was set to ~186 ms, which was adequate for diastole, but may not be optimal for systole. In spite of this, radial data acquisition is less sensitive to motion than Cartesian sampling, and no obvious motion artefacts were observed in the systolic T 1 maps.
Systolic and diastolic T 1 maps were reconstructed for each subject separately in our experiments. Instead of separate reconstruction for each cardiac phase, cardiac-resolved T 1 maps can be reconstructed simultaneously, by which additional redundancy along cardiac phases can be exploited and improved reconstruction performance can be expected at the expenses of memory requirement and computational times. Furthermore, 3D cardiac cine reconstruction from the freerunning sequence should be also possible. This could be done by sorting the respiration motion compensated k-space data with a small temporal window along the cardiac phase (e.g., 50 ms) and a large temporal window along contrast. After binning, the undersampled k-space data of different cardiac phases could be reconstructed simultaneously, e.g., by adding a total variation constraint along the cardiac phase direction 31 to the reconstruction Equation 2. This approach will be investigated in future studies. F I G U R E 8 Representative diastolic and systolic T 1 maps of basal, mid and apical short-axis views using the proposed 3D T 1 mapping sequence from 2 healthy subjects (A, B). Uniform myocardial T 1 distribution can be observed on the T 1 maps, both in diastole and systole. SPGR readout with low flip angle instead of balanced steady-state free precession (bSSFP) readout was used in the proposed sequence. Although bSSFP readout may have SNR benefits, the low flip angle SPGR readout is less sensitive to B 1 and B 0 field inhomogeneities that will influence T 1 estimation and may cause banding artefacts in the image. The proposed technique was investigated at 1.5 T in this study, but has potential to be extended to 3 T. Furthermore, benefits of higher SNR, and therefore shortened scan time or improved spatial resolution can be expected at 3 T.
Several limitations of the proposed technique need to be mentioned. First, the proposed T 1 mapping technique cannot estimate T 1 of flowing blood accurately. This is because the signal model used for T 1 mapping assumes that the tissue of interest is static and experiences all the inversion and excitation pulses in the sequence, which may not be true for flowing blood. When measurement of myocardial extracellular volume is needed, estimation of blood T 1 could be performed using a free-breathing rapid low spatial resolution (enough to only define a region of interest within the blood pool) T 1 mapping sequence acquired in a mid-ventricular slice. Second, in this study the proposed technique was only tested for native myocardial T 1 mapping. However, post-contrast myocardial T 1 mapping using this technique should be also feasible. Furthermore, considering the reduced T 1 after contrast injection, shorter IRTR could be adopted and shorter scan duration can be expected. Third, in the free-running sequence with retrospective binning, although inter-bin respiratory and cardiac motion has been carefully addressed, compared with breath-hold images, there may still be some remaining intra-bin respiratory motion and cardiac motion that result in some blurring of the reconstructed images. Future studies will investigate the performance of the proposed sequence in a cohort of patients with cardiovascular disease. Last, the current implementation of HD-PROST reconstruction is suboptimal and it takes ~3 h to reconstruct a whole heart 3D T 1 map. GPU implementation will be investigated in the future for the non-uniform Fourier Transform 32 and the patch-based denoising process to accelerate the reconstruction.

| CONCLUSIONS
In this study, a free-running 3D myocardial T 1 mapping technique with whole heart coverage and high isotropic spatial resolution is proposed. 3D T 1 mapping was shown to have good accuracy and precision in comparison to reference sequence in phantom and conventional T 1 mapping methods in in vivo experiments. Based on retrospective data selection and combined dictionary-based low-rank inversion and patch-based reconstruction, high resolution 3D T 1 maps could be obtained for different cardiac phases. Future studies will investigate the clinical value of the proposed technique.