Improved parallel magnetic resonance imaging reconstruction with multiple variable density sampling

Generalized auto-calibrating partially parallel acquisitions (GRAPPA) and other parallel Magnetic Resonance Imaging (pMRI) methods restore the unacquired data in k-space by linearly calculating the undersampled data around the missing points. In order to obtain the weight of the linear calculation, a small number of auto-calibration signal (ACS) lines need to be sampled at the center of the k-space. Therefore, the sampling pattern used in this type of method is to full sample data in the middle area and undersample in the outer k-space with nominal reduction factors. In this paper, we propose a novel reconstruction method with a multiple variable density sampling (MVDS) that is different from traditional sampling patterns. Our method can significantly improve the image quality using multiple reduction factors with fewer ACS lines. Specifically, the traditional sampling pattern only uses a single reduction factor to uniformly undersample data in the region outside the ACS, but we use multiple reduction factors. When sampling the k-space data, we keep the ACS lines unchanged, use a smaller reduction factor for undersampling data near the ACS lines and a larger reduction factor for the outermost part of k-space. The error is lower after reconstruction of this region by undersampled data with a smaller reduction factor. The experimental results show that with the same amount of data sampled, using NL-GRAPPA to reconstruct the k-space data sampled by our method can result in lower noise and fewer artifacts than traditional methods. In particular, our method is extremely effective when the number of ACS lines is small.


Results
To demonstrate the effects of the MVDS method proposed in the paper, in vivo data sets were used. Figure 1 shows the standard reference image obtained by using SoS to reconstruct the experimental brain data.
Our method is better in the case of lower ACS. When R nom =4 and ACS = 16, the image reconstructed by NL-GRAPPA with traditional sampled method and its different map were shown in Fig. 2(a,d). Used the VDS method to undersample data, reduce the number of ACS lines from 16 to 12, and add 6 lines of under-sampled data with ORF = 2. The kernel sizes are 15 × 2. The reconstructed image and its different map were shown in Fig. 2(b,e). Artifact power (AP) and Signal-To-Noise Ratio (SNR) values were shown in Table 1. We found that too few ACS lines make the reconstruction effect worse than traditional methods. We could see from the image reconstructed by the data sampled by traditional method already have extremely serious artifacts that were completely unrecognizable. Used proposed method to sample data and increase the number of undersampled data with ORF = 2 and 6 from 0. Figure 2(c,f) were the result of undersampling 20 lines of R = 2 and its difference map. The effect of collecting different numbers of lines with ORF = 2 and 6 on the AP value of the reconstructed image was shown in Table 1 and Fig. 3. It can be seen from Fig. 3 that the AP value of the image shows a downward trend as the number of undersampled data with ORF = 2 on both sides of the central area increases. It indicated that the more data collected near the central area with ORF = 2, the better the quality of the reconstructed image. As can be seen from the table, the AP value after reconstruction with the traditional method is 9.59, and with the MVDS method can reduce the AP to 0.44, the effect was obvious. Figure 2 Fig. 2(a), the artifacts are much reduced. We can see that the MVDS method is obviously better than VDS and traditional methods. When R nom =4 and ACS = 20, the image reconstructed by NL-GRAPPA with traditional sampled method and its different map were shown in Fig. 4(a,d). We found that as the ACS decreases, the reconstructed image artifacts become more serious. Using the VDS method to undersample data, reduces the number of ACS lines from 20 to 16, and adds 6 lines of under-sampled data with ORF = 2. The reconstructed image and its different map were shown in Fig. 4(b,e). Its AP and SNR values were shown in Table 2. We found that due to the reduction of ACS lines, the effect is worse than traditional methods. Sampled data with the MVDS method, reduced the ORF on both sides of the ACS data in the central area, use ORF = 2 to undersample data; increase the ORF at the outermost k-space, use ORF = 6 to undersample data. Figure 4(c,f) were the result of undersampling 20 lines of R = 2. From this we can see that the reconstructed image of Fig. 4(c) was clearly superior to the traditional method.     www.nature.com/scientificreports/ Figure 5 showed with the same sampling lines, the AP of reconstruction results using different kernel sizes. We can see that the quality of MVDS method reconstruction is higher under larger kernel size. When R nom =6 and ACS = 24, another set of eight-channel data 11 reconstructed by NL-GRAPPA with traditional uniform subsampling and MVDS were shown in Fig. 6. The kernel sizes are 15 × 2. The standard reference image was shown in Fig. 6(a). Figure 6(b,d) are the image reconstructed with traditional sampling method and its different map. Using the proposed MVDS method to subsample data, we increased the ORF from 6 to 8 on the outermost k-space and reduced the ORF from 6 to 4 on the both sides of ACS region. The result image reconstructed with MVDS and its different map were shown in Fig. 6(c,e). The AP and SNR values were shown in Table 3. We could see that the image quality of MVDS method was improved compared with the traditional method.
The proposed method was applied to GRAPPA. In addition to NL-GRAPPA, the MVDS method was also effective in GRAPPA. In the case of R = 4 and ACS = 40, the GRAPPA reconstruction images sampled by the three   Figure 7(a,d) were the result image of GRAPPA with traditional sampling method and its difference map. Figure 7(b,c) were the result image with VDS and the proposed method, and Fig. 7(e,f) were their difference maps. The kernel sizes are 15 × 2. The results of sampling number and AP value are shown in Table 4. We still found that the proposed method was better than traditional methods. Compared with other methods, SC-SENSE and JSENSE are well-known image domain reconstruction algorithms. Compared the reconstruction image of proposed method with SC-SENSE and JSENSE, the results were shown in the Fig. 8(a-c). The data in the Table 5 also showed that our image is better than other methods.   SPIRiT is a method that can reconstruct images from arbitrary sampling patterns 30 . SAKE is a calibration data-free parallel imaging reconstruction method 31 . When R = 4 and ACS = 32, the results of SPIRiT and SAKE reconstructed images sampled by traditional uniform subsampling were shown in Fig. 9(a,e), and their different maps were shown in Fig. 9(c,g). The results of SPIRiT and SAKE images sampled by proposed method were shown in Fig. 9(b,f), and their different maps were shown in Fig. 9(d,h). The kernel size in SPIRiT was 5 × 5 and the iterations was 30. The kernel size in SAKE was 6 × 6. We could see that there are still some improvements in the Fig. 9 and Table 6.

