Statistical Noise Analysis in GRAPPA Using a Parametrized Noncentral Chi Approximation Model

Statistical Model in Single- and Multiple-Coil MR Signals

k–space data is typically acquired through repeated application of excitation pulses with a different phase encoding gradient for each readout gradient. Each sampled line of k–space is “frequency encoded,” and the measured signal is uniformly sampled at the desired rate. As a consequence, points in the k–space measured from the MRI scanner are independent samples of an RF signal measured with a receiver coil, for which the noise can be modeled as a complex stationary additive white gaussian noise process, with zero mean and variance ${\mathrm{\sigma }}_K^2$.

For a single–coil acquisition the complex spatial MR data is therefore also modeled as a complex Gaussian process and the magnitude signal M(x) is the Rician distributed envelope of the complex signal (20). In the image background, where the SNR is zero due to the lack of water-proton density in the air, the Rician simplifies to a Rayleigh distribution. In a multiple–coil MR acquisition system, from a merely statistical point of view, the acquired signal in the k–space data for each coil can be modeled as:

with a_l (k) the original signal if no noise is present and $n_l(\mathbf{k};{\mathrm{\sigma }}_K^2)=n_{l_r}(\mathbf{k};{\mathrm{\sigma }}_K^2)+j{n}_{l_i}(\mathbf{k};{\mathrm{\sigma }}_K^2)$ complex uncorrelated Gaussian processes with zero mean and variance ${\mathrm{\sigma }}_K^2$. The noise is initially assumed to be stationary, and therefore parameter ${\mathrm{\sigma }}_K^2$ does not depend on k.

The complex x–space image is obtained via the 2D inverse discrete Fourier Transform of s_l (k) for each slice. Since the discrete Fourier Transform is a linear operator, the noise in the complex signal in the x–space for each coil will also be Gaussian, and it can be expressed as:

where $N_l(\mathbf{x};{\mathrm{\sigma }}_n^2)=N_{l_r}(\mathbf{x};{\mathrm{\sigma }}_n^2)+j{N}_{l_i}(\mathbf{x};{\mathrm{\sigma }}_n^2)$ is a complex Gaussian process with zero mean and variance ${\mathrm{\sigma }}_n^2$. It is easy to prove the relation (21):

where |Ω| is the size of the image in each coil, i.e., the number of points used in the 2D discrete Fourier Transform.

If the k–space is fully sampled, the composite magnitude image may be obtained using methods such as SoS (5,22):

Defining $A_T(\mathbf{x})=\sqrt{{\sum }_{l=1}^L{|A_l(\mathbf{x})|}^2}$, and assuming the noise components to be identically and independently distributed, the envelope of the magnitude signal M_L will follow a non–central Chi (nc-χ) distribution (5) with probability density function (PDF):

with I_L(․) the Lth order modified Bessel function of the first kind and u(․) the Heaviside step function. Eq. 4 reduces to the Rician distribution (20) for L = 1. In the background, this PDF simplifies to a central Chi distribution with PDF:

As a final remark, note that three requirements are needed for the SoS of Gaussian random variables (RV) to be modeled as a stationary nc-χ:

1. The noise is stationary in each of the Complex Gaussian x–space images. If the k–space data is fully sampled, the noise variance will be the same for all the points in the image in both the k–space and x–space domains, and the noise can be assumed to be stationary.

2. The variance of noise ${\mathrm{\sigma }}_n^2$ is the same for each of the coils.

3. No correlation is assumed between the Gaussian RVs.

We will see that these requirements are not fulfilled when GRAPPA is employed. One of the aims of parallel imaging is precisely to accelerate the image acquisition process by sub–sampling data in each coil. From a statistical point of view, such a reconstruction will affect the underlying model and the stationarity of the noise in the reconstructed data.

Statistical model in GRAPPA reconstructed images

The GRAPPA reconstruction strategy estimates the missing lines in a sub–sampled k–space data (14,17,23) acquisition for each coil. While the sampled data $s_l^{𝒮}(\mathbf{k})$ remain the same, the reconstructed lines $s_l^{ℛ}(\mathbf{k})$ are estimated through a linear combination of the existing samples. Weighted data in a neighborhood η(k) around the estimated pixel from several coils is used for such an estimation:

with s_l (k) the complex signal from coil l at point k and ω_m(l, k) the complex reconstruction coefficients for coil l. In self-referenced reconstructions, these coefficients are determined from the low-frequency coordinates of k–space, termed the Auto Calibration Signal (ACS) lines, which are sampled at the Nyquist rate (i.e., unaccelerated). Breuer et al. in (24) pointed out that Eq. 6 can be rewritten using the convolution operator:

where w_m(l, k) is a convolution kernel that can be easily built from the GRAPPA weight set ω_m(l, k). Since a (circular) convolution in the k–space is equivalent to a product into the x-space, we can write (24):

with W_m(l, x) the GRAPPA reconstruction coefficients in the x-space. The variance of noise in the x-space of the sampled data is ${\mathrm{\sigma }}_n^2=\frac{{\mathrm{\sigma }}_K^2}{|\mathrm{\Omega }|·R}$, with R the reduction rate. Note that:

1. As a consequence of the coefficients W_m(l, x) (which have different values for each pixel) the reconstructed signal in each coil $S_l^{ℛ}(\mathbf{x})$ becomes a non-stationary Gaussian process, i.e., a complex distributed Gaussian image with different parameters (mean and variance) in each point of the image:

   $$E\{S_l^{ℛ}(\mathbf{x})\}=A_l^{ℛ}(\mathbf{x})=|\mathrm{\Omega }|\sum _{m=1}^L{A}_m^{𝒮}(\mathbf{x})\times W_m(l,\mathbf{x})$$

   [11]

   $$\text{Var}\{S_l^{ℛ}(\mathbf{x})\}={|\mathrm{\Omega }|}^2[\sum _{m=1}^L{\mathrm{\sigma }}_m^2\times {|W_m(l,\mathbf{x})|}^2+2\sum _{p\ne q}{\mathrm{\sigma }}_{pq}^2ℛe\{W_p(l,\mathbf{x})W_q^{*}(l,\mathbf{x})\}]={\mathrm{\sigma }}_l^{2^{ℛ}}(\mathbf{x}).$$

   [12]

   Therefore, the stationarity requirement in described in section “Statistical Model in single- and multiple-coil MR signals” fails.

2. For the same pixel x, the mean and the variance will also vary along the coil dimension, and therefore the second requirement for the nc-χ model fails.

3. Even assuming that the coils are initially uncorrelated, signals $S_l^{ℛ}(\mathbf{x})$ will be correlated due to the GRAPPA reconstruction. Note that the signal in each coil is a linear combination of the signals at all remaining coils, hence the third requirement also fails.

The composite magnitude image M_L(x) can be obtained using SoS. For convenience, in what follows we will equivalently consider $M_L^2(\mathbf{x})$ instead of M_L(x) (if M_L follows a nc-χ distribution, $M_L^2$ will follow a nc-χ squared).

After GRAPPA reconstruction, the interpolated lines at each coil $S_l^{ℛ}(\mathbf{x})$ can be seen as strongly correlated non-stationary complex Gaussian RVs. Therefore, from a theoretical point of view the magnitude image after SoS, M_L(x), will not strictly be a nc-χ, although it is usually modeled that way in the literature (18).

The SoS image can be expressed as

with

Signal $M_L^2(\mathbf{x})$ can be seen as the SoS of the sum of weighted Normal variables. A deep study of the resultant distribution is done in Appendix A. The mean and the variance (for each x) of the resulting distribution are:

where ${\mathbf{C}}_X^2=\mathbf{W}{\mathbf{\Sigma }}^2{\mathbf{W}}^{*}$ is the covariance matrix of the interpolated data at each spatial location. Accordingly, Σ² is the covariance matrix of the original data at each coil. ‖․‖_F is the Frobenious norm, and A(x) = [A₁, A₂, …, A_L]^T is the noise-free reconstructed signal, from which we define:

Although the resultant distribution is not strictly a nc-χ², our hypothesis is that its behavior will be very similar and could be modeled as such (the error in this approximation will be quantified in the results section). However, the high correlations between the reconstructed signals in each coil should translate in a decrease of the number of Degrees of Freedom (DoF) of the distribution. In Appendix, the method of the moments is used to show that the reduction in the number of DoF translates in an effective number of coils (obviously reduced, since the number of coils L is clearly related to the DoF) and an effective variance of noise:

Finally, note that the origin of the reduced DoF of the nc-χ² model is in the correlation and inhomogeneous variance of the complex gaussians, i.e., in ${\mathbf{C}}_X^2$ not being of the form σ²I with I the identity matrix. With GRAPPA the distortion comes mainly from the interpolation matrix W which is not diagonal. For unaccelerated acquisitions, we may write W = I, so that ${\mathbf{C}}_X^2={\mathbf{\Sigma }}^2$. Even in this case, Σ² is not necessarily diagonal, since a certain correlation does exist between the signals at each channel especially for those systems with a large number of receiving coils L. This means that effective parameters should be considered even if no subsampling is performed, as has been empirically shown by (18). We can summarize our developments as follows:

1. The nc-χ model does not hold for GRAPPA reconstructed data. However, this distribution can be used as a good approximation of the actual one.

2. For this approximation to hold, effective parameters have to be considered which represent an equivalent, non-subsampled configuration with a smaller number of coils (DoF) and, consequently, a greater level of noise.

3. Even when the nc-χ model is feasible, the resulting distribution is non-stationary (the effective parameters are spatially dependent).

Data Sets

To validate the hypothesis posed in the previous section, the following data sets are considered for the experiments, see Fig. 1:

• Data set 1: 1 repetition of eight-channel head coil data acquired on a GE Signa 1.5T EXCITE 11m4 scanner, FGRE Pulse Sequence, TR = 500 ms, TE = 13.8 ms, matrix size = 256 × 256, FOV = 20 × 20 cm, slice thickness = 5 mm.

• Data set 2: the MR brain data-set from (25) (the authors provide their data set on-line¹); eight-channel head array from 3 Tesla GE scanner, fast spoiled gradientecho sequence, TR/TE = 300/10 ms, RBW = 16 kHz, matrix size = 256 × 256, FOV 22 × 22 cm.

• Data set 3: 60 repetitions of the same slice of a phantom, scanned in an 8-channel head coil on a GE Signa 1.5T EXCITE 12m4 scanner with FGRE Pulse Sequence to generate low SNR. matrix size = 128 × 128, TR/TE = 8.6/3.38 ms, FOV21 × 21 cm, slice thickness = 1 mm.

In addition, synthetically generated RV will also be considered.

Synthetic Experiments

First, we numerically validate the nc-χ approximation proposed and measure the error introduced with this approach. To that end synthetically generated RVs are considered in the following experiments and the following distributions are compared:

• A nc-χ² distribution with the effective parameters L_eff and ${\mathrm{\sigma }}_{\mathrm{e}\mathrm{f}\mathrm{f}}^2$.

• A nc-χ² distribution with the original parameters (L = 8).

• A Normal distribution fitted using the mean and the variance in Eqs. 15 and 16.

The validation experiments are as follows:

1. For the sake of illustration, we fit the histogram of the real distribution to the proposed three statistical models. We consider 5 · 10⁴ samples of 8 complex Gaussian RVs with zero mean and unitary variance, N(0, 1). The RVs are combined using 8 × 8 complex random weights (W ~ N(1, 4)) and then combined using the (square) SoS. The experiment is repeated with the same parameters, but considering a signal value of A_i = 3, i = 1,… 8.

2. We will analyze the error committed with the nc-χ² approximation for different configurations of matrix C². To avoid any effect of the signal over the noise, we first assume that μ_i = 0. We define 10⁴ samples of 8 complex Gaussian RVs and a 8 × 8 coefficient matrix W, so that three different scenarios are considered (see Appendix for more details about the relation between matrix C² configuration and the final distribution):

   1. Matrix C² is diagonal (different elements in the diagonal).
   2. Matrix C² is not diagonal and the diagonal elements are greater than the rest of the elements (diagonal dominant).
   3. Matrix C² is totally random.

   Since we are assuming unitary noise variance and no correlation, note that ${\mathbf{C}}^2={\mathbf{C}}_X^2$. These configurations will be measured by two parameters:

   1. The coefficient of variation of the elements of the diagonal: $\mathrm{C}{\mathrm{v}}^2\{C_l^2\}$.
   2. The ratio between the determinant of C² and the trace of the matrix: ‖C²‖/Tr(C²). This is a measure of the weight of the non-diagonal elements over the diagonal. Note that for the first case, this parameter is a constant.

   The RVs are combined following Eq. 14. The analysis in terms of the elements of matrix C² is used to characterize whether a nc-χ² can be accurately fitted in each particular situation.² The following (relative) Mean Squared Error (rMSE) is measured:

   $$\mathrm{r}\mathrm{M}\mathrm{S}\mathrm{E}=\frac{{\int }_{\infty }^{\infty }{|p_O(x)-p_z(x)|}^{2}dx}{{\int }_{\infty }^{\infty }{|p_z(x)|}^{2}dx}$$

   [19]

   where pO(x) is the PDF of the original data (estimated from the histogram of the 10⁴ samples for each pixel) and p_z(x) is the distribution to fit.

