Experiments and Results
We perform three experiments with two objectives: first, quantifying the effect of phase correction on signal debiasing, by processing real data from a HCP dataset; second, assessing the debiasing on diffusion metrics, calculated with DTI and MAP reconstructions, on a digital dataset generated for typical scenarios such as DTI at bvalue 1000, 2000, 3000 s/mm^{2} and DTI and MAP multishell.
In the first experiment, we clustered a HCP dataset to obtain typical signal values at bvalue 1000 and 3000s/mm^{2}. Particularly, for each bvalue we applied kmeans to divide the signal intensities of the DWIs—accounting for all the gradient directions—into four clusters. We used the centroid of each cluster to define respectively background, low, medium, and high mean signal values. Based on these, we created a groundtruth synthetic magnitude image—M_{xy} in Eq. (5)— composed of three circles each containing, from left to right, low, medium and high signal respectively. Outside the circles we added background signal. A synthetic phase was generated and the noisy complex DWI was created. After calculating the average b = 0 signal (S(0)_{avg} = 758 a.u.) in the HCP dataset, noise was added as in Eq. (5) in low SNR regime: SNR_{0} = 10. Figure 1 shows, for each bvalue, the noisy magnitude DWI_{xy} and the estimated phasecorrected real part <(DWIpy) (A set to 0.75 after visual inspection). In addition, an effective SNR map is present along with histograms of the magnitude and phasecorrected real signals for each circle. We conclude that in both cases the phasecorrected real image presents more contrast with the background compared to the magnitude. This is more evident at low SNR values—left circle at b = 1000 s/mm^{2}, left and central circles at b = 3000 s/mm^{2}—that are more likely with high bvalues. The <(DWIpy) shows darker colors, i.e. lower signal intensities, as highlighted by the histograms: the magnitude (green line) has a Rician distribution for low SNRs (typically below SNR = 5) whereas the estimated real part (blue line) always shows a Gaussian distribution, thus including negative signal intensities. We point out that since this is an experiment grounded on real data, the centroid of the clusters—especially at low signal values—are based on Rician data and might overestimate the actual (noise free) ones. This means that the Rician bias in histograms (green line) might be an underestimation of the true one.
In the second experiment, we use the HCP dataset to create a mean groundtruth magnitude DWI, M_{xy}, in order to quantify the Rician bias, i.e. the distance from Gaussianity. We calculate the mean b = 0 image S(0)xy, for a slice of interest, by averaging the 40 nondiffusionweighted images in the dataset. Since the SNR is very high, the averaging procedure is not biased. Then, we select all the DWIs corresponding to b = 1000 s/mm^{2} and perform DTI to obtain the mean diffusivity map, MD_{xy}. At this point, we obtain a groundtruth magnitude DWI at any bvalue by extrapolating with M_{xy} = S(b)_{xy} = S(0)_{xy} exp(—b ? MD_{xy}). Although we assume Gaussian isotropic diffusion, this phantom represents an average description of a
Fig. 1 The signal contrast and distributions for synthetic complex DWIs, at b = 1000,
3000 s/mm^{2}, created from clustering a real HCP dataset (SNR_{0} = 10). In the rectangular frames from left to right: the SNR map, the Rician magnitude (Mg) and the phasecorrected estimated real image (Re). Below, the histograms of the signal intensities corresponding to the circles with low, medium, and high signal/SNR, for Mg (green) and Re (blue). Background SNR: 21.5 for b = 1000 s/mm^{2} and 3.2 for b = 1000 s/mm^{2}. *: s/mm^{2}
magnitude DWI. After generating a synthetic phase image, as described in Sect. 2, we calculate the noisy complex DWI for each bvalue, as in Eq. (5). Figure 2 shows the b = 0 magnitude and the phase used for the phantom (left column). We then calculate the Rician magnitude DWI^y and the phasecorrected real part <(DWI?p (A = 0.75). Additionally, we calculate a magnitude image with Gaussian distributed noise DWI^, by adding Gaussian noise (with the same SNRo) to M_{xy}. This will be used as reference for Gaussianity measures. We generate 1000 noise occurrences and calculate, for each pixel of the images, the signal intensities histogram (as in Fig. 1). For each pixel we generate three histograms: one for the Rician DWI_{xy}, one for the phasecorrected <(DWE(y), and one for the Gaussian DWI!y The hypothesis is that, for each pixel, the phasecorrected signal distribution should be closer to that of the Gaussian magnitude image than the Rician magnitude one. As
Fig. 2 The distance from Gaussianity of complex DWIs obtained by processing a HCP dataset and a synthetic phase image. In the first column the b = 0 magnitude image obtained from real data, and the generated phase. In the columns from the second to the fourth, the distance from the Gaussianity measured with Eq. (6) for the Rician magnitude (first row) and the phasecorrected real part (second row), at different bvalues (columns). Contrarily to the case of the Rician magnitude, the distance from Gaussianity remains visually unchanged as the diffusionweighting increases. *: s/mm^{2}
a distance from Gaussianity, we use the discrete Hellinger distance [18]
where P and Q are two discrete probability distributions, and 0 < H(P, Q) < 1 where 1 means maximum distance. In columns 24, Fig. 2 shows the maps of Hellinger distance from the Gaussian magnitude, for the Rician magnitude (first row) and the phasecorrected real part (second row), at bvalue 1000, 2000 and 3000 s/mm^{2}.
We see that the Rician magnitude shows more bias (higher H), especially in regions where MD is high. As expected, at higher bvalues (from left to right) the signal intensity is lower and the bias occurs in a larger number of pixels. Conversely, the phasecorrected real part does not show a clear change.
In the third experiment, we quantify the bias on the estimated DTI and MAP metrics for the Rician magnitude, and we quantify the debiasing power of phase correction by looking at the change in the distributions of such metrics compared to the Gaussian noise case. We generate complex DWIs with Phantomas [15] as described in Sect. 2. For each gradient direction g = (g_{x}, g_{y}, g_{z}), the synthetic phase image is oriented towards v = (g_{x}, g_{y}) [see Eq. (4)] with constant phase shifts between slices (along the z direction). Figure 3 shows the original noisy data (Rician magnitude, phase, real and imaginary parts) and the one after phase correction for a reference slice (b = 1000 s/mm^{2}). For each SNR_{0} e {10,20,30}
Fig. 3 A slice of data generated with the digital phantom for SNR_{0} = 10, b = 1000 s/mm^{2}. On the left: the original noisy data calculated with a groundtruth magnitude image obtained with Phantomas [15] and a synthetic phase. On the right: the phasecorrected real and imaginary parts (A = 0.75); the signal information is almost entirely contained in the real part, whereas the imaginary part mainly contains Gaussian noise
we generate the Rician magnitude data and, as for the second experiment, the Gaussian DWI image to be used as reference. In this experiment we also investigate the effect of the regularization parameter A of the total variation filtering in Eq. (3). Therefore, for each SNR_{0} we generate six phasecorrected datasets, for A e {0.25,0.5,0.75,1,2,5}. For each combination of SNR_{0} and type of data— Rician, Gaussian and the six phase corrections—we fit DTI with singleshell scheme at bvalue 1000, 2000 and 3000 s/mm^{2}, and DTI and MAP with multishell scheme. For DTI, we calculate the mean diffusivity (MD), the principal diffusivity (PD), and the fractional anisotropy (FA). We calculate qspace metrics based on closed formulas derived for MAP. These are the return to origin (RTOP), axis (RTAP) and plane (RTPP) probabilities [3], the mean squared displacement (MSD) and the qspace inverse variance (QIV) [17]. We create a mask of voxels within fibers, based on the noisefree dataset, by considering only the voxels where RTOP e [0.5e6,0.7e6] (range chosen based on visual inspection). For each value of A and for each calculated DTI and qspace metric, we compute the probability distribution inside the mask. So we do for the Rician and Gaussian data. Figure 4 illustrates the influence of the regularization parameter A (decreasing along the
Fig. 4 Histograms of the principal diffusivity (PD)—for DTI at 1000, 2000 s/mm^{2}—and return to plane probability (RTPP)—for MAP—estimated on Gaussian DWI (“Ga”, red), Rician magnitude (“Mg”, green), and phasecorrected real part (“Re”, blue), SNR_{0} = 10. While the red and green histograms remain unchanged along the rows, the blue histograms change as function of the regularization parameter X [see Eq. (3)]. As the attachment to data decreases (from top to bottom) the phasecorrected histograms overlap more with the red Gaussian ones. *: s/mm^{2}
rows) on the recovered metric’s probability distribution. The figure shows the Gaussian (red line), Rician (green line), and phasecorrected (blue line) probability distributions (SNR_{0} = 10) of PD (DTI at 1000, 2000 s/mm^{2}) and RTPP (MAP). The results confirm the underestimation of PD that increases with the bvalue, i.e. the green histograms are left to the red ones. Consequently, also MD (see Sect. 1) is underestimated [1]. Inverse analogous considerations hold for RTPP. We observe that X has a great influence on the phase correction results. Particularly, a large X implies strong attachment to data, resulting in a poor phasecorrection since the estimated lowfrequency phase is very similar to the original noisy one, /DWI„, & /DWIry in Eq. (2). Indeed, the blue histograms (phasecorrected) in the first row of Fig. 4 almost entirely overlap the green ones (Rician magnitude data). As the attachment to the data decreases (from top to bottom), the blue phasecorrected histograms move towards the (red) Gaussian based distributions, visually reducing the distance from Gaussianity. As in the second experiment, we quantify the distance from Gaussian metrics by using Hellinger’s formula in Eq. (6). Figure 5 illustrates the variation of the H distance for the phasecorrected data as function of X, for each acquisition setup, reconstruction method (DTI, MAP), and diffusion metric. In each image, the dashed lines represent the distance of the metric calculated on Rician magnitude data from the corresponding Gaussian one, whereas the solid lines report the distance from Gaussianity for metrics calculated on phasecorrected data, which varies with X (abscissa). Color codes indicate the SNR_{0} value. We observe that phase correction leads to metric distributions that are closer to the Gaussianity (H distance close to 0) than the Rician magnitude ones, for specific ranges of X. In general, phase correction debiases the metric distributions up to a great extent. The improvement over the Rician magnitude is clearly correlated with the combination of acquisition scheme—especially the maximum bvalue—and SNR_{0} as also indicated by the signal intensities experiments illustrated in Figs. 1 and 2. Indeed, at high bvalues the signal is low—especially along the less restricted diffusion direction—which, in combination with a poor SNR_{0}, causes the effective SNR to fall well below 5 where a Rician distribution diverges from a Gaussian one. Therefore, the best X value (highlighted with a dot in Fig. 5) also depends on these factors. We point out that in some cases, as for DTI at b = 1000 s/mm^{2}, too much filtering (small X) causes the phasecorrected metric distributions to be more distant from Gaussianity compared to those based on the Rician magnitude (dashed lines). The best X also seems to have a dependence on the considered metric. For instance, at SNR_{0} = 30 the best X for RTPP is 0.75 whereas for RTAP is 2. This can be associated to the fact that metrics that are highly related to signal measured along the less restricted diffusion direction (low intensity signal), such as PD and RTPP, benefit more than others of phase correction. See Table 1 for a comprehensive summary.
Fig. 5 Distance from Gaussianity [see Eq. (6)] for DTI and qspace metrics calculated on Rician magnitude data (“Mg”, dashed lines) and on phasecorrected real data (“Re”, solid lines) as function of A. The acquisition scheme changes for DTI in the first four rows. Each image reports colorencoded lines for different values of SNR_{0} €{10, 20, 30g. The minimum distance for each solid line is highlighted by a dot. A lower value signifies more closeness to Gaussianity. *: s/mm^{2}
Table 1 Maximum relative reduction [0, 1] in Я distance after phase correction compared to Rician magnitude (bias reduction)
SNR_{0} 
MD 
PD 
FA 
RTOP 
RTAP 
RTPP 
MSD 
QIV 
(lK,2K,3K,ms) 
(lK,2K,3K,ms) 
(lK,2K,3K,ms) 
ms 
ms 
ms 
ms 
ms 

10 
(0.75,0.84,0.85,0.85) 
(0.76,0.84,0.85,0.84) 
(0.45,0.78,0.86,0.78) 
0.85 
0.68 
0.86 
0.78 
0.84 
20 
(0.62,0.80.0.84.0.81) 
(0.63,0.79,0.86.0.82) 
(0.24,0.72,0.86.0.70) 
0.75 
0.45 
0.85 
0.55 
0.16 
30 
(0.39,0.78,0.80,0.80) 
(0.52,0.80,0.79.0.79) 
(0.13,0.61,0.75.0.62) 
0.68 
0.36 
0.82 
0.56 
0.11 
Values are reported for each acquisition type—b = 1000, 2000, 3000 s/mm^{2} (1K,2K,3K), and multishell (ms)—and SNRq