Discussion
In this paper, we proposed a novel sampling pattern based on an improvement to the traditional sampling pattern. Unlike the traditional sampling pattern, which has only a single reduction factor, we proposed that the sampling pattern includes multiple reduction factors. While keeping the number of ACS unchanged, the MVDS method reduces the value of ORF close to the central k-space by increasing the value of ORF in the outermost k-space, so that the error after reconstruction of the central k-space is smaller. Experimental results prove that the proposed method had resulted in lower noise and fewer artifacts with the same amount of sampled data. At the same time, even when the amount of data sampled by the traditional method is slightly larger, the MVDS method can still maintain the advantages of lower noise and fewer artifacts. The proposed method is more effective at higher reduction factors. At the same time, we only improved the sampling method, and do not involve changes to the reconstruction algorithm. Therefore, it can be directly applied to algorithms similar to the GRAPPA or the algorithms improved from GRAPPA.
The VDS method proposed by Jaeseok increases the number of undersampled data with ORF = 2 by reducing the number of ACS lines in the central area. This method can obviously outperform the traditional method when there are enough ACS lines. But once there are few ACS lines in the traditional method, reducing the number of ACS lines in the VDS method will produce more significant artifacts than the traditional method due to too few ACS lines. Figures 2 and 4(a,b) showed that when there are fewer ACS lines, the image reconstructed by the VDS method is not as good as the traditional method. Figures 2 and 4(a,c) were fully illustrated of the advantages of the MVDS method.
Most of the information of the image is contained in the k-space 29 , so the smaller error in the center of the k-space will enhance image quality. In the method proposed in this paper, the data near the ACS region was sampled with a smaller value of R. The reconstruction results with lower ORF have lower noise 11 . The data sampled by these smaller R and the ACS lines occupy a wider region of the k-space center as compared to traditional methods that only use ACS lines. Thus, the reconstruction data of this part of the k-space region in the proposed method (lower R) contains less noise than the traditional method (higher R). In the reconstructed images of Figs. 2, 4 and 6 and the results shown in Tables 1, 2 Fig. 3, after reducing the ORF value close to the ACS area, the AP of the reconstructed image can be reduced, which means that the quality of the image is improved. When the amount of data undersampled with lower ORF in the central area increases, the AP value shows a downward trend. From the Tables 1 and 2, we can see that when the amount of data sampled by different ORFs are similar, the quality of the reconstructed image will be significantly improved compared to the traditional methods. At the same time, it can be seen from Table 1 that when the numbers of lines with three ORFs are similar, changing the number of lines has less of improvement effect on the results. This may be due to the fact that when there are lower ORF data, the number of corresponding higher ORF data must also increase, which makes the area of higher ORF larger and leads to the large error of this part of the area. This part of the increased error has begun to reduce the quality improvement brought by the increase of lower ORF. Although the outermost region contains less information, the larger error of a large number of higher ORF regions will also reduce the quality of the results in the whole image.
Therefore, we recommend that when sampling data, use the lower ORF in the center of k-space to undersample the appropriate data. For example, sampling a similar number of lines from different ORFs may be effective. We can see from Fig. 5 that MVDS will perform better when the kernel is enlarged, and the reconstruction quality is better than the traditional method using the best kernel size.
In order to keep the collection time consistent with traditional methods, we need to increase the ORF in the outermost region. Since the outermost region of k-space contains less information, increasing the ORF value of this region will not have a significant impact on the results. The proposed method does not reduce the number of ACS lines, so it can still maintain its advantages compared with traditional methods when there are fewer ACS. The method proposed in this paper can effectively reduce the artifacts after reconstruction without prolonging the acquisition time (i.e. the same amount of sampling data).

