Marcelo V. W. Zibetti^{1}, Gabor T. Herman^{2}, and Ravinder R. Regatte^{1}

^{1}Radiology, NYU, New York, NY, United States, ^{2}The Graduate Center, CUNY, New York, NY, United States

Compressed Sensing (CS) and parallel MRI (pMRI) have been successfully applied to accelerate MRI-data acquisition. CS requires incoherence, usually achieved by random undersampling the data, but pMRI does not. Combined, these methods allow even higher acceleration rates. However, it is unknown how the sampling pattern (SP) should be selected. It is also unknown if the SP is dependent on the reconstruction method. Here we demonstrate, using a new algorithm, that the SP can be learned from given data and reconstruction method. Our results show that the learned SP is superior to others such as Poisson disk and variable density.

Here, we propose an algorithm to learn a good SP from given data and reconstruction method, similarly to [6], but with CS-pMRI images. We focus on T1ρ-weighted images of human knee joint when sparse and low rank reconstructions are utilized [7]. Our preliminary results show that the learning-based SP is superior to others such as Poisson disk [8] and variable density [9] for human knee joint.

$$\hat{\mathbf{x}}=\arg\min_\mathbf{x} \left( ||\mathbf{m}- S_ΩFC \mathbf{x}||_2^2+λP(\mathbf{x}) \right)\approx R(\mathbf{m},Ω).$$

Here $$$\mathbf{x}$$$ represents the 2D+time images, of size $$$N_x\times N_y \times N_t$$$ (in our experiments this is $$$128 \times 64 \times 10$$$) which denotes vertical $$$N_x$$$ and horizontal $$$N_y$$$ sizes and time $$$N_t$$$. $$$\mathbf{m}$$$ is the undersampled multicoil k-t-space data. $$$C$$$ denotes the coil sensitivities transform, which maps $$$\mathbf{x}$$$ into multicoil-weighted images of size $$$N_x \times N_y \times N_t \times N_c$$$, with number of coils $$$N_c$$$. $$$F$$$ represents the spatial FFTs, which are $$$N_t \times N_c$$$ repetitions of the 2D-FFT, $$$S_Ω$$$ is the sampling function using SP $$$Ω$$$ (same for all coils) and $$$λ$$$ is the regularization parameter, see details in [7]. The SP contains the k-t-space points to be sampled; it can be displayed as a binary mask with dimensions $$$N_k=N_x \times N_y \times N_t$$$, as shown in Figures 2 to 5. The regularization functions considered are: $$$l_1$$$-norm, as $$$P(x)=||Tx||_1$$$, where $$$T$$$ is the spatiotemporal finite differences (STFD); and low rank (LR), using nuclear-norm of $$$\mathbf{x}$$$ (reshaped as a matrix), given by $$$||\mathbf{x}||_*$$$. We use the iterative algorithm MFISTA-VA [11], to obtain $$$R(\mathbf{m},Ω)$$$, that approximates the above minimization.

Our proposal is similar to [6], but using:

$$\hat{Ω}=\arg\min_Ω\sum_{i=1}^{N_i}||\bar{\mathbf{m}}_{i} –FCR(S_Ω\bar{\mathbf{m}}_{i},Ω)||_2^2 \approx O\left(Ω_{init},R, \{\bar{\mathbf{m}}_{i} \}_{i=1}^{N_i},K,\alpha,L \right).$$

$$$N_i$$$ is the number of images used for learning. We use fully-sampled pMRI data $$$\bar{\mathbf{m}}_{i}$$$, of size $$$N_x \times N_y \times N_t \times N_c$$$ (size $$$128 \times 64 \times 10 \times 15$$$ in the experiments), and $$$O\left(Ω_{init},R, \{\bar{\mathbf{m}}_{i} \}_{i=1}^{N_i},K,\alpha,L \right)$$$ is the algorithm described in Figure 1. Note that the error, denoted by $$$\mathbf{e}_i=[\mathbf{e }_{i,1}^T \ldots \mathbf{e}_{i,N_c}^T ]^T$$$ in the algorithm, is evaluated in the fully-sampled multicoil k-t-space, where $$$\mathbf{e }_{i,c}$$$ is the k-t-space error for the coil $$$c$$$, and $$$i$$$ represent the index of the image in the training set. It is assumed only a fixed fraction of the k-t-space points belongs to $$$Ω$$$, where the acceleration factor (AF) is the ratio of the total number of points to the number of samples. The optimization finds the k-t-space points to be sampled in order to have the approximately best results from the given reconstruction method and dataset.

This combinatorial problem is solved by $$$L$$$ iterations (90 in the experiments) of the proposed algorithm, as described in Figure 1, instead of using slow greedy approaches [6]. At each iteration $$$l$$$, it searches among the k-t-space points for the K elements with the largest mean squared error $$$\boldsymbol{\varepsilon}$$$ and smallest regularized error $$$\mathbf{r}$$$ to decide which points should be included and removed, as described in Figure 1 and illustrated in Figure 2. The parameters used in the experiments are: initial $$$K=$$$410, and $$$\alpha=$$$0.1. The proposed algorithm can be applied to any initial SP. We used $$$N_i=$$$40 for training and a validation set with 12 images that are not in the training set.

