Free-breathing and ungated cardiac cine using navigator-less spiral SToRM

TL;DR

This paper presents a novel kernel low-rank algorithm for free-breathing, ungated cardiac MRI reconstruction from spiral data without navigators, outperforming classical methods in quality and quantitative metrics.

Contribution

It introduces a kernel low-rank matrix completion approach that directly recovers missing k-space data without explicit motion navigation, addressing challenges of motion-induced low-rankness.

Findings

Improved image quality over classical low-rank and XD-GRASP methods.

Quantitative metrics comparable to breath-held cine data.

Validated with simulated and in-vivo data.

Abstract

We introduce a kernel low-rank algorithm to recover free-breathing and ungated dynamic MRI from spiral acquisitions without explicit k-space navigators. It is often challenging for low-rank methods to recover free-breathing and ungated images from undersampled measurements; extensive cardiac and respiratory motion often results in the Casorati matrix not being sufficiently low-rank. Therefore, we exploit the non-linear structure of the dynamic data, which gives the low-rank kernel matrix. Unlike prior work that rely on navigators to estimate the manifold structure, we propose a kernel low-rank matrix completion method to directly fill in the missing k-space data from variable density spiral acquisitions. We validate the proposed scheme using simulated data and in-vivo data. Our results show that the proposed scheme provides improved reconstructions compared to the classical methods such as low-rank and XD-GRASP. The comparison with breath-held cine data shows that the quantitative metrics agree, whereas the image quality is marginally lower.

Figures

Tables

Table 1. TABLE I: Quantitative comparison of the methods on simulated data in Fig. 3 using the signal-to-error ratio (SER), normalized high frequency error (HFEN), structural similarity index (SSIM) and global phase coherence (GPC) metrics. All of these metrics are computed in a square region of interest around the cardiac region. Higher values of the above-mentioned performance metrics correspond to better reconstruction except for the HFEN, where a lower value is better. These comparisons show that the proposed scheme performs better than the other methods except for SToRM: Self-Nav, where dedicated k-space navigators are used to estimate the Laplacian matrix.

Method	SER	SSIM	HFEN	GPC
Low Rank	$16.60\pm 0.90$	$0.80\pm 0.03$	$0.42\pm 0.05$	$115\pm 33$
SToRM: SENSE	$16.47\pm 0.86$	$0.71\pm 0.03$	$0.47\pm 0.04$	$112\pm 60$
Compressed Sensing	$17.90\pm 1.60$	$0.80\pm 0.06$	$0.38\pm 0.07$	$106\pm 29$
kt-SLR	$18.60\pm 0.70$	$0.87\pm 0.02$	$0.25\pm 0.01$	$204\pm 35$
SToRM: Self-Nav	$25.01\pm 1.41$	$0.94\pm 0.02$	$0.128\pm 0.02$	$428\pm 113$
SToRM: Two-Step	$25.77\pm 1.50$	$0.95\pm 0.02$	$0.107\pm 0.02$	$\mathrm{𝟒𝟕𝟓}\pm \mathrm{𝟏𝟏𝟎}$

Table 2. TABLE III: In this experiment, we have fixed the number of iterations=5 and compared all approaches using Signal to error ratio (SER) and reconstruction time. First row corresponds to the reconstructions from the single-step SToRM approach, described in Section III-B , while the two step strategy described in Section III-C is shown in the second row. However, the computational complexity of this approach is 3-4 times higher. In addition, the difference in SER is not significant. In favor of faster experiments, we resort to two-step recovery scheme, where a low-resolution reconstruction is used to estimate the Laplacian matrix as described in Section III-C . We have also compared Laplacian matrices estimated from original images and reconstructed images. As noted in the last row of the Table, we note that the SER of the original images is marginally higher.

Laplacian Estimation Scheme	Iterations	SER	time(minutes)
Low resolution	$5$	$25.4$	$8$
Full resolution	$5$	$25.8$	$33$
Original Images	$0$	$26.1$	$3.2$

Equations

$${\int }_{\mathcal{M}}\Vert \nabla x{\Vert }^2\approx \frac{1}{2}\sum _{i,j=1}^k{\mathbf{w}}_{i,j}\Vert {\mathbf{x}}_i-{\mathbf{x}}_j{\Vert }^2=\text{ trace}(\mathbf{X}\mathbf{L}{\mathbf{X}}^H).$$

$${\mathbf{w}}_{ij}={\mathbf{e}}^{-\frac{\Vert {\mathbf{z}}_i-{\mathbf{z}}_j{\Vert }^2}{{\sigma }^2}}.$$

$${\mathbf{X}}^{\ast }=\arg \underset{\mathbf{X}}{\min }\Vert \mathcal{A}(\mathbf{X})-\mathbf{B}{\Vert }_F^2+\lambda \mathrm{trace}(\mathbf{X}\mathbf{L}{\mathbf{X}}^H).$$

$$\mathcal{S}=\{\mathbf{x}|\psi (\mathbf{x})=0\},$$

$$\varphi (\mathbf{x})=\left[\begin{array}{c}\frac{1}{{\sigma }^2\Vert \mathbf{k}\Vert }\exp \left(j{\mathbf{k}}_1^T\mathbf{x}\right)\\ ⋮\\ \frac{1}{{\sigma }^2\Vert \mathbf{k}\Vert }\exp \left(j{\mathbf{k}}_P^T\mathbf{x}\right),\end{array}\right]$$

$$\Phi (\mathbf{X})=\left[\begin{array}{cccc}\varphi ({\mathbf{x}}_1)& \varphi ({\mathbf{x}}_2)& \dots & \varphi ({\mathbf{x}}_N).\end{array}\right]$$

