High-dimensional imaging has great potential to improve many cardiac applications, but is difficult to perform due to prohibitive data acquisition requirements. For example, first-pass myocardial perfusion quantification may be improved by $$$T_1$$$ mapping¹ (requiring a relaxometry/recovery dimension) and may additionally benefit from analysis at multiple cardiac phases. Here we propose a general framework for high-dimensional cardiac imaging, employing low-rank tensor modeling^2–4 and extending motion-sorted joint reconstruction^5,6 to additional dimensions. We demonstrate a variant of the framework using explicit subspaces (estimated from navigator data and a dictionary of solutions to the Bloch equations, similar to MR fingerprinting⁷) to perform free-breathing, ECG-less $$$T_1$$$ mapping of native myocardium at multiple cardiac phases, as well as ECG-less, cardiac- and time-resolved $$$T_1$$$ mapping during first-pass perfusion.

We model a high-dimensional cardiac image $$$a$$$ as a low-rank tensor $$$\mathcal{A}$$$, partially separable in some combination of space ($$$\mathbf{x}$$$), cardiac phase ($$$t_\mathrm{c}$$$), respiratory phase ($$$t_\mathrm{r}$$$), wall-clock time ($$$t$$$, e.g. to represent contrast agent dynamics), and variable sequence or timing parameters $$$p_1,p_2,\dots,p_N$$$ (e.g., flip angle, recovery time):$$\begin{array}{l}a(\mathbf{x},t_\mathrm{c},t_\mathrm{r},t,p_1,\dots,p_N)\\\quad=\sum_{\ell_\mathrm{x}=1}^{L_\mathrm{x}}\sum_{\ell_\mathrm{c}=1}^{L_\mathrm{c}}\sum_{\ell_\mathrm{r}=1}^{L_\mathrm{r}}\sum_{\ell_\mathrm{t}=1}^{L_\mathrm{t}}\sum_{\ell_\mathrm{1}=1}^{L_\mathrm{1}}\dots\sum_{\ell_\mathrm{N}=1}^{L_\mathrm{N}}\\\quad\qquad{b}_{\ell_\mathrm{x}\ell_\mathrm{c}\ell_\mathrm{r}\ell_\mathrm{t}\ell_1\dots\ell_N}u_{\mathrm{x},\ell_\mathrm{x}}(\mathbf{x})u_{\mathrm{c},\ell_\mathrm{c}}(t_\mathrm{c})u_{\mathrm{r},\ell_\mathrm{r}}(t_\mathrm{r})u_{\mathrm{t},\ell_\mathrm{t}}(t)u_{1,\ell_1}(p_1)\dots u_{N,\ell_N}(p_N).\end{array}\quad(1)$$ The $$$u$$$’s denote basis functions, $$$L$$$’s denote model orders, and $$$b$$$ denotes the core tensor. Equation (1) is more compactly represented in matrix form as$$\mathbf{A}_{(1)}=\mathbf{U}_\mathrm{x}\mathbf{B}_{(1)}(\mathbf{U}_\mathrm{p_N}\otimes\dots\otimes\mathbf{U}_\mathrm{p_1}\otimes\mathbf{U}_\mathrm{t}\otimes\mathbf{U}_\mathrm{r}\otimes\mathbf{U}_\mathrm{c})^T,\quad(2)$$where the columns of each $$$\mathbf{U}$$$ contain basis functions, $$$\otimes$$$ denotes the Kronecker product, and subscript $$$(i)$$$ denotes mode-$$$i$$$ unfolding of the tensor⁸.

Given (1), $$$\mathcal{A}$$$ can be recovered from undersampled measured data $$$\mathbf{y}$$$ via low-rank tensor completion, e.g.,$$\arg\min_\mathcal{A}\|\mathbf{y}-\Omega(\mathbf{FSA}_{(1)})\|_2^2+\lambda\sum_i\|\mathbf{A}_{(i)}\|_*+R(\mathcal{A}),\quad(3)$$with undersampling operator $$$\Omega$$$, Fourier transform $$$\mathbf{F}$$$, coil sensitivity operator $$$\mathbf{S}$$$, and regularization functional $$$R(\cdot)$$$. However, it may instead be preferable to constrain the solution to an explicit subspace, which has many benefits for low-rank imaging⁹. Image reconstruction then reduces to recovery of the $$$L_\mathrm{x}$$$ images in $$$\mathbf{U}_\mathrm{x}$$$:$$\arg\min_{\mathbf{U}_\mathrm{x}}\|\mathbf{y}-\Omega(\mathbf{FSU}_\mathrm{x}\mathbf{\Phi})\|_2^2+R(\mathbf{U}_\mathrm{x}),\quad(4)$$with known $$$\mathbf{\Phi}=\mathbf{B}_{(1)}(\mathbf{U}_\mathrm{p_N}\otimes\dots\otimes\mathbf{U}_\mathrm{p_1}\otimes\mathbf{U}_\mathrm{t}\otimes\mathbf{U}_\mathrm{r}\otimes\mathbf{U}_\mathrm{c})^T$$$.