In Figure 4 the optimized SP for different regularization functions used in $$$R(S_Ω\bar{\mathbf{m}}_{i},Ω)$$$ are shown, one can also note that the optimized SP are different, even with the same initial guess. This difference is also expected since different reconstruction methods recover the missing k-t-space points differently.

In Figure 5, we show some visual results and NRMSE for AF=16. The NRMSE improvements ranged from 10.3% to 38.8% in the validation images.

[1] R. Otazo, D. Kim, L. Axel, and D. K. Sodickson, “Combination of compressed sensing and parallel imaging for highly accelerated first-pass cardiac perfusion MRI,” Magn. Reson. Med., vol. 64, no. 3, pp. 767–776, Sep. 2010.

[2] M. Murphy, M. Alley, J. Demmel, K. Keutzer, S. Vasanawala, and M. Lustig, “Fast $$$\ell_1$$$-SPIRiT Compressed Sensing Parallel Imaging MRI: Scalable Parallel Implementation and Clinically Feasible Runtime,” IEEE Trans. Med. Imaging, vol. 31, no. 6, pp. 1250–1262, Jun. 2012.

[3] M. Lustig, D. L. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging.,” Magn. Reson. Med., vol. 58, no. 6, pp. 1182–95, Dec. 2007.

[4] M. Blaimer, F. Breuer, M. Mueller, R. M. Heidemann, M. A. Griswold, and P. M. Jakob, “SMASH, SENSE, PILS, GRAPPA,” Top. Magn. Reson. Imaging, vol. 15, no. 4, pp. 223–236, Aug. 2004.

[5] K. P. Pruessmann, M. Weiger, M. B. Scheidegger, and P. Boesiger, “SENSE: Sensitivity encoding for fast MRI,” Magn. Reson. Med., vol. 42, no. 5, pp. 952–962, Nov. 1999.

[6] B. Gözcü, R. K. Mahabadi, Y. H. Li, E. Ilicak, T. Çukur, J. Scarlett, and V. Cevher, “Learning-Based Compressive MRI,” IEEE Trans. Med. Imaging, vol. 37, no. 6, pp. 1394–1406, 2018.

[7] M. V. W. Zibetti, A. Sharafi, R. Otazo, and R. R. Regatte, “Accelerating 3D-T 1ρ mapping of cartilage using compressed sensing with different sparse and low rank models,” Magn. Reson. Med., vol. 80, no. 4, pp. 1475–1491, Oct. 2018.

[8] D. Dunbar and G. Humphreys, “A spatial data structure for fast Poisson-disk sample generation,” ACM Trans. Graph., vol. 25, no. 3, p. 503, Jul. 2006.

[9] J. Y. Cheng, T. Zhang, M. T. Alley, M. Lustig, S. S. Vasanawala, and J. M. Pauly, “Variable-Density Radial View-Ordering and Sampling for Time-Optimized 3D Cartesian Imaging,” Proc. ISMRM Work. Data Sampl. Image Reconstr., p. 3, 2013.

[10] M. V. W. Zibetti, A. Sharafi, R. Otazo, and R. R. Regatte, “Compressed sensing acceleration of biexponential 3D‐T1ρ relaxation mapping of knee cartilage,” Magn. Reson. Med., vol. 81, no. 2, pp. 863–880, Feb. 2019.

[11] M. V. W. Zibetti, E. S. Helou, R. R. Regatte, and G. T. Herman, “Monotone FISTA With Variable Acceleration for Compressed Sensing Magnetic Resonance Imaging,” IEEE Trans. Comput. Imaging, vol. 5, no. 1, pp. 109–119, Mar. 2019.

Figure 1: The proposed algorithm for
learning the sampling pattern is described here.

Figure 2: Representative illustration of the
steps of the algorithm. The sampling patterns illustrated here correspond to
the first time point of the k-t data sequence.

Figure 3: Examples of initial sampling pattern
and final sampling pattern with convergence curves for AF=10. Note that the
convergence depends on the initial guess, but the NRMSE of the optimized
sampling patterns are similar and better than the initial patterns. Note that
not necessarily the best initial pattern (variable density) led to the best
optimized pattern. The reconstruction algorithm uses low rank regularization in
this example.

Figure 4: Examples
of the optimized sampling patterns using the same data set, same initial
sampling pattern but comparing two reconstructions: one using sparsity in the
spatiotemporal finite domain (STFD) and other using low rank (LR). Both with
AF=10.

Figure 5:
Example of reconstructed images of the validation set, with their NRMSE, using Poisson
disk and variable density sampling patterns, compared to optimized sampling
pattern (SP). All SPs with AF=16, for CS with spatiotemporal finite domain
(STFD) and low rank (LR). The NRMSE is better with optimized SP in both
methods.