A Divide-and-Conquer Approach toCompressed Sensing MRI

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
Abstract

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.

keywords:
compressed sensing, magnetic resonance imaging, divide-and-conquer
journal: Medical Image Analysis

1 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 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 CS-MRI 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 under-sampled Fourier basis used to directly measure the k-space 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 CS-MRI 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 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 3 (), total variation 3 (); 4 (); 5 (), contourlets 6 (), Walsh 7 () and PCA 8 (); 9 (); 10 (). Patch-based bases include directional wavelets (PBDWS) 30 (); 31 (), a graph-based 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 CS-MRI 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 patch-based directional wavelet (PBDWS) 30 (); 31 (). 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 (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, 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 1. 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 as in the Equation (1), CS-MRI 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 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 II, 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 III. Section IV 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 V.

1.3 Notation

We use the following notation: the fully-sampled 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 fully-sampled. Projecting these into the frequency domain, we represent the corresponding projections as and . The measured k-space 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 under-sampled Fourier basis matrix . We further denote the reconstructed MRI image as and its k-space 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 .

(a) Brain MRI
(b) k-space coefficients
Figure 1: k-space corresponding to a brain MRI.

2 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 1. 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 element denoted as follows,

(2)
(3)

Here, is the fully-sampled k-space data, and is the k-space of the MRI reconstructed from subsamples.

(a) KARE map of PBDWS  
(b) KRRE map of PBDWS  
(c) Redistributed KARE  
(d) Redistributed KRRE  
(e) PBDWS (SSIM = 0.91)
(f) Redistributed (SSIM = 0.97)
Figure 2: CS-MRI reconstruction results using PBDWS and an ad-hoc error redistribution for 2D random sampling. Since the ground truth is not known in real experiments, this motivating example is purely illustrative. (a) The KARE map of the PBDWS reconstruction. (b) The KRRE map of the PBDWS reconstruction. (c) The KARE map of the redistributed k-space error reconstruction. (d) The KRRE map of the redistributed k-space error reconstruction. (e) CS-MRI reconstruction image using PBDWS. (f) CS-MRI reconstruction after error redistribution.

For the same brain MRI, we plot in Figure 2 the reconstructed image and the corresponding KARE and KRRE images using a CS-MRI inversion method called PBDWS 31 () with under-sampling 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 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 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 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 27 () 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 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 non-uniformity in k-space. In 29 (), 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 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 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.111We wish to emphasize here that our proposed method does not simply partition the k-space data into disjoint frequency regions.

3 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 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 CS-MRI framework in Algorithm 1. Here we adopt the HoriVert subspace decomposition for illustration, which we will elaborate in later section.

Figure 3: Flowchart of proposed subspace method.

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 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 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 fully-sampled 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 element-wise matrix multiplication operation in the Fourier domain. Thus Equation (4) is equivalent to

(5)

where is the Fourier transformation of and denotes element-wise multiplication.

Using the under-sampling, a zero-filled partial k-space view, , of ground truth can be written as

(6)

where 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

(7)

Equation (7) 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 of the filters is rewritten as a circulant matrix . (Note that and differ in that the first operates on two-dimensional k-space 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 k-space of Equation (6) is

(10)

We next return to the filter banks considered in Section 3.4.

3.2 Subspace Reconstruction

After dividing k-space 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 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

(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 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 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 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 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 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 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 closed-form solution, but direct computation is infeasible because of the high dimensionality of . However, by transferring the problem to the Fourier domain the reconstructed k-space is calculated element-wise followed by an inverse Fourier transform,

(15)

Here, the division is element-wise, 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 k-space 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 all-ones 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).

(a)
(b)
(c)
(d)
(e)
(f)
Figure 4: The frequency response of the HoriVert redundant filter banks and the Gaussian complementary filter banks. White corresponds to one and black corresponds to zero, with a smooth transition in between. (a)–(f) are respectively the vertical high frequency, vertical low frequency, horizontal high frequency, horizontal low frequency, Gaussian low frequency, Gaussian high frequency filters.

3.4.2 Gaussian subspace decomposition

We also test our DAC method using spatial Gaussian filters. We design the Gaussian low-pass filter, denoted as , with spatial support and unit standard deviation and similar for the corresponding high-pass 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 k-space 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 k-space. Also in Figure 5(a), we observe the magnitudes within the high frequency subspace keep in the same range, which will benefit the k-space relative reconstruction accuracy.

(a) high frequency subspace
(b) low frequency subspace
Figure 5: The high and low frequency subspaces of a brain MRI obtained with the Gaussian filtering decomposition corresponding to Figures 4(e)–(f). The sum of Figure 5(a) and 5(b) is equal to Figure 1(b). However, we note that these are not simply disjoint partitions of the space.

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 CS-MRI 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 right-most 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.

  Input sub-sampled k-space and Fourier basis .
  Select subspace filters .
  Select CS-MRI penalizations .
  for each  do
     Filter using the th filter as in Equation (7).
     Optimize .
  end for
  while not converged do
     Set using Eq. (15).
     Set and normalize to unit L2 length.
  end while
  return  Vectorized reconstructed MRI .
Algorithm 1 Divide-and-conquer CS-MRI

4 Experiments

(a) Cartesian mask
(b) random mask
(c) radial mask
(d) phantom data
(e) T2wBrain slice
(f) T2wBrain slice
(g) T1wBrain
Figure 6: The MRI data and under-sampling masks used in our experiments.
Figure 7: The reconstruction errors results on phantom data using PBDWS and HoriVert PBDWS under three different under-sampling regimes. The first and second row represent the regular PBDWS and HoriVert PBDWS. From left to the right column we show the results with Cartesian mask, random mask and radial mask.

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 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 6(d) and the MRI-acquired complex-valued 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 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 3.2. As performance measures we use the peak signal-to-noise 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 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.

4.2 Illustrative experiment on phantom data

The real-valued phantom image of the size shown in Figure 6(d) is piece-wise 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 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 under-sampled Cartesian mask, under-sampled random mask and under-sampled 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 sub-sampled k-space data is modeled directly or indirectly through different low and high pass filters.

(a) Fully-sampled
(b) FCSA
(c) Gau FCSA
(d) HV FCSA
(e) Fully-sampled
(f) FCSA
(g) Gau FCSA
(h) HV FCSA
(i) Fully-sampled

(j) FCSA
(k) Gau FCSA
(l) HV FCSA
Figure 8: The experiments conducted on T2wBrain slice27 data using FCSA method under Gaussian and HoriVert subspace schemes. The 1D under-sampling Cartesian mask,2D under-sampling random mask and under-sampling radial mask are applied in the first row, second row and third row experiments respectively. The reconstruction details are magnified in the red box.
(a) Fully-sampled
(b) PBDWS
(c) Gau PBDWS
(d) HV PBDWS
(e) Fully-sampled
(f) PBDWS
(g) Gau PBDWS
(h) HV PBDWS
(i) Fully-sampled
(j) PBDWS
(k) Gau PBDWS
(l) HV PBDWS
Figure 9: The experiments conducted on T2wBrain slice27 data using PBDWS method under Gaussian and HoriVert subspace schemes. The 1D under-sampling Cartesian mask,2D under-sampling random mask and under-sampling radial mask are applied in the first row, second row and third row experiments respectively.

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
Table 1: PSNRSSIMHFEN of the reconstruction for the T2wBrain 27slice MRI using various sampling masks and rates. Larger values are better for PSNR and SSIM, smaller values for HFEN. In most cases divide-and-conquer improves the base algorithm.
Figure 10: The error residual images in Figure 8. The first row corresponds to Figure 8(b), Figure 8(c), Figure 8(d). The second row corresponds to Figure 8(f), Figure 8(g), Figure 8(h). The third row corresponds to Figure 8(j), Figure 8(k), Figure 8(l).
Figure 11: The error residual images in Figure 9. The first row corresponds to Figure 9(b), Figure 9(c), Figure 9(d). The second row corresponds to Figure 9(f), Figure 9(g), Figure 9(h). The third row corresponds to Figure 9(j), Figure 9(k), Figure 9(l). The first row corresponds to the error range from to while the second and third from to .
Figure 12: The KRRE maps of the reconstructed images corresponds to the Figure 10. It is noted that the sampled positions in high frequency regions are reconstructed more accurately under the subspace framework.
Figure 13: The KRRE maps of the reconstructed images corresponds to the Figure 11.

We also test our DAC framework on a clinically-obtained 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 complex-valued T2-weighted 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 3-T Siemens Trio Tim MRI scanner using the T2-weighted 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 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.222As mentioned, we set the parameters to their best experimentally-obtained 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 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 12 and Figure 13. 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 2 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 slice of the T2wBrain data are given in Table 1. 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 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.

(a) Fully-sampled
(b) 30% radial mask
(c) 20% random mask
(d) WatMRI
(e) Gau WatMRI
(f) HV WatMRI
(g) PBDWS
(h) Gau PBDWS
(i) HV PBDWS
Figure 14: The experiments conducted on T2wBrain slice7 data using WatMRI with radial under-sampling and PBDWS with random under-sampling under Gaussian subspace schemes.

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.

Figure 15: The error residual images in Figure 14. The first row corresponds to Figure 14(d), Figure 14(e), Figure 14(f). The second row corresponds to Figure 14(g), Figure 14(h), Figure 14(i).

4.4 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 14 (). 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 = , FOV, slices, matrix). We test various CS-MRI algorithms on the T1wBrain data for comparison, including L1-ESPIRiT 43 (), pFISTA 44 (), PANO 15 (), PBDWS 31 (), GBRWT 37 () 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 16 and 17. 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 18. We observe that Gaussian PBDWS achieve the highest PSNR and SSIM value meanwhile the lowest HFEN value.

(a) Fully-sampled
(b) Mask
(c) L1-ESPIRiT
(d) pFISTA
(e) PANO
(f) PBDWS
(g) GRBWT
(h) Gau PBDWS
Figure 16: Experiments conducted on T1wBrain data using various CS-MRI methods with Cartesian under-sampling. We note that L1-ESPIRiT is a parallel imaging CS-MRI algorithm.
Figure 17: The error residual images in Figure 16. We note the proposed DAC framework applied to the PBDWS algorithm using Gaussian subspace decomposition achieves the minimum reconstruction error. The first row corresponds to Figure 16(c), Figure 16(d), Figure 14(e) and The second row corresponds to Figure 16(f), Figure 16(g), Figure 14(h).
Figure 18: The PSNR, SSIM and HFEN index for the experiments conducted on T1wBrain data with Cartesian under-sampling.

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.

(a) PSNR
(b) SSIM
Figure 19: The performance curve as a function of data fidelity regularization parameter . The reconstruction remains stable when the regularization parameters reaches a certain value for both high and low frequency subspace.

5 Discussion

5.1 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 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 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 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
Table 2: The experiments using same parameter setting.

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 fully-sampled 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 fully-sampled 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 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 slice of the T2wBrain image, where we use 2D random mask for under-sampling. We add Gaussian random noise to the k-space 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 k-space. 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.

(a) PSNR
(b) SSIM
(c) PSNR
(d) SSIM
Figure 20: Experiments conducted on T2wBrain slice27 data using FCSA, PBDWS and their DAC counterparts with random under-sampling. Zero-mean Gaussian noise is added with standard deviation ranging from to .

6 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

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, Gradient-based 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 l1-SPIRiT using 3D walsh transform-based 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 low-rank 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 patch-based low-rank 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, Low-rank 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 patch-based 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 non-convex 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 graph-based redundant wavelet transform, Medical Image Analysis 27 (2016) 93–104.
  • (14) S. Ravishankar, Y. Bresler, MR image reconstruction from highly undersampled k-space 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 ADMM-net 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 multi-channel 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 low-rank 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 patch-based 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 low-rank 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 p-norm 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, k-t 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 TVL1-L2 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 shrinkage-thresholding 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 soft-thresholding 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 two-stage 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, High-frequency 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, Non-sparse phantom for compressed sensing mri reconstruction, in: International Society for Magnetic Resonance in Medicine 19th Scientific Meeting-ISMRM, 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
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
133334
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description