One approach to define $$$\mathbf{\Phi}$$$ is to frequently collect a small set of “navigator” $$$\mathbf{k}$$$-space locations (data which are also useful for identifying cardiac/respiratory phases). The $$$(\mathbf{k},t_\mathrm{c},t_\mathrm{r},t,p_1,\dots,p_N)$$$-space tensor $$$\mathcal{D}_\mathrm{nav}$$$—defined only at the navigator $$$\mathbf{k}$$$-space coordinates—will be missing relatively few entries, which are recoverable by small-scale low-rank tensor completion, e.g.,$$\mathcal{\hat{D}}_\mathrm{nav}=\arg\min_{\mathcal{D}_\mathrm{nav}}\|\mathbf{y}_\mathrm{nav}-\Omega(\mathcal{D}_\mathrm{nav})\|_2^2+\lambda\sum_i\|\mathbf{D}_{\mathrm{nav},(i)}\|_*+R(\mathcal{\mathcal{D}_\mathrm{nav}}).\quad(5)$$$$$\mathbf{\Phi}$$$ can then be extracted from $$$\mathcal{\hat{D}}_\mathrm{nav}$$$, e.g., by the truncated HOSVD¹⁰.

For certain applications, we also propose to predetermine some of the matrix subspaces and further constrain (5) so that $$$\mathcal{D}_\mathrm{nav}$$$ (and by extension $$$\mathcal{A}$$$) lives in these subspaces. For example, with inversion recovery fitting ($$$p_1=T_\mathrm{I}$$$), we predetermine the recovery subspace by generating a dictionary of signal curves from the Bloch equations (factoring in the potential $$$B_1$$$ homogeneity and inefficient inversion) and extracting basis functions from this dictionary’s SVD.

All data were acquired on a 3T Siemens Verio scanner using a modified golden-angle radial sequence with parameters in Fig. 1. The $$$0^\circ$$$ radial line was acquired every other $$$\alpha$$$-pulse for$$$\mathbf{y}_\mathrm{nav}$$$. Only 1 minute of data were used for each reconstruction. The full recovery curve was sampled by continually applying alpha pulses in between preparation pulses. We reconstructed the full set of inversion times experienced by the golden angle readouts (345 for IR-FLASH, 71 for SR-FLASH) in order to avoid grouping lines from different recovery times, as frequent collection of center $$$\mathbf{k}$$$-space in radial trajectories temporally blurs image contrast. $$$T_1$$$, amplitude, $$$\alpha$$$, and preparation pulse efficiency were calculated pixel-by-pixel from the reconstructed images.

Figure 2 depicts end-systole and end-diastole images and inversion recovery profiles for Subject 1, as well as native $$$T_1$$$ maps and MOLLI^11,12 reference $$$T_1$$$ maps. Figure 3 depicts the same for Subject 2. Figure 4 shows first-pass perfusion images, saturation recovery curves as a function of wall-clock time, and quantitative $$$R_1$$$ curves for the LV blood pool (at end-diastole for maximum LV volume) and septal myocardial segments (at end-systole for maximum myocardial thickness).

Native $$$T_1$$$ measurements with the proposed framework show good agreement with MOLLI—the absolute difference in mean myocardial $$$T_1$$$s ranged from 25 ms (Subject 2 diastole) to 80 ms (Subject 1 diastole; all other differences were <40ms)—and did so without breathhold- or ECG-reliance. There is no standard method to compare first-pass $$$T_1$$$ maps against, but the measured blood pool and myocardial signal curves obeyed typical patterns regarding contrast agent dynamics.

We have presented a novel framework exploiting the low-rank tensor structure of high-dimensional cardiac images, demonstrating its use for two challenging applications. The framework enhances the practical utility of native $$$T_1$$$ mapping—using efficient, continual acquisition and advanced reconstruction techniques to overcome the practical limitations of ECG and breathholds; and 2) demonstrating time-resolved $$$T_1$$$ mapping during first-pass perfusion, which may allow direct quantification of tissue contrast agent concentration. Future work includes quantitative evaluation and validation of these applications as well as exploration of other high-dimensional applications.