A DivideandConquer Approach to
Compressed Sensing MRI
Abstract
Compressed sensing (CS) theory assures us that we can accurately reconstruct magnetic resonance images using fewer kspace measurements than the Nyquist sampling rate requires. In traditional CSMRI inversion methods, the fact that the energy within the Fourier measurement domain is distributed nonuniformly is often neglected during reconstruction. As a result, more densely sampled lowfrequency information tends to dominate penalization schemes for reconstructing MRI at the expense of highfrequency details. In this paper, we propose a new framework for CSMRI inversion in which we decompose the observed kspace data into “subspaces” via sets of filters in a lossless way, and reconstruct the images in these various spaces individually using offtheshelf algorithms. We then fuse the results to obtain the final reconstruction. In this way we are able to focus reconstruction on frequency information within the entire kspace more equally, preserving both high and low frequency details. We demonstrate that the proposed framework is competitive with stateoftheart methods in CSMRI in terms of quantitative performance, and often improves an algorithm’s results qualitatively compared with its direct application to kspace.
keywords:
compressed sensing, magnetic resonance imaging, divideandconquer1 Introduction
MRI is an important medical imaging technique because of its harmlessness and the high resolution information it measures in soft tissue. The data acquisition process directly measures Fourier coefficients of the object being imaged, which can then be recovered by an inverse Fourier transformation. A significant drawback of MR imaging is that data acquisition is relatively slow, meaning that a patient has to remain still for a long time to avoid producing motion artifacts. This is something especially difficult for children and those who are critically ill. Thus, accelerating the measurement speed while maintaining a diagnosticquality reconstruction is a major challenge in MR imaging.
Compressive sensing (CS) theory 1 (); 2 () has shown that it is possible to accurately reconstruct MR images with significantly fewer measurements than mandated by the Nyquist sampling theorem. The typical approach to CSMRI can be generalized as the optimization problem
(1) 
The goal is to reconstruct the vector , which corresponds to values in the MR image . The matrix is an undersampled Fourier basis used to directly measure the kspace data . In this objective, the first term enforces that the Fourier coefficients of agree with the observations . Since many vectors will satisfy this requirement, the penalty taking parameters searches for an with additional desired properties such as smoothness.
1.1 Related Work
Many CSMRI algorithms have been proposed in the framework of Equation (1). These developments mostly fall into two categories depending on whether they focus on new objective functions or on new efficient optimization algorithms.
Along the first line, new objective functions for CSMRI exploit the sparsity of MR images under different transform domains. Although medical images may not be sparse in the image domain, they can be projected onto a transform basis incoherent with the Fourier basis where they show a high degree of sparsity. Such bases include wavelets 3 (), total variation 3 (); 4 (); 5 (), contourlets 6 (), Walsh 7 () and PCA 8 (); 9 (); 10 (). Patchbased bases include directional wavelets (PBDWS) 30 (); 31 (), a graphbased redundant wavelet transform (GBRWT) 37 (), and dictionary bases constructed in situ using dictionary learning algorithms like KSVD 11 () and BPFA 12 (). As researches in deep learning methods thrive, the popular deep learning architectures are brought into CSMRI 38 (); 39 (); 40 (). In recent years, other approaches have enforced structural sparsity 13 (); 14 (); 41 (), nonlocal priors 15 (); 16 (), and approximations to the desired penalty such as the convex , FOCUSS, (for ), IRLS and smoothing functions 17 (); 18 (); 19 (); 20 ().
One particular method that we highlight and will use is the patchbased directional wavelet (PBDWS) 30 (); 31 (). This CSMRI objective function assumes that patches extracted and vectorized from the reconstructed MRI are sparse in the Haar wavelet domain. A key novelty is that each patch is vectorized in a way that depends on the geometric structure of the signal in that patch, and is chosen such that sparsity is maximized. In this way image details can be preserved better while satisfying the need for sparsity.
Another line of work has been devoted to finding more efficient ways to optimize the various objective functions arising from the framework of Equation (1), such as TVCMRI 21 (), RecPF 22 (), FISTA 23 (), FCSA 24 (), pFISTA 44 () (an algorithmic variation on SparseMRI 3 ()), Bregman 25 () and ADMM 12 (). These methods typically work by representing the objective function in a way that is easier to optimize iteratively.
1.2 Our contribution
In general, kspace data of MRI are not distributed uniformly across all frequency bands in energy and magnitude. As an example, we show the magnitudes of kspace data for an MRI of the brain in Figure 1. As is clearly evident, the energy is concentrated much more in the low frequency region of kspace than the high frequency region. This is a wellknown fact about MRI, which is considered in virtually all compressed sensing frameworks through the design of variable density sampling masks that sample more heavily in the low frequency part of kspace to ensure the basic structure is measured. However, as a result of using a squared error penalty term of the form as in the Equation (1), CSMRI inversion will result in these low frequency measurements dominating the reconstruction of . Since high frequency coefficients encode structural details such as edges and curves, the accumulation of error caused by overlooking high frequency information may not give the detail needed for diagnosis. In this paper we propose a general framework for reconstructing MRI measured with compressed sensing in a way that focuses more equally on all regions of kspace, and can be incorporated into any existing reconstruction algorithm. In particular, our approach breaks down the data fidelity squared error into various parts in a “divideandconquer” (DAC) manner and can incorporate any existing CSMRI inversion algorithm to reconstruct those parts.
The paper is organized as follows: In Section II, we discuss the nonuniformity property of the kspace data and its relationship to the reconstruction accuracy for CSMRI. To address this issue, we propose a threestep divideandconquer (DAC) framework in Section III. Section IV demonstrates the performance of our DAC framework on several MR data with different undersampling masks and ratios. We provide further discussion on parameter setting and noise characteristics in Section V.
1.3 Notation
We use the following notation: the fullysampled or reference MRI is denoted as with the vectorized form . Upper case means the representation in 2D while lower case means a vectorized form of the corresponding 2D representation. The subscript means fullysampled. Projecting these into the frequency domain, we represent the corresponding projections as and . The measured kspace has its unsampled positions padded with zeros, but the vectorized form will only have the sampled locations. The vector can be obtained by multiplying by an undersampled Fourier basis matrix . We further denote the reconstructed MRI image as and its kspace representation as . If we define a filter with limited spatial support, which is usually much smaller than the convolved image, the block Toeplitz matrix of the filter is defined as . Furthermore, the frequency response of filter is defined as .
2 Motivation: Nonuniform kspace density
Kspace data is known to be distributed nonuniformly in energy, as illustrated with the brain MRI in Figure 1. While kspace magnitudes tend to be larger in low frequency bands, much diagnostically important detail information is known to be in the high frequency regions. However, many existing CSMRI methods treat all errors equally, which tends to favor highmagnitude, lowfrequency information in reconstructions at the expense of the details.
To further analyze the phenomenon, we define the following two simple measures: the kspace absolute reconstruction error (KARE) and the kspace relative reconstruction error (KRRE). These are the same size as the kspace, with the element denoted as follows,
(2)  
(3) 
Here, is the fullysampled kspace data, and is the kspace of the MRI reconstructed from subsamples.
For the same brain MRI, we plot in Figure 2 the reconstructed image and the corresponding KARE and KRRE images using a CSMRI inversion method called PBDWS 31 () with undersampling and a 2D random mask. In Figures 2(a) and 2(b), we observe that data in the peripheral high frequency regions suffer almost the same absolute but larger relative errors. As an illustrative experiment under the assumption that we have access to the fullysampled kspace, one can manually redistribute the kspace error from PBDWS more evenly in such a way that the PSNR remains unchanged, as shown in Figures 2(c) and 2(d). However, as shown in Figures 2(e) and 2(f), after redistribution the reconstruction has much better visual quality, and a larger structural similarity measure (SSIM).
This simulation supports the claim that better high frequency kspace reconstruction accuracy, even at the expense of low frequency information, can lead to better visual quality. However, conventional CSMRI methods tend to leave the high frequency information deemphasized in their reconstruction objectives, coupled with heavier subsampling in lowfrequency regions. While in previous work such information is implicitly modeled, however, to our knowledge this problem has not been explicitly addressed for CSMRI inversion problems.
For example, in 27 () the low frequency kspace data is reconstructed first because of its dense distribution, then the reconstructed low frequency portions are padded back to the measurements for the secondstage reconstruction. The performance improvement of this method is limited due to fact that the reconstruction error of the low frequency bands propagate to the later reconstruction. The high and low frequency regions are separated using the rectangle box which can cause Gibbs effects and introduce artifacts. In 4 (); 5 () the horizontal and vertical differential images are reconstructed and then used as gradient constraints, where they mainly focus on the sparse nature of MRI in the gradient domain. In 26 (), a convolutional constraints is proposed, but the work nonuniformity in kspace. In 29 (), a method called HiSub CSMRI formalized a link between the kspace and wavelet domain to apply separate undersampling and reconstruction for high/low frequency kspace data. In the HiSub method, the high/low frequency regions in kspace are defined based on the separation of wavelet subbands; compressed sensing techniques are used for the high frequency region while parallel imaging is used for the lowfrequency region. The HiSub method relies on the specific sampling pattern and is not exclusively based on CSMRI ideas. In 28 (), the local scale mixture model is proposed to decompose the MR images into dual block sparse components: total variation for piecewise smooth parts and wavelet for residuals, but the decomposition only depends on the different priors between the total variation and wavelet in spatial domain.
The nonuniformity in kspace motivates us to reconsider the CSMRI problem using a divideandconquer (DAC) framework that can be implemented using existing inversion algorithms. Our method consists of three steps: subspace decomposition, subspace reconstruction and subspace integration. While the word “subspace” is welldefined for linear algebra, here we mean a specific frequency view into which the kspace measurements are decomposed using standard filtering techniques. This method allows for the algorithms to deal with the high and low frequency kspace data separately to better preserve fine structural details. Although the idea is simple, we note that the proposed subspace method exhibits great potential for recovering fine details in MRI by better preserving the high frequency information possibly improving the diagnosis quality of medical imaging applications.^{1}^{1}1We wish to emphasize here that our proposed method does not simply partition the kspace data into disjoint frequency regions.
3 Divideandconquer Subspace Framework
The proposed subspace method includes three steps: subspace decomposition, subspace reconstruction and subspace integration. We display a flowchart of the method in Figure 3. In this section we first discuss each of the three steps above separately. Then we connect these three steps to the objective function given in Equation (1). At the end of the section we summarize the proposed DAC CSMRI framework in Algorithm 1. Here we adopt the HoriVert subspace decomposition for illustration, which we will elaborate in later section.
3.1 Subspace decomposition
In order to reconstruct the corresponding image in each subspace, we need to first define the subspaces and then obtain the partially observed kspace data in each subspace. We call this process subspace decomposition. According to signal processing theory, a lossless decomposition can be accomplished using filters. In our case, we use a set of linear filters for their simplicity. We let be the impulse responses of the set of chosen filters, where is selected in advance. In principle, a decomposition in the image domain can be formulated as follows,
(4) 
where is the convolved image of the fullysampled in the th subspace and denotes the 2D convolution.
We equivalently work within the Fourier domain and obtain a series of frequency responses denoted as . As is well known, matrix convolution in the image domain can now be converted into an elementwise matrix multiplication operation in the Fourier domain. Thus Equation (4) is equivalent to
(5) 
where is the Fourier transformation of and denotes elementwise multiplication.
Using the undersampling, a zerofilled partial kspace view, , of ground truth can be written as
(6) 
where is the undersampling mask in which the sampled locations contain ones and unmeasured frequencies contain zeros. By simple algebra, the zerofilled, undersampled kspace data in each subspace can now be derived as
(7) 
Equation (7) indicates that the partial kspace data decomposed in each subspace can be derived via the elementwise multiplication between the frequency response of the filter and the original partially observed kspace data, or equivalently between the partially observed frequency response and the complete kspace of the MRI. We will use this latter observation in our reconstruction algorithm.
For the remaining derivation, we work with the problem by converting from matrix elementwise multiplication to matrixvector multiplication. Thus the impulse responses of the filters is rewritten as a circulant matrix . (Note that and differ in that the first operates on twodimensional kspace data while the second is on vectorized images.) An equivalent form of Equation (4) and Equation (5) can now be written as
(8) 
where is the vectorized form of . Similarly, Equation (7) can be written as
(9) 
while the original measured kspace of Equation (6) is
(10) 
We next return to the filter banks considered in Section 3.4.
3.2 Subspace Reconstruction
After dividing kspace into partial frequency views , it is intuitive to exploit this isolated information and reconstruct the corresponding images in each subspace separately. This will ensure that high frequency information, captured in certain by appropriate filter banks, is not sacrificed in favor of the far greater number of high magnitude, low frequency measurements. To this end, we define a subspacespecific optimization problem with an appropriate regularization term to be determined. The reconstruction of each subspace can be formulated by minimizing the following objective function
(11) 
where enforces desired properties of the reconstructed subspace image and is the regularization parameter for the data fidelity term. As with the filters, , these penalty functions can be chosen to be any CSMRI inversion algorithm. We note the new data fidelity squared error is proposed for better preservation for high frequency information, which also distinguish the proposed subspace method with merely regularizing the filtered subspace images with a unified data fidelity term. In this paper we use three recent stateoftheart CSMRI reconstruction methods: FCSA 24 (), WatMRI 14 () and PBDWS 30 (). To summarize, these penalties are the following:
(12) 
All methods use the squared error as a data fidelity term. In FCSA, the objective function is split into TV and wavelet regularization subproblem in an iterative manner, where each subproblem is solved efficiently using proximal gradient descend algorithm. Then the reconstructed MRI is obtained via the weighted average of the solutions of the two subproblems. In this sense it is an algorithmic development of the classic SparseMRI framework 3 (). WatMRI instead imposes wavelet tree group sparsity on the MRI via . The optimization is similar to FCSA, using splitting techniques. Both FCSA and WatMRI are CSMRI methods based on a global sparse regularization with nonadaptive transform bases. Thus they are suitable to stand for the CSMRI methods with fixed transform basis. PBDWS on the other hand is a patchbased method in which the MRI to be reconstructed is divided into patches and wavelets are used in a way that considers the geometric structure of the patch under consideration with the goal of maximizing sparsity. Thus PBDWS is representative of CSMRI algorithms that use an adaptive transform basis, but the proposed framework is not limited to using methods. We choose the three methods for illustrative purposes and because they are high quality algorithms.
We experiment with these three methods to show that our method can provide general improvement to many existing CSMRI models. In our experiments, we apply the same reconstruction algorithm to all subspaces, but with different parameter settings. If the MRI is divided into subspaces via , then the chosen algorithm would be run independently times, once for each subspace. Therefore, our framework increases computation time by a factor of , but the independence of each optimization allows for a straightforward parallelization.
3.3 Subspace Integration
Since the subspace decomposition is a linear decomposition using linear filters, if the decomposition is complementary but not redundant then integrating these results into a final reconstruction can be done by simply adding the images together,
(13) 
We also consider a Tikhonov regularization method by formulating the integration according to the following objective function
(14) 
This objective function admits a closedform solution, but direct computation is infeasible because of the high dimensionality of . However, by transferring the problem to the Fourier domain the reconstructed kspace is calculated elementwise followed by an inverse Fourier transform,
(15) 
Here, the division is elementwise, as is the magnitude of in the denominator. We discuss a method for determining each below in Section 3.5.
3.4 Filter banks for dividing and conquering
We consider two kinds of filter banks in this paper: one based on the horizontal and vertical redundant filter bank (HoriVert), and another based on the Gaussian complementary filter bank (Gaussian).
3.4.1 HoriVert subspace decomposition
Because much of the high frequency details in MRI can be represented as vertical or horizontal edges, we consider a decomposition of kspace into horizontal/vertical high and low frequency subspaces. We adopt the following four filters for decomposition: and for vertical high and low frequency subspaces, and and for horizontal high and low frequency subspace. The frequency responses of these filters satisfy the relationships
(16)  
(17) 
where 1 is the allones matrix. It’s easy to verify that
(18)  
(19) 
Therefore, the proposed filtering scheme is redundant and meets the requirements for completeness, and is thus lossless. We call this proposed decomposition scheme the HoriVert subspace decomposition. We display the frequency responses of these filter banks in Figure 4(a)–(d).
3.4.2 Gaussian subspace decomposition
We also test our DAC method using spatial Gaussian filters. We design the Gaussian lowpass filter, denoted as , with spatial support and unit standard deviation and similar for the corresponding highpass filter, denoted . This gives frequency responses for and for . and have a similar lossless property,
(20) 
We also show the frequency response of the Gaussian complementary filter banks in Figures 4(e)–(f).
To illustrate the proposed decomposition scheme, we show the brain MRI kspace magnitudes of the two subspaces using Gaussian subspace decomposition scheme in Figure 5. The sum of Figure 5(a) and 5(b) is equal to Figure 1(b). However, as Figure 5(a) and 5(b) indicate, these subspaces do not simply correspond to disjoint partitions of kspace. Also in Figure 5(a), we observe the magnitudes within the high frequency subspace keep in the same range, which will benefit the kspace relative reconstruction accuracy.
3.5 An equivalent objective
We next briefly summarize the basic objective function that our DAC algorithm is optimizing. Recall that the typical objective function for CSMRI has the form
(21) 
In the algorithm described above, we modify this to
(22) 
As can be seen, we minimize in two parts. First, we minimize over the two rightmost terms in the “divide” portion of the algorithm. We then minimize over in the first term using the learned . This way, the low and high frequency subspaces can contribute more equally to the reconstruction of .
Finally, the setting of is important when reconstructing the MRI . We found that a uniform setting consistently works well. As another approach, viewing the first term as an augmented Lagrangian that attempts to enforce what is originally a strict equality , we also experiment with maximizing over in an adversarial manner to try and enforce these equalities (subject to ). We use this approach in our experiments. We summarize the entire procedure in Algorithm 1.
4 Experiments
4.1 Experimental setup
For our experiments we adopt three sampling masks 12 (): 1D Cartesian sampling with random phase encodes, 2D random sampling, and pseudo radial sampling. These are shown in Figure 6. The undersampling ratio here means the ratio of kspace data being measured to the total number of fullsampled kspace data, it is negative correlated with reduction factor appears in other CSMRI literatures. We conduct the simulations on the phantom shown in Figure 6(d) and the MRIacquired complexvalued brain images shown in Figures 6(e), 6(f) and 6(g). These images are normalized to have maximum magnitudes of 1. As with other CSMRI methods, the compressed data is acquired by simulating the undersampling of the 2D DFT using the fullysampled MRI.
We compare the Gaussian and HoriVert subspace filtering methods using the three stateoftheart CSMRI methods described in Section 3.2. As performance measures we use the peak signaltonoise ratio (PSNR), the structural similarity index (SSIM), and the high frequency error norm (HFEN) 11 (). The standard PSNR is a function of the MSE, but as we previously indicated, the PSNR measure is not the optimal choice in assessing the quality of an MR image. Therefore we also use SSIM and HFEN. SSIM measures the structural similarity of two images and is more consistent with the evaluation system of the human eye. HFEN has been proposed to evaluate the reconstruction quality of high frequency portions of MRI. In HFEN, the Laplacian of Gaussian (LoG) filter is used to extract the high frequency information within the MRI. HFEN is measured by the norm of the extracted features between the fullysampled image and the reconstructed image.
All the experiments are coded in Matlab (R2014a). Computations are implemented with a Intel Core i5 CPU at 3.20GHz and 8G memory, employing a 64bit Windows 7 operating system. For FCSA, WatMRI and PBDWS, we use the source code available from the authors’ homepage, but we make parameter adjustments obtain the best performances.
4.2 Illustrative experiment on phantom data
The realvalued phantom image of the size shown in Figure 6(d) is piecewise constant and contains various image structures 42 (). The simluation phantom is created via Matlab implementation. Note that there exists rich low contrast, high frequency information in the phantom data. Thus this phantom data is more appropriate to evaluate the performance of the proposed DAC framework compared with conventional SheppLogan phantom. The SheppLogan phantom is extremely sparse under a gradient transform, so most CS algorithms can exactly reconstruct it from very few Fourier samples. To show the advantage of our divideandconquer method, we compare the reconstruction result of the original PBDWS algorithm with HoriVert PBDWS (DAC using HoriVert filters and PBDWS reconstruction) with a undersampled Cartesian mask, undersampled random mask and undersampled radial mask. In Figure 7 we show the error residual images for each reconstruction. As is clear, the proposed DAC method is able to reconstruct the high frequency data more accurately by allowing the PBDWS algorithm to focus on these regions independently from the low frequency information. We again note that the same reconstruction algorithm is being used in both cases; the only difference is whether the subsampled kspace data is modeled directly or indirectly through different low and high pass filters.
4.3 Experiments on T2wBrain data
Mask  FCSA  Gauss FCSA  HoriVert FCSA  PBDWS  Gauss PBDWS  HoriVert PBDWS  