$$\begin{array}{r}{\mathbf{X}}^{\ast }=\arg \underset{\mathbf{X}}{\min }\Vert \mathcal{A}(\mathbf{X})-\mathbf{B}{\Vert }_F^2+{\lambda }_1{\Vert \Phi (\mathbf{X})\Vert }_{\ast }.\end{array}$$

$${\left[\mathcal{K}(\mathbf{X})\right]}_{i,j}=\exp \left(-\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j{\Vert }^2}{2{\sigma }^2}\right).$$

$$\begin{array}{r}{\mathbf{X}}^{(n)}=\arg \underset{\mathbf{X}}{\min }\Vert \mathcal{A}(\mathbf{X})-\mathbf{B}{\Vert }_F^2+{\lambda }_1\text{ trace}(\mathbf{X}{\mathbf{L}}^{(n)}{\mathbf{X}}^H)\end{array}$$

$${\mathbf{L}}^{(n)}={\mathbf{D}}^{(n)}-{\mathbf{W}}^{(n)}.$$

$${\mathbf{W}}^{(n)}=-\frac{1}{{\sigma }^2}\mathcal{K}\left({\mathbf{X}}^{(n-1)}\right)\odot {\left[\mathcal{K}\left({\mathbf{X}}^{(n-1)}\right)+\gamma \mathbf{I}\right]}^{-\frac{1}{2}},$$

$$\begin{array}{r}{\mathbf{X}}_L^{\ast }=\arg \underset{{\mathbf{X}}_L}{\min }\Vert {\mathcal{A}}_L({\mathbf{X}}_L)-{\mathbf{B}}_L{\Vert }_F^2+{\lambda }_1{\Vert \Phi ({\mathbf{X}}_L)\Vert }_{\ast }+\\ {\lambda }_2\underset{\text{trace}({\mathbf{X}}_L{\mathbf{L}}_{\text{tik}}{\mathbf{X}}_L^H)}{\underbrace{\sum _i\Vert {\mathbf{x}}_{(i+1)}-{\mathbf{x}}_i{\Vert }^2}}.\end{array}$$

$$\begin{array}{r}{\mathbf{X}}_L^{(n)}=\arg \underset{{\mathbf{X}}_L}{\min }\Vert {\mathcal{A}}_L({\mathbf{X}}_L)-{\mathbf{B}}_L{\Vert }_F^2+\text{ trace}({\mathbf{X}}_L{\mathbf{L}}_{\text{eq}}^{(n)}{\mathbf{X}}_L^H),\end{array}$$

$$\mathbf{SER}=20{\log }_{10}\frac{||{\mathbf{x}}_{\text{ orig}}|{|}_2}{||{\mathbf{x}}_{\text{orig}}-{\mathbf{x}}_{\text{rec}}|{|}_2},$$

$$\mathbf{HFEN}=20{\log }_{10}\frac{||\text{LoG}({\mathbf{x}}_{\text{orig}})-\text{LoG}({\mathbf{x}}_{\text{rec}})|{|}_2}{||\text{LoG}({\mathbf{x}}_{\text{orig}})|{|}_2},$$

$$\begin{array}{r}{\mathbf{X}}_S^{\ast }=\arg \underset{{\mathbf{X}}_S}{\min }\Vert {\mathcal{A}}_S({\mathbf{X}}_S)-{\mathbf{B}}_S{\Vert }_F^2+\Vert {\mathbf{X}}_S{\Vert }_F^2,\end{array}$$

Full text

I Introduction

Breath-held cine MRI is an integral part of clinical cardiac exams. It is widely used for the anatomical and functional assessment of the heart. Diagnostic cine images require breath holding which results in a long scan time to achieve better spatial and temporal resolution. It is often challenging for children, patients with heart failure and patients with respiratory complications such as chronic obstructive pulmonary disease (COPD) [1]. In addition, multiple breath holds along with intermittent pauses also prolong the scan time, adversely impacting patient comfort and compliance. The scans from different slices may also suffer from inconsistencies between breath-held positions [2]. The acceleration of breath-held cine MRI has been the subject of extensive research in the recent past. Classical approaches include parallel MRI, where the diversity of coil sensitives are exploited to reduce the breath-held duration.Recent approaches further improve the performance by exploiting the structure of x-f space [3], sparsity [4], low-rank property [5, 6], low-rank +sparsity [7], learned dictionaries [8], motion-compensated methods [9], deep learning methods [10] and kernel low-rank methods [11]. When the subjects cannot hold their breath, a standard alternative is real-time imaging, which does not require breath holding or ECG gating. However, these methods have been shown to sacrifice spatial and/or temporal resolution [12, 13]. Another approach is the use of diaphragmatic navigators, which restricts the acquisition to images in the specific respiratory phase [14]. The drawbacks of these schemes include respiratory gating efficiency and variability in the scan time. Several methods that rely on radial acquisitions were introduced in recent years to estimate the cardiac and respiratory phases from the central k-space regions using band-pass filtering [15]. These methods usually require careful selection of receiver coils to obtain self-gating signals, as each coil has different sensitivity to cardiac and respiratory motions. The data is then binned to the respective phases, followed by reconstruction using compressed sensing [15] or low-rank tensor methods [16]. Methods that rely on respiratory motion compensation followed by binning have also been introduced to improve computational efficiency [17]. A challenge with these approaches is the dependency on the phase estimation using band-pass filtering that relies on cardiac and respiratory rates, which may degrade in the presence of irregular respiratory motion or arrhythmia [18, 19]. Since these methods rely on the explicit segmentation of the data into their respective phases, the applicability of these schemes for arrhythmia [20] or for non-cardiac applications (e.g, speech) is not straightforward.