Material and methods
The source data were drawn from http:// bi. tamu. edu/ softw are/ downl oads_ softw are. htm (University of Texas A&M, USA). To test its performance, in vivo data from a healthy volunteer were acquired on a 3 T MR system (GE Healthcare Technologies, Waukesha, WI, USA) using an eight-channel head coil and a 3D spoiled gradientecho (SPGR) pulse sequence (TE = 10 ms, TR = 300 ms, RBW = 16 kHz, FOV = 22 cm × 22 cm, matrix = 256 × 256) 32 . Informed consent was obtained from the volunteer in accordance with the institutional review board policy. All methods were carried out in accordance with relevant guidelines and regulations. All experimental protocols were approved by the institutional review board (IRB) at Hangzhou Dianzi University (IRB-2020002). SMASH is an earlier reconstruction approach based on k-space domain. AUTO-SMASH and VD-AUTO-SMASH were proposed to improve this approach. GRAPPA can be seen as a more generalized implementation of VD-AUTO-SMASH [33][34][35] . The GRAPPA approach needs to acquire a block of extra ACS lines in the central of k-space when undersampling the outer k-space data by some reduction factors (R). The value of the reduction factor R indicates that one line of the data is sampled in every R lines in the outer k-space. These ACS lines are used to determine the complex weights of the linear fit. This estimation process can be simplified as a matrix equation: where A represents the matrix of acquired data of ACS line, b represents the "missing data" but are actually acquired from ACS lines, and x represents the linear weights to be fitted.
The estimation process was shown in the Fig. 10, we used the ACS data around the "missing data" to fit the "missing data" to estimate the weights, and then used theses weights to estimate the truly missing data from the undersampled data.
NL-GRAPPA improves on traditional GRAPPA and introduces extra nonlinear mapping steps. Specifically, the undersampled k-space data A are mapped to a high-dimensional space H through a nonlinear mapping Φ (·): A-> H, and then the unacquired data is reconstructed. The weights of the linear fit can also be estimated from the ACS data. The relationship between the selected kernel and the mapping Φ is as follows: where <,> represents the inner product. There are many different types of kernels 36 . This method chooses the a polynomial kernel 11,37 , it can be expressed as where the parameter value is γ = r = 1 and d = 2.
The formula of nonlinear GRAPPA method using this kernel can be expressed as follows: Figure 10. The basic principle of the GRAPPA algorithm is shown above. The sampled data around a sampled point in the ACS data is used to calculate the weight of these sampled points to the point. Then use the sampled data around the missing point to use the weight to calculate the value of the missing point. ORF = 2 in this case. where S j (k y + r∆k y ,k x ) denotes the missing k-space data at the target coil, S j (k y + tR∆k y ,h∆k x ) denotes the acquired undersampled data, and w i,j (l,t,h)denotes the weights of linear combination. The first-order term is equivalent to the conventional GRAPPA, and the second-order term can represent other nonlinear models. Here R represents the reduction factors, l counts all coils, t and h represents the number of reconstruction blocks, k x and k y represent the coordinates along the frequency-encoding and phase-encoding directions 11 . At the same time, the formula for estimating the fitting weights b = Ax can be transformed to where the matrix Φ(A) is size of M × N k , M represents the number of all matching equations during estimation, N k is the dimension of the new high-dimensional feature space, and a i is the row vector of A.
(4) in the full space is very slow. Because frequency encoding steps are much faster than phase encoding steps in MRI scanners, one way to accelerate MRI data acquisition is to reduce the number of phase encoding lines. This type of undersampled dataset is equal to fold images along phase encoding direction and yield aliasing artifacts along the undersampled direction with a reduced field of view (FOV) as seen as Fig. 11. Fully sampling results in a full FOV image after Fourier transformation, and undersampling yields a reduced FOV image with aliasing artifacts. The traditional sampling pattern is to choose a certain reduction factor for uniform undersampling to acquire data. The sampling pattern of SC-SENSE, GRAPPA, NL-GRAPPA and other methods are similar. While sampling the undersampled data, some ACS data need to be obtained at the Nyquist rate in the central k-space. SC-SENSE uses these ACS data to extract a more accurate sensitivity calibration images. GRAPPA, NL-GRAPPA uses these ACS lines to calculate the coefficients of the linear fit. We use the expected reduction factor (termed nominal reduction factor, R nom ) to represent the outer reduction factor (ORF) used in traditional methods, which means that data is collected once every ORF rows.
Variable density sampling. Park et al. improved the traditional sampling method 26 , using a lower ORF on both sides of the ACS lines for undersampling data. As shown in the Fig. 12(a) is the traditional sampling method and (b) is the method proposed by Park et al., and the R nom =4. In Fig. 12(b), Park et al. increased the undersampled data with ORF = 2 on both sides of the central area by reducing the number of full sampled ACS lines. Note that the full sampled ACS lines as well as the data undersampled with ORF = 2 is larger than the k-space occupied by the ACS lines in (a). And this part of the area occupies most of the central area of k-space. Therefore, the error after reconstruction is small, and the image quality will be improved.
The VDS method proposed by Park et al. keeps the sampling time consistent with the traditional method by reducing the number of ACS lines sampled. The redundant time of sampling ACS data will be reduced to sample the data with ORF = 2. In this way, the error of reconstructed data at the ORF = 2 area becomes smaller. It makes the reconstruction image quality better than the traditional method. However, once the number of ACS lines in the traditional method is very small (for example, ACS = 24, ACS = 16), then reducing the ACS lines will make the results produce extremely significant artifacts. Therefore, we propose a method to keep the number of sampled ACS lines consistent with the traditional method and increase the value of ORF in outermost k-space to reduce the ORF value on both sides of the ACS area. The improved method is applied to NL-GRAPPA, and the effect is much better than the traditional method in the case of lower ACS lines. The MVDS method keeps the number of ACS lines in the center of the k-space unchanged while using VDS in the outer k-space for data sampled. The proposed method was based on the theory that higher reduction factors will cause higher errors and that most of the image information is contained in the central area of k-space. The MVDS method specifically increases the value of the ORF of the data undersampled in the outermost k-space to decrease the value of the ORF undersampled near the ACS lines. The number of ACS lines and the total amount of data sampled remain unchanged. The data undersampled using this method will be more accurate near the center of k-space after reconstruction due to the lower ORF. This method was shown in Fig. 13. Taking R nom =4 as an example, Fig. 13(a) was the traditional sampling method, and Fig. 13(b) showed the MVDS method proposed in this paper. Please note that the number of ACS lines fully sampled by this method is consistent with the traditional method. In the traditional method, we increase the ORF = 4 to R = 6 in the outermost k-space, and reduce the ORF = 4 to ORF = 2 on both sides of the ACS region. In k-space, most of the image information is in the central region. Therefore, a lower ORF value in the central region will make the image error smaller. Unlike the VDS method in Fig. 12(b), the method was to keep the outermost R = 4 unchanged, and reduce the number of ACS lines sampled at the innermost side to reduce the value of ORF on both sides of the lower central area from 4 to 2 to keep the sampling time unchanged. But we reduced the value of ORF on both sides of the central area while we increased the value of ORF in the outermost k-space to keep the sampling time unchanged. The MVDS method does not change the amount of ACS lines, so it can also produce excellent results when the traditional method collects fewer ACS lines. The difference between the MVDS method and the traditional method is that we do not change the number of ACS lines. When using ACS data to calculate weights, more ACS lines make the calculated weights more accurate, and too few ACS lines will produce significant artifacts. Therefore, when the ACS line is sufficient in the VDS method, the reconstruction data error of ORF = 2 area is smaller and the reconstruction result is better. But once the ACS lines used in the traditional method are less, reducing the number of ACS data will produce more serious artifacts to the result. The MVDS method does not change the number of ACS lines, and reduces the ORF value in the central area. This part of the area will have a lower error than the traditional method after restoration. Therefore, the MVDS method can produce excellent results when fewer ACS lines are collected.
The method of restoring the missing k-space data using the data sampled by the MVDS method is similar to the traditional method. As shown in Fig. 14, the first step was to use the corresponding reconstruction method (GRAPPA or NL-GRAPPA, etc.) to estimate the weight of the corresponding ORF, and the second step was to use the weight to calculate the missing data of the corresponding region. In Fig. 5, the thick line was the sampling www.nature.com/scientificreports/ data, and the thin line is the missing value restored by the reconstruction method. After these two steps, the complete k-space data can be obtained. The proposed method in this paper was based on the theory that most image information, including contrast and general shape, was contained in the center of space 29 . We know that a higher acceleration factor will lead to higher noise 11 . Therefore, if undersampling the data in the region near the center with a lower reduction factor, the error in this region will be smaller than traditional after reconstruction. In the traditional method, if the reduction factor is high, the region of reconstruction data except ACS lines will have a large error. The proposed method kept the ACS lines unchanged and sampled data with lower ORF values on both sides of the central area. Therefore, the error in the center of k-space is smaller and the reconstructed image is more accurate. A higher ORF is used for sampling the outermost area of the k-space region to keep the total acquisition time consistent, but it has little influence on the reconstruction images because of the less information contained in this part of the region. The MVDS method is superior to the traditional method when the number of ACS lines is fewer.
Image acquisition and analysis. The full sample k-space data of the eight coils were first to generate the corresponding coil image by performing inverse Fourier transform. We used the sum of squares method to process the images of eight coils and generated the SoS image as a standard image reference. The calculation method is as follows: In the formula, I j represents an image generated by inverse Fourier transform of k-space data of each coil, Coil-Num represents the number of coils, and I sos represents the generated standard image reference.
Each pixel in each reconstruction image is subtracted from the corresponding pixels in the SoS image to obtain equal scale difference maps. These maps can more directly show the severity and distribution of noise and the location of artifacts. In this experiment, the AP and SNR were used to measure the quality of the reconstructed image. The AP is used to describe the error between the reconstructed image and the standard reference image. A higher AP value indicates the reconstructed image has greater noise. The calculation method of AP is as follows: where the I recon represents the reconstructed image, I ref represents the standard reference image, and ROI means the region of interest. In this experiment, we selected the entire image matrix region, which is ROI = 256. It can be seen from the formula that a smaller the value of AP means the reconstructed image is closer to the standard reference image, that is, the higher the reconstruction error. The signal-to-noise ratio SNR of an image is also an indicator to measure the quality of the image. The larger the SNR value indicates the less noise of the reconstructed image. The formula is as follows: In pMRI we use the value of reduction factor to indicate the degree of pMRI being accelerated. The termed nominal reduction factor ( R nom ) represents that sample k-space data every R lines. The larger the value of R nom means the less the total amount of data sampled, and the less time pMRI costs. Some methods such as GRAPPA need to sample extra ACS data in the central k-space, leading to an increase in the total amount of data sampled. These increase data make the actual net reduction factor ( R net ) and R nom different. In the traditional sampling method, the ACS lines are fully sampled in the center of k-space (R = 1), and the R nom is used to sample data in the outer k-space region. Therefore, R net is usually smaller than R nom that we expected. The number of sampled lines and the calculation method of R net are as follows: where N fully represents the total number of lines of k-space data, and N ACSL represents the number of extra ACS lines. SUM represents the amount of data lines sampled.

Conclusion
We proposed an MVDS method that uses multiple reduction factors in the sampling process to improve the reconstruction quality of pMRI. Areas sampled with a factor lower than the nominal reduction factor can have less noise. The proposed method effectively reduces noise and artifacts, especially when there are fewer ACS lines. We apply the proposed method to NL-GRAPPA and GRAPPA, and compare them with traditional sampling methods, SC-SENSE, JSENSE and other methods. The experimental results show that this method has obvious advantages compared with the traditional sampling method and VDS method. This method only improves the sampling pattern, and does not change the algorithm. The sampling pattern introduced here can be used in most existing technologies similar to GRAPPA. At the same time, our method has some shortcomings. For example, it is currently only suitable for k-space based imaging algorithms. Our method needs to be improved overtime with the goal of achieving a more universal method.