Blind compressive sensing dynamic MRI

Abstract

We propose a novel blind compressive sensing (BCS) frame work to recover dynamic magnetic resonance images from undersampled measurements. This scheme models the dynamic signal as a sparse linear combination of temporal basis functions, chosen from a large dictionary. In contrast to classical compressed sensing, the BCS scheme simultaneously estimates the dictionary and the sparse coefficients from the undersampled measurements. Apart from the sparsity of the coefficients, the key difference of the BCS scheme with current low rank methods is the nonorthogonal nature of the dictionary basis functions. Since the number of degrees-of-freedom of the BCS model is smaller than that of the low-rank methods, it provides improved reconstructions at high acceleration rates. We formulate the reconstruction as a constrained optimization problem; the objective function is the linear combination of a data consistency term and sparsity promoting l1 prior of the coefficients. The Frobenius norm dictionary constraint is used to avoid scale ambiguity. We introduce a simple and efficient majorize-minimize algorithm, which decouples the original criterion into three simpler subproblems. An alternating minimization strategy is used, where we cycle through the minimization of three simpler problems. This algorithm is seen to be considerably faster than approaches that alternates between sparse coding and dictionary estimation, as well as the extension of K-SVD dictionary learning scheme. The use of the l1 penalty and Frobenius norm dictionary constraint enables the attenuation of insignificant basis functions compared to the l0 norm and column norm constraint assumed in most dictionary learning algorithms; this is especially important since the number of basis functions that can be reliably estimated is restricted by the available measurements. We also observe that the proposed scheme is more robust to local minima compared to K-SVD method, which relies on greedy sparse coding. Our phase transition experiments demonstrate that the BCS scheme provides much better recovery rates than classical Fourier-based CS schemes, while being only marginally worse than the dictionary aware setting. Since the overhead in additionally estimating the dictionary is low, this method can be very useful in dynamic magnetic resonance imaging applications, where the signal is not sparse in known dictionaries. We demonstrate the utility of the BCS scheme in accelerating contrast enhanced dynamic data. We observe superior reconstruction performance with the BCS scheme in comparison to existing low rank and compressed sensing schemes.

Fig. 1

Comparison of blind compressed sensing (BCS) and blind linear model (BLM) representations of dynamic imaging data: The Casorati form of the dynamic signal Γ is shown in (a). The BLM and BCS decompositions of Γ are respectively shown in (b) and (c). BCS uses a large over-complete dictionary, unlike the orthogonal dictionary with few basis functions in BLM; (R > r). Note that the coefficients/spatial weights in BCS are sparser than that of BLM. The temporal basis functions in the BCS dictionary are representative of specific regions, since they are not constrained to be orthogonal. For example, the 1st, 2nd columns of U_M×R in BCS correspond respectively to the temporal dynamics of the right and left ventricles in this myocardial perfusion data with motion. We observe that only 4-5 coefficients per pixel are sufficient to represent the dataset.

Fig. 2

Comparison of different BCS schemes: In (a), we show the reconstruction error vs reconstruction time for the proposed BCS, alternate BCS, and the greedy BCS schemes. The free parameters of all the schemes were optimized to yield the lowest possible errors, while the dictionary sizes of all methods were fixed to 45 atoms. We plot the reconstruction error as a function of the CPU run time for the different schemes with different dictionary initializations. The proposed BCS and alternating BCS scheme converged to the same solution irrespective of the initialization. However, the proposed scheme is observed to be considerably faster; note that the alternating scheme takes around ten times more time to converge. It is also seen that the greedy BCS scheme converged to different solutions with different initializations, indicating the dependence of these schemes on local minima.

Fig. 3

Model coefficients and dictionary bases. We show few of the estimated spatial coefficients u_i(x) and its corresponding temporal bases υ_i(t) from 7.5 fold undersampled myocardial perfusion MRI data (data in Fig. 2). (a) corresponds to the estimates using the proposed BCS scheme, while (b) is estimated using the greedy BCS scheme. For consistent visualization, we sort the product entries u_i(x)υ_i(t) according to their ℓ₂ norm, and show the first 30 sorted terms. Note that the BCS basis functions are drastically different from exponential basis functions in the Fourier dictionary; they represent temporal characteristics specific to the dataset. It can also be seen that the energy of the basis functions in (a) varies considerably, depending on their relative importance. Since we rely on the ℓ₁ sparsity norm and Frobenius norm dictionary constraint, the representation will adjust the scaling of the dictionary basis functions υ_i(t) such that the ${\Vert \mathrm{U}\Vert }_{{\ell }_1}$ is minimized. Specifically, the ℓ₁ minimization optimization will ensure that basis functions used more frequently are assigned higher energies, while the less significant basis functions are assigned lower energy (see υ₂₅(t) to υ₃₀(t)), hence providing an implicit model order selection. By contrast, the formulation of the greedy BCS scheme involves the setting of ℓ₀ sparsity norm and column norm dictionary constraint; the penalty is only dependent on the sparsity of U. Unlike the proposed scheme, this does not provide an implicit model order selection, resulting in the preservation of noisy basis functions, whose coefficients capture the alias artifacts in the data. This explains the higher errors in the greedy BCS reconstructions in Fig. 2.

