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

Improving GRAPPA reconstruction by frequency discrimination in the ACS lines.

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-s...
3MB Sizes 2 Downloads 7 Views