We recently introduced the smoothness regularization on manifold (SToRM) approach, which enables ungated cardiac cine imaging in the free-breathing mode using radial acquisitions [21, 18]. SToRM algorithm assumes that the images lie on a smooth and low-dimensional manifold, parameterized by a few variables (e.g. cardiac and respiratory phases). We note that the smooth manifold/surface model is a non-linear generalization of the linear subspace/low-rank models. These models represent the dynamic dataset more efficiently as compared to the subspace models, which result in reduced blurring in free-breathing applications with extensive cardiac and respiratory motion [21, 18]. The manifold prior facilitates the implicit sharing of data between images in the dataset that have similar cardiac or respiratory phases, which is an alternative to explicit motion-resolved strategies [15, 16]. While this approach does not perform explicit binning of data as in other studies [15, 16], it still exploits the similarity of images in the time series and can be viewed as a soft-binning strategy; a particular image is not assigned to any phase, but the inter-frame weights indicate the similarity of the image with other images in the time series. Since the framework does not require complex processing steps that assume the periodicity of the cardiac/respiratory motion, it is readily applicable to several dynamic applications, including speech imaging, as shown in previous work [21], and cardiac applications involving arrhythmia. We note that there are similar manifold regularization schemes that also rely on the non-linear structure of the data to recover cardiac MRI data [11, 22]. The main difference of our algorithm from past kernel low-rank methods [11] is that we do not require an explicit evaluation of the image feature maps. The work [22] relies on a sparse optimization scheme to recover the Laplacian matrix from navigator data, which is used to recover the data as in our work [18].

Our previous implementation, which we refer to as SToRM:Self-Nav, as well as [11, 22] relied on explicit radial k-space navigators to estimate the manifold structure. Specifically, a few radial spokes with the same orientations are played out periodically. The manifold Laplacian estimated from the navigator data is used to recover the images [18]. Compared to the short radial readouts, the longer spiral readouts considered in this work offers improved sampling efficiency; this approach enables the acquisition of more k-space samples in a given scan time. These longer readouts along with higher flip angles also offer improved myocardial contrast. However, the direct use of SToRM:Self-Nav in our setting results in a large overhead; this approach would require the acquisition of one navigator readout for every frame (corresponding to 3-4 spiral interleaves). In addition, the navigated approach cannot be readily applied to golden angle sequences implemented on several scanners without dedicated navigators. To minimize the above problems, we generalize the SToRM algorithm to recover free-breathing and ungated cardiac MRI data from a variable-density spiral gradient echo (GRE) acquisition without any navigators. Since the navigators are not available, we use the kernel low-rank penalty in the matrix completion setting to fill the missing entries. To improve computational efficiency, we rely on a two-step approach, where the low-rank matrix completion is first applied to low-resolution data. This step is computationally efficient because the size of the images are small. Once the low-resolution data is obtained, the manifold Laplacian estimated from this data is used to recover the high-resolution images. We rigorously validate the spiral SToRM approach against conventional algorithms as well as breath-held cine, both quantitatively and qualitatively using simulated as well in-vivo multi-slice data.

II Background

The SToRM framework relies on the manifold structure of images in the real-time cardiac MRI. The main focus of this work is to extend the SToRM framework with explicit k-space self-gating navigators to a navigator-free setting, which increases the sampling efficiency. The proposed navigator-free SToRM algorithm is enabled by the variable density sampling offered by the spiral gradient echo (GRE) acquisition. The GRE acquisition is free from banding artifacts and does not require additional frequency scouts, which are needed to minimize banding artifacts in steady state free precession (SSFP) sequences on 3T scanners. In comparison to the navigated radial acquisition scheme in earlier work [21], the spiral acquisition scheme offers higher sampling efficiency and signal-to-noise ratio. In addition, we propose to exploit the central k-space regions that are densely sampled relative to the radial settings, which would eliminate the need for navigator-based acquisition to determine the Laplacian matrix. We now briefly describe the SToRM framework.

II-A Overview of SToRM framework [21]

The SToRM scheme models the images in the dataset as points on a low-dimensional smooth manifold in high-dimensional space, which is equal to the number of pixels in each image. The SToRM algorithm relies on minimizing a smoothness prior on the manifold of images $\mathbf{x}_{i}$ in the time series to exploit this structure:

[TABLE]

Here, $\mathbf{X}$ is the Casorati matrix, whose columns correspond to $\mathbf{x}_{i};i=1,..,k$. The weights $\mathbf{w}_{i,j}$ specify the neighborhood structure on the points/images, or equivalently the similarity between images in the dataset.

In our previous work, we relied on 4 radial navigator spokes with the same orientation that were played out periodically (repeated every 10 spokes) to estimate the weights $\mathbf{w}_{i,j}$ using kernel low-rank estimation [21, 18]. These readouts were referred as k-space navigators. The rest of the radial spokes are played out in the golden angle view ordering. Denoting k-space data from the navigator spokes at the $i^{\rm th}$ image by $\mathbf{z}_{i}$, the weights are estimated from the equation as

[TABLE]

Note that the above choice assigns higher weights to image pairs $\mathbf{x}_{i}$ and $\mathbf{x}_{j}$, if the differences of their k-space navigators specified by $\|\mathbf{z}_{i}-\mathbf{z}_{j}\|^{2}$ indicate that they are neighbors on the manifold. Here, $\sigma$ is a parameter that controls the smoothness of the manifold. $\mathbf{L}=\mathbf{D}-\mathbf{W}$ is the Laplacian matrix in (1). Here, $\mathbf{D}$ is a diagonal matrix with elements defined as $\mathbf{D}_{ii}=\sum_{j}\mathbf{W}_{ij}$.