3. We will test the model fit for different SNR values with a given set of weights. We consider 10⁴ samples of 8 complex Gaussian RVs with unitary variance and mean in the range μ ∈ (0, 3), i.e., ~N(μ, 1). The RVs are combined using 8 × 8 complex random weights (selected in a range similar to an actual GRAPPA reconstruction) and the SoS is considered. We calculate the error of fitting a nc-χ² and a Gaussian for different SNR values, defining the original SNR in each coil as $\mathrm{S}\mathrm{N}\mathrm{R}=\frac{\mathrm{\mu }}{\mathrm{\sigma }}=\mathrm{\mu }$.

Experiment with the MR Data-Sets

In a second stage, GRAPPA coefficients from actual MR acquisitions will be used. In each case, a rectangular 2 k_y - by- 5 k_x GRAPPA convolution kernel η(k) is employed. The following experiments are considered:

1. Data sets 1 and 2: The k–space data are 2× subsampled with 32 ACS lines–skipping every other line except in the central region. From the ACS lines the GRAPPA coefficients ω_m(l, k) are estimated. The k–space domain coefficients are then transformed to x-space, W_m(l, x). The experiments proposed for synthetic data are repeated using the calculated GRAPPA coefficients. For each of the pixels in the image we consider 10⁴ samples of 8 complex Gaussian RVs with unitary variance and mean in the range μ ∈ (0, 4), so that original SNR ∈ (0, 4). To simulate the reconstruction process, for each of the values of x the RVS are combined following Eq. 9, and the composite magnitude image is obtained by SoS. For the sake of simplicity, $M_L^2(\mathbf{x})$ will be used. After reconstruction, 10⁴ samples of the 256 × 256 image are available. The distribution of the data is tested by fitting a theoretical probability distribution to the real data. To that end, the rMSE is measured. The same three distributions used in the synthetic experiments (nc-χ² with effective parameters, nc-χ² with the original parameters and Gaussian) are considered.

2. Data set 1: To further test it, the experiment for Data set 1 is repeated with SNR = 1. A variable number of samples is considered, and a Kolmogorov–Smirnov test at a significance level of α = 0.05 is used to decide whether the sample data follow any of the proposed distributions.

3. Data set 3: The GRAPPA reconstruction coefficients are derived from one of the 60 samples, using 2× subsampling and 32 ACS lines. All the 60 samples are 2× subsampled and reconstructed with the same GRAPPA coefficients and the composite magnitude image is obtained by SoS. Agolden standard (to obtain the A_i values) is created by averaging the complex data in each coil. The original σ_n and σ_K values for each coil are estimated from the background of the fully sampled images (6). Theoretical values of L_eff (x) and ${\mathrm{\sigma }}_{\mathrm{e}\mathrm{f}\mathrm{f}}^2(\mathbf{x})$ are calculated for each point x. A Kolmogorov–Smirnov test at a significance level of α = 0.05 is used to decide whether each x point of the 60 samples follow a nc-χ distributions with parameters L_eff (x) and ${\mathrm{\sigma }}_{\mathrm{e}\mathrm{f}\mathrm{f}}^2(\mathbf{x})$.

Synthetic Experiments

First, we will consider the experiments with synthetically generated weights. Results for the first experiment are shown in Fig. 2. Note that the nc-χ² distribution with effective parameters is the distribution that better follows the actual variations of the data for both experiments. The more asymmetric the distribution, the better is the behavior of the nc-χ² when compared with the Gaussian.

