1 Introduction

Magnetic resonance imaging (MRI) can be used for the comprehensive diagnosis of myocardial infarction. Three important aspects are: Time-resolved cardiac imaging (CINE imaging) allows contractility assessment of the myocardium [2]. Myocardial viability can be assessed with late gadolinium enhancement (LGE) imaging, which visualizes scar tissue [6]. \(T_1\) mapping can be used for the diagnosis of diffuse myocardial fibrosis [1, 11]. In current clinical practice, separate scans are performed for each aspect. For the case of LGE imaging, an additional pre-scan is required for the a-priori selection of a parameter called inversion time (\(T\!I\)), which ensures that the LGE scan suppresses signal from healthy myocardium while highlighting scar tissue (see Fig. 1).

Fig. 1.
figure 1

Example images from the acquired series and corresponding points on the \(T_1\) recovery curve of different tissues. After inversion, longitudinal magnetization starts at the negative maximum value, \(-M_z\), and then recovers back to \(M_z\). Most notably, the 2\({}^\text {nd}\) image from the left shows the optimal \(T\!I\) time for scar visualization, because the myocardium is black and the scar tissue is hyperintense compared to the blood.

Full size image

A recent method [10] allows the acquisition of a series of images from which all of the above information can be extracted from a single measurement. After an inversion pulse, \(4\,\text {s}\) or \(\sim \)4 cardiac cycles of real-time images are acquired. This leads to an image series where both the image contrast and cardiac motion state change over time (see Fig. 1). The contrast change is related to the tissue-dependent time constant for longitudinal magnetization recovery, \(T_1\), and the time between inversion and acquisition of each image, \(T\!I\). Simultaneously acquired ECG data allows the segmentation of the series into cardiac cycles. Detsky et al. [2] also measure viability and contractility in one scan, but use a longer, segmented acquisition, and don’t jointly estimate both parameters.

Related to the three aspects mentioned above, CINE imaging is provided by the last cardiac cycle in the series, where the image contrast has already stabilized. An image with optimal \(T\!I\) for scar visualization can be selected retrospectively in the first cardiac cycle. For \(T_1^*\) mapping, the approach in [10] is to estimate cardiac motion by image registration of the last cardiac cycle. Under the assumption that motion is identical in the first and last cardiac cycles, this estimate is then used to transform the first cardiac cycle into a consistent motion state to perform pixel-wise least-squares regression of measured values to a signal model [8]. Additionally, a joint scar tissue and contractility visualization is generated by animating the image with optimal \(T\!I\) contrast with the motion estimate to show how wall motion abnormalities relate to scar tissue. The limitations of this method are the assumption of identical motion in the first and last cardiac cycles, and the use of just the first cardiac cycle for \(T_1^*\) mapping.

We propose an iterative approach to jointly estimate cardiac motion for all cardiac cycles and utilize all available data for \(T_1^*\) mapping, retaining the advantages of a-posteriori \(T\!I\) selection and combined contractility and scar visualization, while removing the limitations of the previous approach and improving robustness. We tested our method with phantom data and 11 patient data sets.

2 Materials and Methods

