Int J CARS DOI 10.1007/s11548-015-1172-7
ORIGINAL ARTICLE
Improving GRAPPA reconstruction by frequency discrimination in the ACS lines Santiago Aja-Fernández1 · Daniel García Martín1 · Antonio Tristán-Vega1 · Gonzalo Vegas-Sánchez-Ferrero1,2
Received: 21 November 2014 / Accepted: 9 March 2015 © CARS 2015
Abstract Purpose GRAPPA is a well-known parallel imaging method that recovers the MR magnitude image from aliasing by using a weighted interpolation of the data in k-space. To estimate the optimal reconstruction weights, GRAPPA uses a band along the center of the k-space where the signal is sampled at the Nyquist rate, the so-called autocalibrated (ACS) lines. However, while the subsampled lines usually belong to the medium- to high-frequency areas of the spectrum, the ACS lines include the low-frequency areas around the DC component. The use for estimation and reconstruction of areas of the k-space with very different features may negatively affect the final reconstruction quality. We propose a simple, yet powerful method to eliminate reconstruction artifacts, based on the discrimination of the low-frequency spectrum. Methods The proposal to improve the estimation of the weights lays on a proper selection of the coefficients within the ACS lines, which advises discarding those points around the DC component. A simple approach is the elimination of a square window in the center of the k-space, although more developed approaches can be used. Results The method is tested using real multiple-coil MRI acquisitions. We empirically show this approach achieves great enhancement rates, while keeping the same complexity of the original GRAPPA and reducing the g-factor. The reconstruction is even more accurate when combined with other reconstruction methods. Improvement rates of 35 %
B
Santiago Aja-Fernández
[email protected] 1
Laboratorio de Procesado de Imagen (LPI), ETSI Telecomunicación, Universidad de Valladolid, Valladolid, Spain
2
ACIL, Brigham and Women’s Hospital, Harvard Medical School, Boston, MA, USA
are achieved for 32 ACS and acceleration rate of 3. Conclusions The method proposed highly improves the accuracy of the GRAPPA coefficients and therefore the final image reconstruction. The method is fully compatible with the original GRAPPA formulation and with other optimization methods proposed in literature, and it can be easily implemented into the commercial scanning software.
Keywords MRI · GRAPPA · Parallel imaging · Reconstruction · Estimation
Introduction The generalized autocalibrated partially parallel acquisition (GRAPPA) [8] is one of the most prominent parallel imaging reconstruction techniques for subsampled multiple-coil MR data. The reconstruction takes place into the k-space, where the missing lines are interpolated by a weighted combination of the adjacent acquired lines. The weights used in the reconstruction are estimated from the acquired data themselves by solving an overdetermined system of linear equations posed over the so-called autocalibrated (ACS) lines. The ACS lines are a set of lines from the center of the k-space which are sampled at the Nyquist rate (i.e., unaccelerated). As shown in Fig. 1, left, the ACS lines involve the center of the k-space (low frequencies of the image) as well as the high-frequency areas. Since the subsampled areas of the k-space are usually in the medium to high frequencies, the inclusion of the DC and low-frequency components in the estimation may bias the reconstruction weights. In order to improve the reconstruction performance of GRAPPA, many authors have proposed modifications to the original algorithm. Most of the efforts have been oriented
123
Int J CARS
to the reconstruction step rather than the estimation of the interpolation weights, such as using self-calibrating iterative methods [14], banks of filters [4,17]), image support reduction via a high-pass filtering [10] or kernel-based approaches (nonlinear GRAPPA) [3]. However, some authors pointed out the importance of estimation and proper selection of the ACS lines. For instance, in [15] a different sampling pattern of the ACS is proposed, whereas in [11] a method to remove outliers before estimating the weights is suggested. In this paper, we propose a simple but effective method to estimate the GRAPPA weights that highly improves the final quality of the reconstructed image. It is based on the removal of the low-frequency areas of the ACS lines only during the estimation step, to avoid the bias they introduce in the reconstruction weights. The method is fully compatible with the original GRAPPA formulation, and it can be easily implemented into the commercial scanning software. In addition, the method is also compatible with other optimization methods proposed in the literature.
Theory Estimation of the GRAPPA reconstruction weights The original GRAPPA reconstruction strategy estimates the full k-space in each coil from a subsampled k-space [1,8, 9] acquisition. While the sampled data slS (k) remains the same, the reconstructed lines slR (k) are estimated through a linear combination of the existing samples. Weighted data in a neighborhood η(k) around the estimated pixel from several coils are used for such estimation: slR (k) =
L
smS (k − c)ωm (l, c),
(1)
m=1 c∈η(k)
with sl (k) the complex signal from coil l at point k and ωm (l, k) the complex reconstruction coefficients for coil l. These coefficients are determined from the center lines of the k-space, the so-called ACS lines, which are sampled at the Nyquist rate (i.e., unaccelerated) as shown in Fig. 1. k−space
slACS (k) =
L
smACS (k − c)ωm (l, c), l = 1, . . . , L
m=1 c∈η(k)
(2) and fitting the resultant equation with some optimization method. Frequency discrimination of the ACS components Since the estimated weights are the same for each point within the image for every coil (they do not depend on k), there is a major problem with the estimation: While the subsampled lines usually belong to the medium- to high-frequency areas of the spectrum, the ACS lines include the low-frequency areas around the DC component (see Fig. 1 for illustration). The use for estimation and reconstruction of areas of the kspace with very different features may negatively affect the final reconstruction quality. The influence of the different areas of the spectrum over the final image is a very well-known topic in the signal processing field [13]. However, in GRAPPA, the different properties of the spectrum are not taken into account when estimating the reconstruction weights. As a result, there is a bias due to the low-frequency components in the final coefficients that worsens the reconstruction. As an illustration, in Fig. 2, the average of the magnitude of the 32 ACS lines of an eight-coil acquisition is depicted. Note that, clearly, the signal around the DC component has very different properties in shape and magnitude when compared to the medium- and high-frequency areas. Accordingly, the center of the k-space will have most of the energy of the ACS components. In Fig. 3, the ratio between the energy of a centered square window of size N × N with the total ACS energy is depicted for data sets 1 and 2 (later explained in Sect. Materials and methods). Note that a 5 × 5 window represents the 80 % of the energy of the ACS lines, while a 31 × 31 window has more than 95 % of the energy.
High freq. Medium freq.
ACS lines
The weights ωm (l, k) are obtained using Eq. (1) over these ACS lines
Low freq.
Low frequency
} ACS lines
DC comp.
High frequency Fig. 1 Frequency distribution in the k-space of one coil
123
Int J CARS 60
If we assume a symmetrical distribution of the energies of the frequency components around the center of the k-space, the removed window will accomplish all the center of the ACS, leaving some lines in the top and the bottom available for reconstruction. The number of lines will depend on the acceleration rate, R. Thus, for a symmetrical distribution of the energies, we can define the optimal size of N as
50 40 30 20
Nopt = |ACS| − (R + 1)
10 0
50
100
150
200
250
Fig. 2 Average of the magnitude of the 32 ACS lines of a eight-coil acquisition (k-space)
(3)
where |ACS| is the number of ACS lines. However, in practical situations, the energy distribution is not totally symmetrical, and therefore, it is advisable to choose a smaller window:
1
N ≤ |ACS| − (R + 1).
Total Energy Ratio
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2
Data Set 1 Date Set 2
5
10
15
20
25
30
N (Window size)
Fig. 3 Ratio of the energy of a square window centered in the DC component as a function of the size. Data sets 1 and 2 are eight-coil MR acquisitions of the brain from a 1.5T and a 3T scanner, respectively. More details in “Materials and methods” Center of k−space is removed
} ACS lines Data used for estimation
The assumption in Eq. (3) will be tested later on the experiments section. The method proposed is totally compatible with the existing GRAPPA and can be easily included in the scanner software. Regarding the complexity of the proposed methodology, note that the reconstruction is performed in exactly the same way as it is done in the standard GRAPPA. The only appreciable reduction of complexity is obtained in the estimation step, where the size of the systems of equations to estimate the coefficients from the ACS data is reduced as the number of points removed. As a final recall, note that, although the center of the kspace is not used for coefficient estimation, it is acquired and indeed used for reconstruction. Such a simple approach will highly improve the reconstruction, as shown in the following experiments.
Materials and methods
Fig. 4 New estimation area ACSR
Our proposal to improve the estimation of the weights lays on a proper selection of the coefficients within the ACS lines, which advises discarding of those points around the DC component. Although more developed approaches can be made, in this work we simply propose the elimination of a square window in the center of the k-space, as shown in Fig. 4. Thus, the data used for estimation will be a subset of the ACS lines:
For numerical validation of the method, three nonaccelerated acquisitions are considered: – Data set 1 (DS1): an eight-channel head coil data acquired on a GE Signa 1.5T EXCITE scanner, FGRE Pulse Sequence, TR/TE = 500/13.8, matrix size 256 × 256, FOV = 20 × 20cm. – Data set 2 (DS2): MR data set from [12] (the authors provide their data online1 ); eight-channel head array from 3 T GE scanner, FGRE sequence, TR/TE = 300/10, matrix size 256 × 256, FOV = 22 × 22cm. – Data set 3 (DS3): a doped ball phantom, scanned in an eight-channel head coil on a GE Signa 1.5T EXCITE
ACSR = ACS − {Ck [N × N ]} where ACS stands for all ACS points and Ck [N × N ] is a N × N square centered in the DC component. Other shapes may be used; a square was elected for the sake of simplicity.
1
http://www.ece.tamu.edu/~jimji/pulsarweb/index.htm.
123
Int J CARS Table 1 Comparison of removed window size and shape: reconstruction performance for different size N of the Removed Window Ck [N × N ] for two data sets, acceleration factor R = 3, 32 ACS lines and 3 × 2 kernel size.
Size
Data set 1
Data set 2
Square window
Round window
Square Window
Round Window
MSE
SSIM
MSE
SSIM
MSE
SSIM
MSE
SSIM
0
3.9291
0.8949
3.9291
0.8949
5.3963
0.8432
5.3963
0.8432
1
3.9258
0.8950
3.9235
0.8952
5.3749
0.8442
5.3915
0.8435
3
3.6517
0.9077
3.8854
0.8970
4.8786
0.8678
5.3147
0.8470
5
3.3762
0.9200
3.6916
0.9059
4.4947
0.8853
4.9906
0.8622
7
3.1787
0.9283
3.3653
0.9203
4.3540
0.8911
4.5077
0.8843
9
3.0517
0.9334
3.1739
0.9288
4.2193
0.8964
4.3510
0.8912
11
2.9847
0.9362
3.0576
0.9333
4.1173
0.9004
4.2341
0.8960
13
2.8945
0.9401
3.0015
0.9356
4.0164
0.9048
4.1528
0.8992
15
2.8004
0.9443
2.8937
0.9403
3.9126
0.9095
4.0549
0.9035
17
2.7436
0.9466
2.8233
0.9433
3.8266
0.9135
3.9675
0.9073
19
2.7148
0.9478
2.7718
0.9455
3.7549
0.9168
3.8668
0.9119
21
2.6833
0.9491
2.7339
0.9471
3.6689
0.9205
3.8046
0.9146
23
2.6389
0.9511
2.7044
0.9483
3.6015
0.9236
3.7371
0.9177
25
2.6096
0.9525
2.6692
0.9499
3.5278
0.9269
3.6695
0.9208
27
2.5925
0.9538
2.6332
0.9515
3.4716
0.9295
3.5976
0.9241
29
2.5991
0.9549
2.5978
0.9532
3.4619
0.9312
3.5455
0.9271
31
2.8027
0.9532
2.5803
0.9545
3.8033
0.9232
3.5065
0.9295
Square and round windows are considered. Best results for each data set are highlighted
Table 2 Comparison of removed window size and reconstruction kernel size: Reconstruction performance for different size N of the (square) removed window Ck [N × N ] for two data sets, acceleration factor R = 3 and 32 ACS lines; 5 × 2 and 7 × 2 are considered
Size
Data set 1
Data set 2
5×2
7×2
7×2
MSE
SSIM
MSE
SSIM
MSE
SSIM
MSE
SSIM
0
3.3085
0.9213
3.2070
0.9244
4.6976
0.8717
4.5670
0.8741
1
3.3085
0.9213
3.2071
0.9244
4.6913
0.8719
4.5617
0.8743
3
3.1816
0.9270
3.1011
0.9292
4.4325
0.8841
4.3186
0.8859
5
3.0096
0.9346
2.9349
0.9367
4.1616
0.8968
4.0632
0.8983
7
2.8908
0.9395
2.8255
0.9414
3.9980
0.9034
3.8775
0.9070
9
2.7862
0.9438
2.7142
0.9463
3.8954
0.9079
3.7363
0.9142
11
2.7175
0.9465
2.6499
0.9491
3.8022
0.9118
3.6778
0.9171
13
2.6734
0.9483
2.6035
0.9510
3.6925
0.9168
3.5877
0.9214
15
2.5818
0.9522
2.5221
0.9544
3.6090
0.9206
3.4930
0.9256
17
2.5190
0.9548
2.4562
0.9571
3.5427
0.9237
3.4186
0.9287
19
2.4948
0.9557
2.4234
0.9584
3.5096
0.9252
3.3812
0.9305
21
2.4699
0.9567
2.3924
0.9596
3.4649
0.9273
3.3376
0.9327
23
2.4349
0.9582
2.3564
0.9613
3.4130
0.9296
3.2937
0.9348
25
2.4149
0.9593
2.3368
0.9624
3.3528
0.9325
3.2409
0.9374
27
2.4099
0.9600
2.3340
0.9633
3.3095
0.9347
3.2044
0.9393
29
2.4193
0.9608
2.3521
0.9639
3.2972
0.9364
3.1944
0.9408
31
2.5733
0.9594
2.4856
0.9629
3.5217
0.9315
3.4524
0.9356
Best results for each data set are highlighted
123
5×2
Int J CARS
9
3x2 kernel (sq) 5x2 kernel (sq) 7x2 kernel (sq) 3x2 (round)
1.6
3.8
1.55
3.6
8
3.4
7
1.45
MSE
1.5
MSE
MSE
3.2 3 2.8
6 5
2.6
4
1.4 1.35
2.4 0
5
10
15
20
25
2.2
30
3 0
5
10
Window Size
(a)
15
20
25
0
30
5
10
(b)
DS1, 32 ACS, R=2
15
20
25
30
25
30
Window Size
Window Size
(c)
DS1, 32 ACS, R=3
DS1, 32 ACS, R=4
5.5
2.2
14 2.15
5
12
2.1
2
4.5
MSE
MSE
MSE
2.05
4
10 8
1.95
3.5
6
1.9 1.85
0
5
10
15
20
25
3
30
4 0
5
10
20
25
0
30
5
10
Window Size
Window Size
(d)
15
(e)
DS2, 32 ACS, R=2
15
20
Window Size
(f)
DS2, 32 ACS, R=3
DS2, 32 ACS, R=4
Fig. 5 Reconstruction performance (MSE) for different size N of the removed window Ck [N × N ] for different sizes of the reconstruction kernel, 32 ACS lines and different acceleration rates for two data sets (DS1 and DS2)
3.5
MSE
3
MSE
DS1 (3x2) DS1 (5x2) DS2 (3x2) DS2 (5x2) DS3 (3x2) DS3 (5x2)
2.5 2 1.5
9
9
8
8
7
7
6
6
MSE
4
5 4
4
3
3
2
1 0
5
10
2 0
15
5
Window Size
(a)
10
15
0
(b)
MSE, 16 ACS, R=2
MSE, 16 ACS, R=3
0.96
0.96
0.98
0.94
0.94
0.92
0.92
0.95 0.94 0.93
SSIM
SSIM
0.9 0.96
0.88
(d)
10
Window Size SSIM, 16 ACS, R=2
15
15
(c)
MSE, 16 ACS, R=4
(f)
Window Size SSIM, 16 ACS, R=4
0.9 0.88
0.86
0.86
0.84
0.84
0.82 5
10
Window Size
0.99
0
5
Window Size
0.97
SSIM
5
0.82 0
5
(e)
10
Window Size SSIM, 16 ACS, R=3
15
0
5
10
15
Fig. 6 Reconstruction performance (MSE and SSIM) for different size N of the removed window Ck [N × N ] for different sizes of the reconstruction kernel size, 16 ACS lines and different acceleration rates for the three data sets
123
Int J CARS 2.6
DS1 DS2 DS3
2.4 2.3
2.4
g−factor
g−factor
2.2 2.1 2
2.2
2
1.9 1.8
1.8
1.7 1.6
0
5
10
15
20
25
30
1.6
0
5
Window Size (a) 32 ACS lines
10
15
Window Size (b) 16 ACS lines
Fig. 7 Study of the GRAPPA g-factor variations for different size N of the removed window Ck [N × N ] for the three data sets. 16 and 32 ACS lines and acceleration rate R = 2 are considered Table 3 Optimal size N of the removed window Ck as a function of the number of ACS lines and the acceleration factor R
ACS R
32 2
32 3
32 4
24 2
24 3
24 4
16 2
16 3
16 4
DS1
27
27
27
21
19
19
13
11
11
DS2
29
29
27
21
19
19
13
11
11
DS3
21
21
21
15
15
17
11
11
9
DS1 + σ10
31
31
31
23
23
19
15
11
11
DS2 + σ10
31
31
27
23
21
21
13
11
11
DS3 + σ10
27
25
23
19
19
17
13
13
11
N for min{MSE}
Average
27.6
27.3
26.0
20.3
19.3
18.6
13.0
11.3
10.6
Median
28
28
27
21
19
19
13
11
11
12m4 scanner with FGRE Pulse Sequence, matrix size 128 × 128, TR/TE = 8.6/3.38, FOV21 × 21cm. For comparison purposes, the three data sets are normalized so that the composite magnitude signal (CMS) after sum of squares [5] is bounded in [0, 255]. The following experiments were carried out: 1. Test of the reconstruction quality as a function of the removed window size. We simulate accelerated acquisitions with 32 and 16 ACS lines and uniform acceleration factor R = 3 (effective acceleration factors of 2.37 and 2.67, respectively, for the 256 × 256 data sets). Large acceleration rate is considered to highlight the differences in reconstruction. A square window of size N × N is removed from the center of the ACS lines for estimation. N ranges from 1 to 31 (for 32 ACS lines) and from 1 to 15 (for 16 ACS lines). The reconstruction kernel is assumed to be 3 × 2, the standard de facto in GRAPPA acquisitions. The experiment is repeated for a round window of diameter N , different kernel sizes and different acceleration rates. For numerical comparison, two quality indexes are used: the mean squared error (MSE) [7]
123
and the structural similarity index (SSIM) [16]. The first index gives a measure for the global error in the image, while the second compares the structural similarities. For the latter, the closer to one, the better the reconstruction is. The comparison is done over the CMS in the x-space. To avoid any influence of the background over the quality measures, only the signal part of the data is considered for testing by a manual masking of the data. In addition, in order to use a parallel imaging specific metric, the quantitative g-factor for GRAPPA proposed in [2] is also considered, specifically, the g-factor for each coil defined as: SNRfull |W 2 W H |kk k = (4) gk = √ SNRacc R | 2 |kk k · with 2 the correlation matrix of the data and W a matrix derived from the GRAPPA coefficients. For the sake of simplicity, 2 will be assumed to be the identity matrix. To obtain a single g-factor value, the average value along the coils is considered. 2. Validation of the assumption for the optimum size for N proposed in Eq. (3). We simulate acquisitions with 16,
Int J CARS
24 and 32 ACS and acceleration rates of 2, 3 and 4. We will use the MSE of the reconstructed signal as a quality measure. To increase the number of experiments and to test low SNR data, Gaussian noise is added to each of the signals, so that max(SNR) =
max{Sl (x)} = 25 σ
Table 4 Comparison of reconstruction quality of the proposed method (FD-GRAPPA) compared to standard GRAPPA, for different number of ACS lines and different acceleration rates R Acceleration
2
Results Results for the first experiment are given in Table 1 (32 ACS, kernel 3 × 2, square and round windows), Table 2 (32 ACS, kernels 5 × 2 and 7 × 2, square window), Fig. 5 (32 ACS and different parameters) and Fig. 6 (16 ACS). N = 0 denotes the case in which no window is removed from the ACS, i.e., the standard GRAPPA.
N/A
GRAPPA
FD-GRAPPA
SSIM
(5)
3. Test of the reconstruction quality for a fixed window for different acceleration rates and number of ACS lines. The data are subsampled at three different rates, R = 2, R = 3 and R = 4, and 32, 24 and 16 ACS lines are considered. The size of the removed window will be dynamically selected as a function of R and the number of lines accordingly to the results of the previous experiment. For comparison, MSE and SSIM are used. 4. Comparison of the proposed method with other GRAPPA improvements from literature: fast robust (FR) GRAPPA [11], high-pass GRAPPA (HP) and nonlinear (NL) GRAPPA [3]. All the algorithms are similarly built in MATLAB, and MSE and execution time are measured. As an illustration of the possibility of using the proposed method together with existing ones, we will merge it with FR-GRAPPA and HP-GRAPPA. 5. Finally, we test the influence of the SNR level of the ACS data on the reconstruction quality. It has recently been shown in [6] that the presence of noise in the ACS data can have an effect similar to that of a Tikhonov regularization under some SNR conditions. Thus, there are situations in which the reduction of the SNR of the ACS lines during the estimation step yields to an improvement in the quality of the reconstruction. The method here presented eliminates the center of the k-space, precisely those areas with higher SNR, so it implicitly operates over ACS lines with a smaller global SNR. To prove that the gain of the method is not produced by a side effect of the SNR reduction, which may be equivalent to the Tikhonov regularization, the following experiment is conducted: The MSE of the reconstructed data using the proposed approach will be compared to the MSE of the standard GRAPPA algorithm when noise is added to the ACS lines (only for estimation).
ACS
3
4
32
0.9386
0.9799
0.9849
24
0.9177
0.9791
0.9836
16
0.8900
0.9783
0.9821
32
0.9040
0.8949
0.9549
24
0.8720
0.8891
0.9450
16
0.8231
0.8834
0.9284
32
0.8768
0.6755
0.9071
24
0.8338
0.6604
0.8857
16
0.7799
0.6470
0.8443
32
5.3911
1.6061
1.4414
24
6.9047
1.6491
1.5007
16
9.2958
1.6911
1.5971
32
6.8794
3.9291
2.5925
24
8.9175
4.0859
2.8495
16
12.5444
4.2597
3.4628
32
7.8753
9.1302
4.0412
24
10.1017
9.6119
4.6149
16
13.5533
10.1026
5.8176
MSE 2
3
4
The size of the window is dynamically chosen. N/A stands for nonparallel imaging applied (Data set 1)
For the 32 ACS acquisition, see Table 1, note that as N grows the MSE decreases and the SSIM index increases, showing both an improvement in the final quality of the reconstructed image. Note that the MSE reduces up to a 35 % for both data sets. When the window size is too large, close to the number of ACS, we can also detect a slight degradation of the indexes (for DS1), due to an excessive loss of data for estimation, but even then, there is a clear improvement when compared to the original GRAPPA method. In addition, we can see that there is no any particular advantage of using of a round window instead of a square one. When larger GRAPPA reconstruction kernels are considered, see Table 2, the observed behavior is similar; there is an increase in the quality of the reconstructed signal, with slight degradation for the highest window size. Both data sets show a small improvement for larger reconstruction kernels. However, since a larger kernel size goes along with an increase in the complexity of the algorithm (in terms of estimation and reconstruction), in what follows we will use the 3 × 2 kernel size. The experiment is repeated for different acceleration rates and results depicted in Fig. 5 for the sake of comparison. In all cases, square windows outperform round windows, while larger kernels show a smaller error in reconstruction. A similar behavior can be observed for the 16 ACS lines experiment in Fig. 6. In this case, the degradation of the
123
Int J CARS Table 5 Comparison of reconstruction quality of the proposed method (FD-GRAPPA) compared to standard GRAPPA, for different number of ACS lines and different acceleration rates R
Table 6 Comparison of reconstruction quality of the proposed method (FD-GRAPPA) compared to standard GRAPPA, for different number of ACS lines and different acceleration rates R
Acceleration
Acceleration
ACS
N/A
GRAPPA
FD-GRAPPA
SSIM 2
3
4
3
4
GRAPPA
FD-GRAPPA
32
0.9361
0.9666
0.9744
32
0.9012
0.9492
0.9687
24
0.9148
0.9652
0.9720
24
0.8792
0.9450
0.9660
16
0.8845
0.9641
0.9691
16
0.8493
0.9403
0.9624
32
0.8937
0.8432
0.9312
32
0.8836
0.8750
0.9516
24
0.8588
0.8345
0.9127
24
0.8603
0.8635
0.9461
2
16
0.8213
0.8528
0.9341
32
0.8401
0.6278
0.9055
3
16
0.7978
0.8257
0.8872
32
0.8690
0.5463
0.8745
24
0.8292
0.5243
0.8319
24
0.8056
0.6019
0.8924
16
0.7682
0.4927
0.7672
16
0.7563
0.5729
0.8346
32
5.7027
2.1928
1.9449
32
6.5418
2.2255
1.7804
24
7.5406
2.2514
2.0602
24
7.9071
2.3317
1.8660
16
10.7988
2.3118
2.1883
16
10.0525
2.4495
1.9656
32
7.3040
5.3963
3.4619
32
8.3820
3.7654
2.3124
24
9.8880
5.6480
4.0515
24
10.0498
3.9869
2.4743
16
15.0715
5.9588
4.7574
16
12.5803
4.1925
2.7711
32
8.4066
14.6476
5.2268
32
9.4617
8.4333
3.4573
4
MSE 2
3
4
24
11.0553
15.8408
6.5275
24
11.4645
9.0305
3.7595
16
15.8305
17.6657
8.5257
16
14.6068
9.8200
4.9346
The size of the window is dynamically chosen. N/A stands for nonparallel imaging applied (Data set 2)
quality when the removal is excessive is much clearer, due to a smaller size of the estimation area and the smaller amount of data available for estimation. Nevertheless, it also shows a great reduction of the error when the central part of the k-space is removed. An interesting issue can also be noticed when comparing results in Fig. 5 and figure:MSESSIM16. Results for DS1 in Fig. 6b show that the estimation with 16 ACS lines removing a 11 × 11 window is even better than a traditional GRAPPA with 32 ACS acquisition in Fig. 5b. Similar behaviors can be also seen for R = 2 and R = 4. The proper choice of the samples used for estimation increases the accuracy of the reconstruction. The number of ACS lines can be reduced and even the acceleration rate can be increased, getting a higher final effective acceleration rate. Results for the GRAPPA g-factor for R = 2, different window sizes and different number of ACS lines are depicted in Fig. 7 for the three data sets. Note that the g-factor will only depend on the GRAPPA reconstruction coefficients, and not on the quality of the reconstruction itself. The consequence of the removal of the center of the k-space for estimation is a decrease in the g-factor. This will imply an increase in the SNR of the reconstructed signal. Results for the second experiment are given in Table 3. The optimal size for N is selected as
123
N/A
SSIM
MSE 2
ACS
The size of the window is dynamically chosen. N/A stands for nonparallel imaging applied (Data set 3)
Nopt = arg min{MSE}. N
The average and the median for each pair [ACS,R] are considered. The median results for the experiments are totally consistent with Eq. (3). However, there are particular cases in which the maximum lays in smaller sizes. The selected size will depend on several parameters such as the acquisition modality and probably the SNR. Therefore, a calibration step may be necessary. Results for the third experiment are given in Tables 4, 5 and 6. FD-GRAPPA stands for frequency discriminated GRAPPA, the proposed method. The optimal window size in Eq. (3) is used. N/A denotes the case in which the unsampled lines are simply padded with zeros (no reconstruction method is applied). In all the cases, the MSE and SSIM improve when using FD. This is particularly significant for high acceleration rates, like R = 4, where the original GRAPPA shows an smaller SSIM than the non-reconstructed case. The error of the proposed method is significantly smaller than GRAPPA. For the sake of illustration, some visual results for DS1 and DS3 are shown in Figs. 8 and 9. There is a clear improvement in the visual quality of the reconstructed data when FD is used. Note that even for higher acceleration rates (R=4) FDGRAPPA is able to recover images with an acceptable visual quality, when the traditional GRAPPA fails. Similar results
Int J CARS
Fig. 8 Reconstruction of subsampled MRI data from eight coils using GRAPPA and FD-GRAPPA. Data set 1
can be seen in the doped ball phantom, for R = 3. We can appreciate a texture pattern inside the ball in the traditional GRAPPA reconstruction that does not appears in FD. Results for the comparative experiment are given in Table 7 for different data sets and acceleration rates. For FR-GRAPPA, a 10 % outlier rejection is selected (the one that gives the best performance in the three data sets). For HP-GRAPPA, optimal result has been found for c = 10 and w = 2. The time is expressed in seconds and measured over an Intel Core I7 processor with 12 Gb of RAM. From the results, we can observe that again, the use of FD improves the MSE when compared with GRAPPA, but
also note that the acquisition is accelerated, due to the reduction in the number of equations used to estimate the weights. The proposed method even improves FR-GRAPPA and HPGRAPPA for most configurations, while being much faster than the former. In addition, note that, when FD is used together with FD + FR, the MSE is even better and the execution time is reduced. Result of the NL-GRAPPA algorithm is, in all cases, slightly better that the proposed method. However, note that this improvement goes with a huge processing time when compared with FD-GRAPPA (161 vs. 0.24 s for the first experiment). Although reconstruction time may not me an issue for off-line reconstruction, it is for real scanning
123
Int J CARS
Fig. 9 Reconstruction of subsampled MRI data from eight coils using GRAPPA and FD-GRAPPA. Data set 3
applications. If more than 2 min are needed to reconstruct a single slice, a whole volume or a dMRI acquisition will take excessive waiting time. In addition, note that NL fails for data set 3, where the object is uniform and without details. Finally, results for the SNR experiment are shown in Figs. 10 and 11. Results for standard GRAPPA in Fig. 10 show that, even for the optimal SNR reduction, the MSE is higher than in the FD-GRAPPA case. This leads us to conclude that the improvement of the method is not only due to the SNR reduction, but there are some other factors involved;
123
precisely the frequency discrimination proposed. While the effect of the removal of the center of the k-space is an effective reduction of the SNR, the considered values do not change from GRAPPA to FD-GRAPPA. The points used for estimation are those of a lower SNR, but their values are exactly the same than those used for the standard GRAPPA algorithm. Thus, they are still susceptible of further regularization. In Fig. 11, results for this regularization over FD-GRAPPA are shown. Note that the new method can additionally benefit from a reduction of the SNR of the ACS lines, achieving even smaller error rates.
Int J CARS Table 7 Comparison of different GRAPPA improvement algorithms for the three data sets and different subsampling configurations
GRAPPA
FD
FR (10 %)
FD + FR
HP
FD + HP
NL
Data set 1, 32 ACS R = 3 MSE
3.9291
2.5991
3.0227
2.6674
2.5684
2.5931
2.3008
SSIM
0.8949
0.9549
0.9344
0.9569
0.9547
0.9555
0.9637
Time
0.26
0.24
23.61
9.36
0.24
0.24
161.9
Data set 1, 32 ACS R = 4 MSE
10.1026
5.8176
5.7214
6.0203
7.5123
7.7761
8.4592
SSIM
0.6470
0.8389
0.8483
0.8708
0.8569
0.8541
0.8572
Time
0.21
0.18
15.15
16.22
0.18
0.18
104.6
Data set 2, 32 ACS R = 3 MSE
5.3963
3.4619
4.0193
3.4760
3.4941
3.4514
3.1346
SSIM
0.8432
0.9312
0.9043
0.9336
0.9292
0.9325
0.9412
Time
0.26
0.24
23.94
18.32
0.24
0.23
153.3
Data set 2, 32 ACS R = 4 MSE
17.6657
8.5257
7.7726
7.8526
10.1311
10.9858
20.8219
SSIM
0.4927
0.7672
0.8026
0.8265
0.8087
0.8022
0.6736
Time
0.20
0.19
18.12
18.47
0.26
0.25
107.1
Data set 3, 16 ACS R = 2 MSE
2.4495
1.9794
2.4023
1.9872
2.9926
3.6665
4.8660
SSIM
0.9403
0.9624
0.9428
0.9623
0.9497
0.9465
0.9009
Time
0.09
0.09
1.84
0.90
0.05
0.05
16.1
Data set 3, 16 ACS R = 4 MSE
9.8200
5.0585
8.9192
5.8734
7.7021
8.9561
125.4791
SSIM
0.5729
0.8307
0.6360
0.8173
0.8162
0.7979
0.1031
Time
0.11
0.10
9.03
6.27
0.09
0.09
47.4
The time is expressed in seconds, measured over an Intel Core I7 processor with 12 Gb of RAM 4
5.5
6
5.5
5
5
4.5
4.5
5.5
4
MSE
MSE
MSE
MSE
3.5
4
5
3 4.5
3.5
2.5
3.5
3 0
1000
2000
3000
4000
5000
6000
(a)
4
3
0
1000
2000
Noise (σ)
3000
4000
5000
6000
0
2000
Noise (σ)
(b)
D1, 32 acs R=3
4000
6000
8000
0
10000
1000
2000
(c)
D1, 16 acs R=3
3000
4000
5000
6000
Noise (σ)
Noise (σ)
(d)
D2, 32 acs R=3
D2, 16 acs R=3
Fig. 10 Variation of the MSE for a reduction of the SNR of the ACS lines. FD-GRAPPA (dashed) versus standard GRAPPA with noise in the ACS lines (solid) 2.7
3.45
2.6
3.5 3.4
2.5 0
100
200
300
400
500
600
3.2
3.25
3.1
3.2
0
500
Noise (σ)
(a)
D1, 32 acs R=3
3.35
4.6
4.5
3.3
3.3
2.55
4.7
3.4
MSE
MSE
3.6
MSE
2.65
MSE
4.8
3.5
3.7
1000
1500
2000
4.4
4.3 0
500
Noise (σ)
(b)
D1, 16 acs R=3
1000
1500
0
500
(c)
D2, 32 acs R=3
1000
1500
2000
Noise (σ)
Noise (σ)
(d)
D2, 16 acs R=3
Fig. 11 Variation of the MSE for a reduction of the SNR of the ACS lines. FD-GRAPPA (dashed) versus FD-GRAPPA with noise in the ACS lines (solid)
123
Int J CARS
Discussion and conclusions The ACS lines commonly used to estimate the GRAPPA coefficients comprise high-, medium- and also low-frequency samples, whereas the lines to be retrieved by interpolation range from medium to high frequencies. As a result, the estimation of the interpolation weights becomes biased toward the high-energy, low-frequency spectrum. In this paper, we have presented a simple but powerful method based on the elimination of the center samples of the ACS lines for the estimation step. This simple reduction of the data goes together with a highly improvement of the accuracy of the GRAPPA coefficients and therefore the final image reconstruction. At the same time, a reduction in the number of equations will accelerate the estimation. The method shows better results than traditional GRAPPA even with the half of the ACS lines. The improvement in the reconstruction quality makes it feasible using higher acceleration rates. One relevant issue about the method is the calculation of the shape and size of the removal window. While experiments show that the quality improvement is a fact as long as the limitation in Eq. (5) is hold, the optimum value will depend on the area of the body and on the coil configuration. We have provided here a general solution, but we consider it must be further tested for different regions and coil configurations. In addition, the proposal shows best results for the quality indices considered, MSE and SSIM. If other quality measures are considered, the optimal size may vary, but we are confident it will still hold the limitation in Eq. (5). The method presented shows some interesting additional features. First, it is totally compatible with standard GRAPPA; the reconstruction step is exactly the same, while the estimation uses the same methodology with a smaller system of equations. Furthermore, the method proposed is also compatible with other method in the literature, such as Fast Robust GRAPPA or High-pass GRAPPA. Results have shown that these methods improve their quality when used together with FD-GRAPPA. The method can also benefit from regularizations such as the SNR reduction of the ACS lines recently proposed by Ding et al. Finally, the simplicity of the method makes it suitable to be easily incorporated in commercial scanners in which GRAPPA is already implemented. Acknowledgments The authors acknowledge Ministerio de Ciencia e Innovación for grant TEC2013-44194, Institute of Health Carlos III for grant PI11-0149 and Centro de Técnicas Instrumentales, Universidad de Valladolid for MRI acquisitions. Gonzalo Vegas-Sánchez-Ferrero acknowledges Consejería de Educación, Juventud y Deporte of Comunidad de Madrid and the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007– 2013) for REA grant agreement no. 291820.
123
Conflict of interest The authors declare that they have no conflict of interest in relation to this article.
References 1. Blaimer M, Breuer F, Mueller M, Heidemann R, Griswold M, Jakob P (2004) SMASH, SENSE, PILS, GRAPPA: how to choose the optimal method. Top Magn Reson Imaging 15(4):223–236 2. Breuer FA, Kannengiesser SA, Blaimer M, Seiberlich N, Jakob PM, Griswold MA (2009) General formulation for quantitative gfactor calculation in GRAPPA reconstructions. Magn Reson Med 62(3):739–746 3. Chang Y, Liang D, Ying L (2012) Nonlinear GRAPPA: a kernel approach to parallel MRI reconstruction. Magn Reson Med 68(3):730–740 4. Chen Z, Zhang J, Yang R, Kellman P, Johnston LA, Egan GF (2010) IIR GRAPPA for parallel MR image reconstruction. Magn Reson Med 63(2):502–509 5. Constantinides C, Atalar E, McVeigh E (1997) Signal-to-noise measurements in magnitude images from NMR phased arrays. Magn Reson Med 38:852–857 6. Ding Y, Xue H, Ahmad R, Chang Tc, Ting ST, Simonetti OP (2014) Paradoxical effect of the signal-to-noise ratio of GRAPPA calibration lines: a quantitative study. Magn Reson Med. doi:10.1002/ mrm.25385 7. Eskicioglu A, Fisher P (1995) Image quality measures and their performance. IEEE Trans Commun 43(12):2959–2965 8. Griswold MA, Jakob PM, Heidemann RM, Nittka M, Jellus V, Wang J, Kiefer B, Haase A (2002) Generalized autocalibrating partially parallel acquisitions (GRAPPA). Magn Reson Med 47(6):1202–1210 9. Hoge WS, Brooks DH, Madore B, Kyriakos WE (2005) A tour of accelerated parallel MR imaging from a linear systems perspective. Concepts Magn Reson Part A 27A(1):17–37 10. Huang F, Li Y, Vijayakumar S, Hertel S, Duensing GR (2008) Highpass GRAPPA: an image support reduction technique for improved partially parallel imaging. Magn Reson Med 59(3):642–649 11. Huo D, Wilson D (2008) Robust GRAPPA reconstruction and its evaluation with the perceptual difference model. J Magn Reson 27(6):1412–1420 12. Ji JX, Son JB, Rane SD (2007) PULSAR: a matlab toolbox for parallel magnetic resonance imaging using array coils and multiple channel receivers. Concepts in Magn Reson Part B Magn Reson Eng 31B(1):24–36 13. Lim JS (1990) Two dimensional signal and image processing. Prentice Hall, Englewood Cliffs 14. Park S, Park J (2012) Adaptive self-calibrating iterative GRAPPA reconstruction. Magn Reson Med 67(6):1721–1729 15. Wang H, Liang D, King KF, Nagarsekar G, Chang Y, Ying L (2012) Improving GRAPPA using cross-sampled autocalibration data. Magn Reson Med 67(4):1042–1053 16. Wang Z, Bovik AC, Sheikh HR, Simoncelli EP (2004) Image quality assessment: from error visibility to structural similarity. IEEE Trans Image Process 13(4):600–612 17. Wu C, Hu W, Kan R, Yu J, Sun X (2011) An improved GRAPPA image reconstruction algorithm for parallel MRI. In: Control and decision conference (CCDC), 2011 Chinese, pp 4096–4100