A Divide-and-Conquer Approach to Compressed Sensing MRI
Liyan Sun Zhiwen Fan Xinghao Ding Congbo Cai Yue Huang John Paisley








Fujian   Key   Laboratory   of   Sensing   and   Computing   for   Smart   City,  Xiamen   University, Fujian, China

^†Department of Electrical Engineering, Columbia University, New York, NY, USA

Declarations of interest: none



Compressed sensing (CS) theory assures us that we can accurately reconstruct magnetic resonance images using fewer k-space measurements than the Nyquist sampling rate requires. In traditional CS-MRI inversion methods, the fact that the energy within the Fourier measurement domain is distributed non-uniformly is often neglected during reconstruction. As a result, more densely sampled low-frequency information tends to dominate penalization schemes for reconstructing MRI at the expense of high-frequency details. In this paper, we propose a new framework for CS-MRI inversion in which we decompose the observed k-space data into “subspaces” via sets of filters in a lossless way, and reconstruct the images in these various spaces individually using off-the-shelf 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 k-space more equally, preserving both high and low frequency details. We demonstrate that
the proposed framework is competitive with state-of-the-art methods in CS-MRI in terms of quantitative performance, and often improves an algorithm's results qualitatively compared with its direct application to k-space.



compressed sensing, magnetic resonance imaging, divide-and-conquer








§ 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 diagnostic-quality reconstruction is a major challenge in MR imaging.

Compressive sensing (CS) theory <cit.> 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 CS-MRI can be generalized as the optimization problem

    x =  min_x  F_u x - y _2^2 + f_θ(x).

The goal is to reconstruct the vector x ∈ℂ^N, which corresponds to values in the MR image X ∈ℂ^√(N)×√(N).
The matrix F_u∈ℂ^u × N,u < N is an under-sampled Fourier basis used to directly measure the k-space data y∈ℂ^u. In this objective, the first term enforces that the Fourier coefficients of x agree with the observations y. Since many vectors x will satisfy this requirement, the penalty f_θ(x) taking parameters θ searches for an x with additional desired properties such as smoothness.




 §.§ Related Work

Many CS-MRI algorithms have been proposed in the framework of Equation (<ref>). 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 CS-MRI 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 <cit.>, total variation <cit.>, contourlets <cit.>, Walsh <cit.> and PCA <cit.>. Patch-based bases include directional wavelets (PBDWS) <cit.>, a graph-based redundant wavelet transform (GBRWT) <cit.>, and dictionary bases constructed in situ using dictionary learning algorithms like KSVD <cit.> and BPFA <cit.>. As researches in deep learning methods thrive, the popular deep learning architectures are brought into CS-MRI <cit.>. In recent years, other approaches have enforced structural sparsity <cit.>, nonlocal priors <cit.>, and approximations to the desired ℓ_0 penalty such as the convex ℓ_1,
FOCUSS, ℓ_p (for 0<p<1), IRLS and smoothing functions <cit.>.

One particular method that we highlight and will use is the patch-based directional wavelet (PBDWS) <cit.>. This CS-MRI 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 (<ref>), such as TVCMRI <cit.>, RecPF <cit.>, FISTA <cit.>, FCSA <cit.>, pFISTA <cit.> (an algorithmic variation on SparseMRI <cit.>), Bregman <cit.> and ADMM <cit.>. These methods typically work by representing the objective function in a way that is easier to optimize iteratively.



 §.§ Our contribution


In general, k-space data of MRI are not distributed uniformly across all frequency bands in energy and magnitude. As an example, we show the magnitudes of k-space data for an MRI of the brain in Figure <ref>. As is clearly evident, the energy is concentrated much more in the low frequency region of k-space than the high frequency region. This is a well-known 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 k-space to ensure the basic structure is measured. However, as a result of using a squared error penalty term of the form F_ux - y_2^2 as in the Equation (<ref>), CS-MRI inversion will result in these low frequency measurements dominating the reconstruction of x. 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 k-space, 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 “divide-and-conquer” (DAC) manner and can incorporate any existing CS-MRI inversion algorithm to reconstruct those parts.