If the motion between cardiac phases is known, i. e. we have deformation fields \({\varvec{d}}_{i\rightarrow {}j}({\varvec{x}})\) such that an image \(m_i({\varvec{x}})\) in one cardiac phase \(i \in [0, M[\), where M is the total number of images, can be transformed as \(m_i({\varvec{x}} + {\varvec{d}}_{i\rightarrow {}j}({\varvec{x}})\) to match the cardiac phase j in image \(m_j({\varvec{x}})\), then pixel-wise \(T_1^*\) mapping is possible using images transformed to a common reference phase. Conversely, given a \(T_1^*\) mapping of an image series, cardiac motion can be estimated using intensity-based registration, similar to the use in [12] for imperfect breath-hold correction.

We propose a fixed-point iteration scheme which alternately improves \(T_1^*\) mapping results and motion estimates:

Fig. 2.
figure 2

Overview of the iterative estimation: The measured image series is transformed to a common cardiac phase using the current motion estimate (a), followed by pixel-wise \(T_1^*\) mapping (b). From the \(T_1^*\) map, synthetic images at the same \(T\!I\) times as the measured images are generated (c), which are then transformed to the different cardiac phases, followed by pair-wise registration of the measured image series and current synthetic image series to update the current motion estimate (d).

Full size image
  1. 1.

    Compute initial motion estimate from last cardiac cycle (see Sect. 2.1)

  2. 2.

    Iteratively update \(T_1^*\) map and motion parameter estimates (see Fig. 2)

    1. (a)

      Transform the originally acquired image series to a single cardiac phase using the current motion estimate

    2. (b)

      Compute pixel-wise \(T_1^*\) parameter estimates (see Sect. 2.2)

    3. (c)

      Using the resulting parameters, compute synthetic images corresponding to the inversion times of the acquired image sequence

    4. (d)

      Update the motion estimate by performing image registration between the original and synthetic images (see Sect. 2.3)

  3. 3.

    Compute the true \(T_1\) map from the final \(T_1^*\) map (see Sect. 2.2)

2.1 Initial Motion Estimation

The first motion estimate is computed with intensity-based registration of the last cardiac cycle of the acquisition. The first frame of the last cycle \(m_f({\varvec{x}})\) is selected as the reference frame, so we obtain deformation fields \({\varvec{d}}_{i\rightarrow {}f}({\varvec{x}})\) for \(i \in [f, M[\). Non-rigid registration based on the elastix toolbox [5] with these components is used: the transform model is uniform cubic B-splines with a \(16\,\text {mm}\) isotropic control point spacing, the similarity metric is cross-correlation, optimized with the L-BFGS algorithm on a 3-level multi-resolution pyramid.

Deformation fields for the other cardiac cycles are then determined by linear interpolation of those from the last cardiac cycle, as we no longer assume that cardiac cycles are of the same length. For the initial motion estimation, we assume that there is no motion between the first images of each cardiac cycle. Thus, a deformation field \({\varvec{d}}_{i\rightarrow {}f_i}({\varvec{x}})\), where \(f_i\) is the first frame of the cardiac cycle containing image i, can also be used as a deformation field \({\varvec{d}}_{i\rightarrow {}0}({\varvec{x}})\) to transform image i to the cardiac state of the very first frame. This allows us to transform all images \(i \in [0, M[\) into a reference cardiac phase for \(T_1^*\) mapping.

2.2 \(T_1^*\) Mapping and \(T_1\) Correction

For \(T_1^*\) mapping, the current motion estimate \({\varvec{d}}_{i\rightarrow {}0}({\varvec{x}})\) (\(i \in [0, M[\)) is used to transform all images \(m_i({\varvec{x}})\) to the cardiac phase of the first image:

$$\begin{aligned} m_i'({\varvec{x}}) = m_i({\varvec{x}} + {\varvec{d}}_{i\rightarrow {}0}({\varvec{x}})). \end{aligned}$$
(1)

The \(m_i'({\varvec{x}})\) now only differ in their contrast, which is dependent on \({T\!I}_i\), the inversion time of image i, as well as tissue-dependent parameters \(a({\varvec{x}})\), \(b({\varvec{x}})\) and \( T_1^*({\varvec{x}})\), which can be determined by a pixel-wise non-linear least-squares fit [8]:

$$\begin{aligned} \min _{a({\varvec{x}}),\ b({\varvec{x}}),\ T_1^*({\varvec{x}})} \textstyle \sum _{i=0}^{M-1} ( m_i'({\varvec{x}}) - | a({\varvec{x}}) + b({\varvec{x}}) \cdot \exp ( -\tfrac{ {T\!I}_i }{ T_1^*({\varvec{x}}) } ) | )^2. \end{aligned}$$
(2)

Unlike [10], all cardiac cycles are considered for the fit, not just the first one, to increase robustness. The residual error of this fit is due to noisy measurements and errors due to misalignment, which can be reduced by improving motion estimates in each iteration. As the MRI measurement influences the longitudinal magnetization recovery, this does not yield the true \(T_1\) value, but an “apparent” \(T_1\) called \(T_1^*\). A correction from \(T_1^*\) to \(T_1\) maps can then be computed [7, 8].

2.3 Motion Estimate Updating

The current set of parameters \(a({\varvec{x}})\), \(b({\varvec{x}})\) and \(T_1^*({\varvec{x}})\) can be used to generate a synthetic image series \(s_i'({\varvec{x}})\) in the reference cardiac phase as

$$\begin{aligned} s_i'({\varvec{x}}) = | a({\varvec{x}}) + b({\varvec{x}}) \cdot \exp ( -\tfrac{ T\!I_i }{ T_1^*({\varvec{x}}) } ) | \end{aligned}$$
(3)

and transformed to the cardiac phase it was acquired in:

$$\begin{aligned} s_i({\varvec{x}}) = s_i'({\varvec{x}} + {\varvec{d}}_{i\rightarrow {}0}^{-1}({\varvec{x}})). \end{aligned}$$
(4)

Pairwise registration between originally measured images \(m_i({\varvec{x}})\) and corresponding synthetic images \(s_i({\varvec{x}})\) then yields the error \({\varvec{e}}_{m_i\rightarrow {}s_i}({\varvec{x}})\) of the current motion estimate \({\varvec{d}}_{i\rightarrow {}0}({\varvec{x}})\), so the motion estimate is updated by concatenation:

$$\begin{aligned} {\varvec{d}}_{i\rightarrow {}0}^\text {(new)} = {\varvec{d}}_{i\rightarrow {}0}^\text {(old)} \circ {\varvec{e}}_{m_i\rightarrow {}s_i}. \end{aligned}$$
(5)

As this step considers each cardiac cycle individually, errors of the initial motion estimate will be corrected at this stage.

2.4 Experiments

Phantom data. A cardiac \(T_1\) map, acquired from a patient with a dedicated \(T_1\) mapping protocol, was used to generate different phantom data sets for the same \(T\!I\) times a real acquisition would use. This allowed to evaluate the algorithm against ground truth motion and \(T_1\) values. Phantom deformation fields mimicking 3 cardiac cycles were generated to simulate contraction of the heart. The resulting magnitude images were distorted with Rician noise [4], the appropriate distribution for magnitude MRI images. The noise parameter was chosen as \(5\,\%\) of the maximum magnitude. Figure 3 shows examples of the phantom image series.

Fig. 3.
figure 3

Phantom images for different \(T\!I\) times and cardiac states.

Full size image

Different datasets for the following 3 scenarios were simulated, with the last two simulating arrhythmia: (i) Same cardiac motion in all cardiac cycles, the ideal case assumed in method [10], (ii) Deviation of the magnitude of cardiac motion, (iii) Deviation of the position of end-systole within the cardiac cycle. For the quantitative evaluation of the \(T_1\) mapping step, the root-mean-squared error (RMSE) of estimated \(T_1\) maps to ground truth and the residual of the least-squares parameter estimation were computed. To evaluate the quality of the motion estimate, the blood pool of the heart was manually segmented in one cardiac phase. The resulting mask was transformed to all other cardiac phases using the ground truth deformation fields as well as the estimated deformation field. Finally, they were compared using Dice [3] and Hausdorff [9] coefficients.

Patient data. Data acquisition as described in [10] was performed in 11 patients (age \(63 \pm 10\), 4 female) with known or suspected status post myocardial infarction on a \(3\,\text {T}\) clinical scanner (MAGNETOM Skyra, Siemens Healthcare, Erlangen, Germany). Prototype sequence parameters include: \(4\,\text {s}\) breath-hold acquisition, \(192 \times 150\) matrix size, \((1.9\,\text {mm})^2\) in-plane, \(8\,\text {mm}\) slice, \(33\,\text {ms}\) temporal resolution.

For \(T_1\) mapping evaluation, the residual was computed. To evaluate the motion estimates, the blood pool of the heart was manually segmented for each end-diastolic (ED) and end-systolic (ES) frame. The last ED frame is also the reference frame for registration. The resulting mask for the reference frame was transformed to all other segmented frames using the motion estimates, and compared to the manual segmentations using the Dice and Hausdorff coefficients.

3 Results and Discussion

Phantom data. Figure 4 plots the Dice coefficients for data sets (i)–(iii) for all images in the series without and with initial motion correction (corresponding to method [10]), and after 1 and 2 fixed-point iterations, after which the mean magnitude of motion estimate updates dropped below \(0.2\,\text {px}\) and no more substantial changes were observed. Table 1 lists \(T_1\) mapping RMSE, residual, and the mean Dice and Hausdorff coefficients over all images. The mean Pearson correlation between \(T_1\) mapping RMSE and residual over all data sets is \(0.98 \pm 0.01\).

Fig. 4.
figure 4

Plots of Dice coefficients for all images in all phantom data sets with the same motion in all cardiac cycles (left), different motion strength in each cardiac cycle (middle), and different positions of end-systole (right). Results without motion correction represent the baseline, followed by initial motion correction using motion information from the last cardiac cycle, and two iterations of the joint \(T_1^*\) and motion estimation.

Full size image
Table 1. Quantitative results for phantom data sets for different stages of motion correction. For the evaluation of the \(T_1\) maps, the RMSE compared to ground truth and residual error of least-squares estimation are given. For the evaluation of the motion estimates, mean Dice and Hausdorff coefficients of all images are given
Full size table

For phantom data set (i), the prerequisites of method [10] are met and fixed-point iteration shows no benefit. For the cases of simulated cardiac arrhythmia (ii) and (iii), our method improves the \(T_1\) mapping RMSE and residual as well as the Dice and Hausdorff coefficients compared to [10], in terms of both mean and standard deviation. The non-zero RMSE using ground truth motion is due to interpolation. The correlation between \(T_1\) mapping RMSE and residual suggests that we can use the latter as a surrogate to evaluate the quality of \(T_1\) maps for patient data, where no ground truth exists.

Patient data. Table 2 lists mean \(T_1^*\) mapping residuals over all patients and mean ED and ES Dice and Hausdorff coefficients without and with initial motion correction (as in method [10]), and after 1 and 5 fixed-point iterations, after which the mean magnitude of motion estimate updates dropped below \(0.4\,\text {px}\) and no more substantial changes were observed. One patient data set contained a premature ventricular contraction (PVC), where the Dice coefficient of ED phase between normal systole and PVC systole using initial motion correction was \(86\,\%\) and after 5 iterations was \(94\,\%\). Our websiteFootnote 1 contains visual results.

Table 2. Quantitative results of patient data evaluation for different stages of motion correction. For the evaluation of the \(T_1\) maps, the residual error of least-squares estimation is given. For the quality of the motion estimates, Dice and Hausdorff coefficients of all ED and ES images are given
Full size table

For patient data, fixed-point iteration improves all quality measures compared to [10], most notably also the standard deviations. The appearance of papillary muscles near the endocardial boundary rendered blood pool segmentation in ES ambiguous, which could explain the comparatively lower ES Dice coefficients. The assumption that there is no motion between the first frames of each cardiac cycle (which are ED frames) is shown to be reasonable by the ED Dice coefficient without motion correction of \(96\,\%\). The results for ED Dice coefficients without and with initial motion correction are necessarily identical due to this assumption. However, cardiac arrhythmia violates the assumption, as demonstrated by the PVC case. Still, the \(T_1^*\) mapping step is robust to some images being misaligned, with only 3 parameters fitted to \(\sim 100\) observations. Thus, fixed-point iteration is able to correct the misalignment and substantially raise the Dice coefficient closer to the mean level for ED. While its application to 2-D imaging is susceptible to through-plane motion artifacts, our proposed algorithm is equally applicable to 3-D data, which would eliminate this problem.

4 Conclusion

Our proposed method provides joint contractility and scar tissue visualization as well as \(T_1\) mapping for the comprehensive assessment of myocardial infarction patients. Evaluations on phantom and patient data show substantial improvements to the state of the art, most notably robustness to cardiac arrhythmia, which is common in this patient population.