Once $\mathbf{L}$ is available, SToRM performs the joint recovery of the images in the dataset by solving the following problem:

[TABLE]

Here $\mathcal{A}$ is the measurement operator that accounts for the multichannel spiral sampling of the columns of $\mathbf{X}$, which are the image frames.

II-B Bandlimited SToRM model

We consider a bandlimited surface model, where the images $\mathbf{x}_{1},..\mathbf{x}_{N}$ are modeled as high-dimensional points on a smooth surface [23, 18]. We model the surface as the zero level-set of a band-limited function $\psi(\mathbf{x})$:

[TABLE]

where $\psi(x)$ are linear combination of exponentials, whose frequencies are supported at $\mathbf{k}_{1},\ldots,\mathbf{k}_{P}$ on a discrete lattice.

With the bandlimited assumption, we have shown that the exponential feature maps $\Phi(\mathbf{x})$ of the images specified by

[TABLE]

live in a low-dimensional subspace. This implies that the feature matrix

[TABLE]

is low-rank [18]. If the rank is $r$, we can find an orthonormal basis $\mathbf{Q}\in\mathbb{C}^{(n-r)\times n}$ of the null-space such that $\Phi(\mathbf{X})\mathbf{Q}=0$.

Similar to PSF methods [5], we have estimated the above subspace from the inverse Fourier transform of the navigator readouts, denoted by $\mathbf{Z}$. Note that in the radial setting $\mathbf{Z}$ corresponds to projections of the $\mathbf{X}$ along specific orientations. We assume that the null-space $\mathbf{Q}$ of the feature matrix $\Phi(\mathbf{X})$ to be the same as the null-space of $\Phi(\mathbf{X})$. We estimate $\mathbf{Q}$ by picking the $n-r$ lowest singular vectors of the kernel matrix $\mathcal{K}(\mathbf{Z})=\Phi(\mathbf{Z})^{H}\Phi(\mathbf{Z})$.

III Proposed approach

The above approach works well with radial k-space navigator lines, which are fully sampled along the readout direction. In addition 2-4 lines are often needed to reliably estimate the null-space, which reduces the sampling efficiency. In this work, we propose to use a spiral trajectory which provides improved sampling efficiency than radial acquisition. We also eliminate the need for k-space navigators to further improve the sampling efficiency.

III-A Kernel low-rank matrix completion for spiral cine data

We propose to recover the images $\mathbf{x}_{1},..\mathbf{x}_{N}$ from their undersampled measurements by relying on kernel low-rank matrix completion:

[TABLE]

Due to extensive cardiac and respiratory motion, the matrix $\mathbf{X}$ may have a higher rank since the images may not lie on a subspace with small dimension. However, the images may lie on a smooth surface, resulting in the feature matrix being low-rank. The second term in (7) is the nuclear norm of the non-linear features of the images $\mathbf{x}_{i}$, which promotes the low-rank nature of $\Phi(\mathbf{X})$. This prior forces the feature maps $\Phi(\mathbf{x}_{i})$ to a subspace, which is equivalent to encouraging the images $\mathbf{x}$ to lie on smooth surface specified by $S$. This approach is a non-linear generalization to classical low-rank/subspace models [6], which are widely used in dynamic imaging. While this approach was used in the denoising setting [18], the utility of this scheme in completing dynamic MRI datasets has not been reported. Results are available for polynomial varieties [24], but their utility in medical imaging have not been explored.

III-B Iterative reweighted algorithm for matrix completion

The direct implementation of eq. (7) would require the non-linear mapping between the images $\mathbf{x}_{i}$ and their features $\phi(\mathbf{x}_{i})$, as well as their inverse [11]. However, this approach is computationally infeasible in our setting since the dimension of the feature matrix $\Phi(\mathbf{X})$ is too large. Therefore, we use an algorithm that relies on the Gram matrix of $\Phi(\mathbf{X})$, denoted by $\mathcal{K}(\mathbf{X})=\Phi(\mathbf{X})^{H}\Phi(\mathbf{X})$, which is referred as the kernel matrix. For the specific choice of exponential maps as in (5), the entries of $\mathcal{K}(\mathbf{X})$ can be computed directly as:

[TABLE]

without requiring the evaluation of the features $\Phi(\mathbf{x}_{i});i=1,..,N$. This approach is widely known as the kernel trick in machine learning [25]. Specifically, we use the iterative reweighted least squares algorithm with gradient linearization [18] to obtain an alternating algorithm to solve (7). This algorithm alternates between

[TABLE]

and update of the matrix $\mathbf{L}^{(n)}$:

[TABLE]

Here, the weight matrix at the $n^{\textrm{th}}$ iteration is specified by

[TABLE]

and $\mathbf{D}^{(n)}$ is the diagonal matrix with diagonal entries $\mathbf{D}^{(n)}_{ii}=\sum_{j}\mathbf{W}^{(n)}_{i,j}$.

Note that the above approach aims to recover a dataset with approximately 400 images. Hence, the computation cost associated with this scheme is high, especially in the high spatial resolution settings. In later section, we have compared our approach with full resolution Laplacian estimation approach to show the benefit of our scheme.

III-C Improve computational efficiency using SToRM: Two-Step

We rely on a two step approach to keep the computational complexity of the algorithm minimal. We call this approach as SToRM: Two-Step. First we propose to recover the low-resolution images from the central k-space samples as shown in Fig. 1. Note that the low-resolution images are still undersampled. We use the formulation in (7) with the addition of a Tikhnonov temporal prior:

