Introduction

Parallel imaging^1-4 is a widely used technique to accelerate scans, for which SENSE³-based reconstruction approaches require explicit knowledge of coil sensitivity maps. Dynamic MRI reconstructions often use a single set of sensitivity maps across time, which is estimated either from a calibration scan or a temporal average of the k-space data⁵. However, these static maps suffer from inaccuracies since sensitivities change with respiratory and other subject motion, as shown in Figure 1.

Although previous works^6,7 have calculated sensitivities frame-by-frame, this approach is challenging for iterative reconstructions with temporal constraints (such as total variation or low-rank regularizations). These reconstructions do not recover images frame-by-frame, but instead, continuously iterate through sets of frames to enforce temporal constraints. Hence, sensitivities cannot be computed on-the-fly and should be stored prior to the reconstruction. However, frame-wise sensitivities occupy a large amount of memory (see Figure 2), which often prohibits their use, especially on memory-limited GPUs. Therefore, we propose Compact Maps, a method that solves for a compact representation of time-resolved sensitivities using a low-dimensional approach.

Method

Formulation
Since physical sensitivities are smooth in magnitude and phase, they have a compact representation in the spatial Fourier domain. To enforce this prior knowledge, we represent the maps as compact k-space convolution kernels. While maps change over time with motion, that change is smooth, and can be represented by a lower-dimensional temporal basis. Therefore, we define $$$(\hat{φ}_1,\hat{φ}_2,...\hat{φ}_N)$$$ as a temporal basis for the map kernels and $$$(\alpha_1,\alpha_2,...,\alpha_N)$$$ as corresponding basis coefficients (see Figure 3a, where $$$N$$$=3). To estimate these quantities, we perform a calibration using only central k-space data. The calibration jointly solves for the basis terms and a low-resolution image using an iterative SENSE⁸-like formulation, as shown below:$$(1)\hspace{0.5cm}\min_{m,\alpha,\hat{φ}}\sum_{i=1}^Tf(m,\hat{φ},\alpha)+R$$$$\textrm{where}\hspace{0.5cm}f(m,\hat{φ},\alpha)=\frac{1}{2}\|PF\{m_i\cdot(F^{-1}Z\sum_{j=1}^N\hat{φ}_j\alpha_{ij}^T)\}-y_i\|_2^2$$
Here, $$$m_i$$$ is the low-resolution image, $$$y_i$$$ is the low-resolution k-space data, $$$T$$$ is the number of time frames, $$$F$$$ is the 2D Fourier Transform, $$$Z$$$ is the zero-padding operation, $$$P$$$ is the projection onto the sampling pattern, and $$$R$$$ stands for the regularization terms, described further below.

Optimization
Although (1) is a non-convex problem, it is convex in each of $$$m$$$, $$$\hat{φ}$$$, and $$$\alpha$$$ when holding the other two quantities fixed. Therefore, we solve this problem with cyclic block coordinate descent, where each variable is updated individually by solving its corresponding convex subproblem. For large problems, we use stochastic optimization by randomly choosing a small batch of frames for each cycle (see Figure 3b).

The three subproblems for each variable update are as follows, where superscripts denote the cycle number:$$(2)\hspace{0.5cm}\alpha^{k+1}=\textrm{arg}\min_{\alpha}\sum_{i\in{batch}}f(m^k,\hat{φ}^k,\alpha)+\frac{\lambda_1}{2}\|D_t\alpha\|_2^2$$$$(3)\hspace{0.5cm}\hat{φ}^{k+1}=\textrm{arg}\min_{\hat{φ}}\sum_{i\in{batch}}f(m^k,\hat{φ},\alpha^{k+1})+\frac{\lambda_2}{2}\|W\hat{φ}\|_2^2+\beta^k\|\hat{φ}-\hat{φ}^k\|_2^2$$$$(4)\hspace{0.5cm}m^{k+1}=\textrm{arg}\min_m\sum_{i\in{batch}}f(m,\hat{φ}^{k+1},\alpha^{k+1})+\frac{\lambda_3}{2}\|m\|_2^2$$
In (2), $$$D_t$$$ represents a finite difference operator across time, which is used to impose temporal smoothness of the coefficients. Hence, the batches used for stochastic optimization must contain consecutive time frames. In (3), $$$W$$$ is a weighting matrix that enforces spatial smoothness by penalizing higher spatial frequencies in the basis components, inspired by the approach in ENLIVE⁹. The last term in (3) is a proximal point term¹⁰, which penalizes the difference between the updated variable value and its previous value to reduce noisy stochastic gradients. $$$\beta_k$$$ is increased for each epoch to gradually decrease the learning rate.

We use 3 iterations of conjugate gradient descent for each subproblem. Every time $$$\hat{φ}$$$ is updated, the Gram-Schmidt procedure is performed to orthonormalize the basis. We initialize $$$m$$$ with 1s, $$$\hat{φ}$$$ as 0s, and $$$\alpha$$$ as 1s. The optimization is warm-started by first using an 8x8 kernel size and $$$N$$$=1, then incrementally ramping these quantities up to the desired values with each warm start step. The described optimization was implemented in Python using the SigPy¹¹ package.

Evaluation
With IRB approval, we acquired a 2D axial bSSFP dataset from a healthy volunteer on a GE 3T MR750W (3.5ms TR, 1.4x1.4mm spatial resolution, 0.5s temporal resolution). The dataset had 150 time frames, 18 channels, and 192x154 image matrix size. Retrospective undersampling of R=3 was performed by uniformly discarding phase-encode lines outside a calibration region of 20 lines. For the optimization, we used $$$N$$$=4 basis components, 24x24 kernel size, batch size of 5 frames, and 72x72 of central k-space data. The stochastic optimization was run for 20 epochs, after 5 epochs for each of three warm start steps. On a Nvidia Titan-Xp GPU, the calibration pipeline took ~500MB GPU memory and ~2min runtime.

Results

To evaluate the quality of Compact Maps, we compare them against time-averaged and frame-wise maps calculated using ESPIRiT⁴. In Figure 4, we use the calibration error metric¹², which is the residual after projecting fully-sampled data onto the span of normalized sensitivities. Compact Maps achieve much lower calibration error than do time-averaged maps and perform similarly to frame-wise maps.

For Figure 5, images are reconstructed using iterative SENSE with l₂-regularization of λ=0.001. For time-averaged maps, residual aliasing artifacts are visible in the reconstructed images and in the difference with the reference, but are significantly reduced in those for Compact Maps. Compact Maps also result in a lower normalized root-mean squared error (NRMSE).

Conclusion

We have introduced a method that solves for a compact representation of time-resolved coil sensitivity maps to dramatically reduce their size in memory. Compact Maps can enable improvements in image quality for dynamic MRI reconstructions.

Acknowledgements

We thank the following funding sources: NIH R01EB009690, NIH R01HL136965, NIH R01EB026136, NIH U01 EB025162, and GE Healthcare.