Valery Vishnevskiy^{1} and Sebastian Kozerke^{1}

Adaptive spatiotemporal dictionaries offer improved reconstruction accuracy for dynamic cardiac MRI. However, most modern methods perform local encoding of image patches treating them independently. In order to increase reconstruction quality, we present a convex model that allows global control of encoding sparsity. The proposed method has a single tunable parameter and delivers 9% peak signal-to-noise ratio improvement of reconstruction compared to the state-of-the-art dictionary-based approach. Moreover, the implemented numerical scheme allowed 3-fold reconstruction time reduction.

Dynamic 2D cardiac MRI datasets were acquired on 1.5 and 3T systems (Philips Healthcare, Best, The Netherlands) with voxel size of 1.4x1.4x8 mm^{3 }and 25 dynamics, yielding 375 short-axis scans of size $$$N_d=128\times128\times25$$$ pixels. Following Caballero et al. [3], we use spatiotemporal patches of $$$4\times4\times4=N_p$$$ size and set the size of the dictionary to $$$N_a=300$$$. Since the dictionary is learned on fully sampled data, we used an efficient implementation of the generic online dictionary learning algorithm [4] provided by the Sklearn package [5], which iteratively solves the following optimization problem:

$$\min_{D\in\mathbb{R}^{N_p\times N_a}, U\in\mathbb{R}^{N_a\times N_s}} \frac12 \|DU − Y\|_{2,2}^2 + \alpha\|U\|_{1,1},\qquad\text{s.t.}\quad\|D_k\|_2=1,\quad k=1,\dots,N_a.\qquad(1)$$

Where $$$Y$$$ are training patches extracted from fully sampled data, $$$U$$$ is the matrix of patch encodings, $$$D$$$ is the learned dictionary, and the regularization parameter $$$\alpha=0.4$$$. We use the overcomplete discrete cosine transform as an initialization for $$$D$$$ and conduct 1000 iterations with 2000 patches in each batch. The resulting dictionary is shown in Figure 1. During reconstruction, the sparsity measure was relaxed to convex L1 norm, yielding the following constrained optimization problem:

$$\min_{U, X} \|U\|_{1,1} + \mathcal{I}[\|\mathcal{P}(X) - DU\|_{2,2} \leq \varepsilon],\quad\text{s.t.} \quad M_iFX_i=b_i, \quad i=1,\dots,N_t.\qquad(2)$$

Here indicator $$$\mathcal{I}[condition]=0$$$ when $$$condition$$$ is fulfilled and is $$$+\infty$$$ otherwise, $$$M_i$$$ and $$$b_i$$$ are undersampling matrices and k-spaces, respectively. The orthogonal linear operator $$$\mathcal{P}$$$ implements patch extraction (with circulant wrapping) from the given 2D+t image, while its adjoint $$$\mathcal{P}^T$$$ reconstructs an image from a patch representation by averaging of overlaps. Parameter $$$\varepsilon$$$ limits the global residual between the original signal and its sparse representation. Low values of $$$\varepsilon$$$ yield a sparser solution, while high values of $$$\varepsilon=10^{−5}N_d$$$ allow to account for noise in the k-space measurements $$$b_i$$$. To solve the problem (2) we cast it to the equivalent optimization problem

$$\min_{U,X,Z,T} \|Z\|_{1,1} + \mathcal{I}[\|T\|_{2,2} \leq \varepsilon],\quad\text{s.t.} \; M_iFX_i=b_i,\,i=1,\dots,N_t,\quad U = Z,\quad\mathcal{P}(X)-DU=T.\qquad(3)$$

Problem (3) can be solved efficiently with the ADMM algorithm [6] as described in Figure 2. We will further refer to the proposed reconstruction method as GSDR (Globally Sparse Dictionary based Reconstruction).

[1] Lustig M, Donoho D, Pauly JM. Sparse MRI: The application of compressed sensing for rapid MR imaging. Magnetic resonance in medicine. 2007;58(6):1182-95.

[2] Schmidt J, Kozerke S. Dynamic cardiac MR image reconstruction models using machine learning on large training data sets. ISMRM 2017.

[3] Caballero J, Price AN, Rueckert D, Hajnal JV. Dictionary learning and time sparsity for dynamic MR data reconstruction. IEEE transactions on medical imaging. 2014;33(4):979-94.

[4] Mairal J, Bach F, Ponce J, Sapiro G. Online dictionary learning for sparse coding. InProceedings of the 26th annual international conference on machine learning 2009. (pp. 689-696). ACM.

[5] Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, Blondel M, Prettenhofer P, Weiss R, Dubourg V, Vanderplas J. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research. 2011;12:2825-30.

[6] Boyd S, Parikh N, Chu E, Peleato B, Eckstein J. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends® in Machine Learning. 2011;3(1):1-22.

[7] Parikh N, Boyd S. Proximal algorithms. Foundations and Trends® in Optimization. 2014 Jan 13;1(3):127-239.

Figure 1: Illustration of the learned spatiotemporal dictionary for dynamic cardiac MR data using spatiotemporal patches of 4x4x4 size. Each temporal slice is mosaiced in the corresponding image.

Figure 2:ADMM algorithm for globally sparse GSDR MRI reconstruction. The $$$X$$$-update step involves trivial inversion of a diagonal matrix and weighting of the measured samples and current estimate in Fourier domain for every dynamic. The $$$Z$$$-update step is implemented with element-wise soft thresholding [7]. The $$$U$$$-update step involves solving a set of linear equations, which can be carried out efficiently by caching the Cholesky factorization of the left-hand side. The $$$T$$$-update step is implemented by proximal mapping onto L2-norm ball [7]. $$$S_X$$$, $$$S_Z$$$, $$$S_Y$$$ are the scaled dual variables.

Figure 3: Peak signal to noise ratio of reconstructions of the DLTG and the proposed GSDR algorithms for different undersampling factors. PSNR is calculated and averaged for 4 images that were scaled to [0, 1] interval.

Figure 4: Illustration of image reconstructions for the undersampling factor of 0.15 (acceleration x6.7). Column (a) shows fully sampled ground-truth (upper row) and temporal profile (bottom row) at the green dashed line. (b) and (c) show images reconstructed with DLTG and GSDR, and their respective errors in (d) and (e).

Figure 5: Local code length maps for DLTG and GSDR for reconstructions with an undersampling factor 6.7. Each point on the maps shows the number of dictionary atoms that were used to encode patch centered at the corresponding point.