3531

Patch-Tensor Low-n-Rank Reconstruction for Oscillating Steady State fMRI Acceleration

Shouchang Guo¹ and Douglas C. Noll²

¹Electrical and Computer Engineering, University of Michigan, Ann Arbor, MI, United States, ²Biomedical Engineering, University of Michigan, Ann Arbor, MI, United States

Synopsis

Oscillating steady-state imaging is a new acquisition method for T2*-weighted functional MRI that offers very high SNR, but longer acquisition times. The oscillations are highly reproducible, which make low-rank models suitable. In this work, a sparse sampling scheme combined with a patch-based low-rank tensor reconstruction is introduced to speed the image acquisitions. The low-n-rank algorithm was applied to oscillating steady state data to demonstrate the utility of this approach for functional MRI, demonstrating a 17-fold speed up with error levels less than 3%.

Introduction

The goals of functional MRI (fMRI) acquisition include high spatial and temporal resolution while maintaining a high signal to noise ratio (SNR). Oscillating steady state (OSS) acquisition is a new approach for fMRI acquisition that offers higher SNR but does so at the expense of imaging time. In this work, we propose to reduce imaging time through the use of undersampling in k-space and a novel reconstruction approach based on a well-generalizable patch-tensor low-n-rank model almost fully recovers OSS images from very few k-space samples. Combined with parallel imaging and non-uniform Fourier transform (NUFFT), the model uses a practical spiral random sampling pattern and can achieve as much as a 17-fold acceleration.

Methods

OSS imaging is a fundamentally new fMRI method similar to balanced steady-state free precession in the sense that multiple images are collected with RF cycling, except that it uses quadratic phase increments. Depending on parameters, the resultant signals are T2*-weighted and suitable for fMRI, but oscillate with a periodicity dictated by the phase progression.

There is past work on the use of low-rank (LR) models or LR and sparse (e.g. LR+S) models for fMRI ¹. It uses space and time as the two dimension in the LR model. In the current work, there are two time dimensions, the oscillating or fast time dimension, and the time series or slow time dimension, which is a good match for patch-based low-rank tensors. We proposed to process OSS fMRI data as low-n-rank tensors, which imposes low-rank on all the tensor unfoldings ². Since local information is more likely to be low-rank, the original tensor is partitioned into patch-tensors, which are further reshaped to satisfy the low-n-rank assumption. Using this non-direct low-n-rank structure of the original data greatly improves model adaptivity and potentially benefits random sampling pattern selection.

The convex optimization problem based on patch-tensor low-n-rank is expressed as

$$\text{minimize }\sum_{\mathcal{P^*}} \sum_{i = 1}^{N}\lVert\mathcal{P}(\mathbf{X})_{(i)}\rVert_* + \frac{\lambda}{2}\lVert\mathcal{A}(\mathbf{X})-b\rVert_2^2 ,$$

where $$$\mathbf{X}\in\mathbb{R}^{I_x\times I_y\times I_t}$$$ represents the OSS fMRI images to be recovered, $$$\mathcal{P}(\mathbf{\cdot})$$$ divides the original tensor into patches and reshapes them to low-n-rank tensors. $$$\mathcal{P}(\mathbf{X})\in\mathbb{R}^{I_{xy}\times I_{ft} \times I_{st}}$$$, where $$$I_{xy}$$$ is a vectorized spatial dimension, $$$I_{ft}$$$ is the oscillating acquisition dimension and equals the number of TR's for the oscillation, $$$I_{st}$$$ denotes a temporal block of the functional signal. $$$\mathcal{P}(\mathbf{X})_{(i)}$$$ is the mode-$$$n$$$ unfolding of $$$\mathcal{P}(\mathbf{X})$$$. $$$\mathcal{A}(\mathbf{\cdot})$$$ can be composed of coil sensitivities and NUFFT. $$$b$$$ contains undersampled k-space measurements. We used two kinds of retrospective sampling patterns. The main pattern is given in Fig. 1 (a). It takes an undersampled variable-density spiral for each temporal frame and rotates with the golden-angle between frames. The other is Poisson-disk sampling without multi-coil to test the model capacity. The patch-tensors are treated as independent modules, and the updating algorithm derived from Alternating Direction Methods of Multipliers (ADMM) ³ is given in Fig. 2.