[TABLE]

Here, $\mathcal{A}_{L}$ and $\mathbf{B}_{L}$ are the forward model and the measured multichannel k-space data corresponding the the central k-space regions. By exploiting the similarity of adjacent temporal neighbors, this approach is expected to further improve performance over (7). The smaller size of the images translates to a faster algorithm. Note that the Tikhonov prior can be rewritten as $\mathrm{{\rm trace}}(\mathbf{X}_{L}{\mathbf{L}_{\rm{tik}}}\mathbf{X}_{L}^{H})$, where $\mathbf{L}_{\rm{tik}}$ is the matrix with block diagonal matrix with entries as $[1,-2,1]$. Once the above algorithm converges, the estimated Laplacian matrix is then used to recover the high-resolution image frames from their undersampled measurements by solving (9) as shown in Fig. 1. The iterative algorithm for solving this cost function is similar to the one in Section III-B, where the step (9) is modified as

[TABLE]

where $\mathbf{L}_{\textrm{eq}}^{(n)}=\lambda_{1}~{}\mathrm{\mathbf{L}^{(n)}}+\lambda_{2}~{}\mathrm{\mathbf{L}_{\rm{tik}}}$. Since we do not alternate between the Laplacian update and the image update in the high-resolution setting, we obtain a fast algorithm. Once the above iterative algorithm converges, we use the $\mathbf{L}_{\rm eq}$ matrix to recover the high-resolution data by solving the quadratic optimization scheme specified by (9). This two-step approach is illustrated in Fig. 1.

IV Experimental details

IV-A Datasets

We use the following datasets for the experimental evaluation of the proposed algorithm:

Simulated Dataset: A retrospective ECG-gated, breath-held cardiac MRI is used to create simulated ungated, free-breathing data, as described in by Zhao et al.[5]. The ground truth breath-held SSFP dataset is warped in space and time to mimic respiratory motion and temporally varying heart rate. Please see Zhao et al. [5] for details. The deformed datasets are combined to form an image sequence with multiple cardiac cycles. This free-breathing dataset has a reasonable amount of inter-frame motion due to respiratory dynamics. The dataset has 200 phase encodings, 256 samples per readout, and 256 temporal frames. FOV= $273$mm $\times 350$mm, spatial resolution$=1.36$mm $\times 1.36$mm and TR$=3$ms. For simulated data, variable density spirals are used with 12 spirals per frame. This simulated dataset enables the quantitative comparison of methods, especially in the free-breathing setting where ground truth is not available.

In-vivo Datasets: Six single slice cardiac data and five whole-heart multi-slice datasets were collected in the free-breathing mode using a golden angle spiral trajectory. We compare the proposed scheme with breath-held whole heart bSSFP Cartesian acquisition on five datasets with 10-13 slices to cover the whole heart. Acquisition parameters: TR/TE= 3.1ms-3.4ms/1.18ms-1.28ms, flip angle= 31-39 degrees. All acquisitions except one were performed on a 3T scanner (MAGNETOM Prisma, Siemens Healthineers, Erlangen, Germany). One dataset was acquired on the GE 3T scanner. Image datasets were acquired using the standard body phased-array RF coil. Subjects included three females (age: 25-27) and eight males (age: 20-30) with short-axis view cine data. The institutional review board at the local institution approved all the in-vivo acquisitions, and written consent was obtained from all subjects. The sequence parameters were: TR/TE= 7.8 ms/1 ms, FOV= 320 mm, Base resolution= 256, Bandwidth= 390 Hz/pixel, flip angle= 15 degrees, slice thickness= 8 mm. Dual-density spirals were generated using a Fermi function with a k-space density of 0.2x Nyquist for the first 20% of the trajectory and an ending density of 0.02x Nyquist [26]. The spirals were continuously acquired with rotation of the trajectory by the golden angle between spirals. Off-resonance effects were minimized by using a short spiral readout duration (5 ms) and by using the vendor-provided cardiac shim routine. Post-acquisition, five spirals per frame were binned to obtain the temporal resolution of 40 ms. For multi-slice data, we have collected 10-13 slices from apex to base to cover the whole heart. The parameters of our reconstruction algorithm were manually optimized on one dataset and kept fixed for rest of the datasets: $\lambda_{1}$=0.01, $\lambda_{2}$=1e-5, $\sigma=4.5$, and $\lambda$=0.025. $\lambda_{1}$ and $\lambda_{2}$ involve the trade off between blurring and aliasing artifacts. Here, $\sigma$ is the width of the temporal kernel. We notice that the results are not too sensitive to $\sigma$; kernel low-rank regularization makes the algorithm relatively insensitive to kernel width. The above mentioned, $\gamma$ $(=100)$ is a smoothing parameter that is decreasing in each iteration. In order to have a long-run stability of the algorithm, we decrease until it approaches to pre-determined minimum value. Further details are mentioned in Ongie et.al. [27]

IV-B Imaging Experiments

All the results were generated using a single node of a high-performance Argon Cluster at the University of Iowa, equipped with an Intel Xeon CPU with 28 Cores at 2.40 GHz with 128 GB of memory running on Red Hat Linux MATLAB R2016b. The reconstruction time of the proposed method was between 8 to 10 minutes (400 time frames). This reconstruction is for single slice reconstruction. Low resolution stage takes around 65% of the time and rest is spent in the final reconstruction, as shown in the Table. III.

