We describe the T1 mapping sequence parameters, T1 map calculation methods, and their implementations and evaluations. Figure 1 shows a custom user interface tool for the study, which is available at https://sites.google.com/site/yoonckim1/software/t1_map_compare. Source code is available at https://github.com/prime52/fT1fit.
Data acquisition
Cardiac MRI scans were performed on a 1.5 T scanner (Siemens Avanto, Erlangen, Germany). Clinical MR examinations were approved by our institutional review board, and informed consent was obtained from the subjects prior to MRI scans. Subjects with suspected cardiovascular diseases were enrolled between May 2014 and May 2015. Subjects who took cardiac pre- and post-contrast T1 mapping scan exams were included for the present study. Subjects with inadequate image quality due to severe motion were excluded. A total of 30 subjects were considered for our study. They consisted of 9 HCM patients, 10 cardiac amyloidosis patients, 7 coronary artery disease patients, and 4 healthy volunteers.
The modified Look-Locker inversion recovery (MOLLI) sequence was used for cardiac T1 mapping [19, 20]. Imaging parameters were slice thickness = 8 mm, echo time (TE) = 1.01 ms, spacing between slices = 20 mm, the number of phase encoding steps = 104, pixel bandwidth = 1085 Hz, acquisition matrix = 192 × 120, image matrix = 384 × 288, pixel spacing = 0.9375 mm × 0.9375 mm, and field of view (FOV) = 360 mm × 270 mm. The MOLLI protocols used were different for the pre-contrast and post-contrast T1 mapping. The MOLLI 5(3)3 protocol used for pre-contrast T1 mapping consisted of 5 inversion time (TI) image acquisitions after the first inversion pulse, a three-heartbeat pause for the recovery of the longitudinal magnetization, and 3 TI image acquisitions after the second inversion pulse. The MOLLI 4(1)3(1)2 protocol used for post-contrast T1 mapping consisted of 4 TI image acquisitions after the first inversion pulse, a one-heartbeat pause, 3 TI image acquisitions after the second inversion, a one-heartbeat pause, and 2 TI image acquisitions after the third inversion. Since the post-contrast T1 relaxation time is approximately less than 500 ms, which is much shorter than the pre-contrast T1 (~ 950 ms for myocardium and ~ 1500 ms for blood), the 3 inversion pulses used in the 4(1)3(1)2 protocol enables more adequate sampling of the early part of T1 relaxation than the 2 inversion pulses used in the 5(3)3 protocol. Instead of 5 images after the first inversion, 4 images were acquired in the post-contrast T1 mapping because for the short T1 (< 500 ms) relaxation the image from the fifth heartbeat appears similar to the image from the fourth heart beat due to fast T1 recovery and thus is unnecessary [21]. The TI images were acquired in a diastolic cardiac phase. All of the TI images were aligned using the motion correction algorithm [22]. The TI images were exported as digital imaging and communications in medicine (DICOM) files for T1 map calculation.
T1 Map calculation
T1 map calculation was performed pixel-by-pixel. In general, the TI images can be available as either complex-valued or magnitude-valued. In the present study, the DICOM dataset was available as magnitude TI images, and most cardiac MRI scanner systems in hospitals save the DICOM dataset in the magnitude scale by default. At each voxel’s location (x, y), the signal intensity \(S(t)\) can be modeled as the following T1 relaxation curve.
$$S\left(t\right)=a\left(1-{e}^{-bt}\right)+c$$
(1)
For a set of inversion times t = [\({TI}_{1}\), \({TI}_{2}\), …, \({TI}_{N}\)] and a set of corresponding signals [\({S(TI}_{1}),{S(TI}_{2}), \dots , {S(TI}_{N})\)], there are N equations and 3 unknowns, which are \(a\), b, and c. In the present study, N was 8 for the pre-contrast data, while N was 9 for the post-contrast data. Notably, b is the reciprocal of T1* (i.e., the apparent T1). The objective function \(f\) to be minimized was then given as the sum of squares of the difference between the relaxation model and the data.
$$f(a,b,c)=\sum_{i=1}^{N}{\left(S\left({TI}_{i}\right)-a(1-{e}^{-b{TI}_{i}}\right)-c)}^{2}$$
(2)
This is a nonlinear least squares problem, and the Levenberg–Marquardt (LM) algorithm [23] can be used to iteratively estimate the parameters \(a\), b, and c. We empirically selected the initial values for the parameter set (\(a\), b, c) = (350, 0.001, − 150) for the pre-contrast data and (350, 0.005, − 150) for the post-contrast data. Since our dataset was in the magnitude scale and \(S(t)\) of Eq. (1) can cover the negative-value range, we incrementally flipped the polarity of \(S({TI}_{i})\) to find the best fit of Eq. (1) [24]. Notably, the estimated T1* does not incorporate the effect of the tip down of the spin magnetization for every repetition time. Hence, the Look-Locker correction was applied in the following way to arrive at the corrected T1 value:
$$T1={T1}^{*}\left(\frac{a}{a+c}-1\right)$$
(3)
Meanwhile, the RD-NLS (or RD) method expands the objective function Eq. (2) and finds the optimal estimates of the parameters separately (refer to Appendix B in Barral et al. [8] for detailed mathematical derivations). The function to be optimized for T1 is independent of the other two parameters. Hence, the decoupling leads to a one-dimensional search problem for T1 estimation.
Python and C++ implementation
For the LM method, we implemented the T1 map calculation in Python. To solve the non-linear least squares problem, we used the scipy.optimize.curve_fit() function, which is based on the Levenberg–Marquardt algorithm. The method is referred to as LM_python. For the RD method, we converted the original MATLAB script written by Barral et al. (code available at http://www-mrsrl.stanford.edu/~jbarral/t1map.html) to Python. We noted that the original RD algorithm failed to work in some cases where the two candidate TIs produced unsatisfactory fitting results. Hence, we modified the RD algorithm by introducing three candidate TIs. This helped remove the unsatisfactory fitting and improved the accuracy (Fig. 2). The method is referred to as RD_python.
We also implemented both the LM- and RD-based T1 map calculations in C++ and used pybind11 [13] to make the compiled C++ code compatible with the Python environment. The C++ implementation was performed on a Windows OS, Microsoft Visual Studio 2017 platform. For the LM-based T1 parameter estimation, we used the solve_least_squares_lm function of the Dlib library [25]. In addition, for the multi-core implementation, we used the OpenMP library [26] for the parallelization of the ‘for’ loop in the pixel-wise T1 map calculation. We chose the static schedule option for parallelization. For evaluation, the methods are referred to as LM_C++_single-core, RD_C++_single-core, LM_C++_multi-core, and RD_C++_multi-core.
To ensure that the curve fitting is correctly performed at each pixel location and the samples are placed close to the T1 relaxation curve, we developed a graphical user interface (GUI) that enables a user to load the DICOM TI image data, select a method for T1 map calculation, and mouse-click a pixel location for displaying its curve fitting result. The GUI also performs the calculation and display of the pre- and post-contrast T1 maps as well as the ECV map (Fig. 1). The GUI was implemented using the PyQT library [27] on a 64-bit Windows PC.
Evaluation
We evaluated the following methods in terms of speed on a Windows PC (AMD Ryzen 7 1800X Eight-Core Processor and 16.0 GB RAM): MRmap, LM_python, RD_python, LM_C++_single-core, RD_C++_single-core, LM_C++_multi-core, and RD_C++_multi-core.
For the evaluation, we used the MRmap software [9], which is publicly available for download at https://sourceforge.net/projects/mrmap/. We chose the following options for T1 map calculation: Limits of T1, T2, and Noise for pre-contrast were set to 3000, 350, and 0, respectively. Limits of T1, T2, and Noise for post-contrast were set to 1500, 350, and 0, respectively. Specifically, the setting of the Noise value had to be consistent, since the choice of the value significantly affected the computation speed. Registration was set to None. Process was set to “T1 mapping—MOLLI,” and Correction was set to “Look-Locker.”
For the evaluation of accuracy in T1 value, we manually drew a region of interest (ROI) in the myocardium in either a pre- or a post-contrast T1 map of each subject and used the same ROI mask for the T1 maps estimated using MRmap and the RD_C++_multi-core method. Mean T1 values were calculated within the myocardial ROIs. Bland–Altman analysis was performed by computing the mean difference and 95% limits of agreement between the two T1 measurements.