The paper is organized as follows: In Section 2, we discuss the non-uniformity property of the k-space data and its relationship to the reconstruction accuracy for CS-MRI. To address this issue, we propose a three-step divide-and-conquer (DAC) framework in Section 3. Section 4 demonstrates the performance of our DAC framework on several MR data with different under-sampling masks and ratios. We provide further discussion on parameter setting and noise characteristics in Section 5.



 §.§ Notation


We use the following notation: the fully-sampled or reference MRI is denoted as X_f ∈ℂ^√(N)×√(N) with the vectorized form x_f ∈ℂ^N. Upper case means the representation in 2D while lower case means a vectorized form of the corresponding 2D representation. The subscript f means fully-sampled. Projecting these into the frequency domain, we represent the corresponding projections as Y_f ∈ℂ^√(N)×√(N) and y_f ∈ℂ^N. The measured k-space Y ∈ℂ^√(N)×√(N) has its unsampled positions padded with zeros, but the vectorized form y ∈ℂ^u will only have the sampled locations. The vector y can be obtained by multiplying x_f by an under-sampled Fourier basis matrix F_u. We further denote the reconstructed MRI image as X_r∈ℂ^√(N)×√(N) and its k-space representation as Y_r∈ℂ^√(N)×√(N). If we define a filter H with limited spatial support, which is usually much smaller than the convolved √(N)×√(N) image, the block Toeplitz matrix of the filter H is defined as H∈ℝ^√(N)×√(N). Furthermore, the frequency response of filter h is defined as H∈ℝ^√(N)×√(N).





§ MOTIVATION: NONUNIFORM K-SPACE DENSITY



K-space data is known to be distributed non-uniformly in energy, as illustrated with the brain MRI in Figure <ref>. While k-space 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 CS-MRI methods treat all errors equally, which tends to favor high-magnitude, low-frequency information in reconstructions at the expense of the details.