Coil selection and compression: We acquired the dataset using 34 coils. However, we excluded the coils with low sensitivities in the region/slice of interest. We used an automatic algorithm to pre-select the 10 best coil images that provided the best signal to noise ratio in the heart region; we observed that removing the unreliable coils resulted in improved reconstructions [26]. This algorithm binned the k-space data from several images to recover the low-resolution coil images. We then used PCA-based coil combination using SVD such that the approximation error was $<5\%$. In most cases, we noted that 5-6 coils were sufficient to bring the approximation error to $<5\%$. The coil sensitivity maps were estimated from these coil-combined virtual channels using the method designed by Walsh et al. [28] and assumed to be constant over time. Our experiments (not included in the paper) show that this coil combination has minimal impact on image quality. The main motivation for the combination was to reduce the memory requirement so that it fit on our GPU device, which significantly reduced the computational complexity.

Performance Metrics: We used four quantitative metrics to compare our method against the existing schemes:

• •

Signal-to-Error Ratio (SER):

[TABLE]

where $||\cdot||_{2}$ donates the $\ell_{2}$ norm, and $\mathbf{x}_{orig}$ and $\mathbf{x}_{rec}$ denote the original and the reconstructed images, respectively.

• •

Normalized High Frequency Error (HFEN) [29]: This measures the quality of fine features, edges, and spatial blurring in the images and is defined as:

[TABLE]

where LoG is a Laplacian of Gaussian filter that captures edges. We use the same filter specifications as in Ravishankar et al. [29]: kernel size of 15 $\times$ 15 pixels, with a standard deviation of 1.5.

• •

The Structural Similarity index (SSIM) is a perceptual metric introduced by Wang et al. [30]. We used the toolbox introduced by Wang et al. [30]: with default contrast values, Gaussian kernel size of 11 $\times$ 11 pixels with a standard deviation of 1.5 pixels.

• •

Global phase coherence (GPC) index [31] provides a measure of image sharpness by estimating the volume of all possible phase functions associated with the measured modulus, which produces images that are not less likely than the original image. The likelihood is measured with the total variation implicit prior, and is numerically evaluated using a Monte-Carlo simulation. We used the toolbox introduced by Blanchet et al. [31] to compute this index for our images.

Higher values of the above-mentioned performance metrics correspond to better reconstruction, except for the HFEN, where a lower value is better.

IV-C Algorithms for comparisons

We have used both simulated and in-vivo data to compare the following algorithms:

• •

SToRM: Two-Step (Proposed): The manifold Laplacian is estimated iteratively by alternating between the estimation of the Laplacian matrix and the update of the images on the low-resolution data. Once the Laplacian is obtained, the high-resolution images are recovered by solving (9) using all the k-space samples.

• •

SToRM: Single-step (Full resolution): we alternate between (9) and (10) on the high-resolution data. This approach may offer improved quality than SToRM:Two-Step since the high-resolution details can potentially yield improved estimation of the Laplacian, and consequently improved results. However, this approach is associated with higher computational complexity.

• •

SToRM:Original: We estimate the Laplacian matrix from ground truth data. This approach is only possible in the simulated setting, and provides an upper bound for the image quality.

• •

SToRM: SENSE: In this method, we estimate the Laplacian matrix from the CG-SENSE reconstructions. Equation (12) is modified as:

[TABLE]

The manifold Laplacian is recovered from $\mathbf{X}_{S}$. This Laplacian matrix is then used to recover the high-resolution images by solving (9) from all the k-space samples.

• •

SToRM: Self-Nav [18]: The manifold Laplacian is recovered from the self-gating navigators acquired in k-space, followed by (9) using all the k-space samples.

• •

XD-GRASP [15]: This self-gated strategy estimates the cardiac and respiratory phases from the center sample of k-space regions [15]. It estimates the cardiac and respiratory signals by filtering the central regions with different band-pass filters, each corresponding to the cardiac and respiration frequencies. We used the author-provided MATLAB code for XD-GRASP implementation [15].

• •

k-t SLR [6]: The image time series is recovered by Schatten p-norm (p=0.5) and total variation regularization minimization. We have tuned the sparsity and low rank regularization parameters to get the optimum results on our dataset.

• •

Low-Rank [5, 6]: The image time series is recovered by nuclear norm minimization. The nuclear norm minimization approach models the images as points living on a subspace and we are setting the sparsity regularization parameter to zero.

• •

Compressed sensing [4]: The image time series is recovered by $l_{1}$ sparsity regularization with total variation.

For image quality and ejection fraction comparison, all data sets were assessed by a board certified radiologist and a cardiovascular imager in a blinded manner. Image quality was evaluated on a 5-point score ranging from 1 (poor and not acceptable clinically) to 5 (excellent clinically). The ejection fraction comparison between breath-held and SToRM: Two-Step results, was performed using a 2-way analysis of variance (ANOVA) analysis. Statistical analysis was performed using SAS software (version $9.4$; SAS Institute Inc., Cary, NC).

V RESULTS

The illustrations in Fig. 1 show the benefit of the iterative strategy in SToRM:Two-Step. The first iteration shows the SToRM: SENSE, where the Laplacian is estimated from SENSE reconstructions using the central k-space regions. Full-resolution reconstruction is obtained using this Laplacian and all of the k-space data. We observe that this results in residual aliasing artifacts. With the SENSE reconstruction and the Laplacian estimated from it as the initial guess, the kernel low-rank algorithm is run with different numbers of iterations, as shown in Fig. 1 on the central k-space regions. The Laplacian matrix estimated from these iterations, was used to obtain the high resolution reconstruction. We note that the image quality improves significantly with iterations. In particular, the Laplacian estimated from the fifth iteration yields improved reconstructions with reduced artifacts in the liver regions and minimal myocardial blurring.