Fig. 4

Blind CS model dependence on the regularization parameter and the dictionary size: (a) shows the reconstruction error (ζ) as a function of different λ in the BCS model. (b) and (c) respectively show the reconstruction error (ζ) and the average number of non zero model coefficients of the BCS and the BLM schemes as a function of the number of bases in the respective models. As depicted in (a), we optimize our choice of λ such that the error between the fully sampled data and the reconstruction is minimal. From (b), we observe that the BCS reconstruction error reduces with the dictionary size and hits a plateau after a size of 20 basis functions. This is in sharp contrast with the BLM scheme where the reconstructions errors increase when the basis functions are increased. The average number of BCS model coefficients unlike the BLM has a non-linear relation with the dictionary size reaching saturation to a number of 4-4.5. The plots in (b) and (c) depict that the BCS scheme is insensitive to dictionary size as long as a reasonable size (atleast 20 in this case) is chosen. We chose a dictionary size of 45 bases in the experiments considered in this paper.

Fig. 5

The numerical phantoms Γ^j, which are used in the simulation study in figure 6. Here j is the number of non zero coefficients (sparsity levels) at each pixel. The top and bottom rows respectively show one spatial frame and the image time profile through the dotted white line. Note that the sparse decomposition provides considerable temporal detail even for a sparsity of one. This is possible since different temporal basis functions are active at each pixel.

Fig. 6

Phase transition behavior of various reconstruction schemes: Top row: Normalized reconstruction error ζ is shown at different acceleration factors (or equivalently different number of radial rays in each frame) for different values of j. Bottom row: ζ thresholded at 1 percent error; black represents 100 percent recovery. We study the ability of the algorithms to reliably recover each of the data sets Γ^j from different number of radial samples in kspace. The Γ^j, shown in Fig. 5 are the j sparse approximations of a myocardial perfusion MRI dataset with motion. As expected, the number of lines required to recover the dataset increases with the sparsity. The blind CS scheme outperformed the compressed sensing scheme considerably. The learned dictionary aware scheme yielded the best recovery rates. However due to a small over head in estimating the dictionary, the dictionary unaware (blind CS) scheme was only marginally worse than the dictionary aware scheme.

Fig. 7

Comparison of the proposed scheme with different methods on a retrospectively downsampled Cartesian myocardial perfusion data set with motion at 7.5 fold acceleration: A radial trajectory is used for downsampling. The trajectory for one frame is shown in (i). The trajectory is rotated by random shifts in each time frame. Reconstructions using different algorithms, along with the fully sampled data are shown in (i) to (v). (a-b), (c), (d-e), (f) respectively show few spatial frames, image time profile, corresponding error images, error in image time profile. The image time profile in (c) is through the dotted line in (i.b). The ripples in (i.c) correspond to the motion due to inconsistent gating and/or breathing. The location of the spatial frames along time is marked by the dotted lines in (i.c). We observe the BCS scheme to be robust to spatio-temporal blurring, compared to the low rank model; eg: see the white arrows, where the details of the papillary muscles are blurred in the Schatten p-norm reconstruction while maintained well with BCS. This is depicted in the error images as well, where BCS has diffused errors, while the low rank scheme (iii) have structured errors corresponding to the anatomy of the heart. The BCS scheme was also robust to the compromises observed with the CS scheme; the latter was sensitive to breathing motion as depicted by the arrows in iv.

Fig. 8

Comparisons of the different reconstructions schemes on a brain perfusion MRI dataset. The fully sampled data in (a) is retrospectively undersampled at a high acceleration of 10.66. The radial sampling mask for one frame is shown in (a), subsequent frames had the mask rotated by random angles. We show a spatial frame, the image time series, and the corresponding error images for all the reconstruction schemes. Note from (b,c), the low rank and CS schemes have artifacts in the form of spatiotemporal blur; the various fine features are blurred (see arrows). In contrast, the BCS scheme had crisper features, and superior spatiotemporal fidelity. The reconstruction error and the HFEN error numbers were also considerably less with the BCS scheme.

Fig. 9

Comparisons of different reconstruction schemes on a stress myocardial perfusion MRI dataset with breathing motion: Retrospective sampling was considered by picking 24 radial rays/frame from the acquired 72 ray data; the rays closest to the golden ratio pattern was chosen. Few spatial frames, the corresponding image time profile, error frames, and error in image time profile are shown for all the schemes. We specifically observe loss of important borders and temporal blur with the low rank and CS schemes while the blind CS reconstructions have crisper borders and better temporal fidelity. Also note from the columns d,e,f that the errors in the BCS scheme are less concentrated at the edges, compared to the other methods. This indicates that the edge details and temporal dynamics are better preserved in the BCS reconstructions.