Cartesian  25  28.700.8501.476  28.450.8501.559  28.930.8691.452  34.120.9240.824  34.520.9440.833  34.440.9500.741 
30  31.710.8931.042  32.540.9150.964  32.480.9180.968  36.790.9400.558  37.810.9650.544  37.190.9680.483  
35  32.160.9030.944  33.080.9260.888  33.060.9290.884  37.580.9450.511  38.670.9710.487  38.180.9730.435  
40  33.330.9140.758  35.030.9420.671  34.800.9410.668  39.050.9510.420  40.430.9760.395  39.710.9780.355  
Random  15  31.040.8960.676  32.200.9220.608  32.390.9270.596  34.720.9330.462  35.940.9600.390  35.370.9610.365 
20  32.550.9160.562  34.070.9440.429  34.240.9450.435  36.760.9500.340  38.280.9710.281  37.580.9720.261  
25  33.220.9260.530  35.230.9520.360  35.190.9500.385  38.030.9510.291  40.010.9780.226  39.060.9790.215  
30  34.240.9380.485  36.860.9630.279  36.580.9590.320  39.770.9570.230  42.400.9840.162  41.150.9840.155  
Radial  15  30.210.8791.057  30.190.8781.100  30.790.9011.002  33.890.9250.667  34.220.9470.653  34.140.9530.577 
20  31.900.9060.775  32.290.9080.742  32.920.9250.678  35.910.9400.492  36.880.9650.440  36.470.9680.399  
25  33.200.9260.617  34.240.9300.530  34.620.9400.501  37.850.9500.360  39.370.9760.305  38.690.9770.276  
30  34.090.9350.546  35.750.9430.413  35.770.9470.424  39.240.9550.290  41.420.9810.233  40.360.9820.212 
We also test our DAC framework on a clinicallyobtained brain MRI also experimented with in 6 (), 12 (), 15 (), 30 () and 31 (). In particular, we use the and slices (named in the acquisition order) of a complexvalued T2weighted brain MRI (size ) volume data called T2wBrain (“slice7” and “slice27” respectively, as shown in Figure 6), which is 2D acquired with 32 coils from a healthy volunteer by a 3T Siemens Trio Tim MRI scanner using the T2weighted turbo spin echo sequence (TR/TE=6100/99 , field of view, slice thickness). We do SENSE reconstruction as the parallel imaging technique with reduction factor 1 to compose full kspace of gold standard images. The full kspace data will be used for emulate singlechannel MRI.
We first test the T2wBrain slice27 data using Gaussian and HoriVert versions of FCSA and PBDWS.^{2}^{2}2As mentioned, we set the parameters to their best experimentallyobtained values according to SSIM index. For FCSA, these were , and for the the high frequency subspace reconstruction and , and for the low frequency subspace reconstruction. For PBDWS we set the data fidelity parameter for each subspace. We apply the same parameter setting for all the tested data. We detail the parameter setting in later discussion section. We show the reconstruction results for FCSA in Figure 8. The MR image details in the red box are magnified for better comparison. As is evident, for these sampling rates and patterns, the reconstructed images of FCSA have severe jagged visual effects because of poor high frequency reconstruction. These details are clearer using the proposed subspace method.
We also test with PBDWS and its DAC Gaussian and HoriVert extensions in Figure 9. Although PBDWS outperforms FCSA, high frequency degradation is similarly observable. With Gaussian and HoriVert PBDWS these fine structures are recovered better. We show the corresponding error residual images for both the FCSA and PBDWS based DAC framework in Figure 10 and Figure 11 respectively. The left column is the results directly using the CSMRI method, the middle and right column correspond to these results with the Gaussian and HoriVert based DAC framework. The first row, second row and the third row corresponds to the Cartesian, random and radial undersampling masks. The proposed subspace method shows smaller reconstruction errors in the structural details compared with the direct application of the same algorithms.
For further illustration, we plot the KRRE maps of the reconstructions in Figure 12 and Figure 13. As can be seen, for our DAC method the high frequency regions of reconstructed kspace suffers less from errors than the direct method. This helps confirm our claim in Section 2 that isolating frequency content into subspaces for independent reconstruction allows for a more uniform reconstruction of kspace than the common squared error penalty.
The quantitative results for the slice of the T2wBrain data are given in Table 1. As is clear, CSMRI methods like FCSA and PBDWS can be significantly improved using the proposed divideandconquer method, which is consistent with the previous subjective evaluation. One interesting phenomenon is that, while the KRRE of Gaussian PBDWS is worse in high frequency regions than HoriVert PBDWS, the PSNR of Gaussian PBDWS is better than HoriVert PBDWS, while the SSIM and HFEN evaluation gives the opposite conclusion. This helps confirms the claim in Section 2 that the PSNR index does not provide a completely convincing measure of reconstruction quality in terms of detail recovery. SSIM and HFEN are designed to measure this, and these quantitative measures are more in agreement with the shown KRRE maps and subjective evaluation. From Figure 11, we observe here that the HoriVert subspace DAC framework slightly outperforms the Gaussian counterparts in visual performance.
Finally, we also conduct experiments on the T2wBrain slice7 using Gaussian WatMRI and Gaussian PBDWS, as shown in the Figure 14. We also give the error residual images in Figure 15. Again, Gaussian WatMRI and Gaussian PBDWS achieve better performance than their standard counterparts, WatMRI and PBDWS.
4.4 Experiments on T1wBrain data
In addition to the complexvalued invivo T2weigthed brain MRI data, we also test on a complexvalued invivo T1weighted brain MRI image called T1wBrain to validate the proposed framework on different MRI modalities 14 (). The T1wBrain image is an axial brain image from a 3T commercial scanner (GE Healthcare, Waukesha, WI) with an eightchannel head coil (In Vivo, Gainesville, FL) using a twodimensional T1weighted spin echo protocol (TE/TR = , FOV, slices, matrix). We test various CSMRI algorithms on the T1wBrain data for comparison, including L1ESPIRiT 43 (), pFISTA 44 (), PANO 15 (), PBDWS 31 (), GBRWT 37 () and the proposed DAC Gaussian PBDWS. Note that L1ESPIRiT uses the parallel imaging technique, while we use a strategy similar to the T2wBrain data to emulate the single coil imaging for other algorithms. We have adjusted the parameters of these algorithms to their best performance in PSNR.
We give the reconstruction results and corresponding residual error images in Figures 16 and 17. We see that structural information is preserved better under the DAC Gaussian PBDWS compared with other CSMRI methods. To quantitatively assess these CSMRI methods, we show the PSNR, SSIM and HFEN indexes in Figure 18. We observe that Gaussian PBDWS achieve the highest PSNR and SSIM value meanwhile the lowest HFEN value.
For the Gaussian and HoriVert subspace methods, the computational time required is roughly two and four times greater than the corresponding regular methods because there are two and four corresponding optimizations, respectively, rather than one. This constituted the most computationally intensive part of the proposed DAC framework, but we observe it is easily parallelizeable. For the subspace decomposition and subspace integration steps, the matrix Hadamard multiplication is computationally efficient.
5 Discussion
5.1 Discussion on parameter setting
In the proposed divideandconquer framework, parameters requiring tuning are in both the subspace reconstruction and subspace integration stages. For subspace reconstruction, the number of parameters to be tuned only depends on the specific base algorithm adopted. If we divide kspace into subspaces, and the number of parameter for single subspace reconstruction is , the total number of parameters for subspace reconstruction is . Therefore, if the chosen subspace reconstruction algorithm is robust to variations in parameter setting, the DAC extension of that algorithm will also be robust. For example, in the Gaussian PBDWS method the parameter to be adjusted is the data fidelity parameter . For high and low frequency subspace reconstruction, we adjust the data fidelity regularization parameter low and high ranging from to , and we plot the performance curve in PSNR and SSIM in Figure 19. The experiment is conducted using PBDWS on T2wBrain slice with 2D mask.
We note the PSNR and SSIM index for both high and low frequency subspace reconstruction reach the optimal around 1e4 and above, meaning for Gaussian PBDWS, when the regularization parameters exceeds 1e4, the method is not susceptible to them. We can choose an arbitrary regularization value greater than 1e4 for any data used by our DAC framework with PBDWS even if we have no access to the fullysampled kspace data in real application scenarios. Hence we choose the regularization parameter in PBDWS used for subspace reconstruction to be 1e6, which is also recommended in the original paper in PBDWS.
For high and low frequency subspace reconstruction using FCSA/WatMRI, the regularization parameters for both the TV and wavelet terms can influence the subspace reconstruction, but we note that even when we use no parameter tuning in the subspace reconstruction phase, meaning the parameter setting is kept the same as the regular FCSA and WatMRI, the proposed method still outperforms these original methods by a considerable margin in Table 2. This shows the improvement of the proposed DAC framework is not simply a result of parameter tuning.
PBDW  Gaussian FCSA  HoriVert FCSA  