We first compare SToRM:Two-Step, and SToRM: Self Nav with breath-held cine in Fig. 2 to determine if the performance of the matrix completion scheme is comparable to the setting with navigators and breath-held data. The data was acquired with navigated acquisition, where one spiral acquisition was repeated after every 5 readouts. The navigators were not included in the estimation of the Laplacian matrix in the SToRM:Two-Step. We manually identified a cardiac cycle from the SToRM reconstructions, which closely matched in the end inspiration phase in which breath-held data was acquired. Three frames (end-diastole, mid-frame, end-systole) from the image series are shown. We note that the proposed scheme provides similar visual quality to the breath-held acquisitions. The experiments also show that the visual quality of the SToRM:Two-Step scheme is quite comparable to that of StoRM: Self Nav.

We compare the proposed scheme against competing methods on a numerical phantom in Fig. 3 and Table I. Fig. 3 shows the visual comparison between the proposed SToRM:Two-Step, low-rank, SToRM: SENSE, compressed sensing (CS), kt-SLR and SToRM: Self-Nav methods. We observe that the proposed scheme significantly reduces the spatial and temporal blurring compared to the low-rank scheme CS and kt-SLR, which demonstrates the ability of the kernel low-rank algorithm in capturing non-linear redundancies. These visual observations are also confirmed by the quantitative results in Table I. We have used four metrics (SER, SSIM, HFEN, GPC) to evaluate the performance. We optimized the parameters for one dataset. We observe that the performance of the SToRM:Two-Step scheme is significantly better than other methods, while it is marginally higher than SToRM:Self-Nav.

We have also investigated the impact of different Laplacian estimation strategies as shown in the supplementary figure S1. Results show that the image quality is not different notably, however, we get lower computational complexity with the low-resolution approach, as mentioned in the section III-C.

We qualitatively compare the proposed algorithm with competing methods on six single-slice experimental datasets in Fig 4-6. The visual comparison against the low-rank and the SToRM:SENSE algorithm on one of the datasets is shown in Fig 4. We have picked three frames (end-diastole, mid-frame, end-systole) from the image series to show the spatial quality of the proposed method as compared to the low-rank method and SToRM: SENSE. The fourth column of Fig. 4 shows the temporal profiles. We observe that the proposed scheme provides improved visual quality with reduced spatial and temporal blurring, when compared to the low-rank method and SToRM: SENSE scheme. The comparison of the proposed scheme against the XD-GRASP algorithm on the same dataset is shown in Fig. 5. Since XD-GRASP and the proposed scheme use different reconstruction methodologies, we illustrate the results in two ways. In Fig. 5, we display the reconstructed images in the different cardiac and respiratory phases in the bottom rows, which are identified by XD-GRASP. We create a cine movie by picking each image in the time series from the XD-GRASP reconstructions depending on the specified cardiac and respiratory phase; the time profile in the top row corresponds to a cut along the myocardium identified by the blue line in one of the images. We also display the rows of the weight matrix $\mathbf{W}$ corresponding to two frames, identified by the yellow and red arrows in the top right columns. As discussed earlier, the weights indicate similarity of the specific frame with other frames in the dataset. The top row in Fig. 5.a shows the temporal profile of XD-GRASP, while the corresponding temporal profiles of the proposed scheme are shown in Fig. 5.b. We note that the cardiac and respiratory phases identified by XD-GRASP are roughly in agreement with the motion patterns in the temporal profiles of the proposed scheme in Fig. 5.b. However, the motion patterns in the temporal profiles in XD-GRASP appear attenuated. We also observe sharp transitions in contrast between frames from different cardiac/respiratory phases. The top rows of Fig. 5.(b) show the weights corresponding to two frames in end inspiration and end expiration, respectively. Note that the weight patterns agree reasonably well visually with the identified cardiac respiratory phases. The weights indicate soft-binning of the phases offered by the proposed scheme. The bottom rows of Fig. 5 show the reconstructed images arranged in the cardiac and respiratory phases, which were identified using the self-gating strategy in XD-GRASP. We note that similar binning can be performed using the eigenvectors of the Laplacian matrix as shown by Poddar et al. [18]. Fig. 5.a shows the recovered images using XD-GRASP, while Fig. 5.b corresponds to the proposed method. We note that some of the phase images are blurred in the XD-GRASP reconstructions. These phases correspond to the poorly sampled cardiac and respiratory phases. The soft-binning offered by the weighting strategy allows for more data-sharing between the phases, resulting in reduced myocardial blurring and improved fidelity of the temporal profiles.

In Fig. 6, we show the image quality scores of different methods, which are rated by two experts in a blinded fashion. These results show that the image quality of SToRM:Two-Step is better than competing methods such as XD-GRASP and low-rank. We compare the scores using ANOVA with a p value of $p=0.0001$, which show that SToRM:Two-Step and XD-GRASP are significantly different from low-rank and SToRM:SENSE, while the improvement in quality score of SToRM:Two-Step over XD-GRASP is not statistically different. However, both the raters consistently scored SToRM:Two-Step over XD-GRASP, except one expert on one of the datasets.

We compare the whole-heart breath-held SSFP data and SToRM:Two-Step on five whole heart datasets in Fig. 7 and Table II. Fig. 7.(i) shows ten slices of SSFP based breath-held data, from base to apex of the heart, while (ii) shows the corresponding free-breathing SToRM:Two-Step reconstructions. The image quality scores as well as ejection fraction calculations, by two experts in a blinded fashion, are reported in Table II. The quality scores in Table II shows that the image quality of the proposed free-breathing strategy is slightly lower than breath-held SSFP. An ANOVA procedure with $p=0.07$ revealed that two scores are not significantly different. The ejection fraction estimated from the proposed SToRM:Two-Step scheme closely matches the breath-held SSFP method. A blocked ANOVA test showed no significant difference between breath-held and SToRM:Two-Step results ($p=0.93$), showing the accuracy of LVEF using our proposed method.