To further analyze the phenomenon, we define the following two simple measures: the k-space absolute reconstruction error (KARE) and the k-space relative reconstruction error (KRRE). These are the same size as the k-space, with the (i,j)^th element denoted as follows,

    ( )_ij = |( Y_r )_ij - (Y_f )_ij|
       ( )_ij = |( Y_r )_ij - (Y_f )_ij|/| ( Y_f )_ij|.
    
    Here, Y_f∈ℂ^√(N)×√(N) is the fully-sampled k-space data, and Y_r∈ℂ^√(N)×√(N) is the k-space of the MRI reconstructed from subsamples.
    
    
    
    For the same brain MRI, we plot in Figure <ref> the reconstructed image and the corresponding KARE and KRRE images using a CS-MRI inversion method called PBDWS <cit.> with 10% under-sampling and a 2D random mask. In Figures <ref> and <ref>, 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 fully-sampled k-space, one can manually redistribute the k-space error from PBDWS more evenly in such a way that the PSNR remains unchanged, as shown in Figures <ref> and <ref>. However, as shown in Figures <ref> and <ref>, 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 k-space reconstruction accuracy, even at the expense of low frequency information, can lead to better visual quality. However, conventional CS-MRI methods tend to leave the high frequency information deemphasized in their reconstruction objectives, coupled with heavier sub-sampling in low-frequency regions.
    While in previous work such information is implicitly modeled, however, to our knowledge this problem has not been explicitly addressed for CS-MRI inversion problems.
    
    For example, in <cit.> the low frequency k-space data is reconstructed first because of its dense distribution, then the reconstructed low frequency portions are padded back to the measurements for the second-stage 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 <cit.> 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 <cit.>, a convolutional constraints is proposed, but the work non-uniformity in k-space. In <cit.>, a method called HiSub CS-MRI formalized a link between the k-space and wavelet domain to apply separate under-sampling and reconstruction for high/low frequency k-space data. In the HiSub method,
     the high/low frequency regions in k-space are defined based on the separation of wavelet sub-bands; compressed sensing techniques are used for the high frequency region while parallel imaging is used for the low-frequency region. The HiSub method relies on the specific sampling pattern and is not exclusively based on CS-MRI ideas.
    In <cit.>, 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 non-uniformity in k-space motivates us to reconsider the CS-MRI problem using a divide-and-conquer (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 well-defined for linear algebra, here we mean a specific frequency view into which the k-space measurements are decomposed using standard filtering techniques. This method allows for the algorithms to deal with the high and low frequency k-space 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.[We wish to emphasize here that our proposed method does not simply partition the k-space data into disjoint frequency
    regions.]
    
    § DIVIDE-AND-CONQUER 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 <ref>. 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 (<ref>). At the end of the section we summarize the proposed DAC CS-MRI framework in Algorithm <ref>. Here we adopt the HoriVert subspace decomposition for illustration, which we will elaborate in later section.
    
    
    
     §.§ 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 k-space 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 {H_i}_i = 1^S be the impulse responses of the set of chosen filters, where S is selected in advance. In principle, a decomposition in the image domain can be formulated as follows,
    
        X_f_i = X_f * H_i,
    
    where X_f_i∈ℂ^√(N)×√(N) is the convolved image of the fully-sampled in the ith subspace and * denotes the 2D convolution.
    
    We equivalently work within the Fourier domain and obtain a series of frequency responses denoted as {H_i∈ℂ^√(N)×√(N)}_i = 1^S. As is well known, matrix convolution in the image domain can now be converted into an element-wise matrix multiplication operation in the Fourier domain. Thus Equation (<ref>) is equivalent to
    
        Y_f_i = Y_f ⊙H_i,
    
    where Y_f_i∈ℂ^√(N)×√(N) is the Fourier transformation of X_f_i and ⊙ denotes element-wise multiplication.
    
    Using the under-sampling, a zero-filled partial k-space view, Y ∈ℂ^√(N)×√(N), of ground truth X can be written as
    
        Y = M ⊙ Y_f,
    
    where M ∈ℝ^√(N)×√(N) is the under-sampling mask in which the sampled locations contain ones and unmeasured frequencies contain zeros. By simple algebra, the zero-filled, under-sampled k-space data in each subspace can now be derived as
    
        [ Y_i = M ⊙Y_f_i = M ⊙ Y_f ⊙H_i =  Y ⊙H_i. ]
    
    Equation (<ref>) indicates that the partial k-space data decomposed in each subspace can be derived via the element-wise multiplication between the frequency response of the filter and the original partially observed k-space data, or equivalently between the partially observed frequency response and the complete k-space 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 element-wise multiplication to matrix-vector multiplication. Thus the impulse responses {H_i}_i = 1^S of the filters is rewritten as a circulant matrix {H_i∈ℝ^ N  ×  N}_i = 1^S. (Note that H_i and H_i differ in that the first operates on two-dimensional k-space data while the second is on vectorized images.)
    An equivalent form of Equation (<ref>) and Equation (<ref>) can now be written as
    
        x_f_i = H_i x_f,
    
    where x_f_i∈ℂ^N is the vectorized form of X_f_i. Similarly, Equation (<ref>) can be written as
    
        y_i = F_ux_f_i,
    
    while the original measured k-space of Equation (<ref>) is
    
        y = F_ux_f.
    
    We next return to the filter banks considered in Section <ref>.
    
    
    
     §.§ Subspace Reconstruction
    
    
    After dividing k-space into partial frequency views y_i, 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 y_i 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 subspace-specific optimization problem with an appropriate regularization term to be determined. The reconstruction of each subspace can be formulated by minimizing the following objective function
    
        min_x_iμ_i /2F_ux_i - y_i_2^2 + ρ _i(x_i),
    
    
    where ρ_i( x_i) enforces desired properties of the reconstructed subspace image and μ_i is the regularization parameter for the data fidelity term. As with the filters, H_i, these penalty functions ρ_i can be chosen to be any CS-MRI 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 state-of-the-art CS-MRI reconstruction methods: FCSA <cit.>, WatMRI <cit.> and PBDWS <cit.>. To summarize, these penalties are the following:
    
        ρ_i(x_i) = {[    αTV(x_i) + βΦ x_i_1                   FCSA; αTV(x_i) + βΦ x_i_tree                 WatMRI;              B_w x_i_0                  PBDWS ].
    
    
    All methods use the squared error as a data fidelity term. In FCSA, the objective function is split into TV and wavelet L_1 regularization sub-problem in an iterative manner, where each sub-problem is solved efficiently using proximal gradient descend algorithm. Then the reconstructed MRI is obtained via the weighted average of the solutions of the two sub-problems. In this sense it is an algorithmic development of the classic SparseMRI framework <cit.>. 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 CS-MRI methods based on a global sparse regularization with non-adaptive transform bases. Thus they are suitable to stand for the CS-MRI methods with fixed transform basis. PBDWS on the other hand is a patch-based 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 CS-MRI 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 CS-MRI models. In our experiments, we apply the same reconstruction algorithm to all subspaces, but with different parameter settings. If the MRI is divided into n subspaces via H_1,…,H_n, then the chosen algorithm would be run independently n times, once for each subspace. Therefore, our framework increases computation time by a factor of n, but the independence of each optimization allows for a straightforward parallelization.
    
    
     §.§ 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,
    
        x = ∑_i = 1^N x_i.
    
    We also consider a Tikhonov regularization method by formulating the integration according to the following objective function
    
        x = min_x ∑_i = 1^n λ _ix_i - H_ix_2^2.
    
    
    This objective function admits a closed-form solution, but direct computation is infeasible because of the high dimensionality of H_i.  However, by transferring the problem to the Fourier domain the reconstructed k-space is calculated element-wise followed by an inverse Fourier transform,
    
        X = F^H (∑_i λ_i H_i^H ⊙X_i)/(∑_i λ_iH_i _2^2).
    
    Here, the division is element-wise, as is the magnitude of H_i in the denominator. We discuss a method for determining each λ_i below in Section <ref>.
    
    
    
    
     §.§ 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).
    
    
    
      §.§.§ 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 k-space into horizontal/vertical high and low frequency subspaces. We adopt the following four filters for decomposition: h_1 = [-0.5,0.5] and h_3 = [0.5,0.5] for vertical high and low frequency subspaces, and h_2 = [-0.5,0.5]^T and h_4 = [0.5,0.5]^T for horizontal high and low frequency subspace. The frequency responses of these filters satisfy the relationships
    
        H_1 + H_3 = 1
        H_2 + H_4 = 1
    
    where 1 is the all-ones matrix. It's easy to verify that
    
        Y_1 + Y_3 = Y ⊙H_1 + Y ⊙H_3 = Y( H_1 + H_3) = Y
        Y_2 + Y_4 = Y ⊙H_2 + Y ⊙H_4 = Y( H_2 + H_4) = Y
    
    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 <ref>(a)–(d).
    
    
    
      §.§.§ Gaussian subspace decomposition
    
    
    We also test our DAC method using spatial Gaussian filters. We design the Gaussian low-pass filter, denoted as g_lp, with 5×5 spatial support and unit standard deviation and similar for the corresponding high-pass filter, denoted g_hp. This gives frequency responses G_lp for g_lp and G_hp for g_hp. G_lp and G_hp have a similar lossless property,
    
        G_lp + G_hp = 1.
    
    We also show the frequency response of the Gaussian complementary filter banks in Figures <ref>(e)–(f).
    
    To illustrate the proposed decomposition scheme, we show the brain MRI k-space magnitudes of the two subspaces using Gaussian subspace decomposition scheme in Figure <ref>. The sum of Figure <ref>(a) and <ref>(b) is equal to Figure <ref>. However, as Figure <ref>(a) and <ref>(b) indicate, these subspaces do not simply correspond to disjoint partitions of k-space. Also in Figure <ref>(a), we observe the magnitudes within the high frequency subspace keep in the same range, which will benefit the k-space relative reconstruction accuracy.
    
    
    
     §.§ An equivalent objective
    
    We next briefly summarize the basic objective function that our DAC algorithm is optimizing. Recall that the typical objective function for CS-MRI has the form
    
        x = min_x y - F_u x_2^2 + ρ(x).
    
    In the algorithm described above, we modify this to
    
        x = min_x ∑_i=1^n λ_iH_i x-x_i^2_2 +  min_x_iy_i - F_u x_i^2_2 + ρ_i(x_i).
    
    As can be seen, we minimize in two parts. First, we minimize x_i over the two right-most terms in the “divide” portion of the algorithm. We then minimize over x in the first term using the learned x_i. This way, the low and high frequency subspaces can contribute more equally to the reconstruction of x.
    
    Finally, the setting of λ_i is important when reconstructing the MRI x. 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 x_i = H_i x, we also experiment with maximizing over λ_i in an adversarial manner to try and enforce these equalities (subject to λ_2 = 1). We use this approach in our experiments. We summarize the entire procedure in Algorithm <ref>.
    
    
    
    § EXPERIMENTS
    
    
     §.§ Experimental setup
    
    
    For our experiments we adopt three sampling masks <cit.>: 1D Cartesian sampling with random phase encodes, 2D random sampling, and pseudo radial sampling. These are shown in Figure  <ref>. The under-sampling ratio here means the ratio of k-space data being measured to the total number of full-sampled k-space data, it is negative correlated with reduction factor appears in other CS-MRI literatures. We conduct the simulations on the phantom shown in Figure <ref> and the MRI-acquired complex-valued brain images shown in Figures <ref>, <ref> and <ref>. These images are normalized to have maximum magnitudes of 1. As with other CS-MRI methods, the compressed data is acquired by simulating the under-sampling of the 2D DFT using the fully-sampled MRI.
    
    We compare the Gaussian and HoriVert subspace filtering methods using the three state-of-the-art CS-MRI methods described in Section <ref>. As performance measures we use the peak signal-to-noise ratio (PSNR), the structural similarity index (SSIM), and the high frequency error norm (HFEN) <cit.>. 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 l_2 norm of the extracted features between the fully-sampled 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 64-bit 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.
    
    
    
     §.§ Illustrative experiment on phantom data
    
    
    The real-valued phantom image of the size 256×256 shown in Figure <ref> is piece-wise constant and contains various image structures <cit.>. 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 Shepp-Logan phantom. The Shepp-Logan 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 divide-and-conquer method, we compare the reconstruction result of the original PBDWS algorithm with HoriVert PBDWS (DAC using HoriVert filters and PBDWS reconstruction) with a 40% under-sampled Cartesian mask, 30% under-sampled random mask and 30% under-sampled radial mask.
    In Figure <ref> 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 sub-sampled k-space data is modeled directly or indirectly through different low and high pass filters.
    
    
    
     §.§ Experiments on T2wBrain data
    
    
    We also test our DAC framework on a clinically-obtained brain MRI also experimented with in <cit.>, <cit.>, <cit.>, <cit.> and <cit.>. In particular, we use the 7^th and 27^th slices (named in the acquisition order) of a complex-valued T2-weighted brain MRI (size 256×256) volume data called T2wBrain (“slice7” and “slice27” respectively, as shown in Figure <ref>), which is 2D acquired with 32 coils from a healthy volunteer by a 3-T Siemens Trio Tim MRI scanner using the T2-weighted turbo spin echo sequence (TR/TE=6100/99 ms, 220×220mm^2 field of view, 3 mm slice thickness). We do SENSE reconstruction as the parallel imaging technique with reduction factor 1 to compose full k-space of gold standard images. The full k-space data will be used for emulate single-channel MRI.
    
    We first test the T2wBrain slice27 data using Gaussian and HoriVert versions of FCSA and PBDWS.[As mentioned, we set the parameters to their best experimentally-obtained values according to SSIM index. For FCSA, these were μ = 2, α = 0.003 and β = 0.001 for the the high frequency subspace reconstruction and μ = 2, α = 0.002 and β = 0.002 for the low frequency subspace reconstruction. For PBDWS we set the data fidelity parameter μ=1e6 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 <ref>. 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 <ref>. 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 <ref> and Figure <ref> respectively. The left column is the results directly using the CS-MRI 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 under-sampling 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 <ref> and Figure <ref>. As can be seen, for our DAC method the high frequency regions of reconstructed k-space suffers less from errors than the direct method. This helps confirm our claim in Section <ref> that isolating frequency content into subspaces for independent reconstruction allows for a more uniform reconstruction of k-space than the common squared error penalty.
    
    The quantitative results for the 27^th slice of the T2wBrain data are given in Table <ref>. As is clear, CS-MRI methods like FCSA and PBDWS can be significantly improved using the proposed divide-and-conquer 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 <ref> 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 <ref>, 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 <ref>. We also give the error residual images in Figure <ref>. Again, Gaussian WatMRI and Gaussian PBDWS achieve better performance than their standard counterparts, WatMRI and PBDWS.
    
    
    
     §.§ Experiments on T1wBrain data
    
    
    In addition to the complex-valued in-vivo T2-weigthed brain MRI data, we also test on a complex-valued in-vivo T1-weighted brain MRI image called T1wBrain to validate the proposed framework on different MRI modalities <cit.>. The T1wBrain image is an axial brain image from a 3T commercial scanner (GE Healthcare, Waukesha, WI) with an eight-channel head coil (In Vivo, Gainesville, FL) using a two-dimensional T1-weighted spin echo protocol (TE/TR = 11/700 ms, 22 cm FOV, 10 slices, 256 × 256 matrix). We test various CS-MRI algorithms on the T1wBrain data for comparison, including L1-ESPIRiT <cit.>, pFISTA <cit.>, PANO <cit.>, PBDWS <cit.>, GBRWT <cit.> and the proposed DAC Gaussian PBDWS. Note that L1-ESPIRiT 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 <ref> and <ref>. We see that structural information is preserved better under the DAC Gaussian PBDWS compared with other CS-MRI methods. To quantitatively assess these CS-MRI methods, we show the PSNR, SSIM and HFEN indexes in Figure <ref>. 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.
    
    
    
    
    § DISCUSSION
    
    
     §.§ Discussion on parameter setting
    
    
    In the proposed divide-and-conquer 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 k-space into N subspaces, and the number of parameter for single subspace reconstruction is P, the total number of parameters for subspace reconstruction is NP. 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 μ_i. For high and low frequency subspace reconstruction, we adjust the data fidelity regularization parameter low and high ranging from 1 to 1e10, and we plot the performance curve in PSNR and SSIM in Figure <ref>. The experiment is conducted using PBDWS on
    T2wBrain 27^th slice with 2D 15% 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 fully-sampled k-space 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 <ref>. This shows the improvement of the proposed DAC framework is not simply a result of parameter tuning.
    
    
    
    
    
    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.
    
    
    
     §.§ 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 fully-sampled MRI image is the ratio between the mean of the magnitudes and the noise standard deviation estimated from the background. Taking the 27^th slice of T2wBrain data for example, the SNR of the fully-sampled is 37.97. 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 CS-MRI algorithm, we can denoise the high frequency subspace MR images while retaining image structures because CS-MRI 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 under-sampled k-space. We conduct experiments on the 27^th slice of the T2wBrain image, where we use 25% 2D random mask for under-sampling. We add Gaussian random noise to the k-space with various standard deviations from 0.01 to 0.05 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 <ref>. The experiments show that the proposed DAC framework is robust to noise contamination in k-space. From Figures <ref> and <ref>, 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.
    
    
    
    § CONCLUSION
    
    
    Based on the common observation that the energy and sparsity of k-space is non-uniformly distributed, we propose a divide-and-conquer framework for CS-MRI inversion. We first apply a series of linear filters to decompose the subsampled k-space measurements into separate frequency views called subspaces. For this we use two filtering schemes called HoriVert decomposition and Gaussian decomposition based on the linear-vertical and Gaussian filters. We then reconstruct the corresponding MRI in each subspace independently using any off-the-shelf CS-MRI 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 complex-valued T2wBrain and T1wBrain MRI data show that the proposed subspace method can improve the performance of existing state-of-the-art CS-MRI methods considerably. We also observe that the proposed method has potential for recovering finer high-frequency details for diagnosis, which may improve the reliability and effectiveness of CS-MRI.
    
    
    
    § 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
