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 b-value 1000, 2000, 3000 s/mm^{2} and DTI and MAP multi-shell.

In the first experiment, we clustered a HCP dataset to obtain typical signal values at b-value 1000 and 3000s/mm^{2}. Particularly, for each b-value we applied k-means 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 ground-truth 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 b-value, the noisy magnitude |DWI|_{xy} and the estimated phase-corrected 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 phase-corrected real signals for each circle. We conclude that in both cases the phase-corrected 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 b-values. 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 ground-truth 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 non-diffusion-weighted 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 ground-truth magnitude DWI at any b-value 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 phase-corrected 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 b-value, 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 phase-corrected 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 phase-corrected <(DWE(y), and one for the Gaussian |DWI|!y The hypothesis is that, for each pixel, the phase-corrected 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 phase-corrected real part (second row), at different b-values (columns). Contrarily to the case of the Rician magnitude, the distance from Gaussianity remains visually unchanged as the diffusion-weighting 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 2-4, Fig. 2 shows the maps of Hellinger distance from the Gaussian magnitude, for the Rician magnitude (first row) and the phase-corrected real part (second row), at b-value 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 b-values (from left to right) the signal intensity is lower and the bias occurs in a larger number of pixels. Conversely, the phase-corrected 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 ground-truth magnitude image obtained with Phantomas [15] and a synthetic phase. On the right: the phase-corrected 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 phase-corrected 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 single-shell scheme at b-value 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 q-space 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 q-space inverse variance (QIV) [17]. We create a mask of voxels within fibers, based on the noise-free 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 q-space 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 phase-corrected 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 phase-corrected 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 phase-corrected (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 b-value, 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 phase-correction since the estimated low-frequency phase is very similar to the original noisy one, /DWI„, &/DWIry in Eq. (2). Indeed, the blue histograms (phase-corrected) 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 phase-corrected 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 phase-corrected 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 phase-corrected 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 b-value—and SNR_{0} as also indicated by the signal intensities experiments illustrated in Figs. 1 and 2. Indeed, at high b-values 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 phase-corrected 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 q-space metrics calculated on Rician magnitude data (“Mg”, dashed lines) and on phase-corrected real data (“Re”, solid lines) as function of A. The acquisition scheme changes for DTI in the first four rows. Each image reports color-encoded 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 multi-shell (ms)—and SNRq