VI DISCUSSION

We introduced an iterative spiral-SToRM framework for the recovery of free-breathing and ungated cardiac images from 2-D spiral acquisition. The framework assumes the images to be on a smooth surface in high dimensions and relies on a kernel low-rank prior to recover the dataset. The main difference of our scheme from our prior work [21] is the use of an iterative kernel low-rank matrix completion algorithm and we are using $l_{2}$ smoothness regularization as it is giving us better results as compared to the $l_{1}$ smoothness regularization. The proposed approach eliminates the need for explicit k-space navigators and relies on variable-density spiral acquisitions, where the central k-space regions are acquired with higher sampling density. By eliminating the need for navigators, the proposed scheme improves sampling efficiency and hence image quality. The use of spiral trajectory improved sampling efficiency and contrast. Specifically, the increased TR and increased flip angle offer improved contrast compared to the low TR radial acquisitions. To improve computational efficiency, we rely on a two-step strategy. In the first step, we estimate low-resolution reconstructions as well as the Laplacian matrix from the central k-space region using a kernel low-rank optimization scheme. Once the Laplacian matrix is estimated, we solve for the high-resolution image from the entire k-space data using the manifold Laplacian. The benefit of using low resolution approach is shown in the Table III and in the supplementary figure. We have compared low and high resolution approaches using SER and processing time. Low resolution reconstruction reduces the computational complexity significantly. We also approximate the Laplacian using a few basis functions, which reduces the computational complexity and memory demand of the algorithm by an order of magnitude. We observe that the SToRM: Two-Step approach recovers 2D cine images with reduced spatial and temporal blurring in a short free-breathing self-gated acquisition, compared to low-rank and explicit binning strategies. We have also compared Laplacian matrices estimated from original images and reconstructed images. The SER obtained from original images $=26.1$, whereas SER from the reconstructed images is $25.8$

The gradient echo (GRE) acquisition schemes have few advantages for simplifying 3T cine imaging, even though SSFP sequences are typically used for cine imaging. The longer repetition time (TR) in the spiral trajectory provides inflow-enhancement of the LV blood pool; the resulting contrast is similar to the Cartesian SSFP imaging as compared to the shorter TR Cartesian GRE imaging. Furthermore, the spoiled GRE-based approach used for the acquisition is robust to banding artifacts, which SSFP methods are vulnerable to, without any frequency scout requirement. In addition, GRE schemes are less sensitive to eddy current artifacts caused by the large angular increment of the golden angle ordering [32]. With 16 seconds per slice, the whole heart is imaged in 3 minutes. The difference in blood pool-myocardium contrast between diastole and systole seen in SToRM:Two-Step are due to inflow effects associated with GRE acquisitions.

Our experiments in Fig. 5 show that the proposed scheme provides less-blurred reconstructions compared to XD-GRASP. As discussed previously, XD-GRASP relies on binning each image to appropriate cardiac/respiratory phases. We note that the time duration of heart in all of the cardiac/respiratory phases is not equal, with some phases (e.g. mid-systolic and inspiration phases ) having significantly fewer spokes than others. The recovery of these images from a few k-space points is significantly more challenging in XD-GRASP, which results in the residual blurring. By contrast, the SToRM strategy relies on a soft-gating strategy with no explicit binning. Our experiments in Fig. 5 show that this approach is more robust to residual blurring.

Our quantitative and qualitative experiments show that the SToRM:Two-Step scheme provides reconstructions that are marginally better in quality than SToRM:Self-Nav in the simulation settings. While both methods rely on the kernel low-rank algorithm, SToRM:Two-Step scheme does not require specialized k-space navigators to estimate the manifold Laplacian. The major benefit of this scheme is its application to golden angle radial or spiral sequences already in place on several scanners. In this work, we have not compared the two approaches on experimental MR data. Specifically, the experimental spiral data was acquired without any navigators, which makes it difficult to realize SToRM:Self-Nav. Thus, based on the current experiments, we cannot conclude that SToRM: Two-Step is superior to SToRM: Self-Nav; we plan to pursue these comparisons in our future work.

The comparisons with breath-held SSFP data in Fig. 7 and Table II show that the image quality and ejection fraction measures of the proposed scheme are comparable to breath-held acquisitions. However, the two sequences differ significantly in several aspects, including image contrast and bias-fields. We notice minor blurring of the myocardium in the slices close to the apex. This problem may be due to the inaccurate estimation of cardiac motion patterns from the navigators. Specifically, when the heart occupies a tiny fraction of the field of view as with apical slices, the in-plane navigators may not be sensitive to cardiac motion. In the future, we propose to extend SToRM with 3D navigators to minimize the above problem.

VII Conclusion

In this paper, we have proposed an iterative SToRM algorithm (SToRM:Two-Step) for the recovery of free-breathing and ungated cardiac MR images using spiral acquisitions with no k-space navigators. Our experiments show that the proposed scheme offers better performance compare to our previous STORM:Self-Nav method, which requires k-space navigators. Qualitative comparisons by experts also show that the proposed scheme provides less spatial and temporal blurring compared to low-rank methods, which do not require explicit binning to cardiac/respiratory phases, and XD-GRASP, which bins the data. Our preliminary study also shows that the proposed framework provides ejection fraction measures, which are statistically equivalent to the breath-held measurements. The MATLAB code for the SToRM-Iterative scheme can be downloaded from

https://github.com/ahaseebahmed/SpiralSToRM-Iterative.