All the data were collected by the OSS acquisition, the parameters includes TR = 15 ms, $$$I_{ft} = 10$$$, FOV = 22 cm, and slice thickness of 4 mm. The human experiments were performed with a finger tapping task repeated 4 cycles, and the activation maps were generated by correlation with the reference waveform.

Results

The model achieved at least 7-fold acceleration with a single-coil acquisition and Poisson disk sampling (Fig. 3). It also helped on sparse noise reduction of the functional time course. Thought the use of multi-coil and the proposed spiral sampling pattern, the original phantom images were recovered with less than 2% error given 5.83% of the k-space samples (Fig. 4). In the case of multi-coil functional human data, the error was less than 3% and the activation map was also recovered (Fig. 5).

Conclusion

A novel fMRI reconstruction model based on patch-tensor low-n-rank is proposed. It uses a spiral random sampling pattern that is easy to implement and can recover OSS images with less than 6% of the original k-space samples. By fully exploring local spatial-temporal relationships as demonstrated, the algorithm shows a strong capacity for data completion. It can be generalized to other image reconstruction problems as long as the non-direct patch-tensor low-n-rank assumption is satisfied, which is easier to satisfy and more flexible than assumptions enforced on the whole dataset.

Acknowledgements

This work was supported by NIH Grant R01EB023618.

References

1. Chiew, M., Graedel, N. N., McNab, J. A., Smith, S. M. and Miller, K. L. (2016), Accelerating functional MRI using fixed-rank approximations and radial-cartesian sampling. Magn. Reson. Med., 76: 1825–1836.

2. T. G. Kolda and B. W. Bader, "Tensor decompositions and applications", SIAM Review, 2009, pp. 455-500.

3. S. Boyd, N. Parikh, E. Chu, et al, “Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers,” Foundations and Trends in Machine Learning, vol. 3, no. 1, pp. 1-122, 2010.

Figures

Figure 1. The model uses very limited k-space samples compared to the fully sampled condition (a) and recovers the oscillating steady-state images (b) by enforcing low-rank on all three unfoldings of the reshaped tensor (c), which has two time dimensions and one spatial dimension. Fig. 1. (a) also shows how the variable-density, undersampled spiral rotates with the golden angle between frames.

Figure 2. The main steps of the algorithm are given. The parameters $$$k$$$, $$$\lambda$$$, $$$\beta$$$ were tuned empirically. $$$N = 3$$$ is the number of unfoldings or the dimension of the patch-tensor. The updates treat patch-tensors independently, so the computation can be easily parallelized.

Figure 3. The Poisson disk sampling was performed with a densely sampled core (upper right). The chosen patch-tensor size is 64 x 10 x 66 (space x fast time x slow time). The activation map obtained from the recovered images looks very similar to the original fully sampled result (left column). The time course recovered has a lower spike (green arrow) compared to the fully sampled, which indicates that the patch-tensor low-n-rank model can remove sparse noise without sacrificing accuracy.

Figure 4. The phantom data was acquired with a 32-channel coil and the proposed spiral undersampling. Patch-tensor size 64 x 10 x 40 (space x fast time x slow time) was chosen. The reconstruction result using the proposed algorithm is compared to SENSE reconstruction with a roughness penalty. The model can almost fully recover the original images.

Figure 5. This set of human data was collected with multi-coil and processed with the spiral sampling. Patch-tensor size 16 x 10 x 44 (space x fast time x slow time) was selected. The decreased spatial dimension often leads to better reconstruction performance. The activation maps (lower panel) used a contiguity threshold of 2, and proposed algorithm presents a much better activation map than the SENSE reconstruction.

Proc. Intl. Soc. Mag. Reson. Med. 26 (2018)

3531