PSNR (dB)  33.78  34.74  34.20 
SSIM  0.914  0.922  0.915 
As for the parameter setting in subspace integration, we find that model performance is already good by setting all the subspace integration parameters the same. These parameter can also be estimated via the proposed scheme in Algorithm 1 using the augmented Lagrangian method. In this way the strict equality of the subspace decomposition holds.
5.2 Discussion on noisy environments
During the acquisition of MRI measurements, the contamination brought by noise is inevitable. Usually the SNR of the magnitude of a fullysampled MRI image is the ratio between the mean of the magnitudes and the noise standard deviation estimated from the background. Taking the slice of T2wBrain data for example, the SNR of the fullysampled is . The SNR index for the high frequency subspace decreases because the noise is amplified. As discussed in Section II, the magnitudes of the high frequency subspace MRI are small yet important, because it contains structural information and fine details. With an efficient CSMRI algorithm, we can denoise the high frequency subspace MR images while retaining image structures because CSMRI algorithms can benefit from high sparsity.
We have experimented with simulated noisy environments to evaluate the performance of the proposed DAC framework by adding the Gaussian noise into the undersampled kspace. We conduct experiments on the slice of the T2wBrain image, where we use 2D random mask for undersampling. We add Gaussian random noise to the kspace with various standard deviations from to to evaluate its robustness to noise using both FCSA and PBDWS methods. We plot the performance curve for PSNR and SSIM with respect to different noise level in Figure 20. The experiments show that the proposed DAC framework is robust to noise contamination in kspace. From Figures 20(b) and 20(d), we observe that the proposed DAC framework also outperforms in SSIM, meaning the high frequency information still better reconstructed in the presence of noise, despite its larger relative magnitude in this region. We also observe the margin to which the DAC frameworks outperform the regular counterparts increases as the noise level goes up.
6 Conclusion
Based on the common observation that the energy and sparsity of kspace is nonuniformly distributed, we propose a divideandconquer framework for CSMRI inversion. We first apply a series of linear filters to decompose the subsampled kspace measurements into separate frequency views called subspaces. For this we use two filtering schemes called HoriVert decomposition and Gaussian decomposition based on the linearvertical and Gaussian filters. We then reconstruct the corresponding MRI in each subspace independently using any offtheshelf CSMRI inversion algorithm. We obtain the final reconstructed MRI by integrating all the reconstructed subspace images using Tikhonov regularization.
The experimental results on simulated phantom data and acquired complexvalued T2wBrain and T1wBrain MRI data show that the proposed subspace method can improve the performance of existing stateoftheart CSMRI methods considerably. We also observe that the proposed method has potential for recovering finer highfrequency details for diagnosis, which may improve the reliability and effectiveness of CSMRI.
Acknowledgment
This work was supported in part by the National Natural Science Foundation of China under Grants 61571382, 81671766, 61571005, 81671674, U1605252, 61671309 in part by the Guangdong Natural Science Foundation under Grant 2015A030313007, in part by the Fundamental Research Funds for the Central Universities under Grant 20720160075, 20720180059, in part by the National Natural Science Foundation of Fujian Province, China under Grant 2017J01126.
References
References
 (1) E. J. Candès, J. Romberg, T. Tao, Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information, IEEE Transactions on Information Theory 52 (2) (2006) 489–509.
 (2) D. L. Donoho, Compressed sensing, IEEE Transactions on Information Theory 52 (4) (2006) 1289–1306.
 (3) M. Lustig, D. Donoho, J. M. Pauly, Sparse MRI: The application of compressed sensing for rapid MR imaging, Magnetic Resonance in Medicine 58 (6) (2007) 1182–1195.
 (4) V. M. Patel, R. Maleh, A. C. Gilbert, R. Chellappa, Gradientbased image recovery methods from incomplete Fourier measurements, IEEE Transactions on Image Processing 21 (1) (2012) 94–105.
 (5) Q. Liu, S. Wang, L. Ying, X. Peng, Y. Zhu, D. Liang, Adaptive dictionary learning in sparse gradient domain for image recovery, IEEE Transactions on Image Processing 22 (12) (2013) 4652–4663.
 (6) X. Qu, W. Zhang, D. Guo, C. Cai, S. Cai, Z. Chen, Iterative thresholding compressed sensing MRI based on contourlet transform, Inverse Problems in Science and Engineering 18 (6) (2010) 737–758.
 (7) Z. Feng, F. Liu, M. Jiang, S. Crozier, H. Guo, Y. Wang, Improved l1SPIRiT using 3D walsh transformbased sparsity basis, Magnetic Resonance Imaging 32 (7) (2014) 924–933.
 (8) S. G. Lingala, Y. Hu, E. DiBella, M. Jacob, Accelerated dynamic MRI exploiting sparsity and lowrank structure: kt SLR, IEEE Transactions on Medical Imaging 30 (5) (2011) 1042–1054.
 (9) H. Yoon, K. S. Kim, D. Kim, Y. Bresler, J. C. Ye, Motion adaptive patchbased lowrank approach for compressed sensing cardiac cine MRI, IEEE Transactions on Medical Imaging 33 (11) (2014) 2069–2085.
 (10) R. Otazo, E. Candès, D. K. Sodickson, Lowrank plus sparse matrix decomposition for accelerated dynamic MRI with separation of background and dynamic components, Magnetic Resonance in Medicine 73 (3) (2015) 1125–1136.
 (11) X. Qu, D. Guo, B. Ning, Y. Hou, Y. Lin, S. Cai, Z. Chen, Undersampled MRI reconstruction with patchbased directional wavelets, Magnetic Resonance Imaging 30 (7) (2012) 964–977.
 (12) B. Ning, X. Qu, D. Guo, C. Hu, Z. Chen, Magnetic resonance image reconstruction using trained geometric directions in 2D redundant wavelets domain and nonconvex optimization, Magnetic Resonance Imaging 31 (9) (2013) 1611–1622.
 (13) Z. Lai, X. Qu, Y. Liu, D. Guo, J. Ye, Z. Zhan, Z. Chen, Image reconstruction of compressed sensing MRI using graphbased redundant wavelet transform, Medical Image Analysis 27 (2016) 93–104.
 (14) S. Ravishankar, Y. Bresler, MR image reconstruction from highly undersampled kspace data by dictionary learning, IEEE Transactions on Medical Imaging 30 (5) (2011) 1028–1041.
 (15) Y. Huang, J. Paisley, Q. Lin, X. Ding, X. Fu, X.P. Zhang, Bayesian nonparametric dictionary learning for compressed sensing MRI, IEEE Transactions on Image Processing 23 (12) (2014) 5007–5019.
 (16) S. Wang, Z. Su, L. Ying, X. Peng, S. Zhu, F. Liang, D. Feng, D. Liang, Accelerating magnetic resonance imaging via deep learning, in: 2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI), IEEE, 2016, pp. 514–517.
 (17) J. Schlemper, J. Caballero, J. V. Hajnal, A. Price, D. Rueckert, A deep cascade of convolutional neural networks for MR image reconstruction, in: International Conference on Information Processing in Medical Imaging, Springer, 2017, pp. 647–658.
 (18) J. Sun, H. Li, Z. Xu, et al., Deep ADMMnet for compressive sensing MRI, in: Advances in Neural Information Processing Systems, 2016, pp. 10–18.
 (19) C. Chen, Y. Li, J. Huang, Forest sparsity for multichannel compressive sensing, IEEE Transactions on Signal Processing 62 (11) (2014) 2803–2813.
 (20) C. Chen, J. Huang, The benefit of tree sparsity in accelerated MRI, Medical Image Analysis 18 (6) (2014) 834–842.
 (21) G. Ongie, M. Jacob, A fast algorithm for structured lowrank matrix recovery with applications to undersampled MRI reconstruction, in: 2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI), IEEE, 2016, pp. 522–525.
 (22) X. Qu, Y. Hou, F. Lam, D. Guo, J. Zhong, Z. Chen, Magnetic resonance image reconstruction from undersampled measurements using a patchbased nonlocal operator, Medical Image Analysis 18 (6) (2014) 843–856.
 (23) W. Dong, G. Shi, X. Li, Y. Ma, F. Huang, Compressive sensing via nonlocal lowrank regularization, IEEE Transactions on Image Processing 23 (8) (2014) 3618–3632.
 (24) J. Trzasko, A. Manduca, Highly undersampled magnetic resonance image reconstruction via homotopic minimization, IEEE Transactions on Medical Imaging 28 (1) (2009) 106–121.
 (25) A. Majumdar, R. K. Ward, An algorithm for sparse MRI reconstruction by Schatten pnorm minimization, Magnetic Resonance Imaging 29 (3) (2011) 408–417.
 (26) J. C. Ye, S. Tak, Y. Han, H. W. Park, Projection reconstruction MR imaging using FOCUSS, Magnetic Resonance in Medicine 57 (4) (2007) 764–775.
 (27) H. Jung, K. Sung, K. S. Nayak, E. Y. Kim, J. C. Ye, kt FOCUSS: A general compressed sensing framework for high resolution dynamic MRI, Magnetic Resonance in Medicine 61 (1) (2009) 103–116.
 (28) S. Ma, W. Yin, Y. Zhang, A. Chakraborty, An efficient algorithm for compressed MR imaging using total variation and wavelets, in: 2008 IEEE Conference on Computer Vision and Pattern Recognition, 2008, pp. 1–8. doi:10.1109/CVPR.2008.4587391.
 (29) J. Yang, Y. Zhang, W. Yin, A fast alternating direction method for TVL1L2 signal reconstruction from partial Fourier data, IEEE Journal of Selected Topics in Signal Processing 4 (2) (2010) 288–297.
 (30) A. Beck, M. Teboulle, A fast iterative shrinkagethresholding algorithm for linear inverse problems, SIAM Journal on Imaging Sciences 2 (1) (2009) 183–202.
 (31) J. Huang, S. Zhang, D. Metaxas, Efficient MR image reconstruction for compressed MR imaging, Medical Image Analysis 15 (5) (2011) 670–679.
 (32) Y. Liu, Z. Zhan, J.F. Cai, D. Guo, Z. Chen, X. Qu, Projected iterative softthresholding algorithm for tight frames in compressed sensing magnetic resonance imaging, IEEE Transactions on Medical Imaging 35 (9) (2016) 2130–2140.
 (33) X. Ye, Y. Chen, F. Huang, Computational acceleration for MR image reconstruction in partially parallel imaging, IEEE Transactions on Medical Imaging 30 (5) (2011) 1055–1063.
 (34) Y. Yang, F. Liu, W. Xu, S. Crozier, Compressed sensing MRI via twostage reconstruction, IEEE Transactions on Biomedical Engineering 62 (1) (2015) 110–118.
 (35) X. Peng, D. Liang, MR image reconstruction with convolutional characteristic constraint (CoCCo), IEEE Signal Processing Letters 22 (8) (2015) 1184–1188.
 (36) K. Sung, B. A. Hargreaves, Highfrequency subband compressed sensing MRI using quadruplet sampling, Magnetic Resonance in Medicine 70 (5) (2013) 1306–1318.
 (37) S. Park, J. Park, Compressed sensing MRI exploiting complementary dual decomposition, Medical Image Analysis 18 (3) (2014) 472–486.
 (38) D. Smith, E. Welch, Nonsparse phantom for compressed sensing mri reconstruction, in: International Society for Magnetic Resonance in Medicine 19th Scientific MeetingISMRM, Vol. 11, 2011, p. 2845.

(39)
M. Uecker, P. Lai, M. J. Murphy, P. Virtue, M. Elad, J. M. Pauly, S. S.
Vasanawala, M. Lustig,
ESPIRiT—an eigenvalue approach
to autocalibrating parallel MRI: Where SENSE meets GRAPPA, Magnetic
Resonance in Medicine 71 (3) (2014) 990–1001.
doi:10.1002/mrm.24751.
URL http://dx.doi.org/10.1002/mrm.24751