The second experiment measures the error of the different approximations as a function of the configuration of matrix C², i.e., the $C_{pq}^2(\mathbf{x})$ coefficients. Results are depicted in Fig. 3. In all the cases, the nc-χ² distribution with effective values is the one that shows a lower rMSE for a wider range of variation of parameters. In all the cases, the error of approximating the actual distribution by a nc-χ² is always smaller than 0.1, usually around 0.05. In Fig. 3a (diagonal case) we can see that as the values in the diagonal are more different, the Gaussian and the nc-χ² with original parameters produce a greater error. However, the approximation with effective parameters remains constant. Fig. 3(b–c) shows the results for the dominant diagonal, and Fig. 3(d–e) the results for the random C². An interesting effect in these two experiments is that for the original nc-χ² and the Gaussian approaches the error increases when the ratio ‖C²‖/Tr(C²) grows. It means that, as the diagonal becomes less significant, these two approaches are not a good representation of the data. The nc-χ² approximation with effective parameters due precisely to the fine tune of these effective parameters, is able to properly fit the actual distribution.

The evolution of the effective number of coils for the experiment in Fig 3(d–e) is reproduced in Fig. 4. As the ratio ‖C²‖/Tr(C²) grows, the correlation between variables also grows. As a consequence, the number of DoF of the system will decrease together with the effective number of coils, L_eff. This reduction in L_eff implies an increase in ${\mathrm{\sigma }}_{\mathrm{e}\mathrm{f}\mathrm{f}}^2$ and consequently a decrease in the SNR.

The third synthetic experiment deals with the error of the approximation for different SNR configurations. Results are depicted in Fig. 5. For large SNR values (SNR > 2), the nc-χ² and the Gaussian approximation converges, as expected.

Experiment with the MR Data-Sets

The next step is to analyze the results generated by estimated GRAPPA coefficients. Results for the data fit for different distributions are shown in Fig. 6 for the data set 1 and in Fig. 7 for the data set 2. In both cases, the rMSE for each of the points of the final image is shown. For better illustration, the average of the rMSE in the whole image is depicted in Fig. 8. Regardless of the SNR level, the error committed when approximating the distribution by a nc-χ² with effective parameters is always lower than 6%, while the other two models (Gaussian and nc-χ²) show greater values. What is more, the rMSE for nc-χ² is quite independent of the SNR value. As expected, the Gaussian model is a good approximation when SNR increases. On the other hand, if nc-χ² is considered but the effective parameters are not used, the error even increases together with the SNR in certain areas of the image.

From these results, we can conclude that the nc-χ model (with effective parameters) proposed gives a small error when considering the GRAPPA reconstruction coefficients from actual parallel acquisitions.

It is also interesting analyzing the distribution of the effective parameters along the image. Note that these parameters only depend on the GRAPPA coefficients. The effective number of coils (L_eff) is shown in Fig. 9 and the effective variance of noise $({\mathrm{\sigma }}_{\mathrm{e}\mathrm{f}\mathrm{f}}^2)$ is shown in Fig. 10. Note that those cases with low SNR also show low levels of the effective number of coils. Note also that the distribution of L_eff correlates with the error committed when considering nc-χ² without effective parameters. The points with lower rMSE are logically those with greatest values of L_eff.

The variance of noise, on the other hand, shows a great variation along the image, which confirms the non-stationarity assumption previously done. Values of σ_n range between 2 and 6. However, note that the variation of this parameter across the image is soft. So, in some applications (like filtering methods) local stationarity may be assumed if needed.

In Fig. 11, the acceptance rate by the Kolmogorov-Smirnov test for different set sizes is depicted. The acceptance rate for the nc-χ² model (with effective parameters) is over the 90% even for large population sizes. The larger the size of the sample, the easier for the test to reject the null hypothesis. For the Gaussian and nc-χ² models, the rate goes down very fast, so that they can only be accepted for very small populations. It may be concluded that the effective nc-χ² model accurately fits data generated following the GRAPPA reconstruction algorithm.

Finally, the Kolmogorov-Smirnov test is run over the 60 GRAPPA reconstructed samples of the MR phantom (data set 3). The 90.13% of the points in the image passed the test, which correlates well with the results of the previous experiment, and confirms the nc-χ (with effective parameters) as a suitable distribution to accurately model data after GRAPPA reconstruction and SoS.