Noise-Corrected, Exponentially Weighted, Diffusion-Weighted MRI (niceDWI) Improves Image Signal Uniformity in Whole-Body Imaging of Metastatic Prostate Cancer

Purpose: To characterize the voxel-wise uncertainties of Apparent Diffusion Coefficient (ADC) estimation from whole-body diffusion-weighted imaging (WBDWI). This enables the calculation of a new parametric map based on estimates of ADC and ADC uncertainty to improve WBDWI imaging standardization and interpretation: NoIse-Corrected Exponentially-weighted diffusion-weighted MRI (niceDWI). Methods: Three approaches to the joint modeling of voxel-wise ADC and ADC uncertainty (σADC) are evaluated: (i) direct weighted least squares (DWLS), (ii) iterative linear-weighted least-squares (IWLS), and (iii) smoothed IWLS (SIWLS). The statistical properties of these approaches in terms of ADC/σADC accuracy and precision is compared using Monte Carlo simulations. Our proposed post-processing methodology (niceDWI) is evaluated using an ice-water phantom, by comparing the contrast-to-noise ratio (CNR) with conventional exponentially-weighted DWI. We present the clinical feasibility of niceDWI in a pilot cohort of 16 patients with metastatic prostate cancer. Results: The statistical properties of ADC and σADC conformed closely to the theoretical predictions for DWLS, IWLS, and SIWLS fitting routines (a minor bias in parameter estimation is observed with DWLS). Ice-water phantom experiments demonstrated that a range of CNR could be generated using the niceDWI approach, and could improve CNR compared to conventional methods. We successfully implemented the niceDWI technique in our patient cohort, which visually improved the in-plane bias field compared with conventional WBDWI. Conclusions: Measurement of the statistical uncertainty in ADC estimation provides a practical way to standardize WBDWI across different scanners, by providing quantitative image signals that improve its reliability. Our proposed method can overcome inter-scanner and intra-scanner WBDWI signal variations that can confound image interpretation.


INTRODUCTION
Whole-body diffusion-weighted MRI (WBDWI) is fast gaining acceptance as a powerful tool for diagnosing, staging, and assessing the response of myeloma (1), lymphoma (2), and metastatic prostate (3,4) and breast (5) cancers to systemic treatments. WBDWI provides high contrast between disease and healthy tissue, without the need for exogenous contrast agent injection. The technique has been shown to be particularly helpful in defining the extent of metastatic bone disease, allowing quick assessment of cancer spread "at a glance." By acquiring images with at least two diffusion weightings (b-values), WBDWI also enables calculation of the apparent diffusion coefficient (ADC) of tissues, a putative response imaging biomarker that reflects tumor cellularity (6); ADC increase is observed in effective treatments and early patient response to therapies (7). However, standardization of this technique is complicated due to the number of parameters that require optimization for each scanner (8).
Computed DWI (cDWI) (9) is a post-processing MRI technique that combines voxel-wise estimates of the MR signal at b = 0 s/mm 2 (S0) and ADC to synthetically generate higher bvalue images (assuming a monoexponential relationship between image signal and ADC). Synthesized b-values are typically higher than those that can be directly acquired on MRI scanners due to constraints in image signal-to-noise ratio (SNR) and/or measurement time; such images help to overcome the "T2 shinethrough effect" that can lead to misinterpretation of metastatic disease (10), and have improved SNR compared with directly acquired high b-value images. By optimizing the signal from diseased tissues and minimizing the signal from normal tissues, cDWI has been applied as a pre-processing step in semiautomatic segmentation of disease in WBDWI studies, providing quantification of whole-body volumetric tumor burden (tDV) and global disease ADC (gADC) as biomarkers for treatment response (11,12).
Estimated cDWI images are affected by the T1-relaxivity, T2-relaxivity, and proton density of the imaged tissue, such that the resulting contrast on cDWI may not wholly depend on differences in ADC. Moreover, image contrast is heavily influenced by coil sensitivity, making images susceptible to (i) signal inhomogeneities across the imaging field of view (a "bias field"), and (ii) non-uniform signal between anatomical acquisition stations on WBDWI (3). Such signal inhomogeneities hinder the development of reproducible disease segmentation techniques in WBDWI. A previous attempt to provide images with pure diffusion-weighted contrast was exponentially weighted DWI (eDWI) (13). Using eDWI, S0 is set to a constant value across the entire field of view so that the synthetic higher b-value images are obtained using only ADC generated contrast. However, this technique is suboptimal due to inherently low image contrast-to-noise ratio (CNR).
In this paper, we propose NoIse-Corrected, Exponentiallyweighted DWI (niceDWI) as a technique to improve the limitations of previously described methods. By combining voxel-wise measurements of ADC with estimates of its statistical uncertainty, σ ADC (calculated through linear weighted least-squares fitting), it is possible to apply a noise reduction weighting to conventional eDWI images and thus improve the CNR of these images. Furthermore, we hypothesize that niceDWI is less susceptible to the bias fields and inter-station signal inhomogeneities that are observed in conventional cDWI. This technique can potentially normalize the signal within and between WBMRI acquisitions, thus facilitating comparison of longitudinal WBMRI studies within and between institutions. This may lead to improved visual assessment of disease and reader confidence when using WBDWI. We also derive and assess the statistical properties of three methods for calculating σ ADC . The first of these, Direct Weighted Linear Least-Squares Estimation (DWLS), requires the acquisition of 3 or more images at each b-value measured (≥ 2 unique b-values), whilst the second method, Iterative Weighted Linear Least-Squares Estimation (IWLS), places the slightly weaker restriction that the total number of (not necessarily unique) b-values measured is 3 or more. Our third methodology, Smoothed Iterative Weighted Linear Least-Squares Estimation (SIWLS), improves precision of σ ADC from IWLS in real-world WBDWI protocols where it is typical to acquire only 3 independent b-values (the lower bound for IWLS).

THEORY
The concept of niceDWI is illustrated in Figure 1. Voxelwise estimates of ADC and the uncertainty in its estimation (σ ADC -characterizing the standard error of the ADC value) are calculated using weighted linear least-squares regression (WLS) to produce maps of both parameters. We combine the resultant quantitative maps to generate a novel computed image in postprocessing, S nc , according to the model: where a c ∈ R + and b c ∈ R + can be reduced/increased in realtime to provide weaker/stronger σ ADC and ADC weighting in the computed image, respectively. The resulting signal intensity remains purely quantitative in contrast to conventional high bvalue DWI where image intensity is influenced a multitude of factors including T1/T2-weighting, RF-receive coil sensitivity, and proton density. Images may be displayed as a volumetric 3D maximum intensity projection (MIP) that enables visualization of bone and soft tissue disease "at-a-glance, " where hot-spots reveal regions with (i) low ADC and (ii) high precision of the ADC estimate. In the remainder of this section we describe our approach to joint estimation of ADC and σ ADC from noisy DWI data, and then derive models for the expected image noise in calculated σ ADC maps.

Estimation of ADC and σ ADC
The conventional monoexponential model for DWI is given by where S 0 represents the signal-intensity in the absence of any diffusion-weighting (b = 0 s/mm 2 ), and ν i ∼ N (0, σ ν ) is It should be noted that we display axial diffusion weighted imaged images (S nc ) using a grayscale colormap, where regions of low ADC and low σ ADC are displayed as bright white regions against a dark background, whereas volume-rendered MIPs are displayed using an inverse grayscale where regions of low ADC and low σ ADC appear dark against a white background. This is common practice in whole-body DWI applications.
additive homoscedastic noise sampled from a univariate, zeromean normal distribution with standard deviation σ ν . ADC estimation is commonly performed through linear least-squares regression (LLS) following linearization of this function by a log-transform: The magnitude of imaging noise following log-transformation becomes dependent on the acquisition b-value (b i ), ADC, and S 0 according to (by error propagation): Such heteroscedastic noise warrants the use of a weighted linearleast squares optimization approach (WLS) for ADC estimation, as has been previously explored (14). WLS is also convenient for voxel-wise estimation of the ADC uncertainty, σ ADC . In this article we explore three potential strategies for WLS fitting of these parameters: Direct Weighted linear Least-Squares (DWLS), Iterative Weighted linear Least-Squares (IWLS), and smoothed IWLS (SIWLS). All algorithms were written in Python using the NumPy Einstein summation convention routines, and are provided in the Supplementary Data.

Direct Weighted Linear Least-Squares Estimation (DWLS)
In most conventional WBDWI studies, multiple acquisitions are performed at each b-value, and then averaged to increase the signal-to-noise ratio (SNR) in the final trace-weighted image. Vendors typically return only these averaged images to reduce the large storage requirements for retaining individual image excitations. However, individual acquisition images provide direct calculation of weights for use in WLS estimation of ADC and σ ADC , as demonstrated in Figure 1. Consider that M repeat excitations are acquired for each of N different b-values. A weight for each b-value b i at each pixel location can be derived as the inverse variance of the M data acquired at b i : where y i represents the average log-signal at a particular voxel location over the repeat acquisitions. DWLS proceeds by converting these weights into a diagonal matrix, W, and combining this with matrices representing the known b-values, B, and measured log-signals, Y: (where w ij = w ik ). We then calculate the vector of parameter estimates α = ADC, ln S 0 ⊺ , and the covariance matrix for these parameters, α , by: ADC and ADC uncertainty estimation are then derived from the first element of the parameter and parameter covariance matrices, respectively:

Iterative Weighted Linear Least-Squares Estimation (IWLS)
Our second approach is an iterative solution for joint estimation of α and weight matrix W, with subsequent calculation of α . This method provides the advantage that weights do not need to be calculated a-priori, nor are repeat measurements required at each b-value. Using the same matrix notation outlined in Equation (6), the following algorithm is performed for each voxel at spatial location (x, y) within the image: Return: α t , α where N t is the maximum number of iterations, and ǫ is a convergence tolerance (in this article we set default values of N t = 100 and ε = 10 −5 ). The ADC and its uncertainty can be calculated from the resulting estimates of α t and α using Equation (8) accordingly.

Smoothed Iterative Weighted Linear Least-Squares Estimation (SIWLS)
Calculation of σ 2 ν and α using IWLS is only possible if NM ≥ 3 and N ≥ 2. However, for low numbers of M and N, estimates of σ ADC become especially noisy (see section 2.2). We therefore propose an optional modification to the above algorithm where we replace the estimate of the weighted data variance, σ 2 ν , with a smoothed version over the entire spatial field before being used to calculate α : where ⊗ is the convolution operation performed over spatial field (x, y), and ρ is some smoothing kernel with width chosen by the user (with ∞ −∞ ρ(x, y, )dxdy = 1). In this article we use a box kernel, ρ(x, y, ) = rect(x/ , y/ )/ 2 , due to an advantage in the statistical properties of this filter, but in principle any linear or non-linear smoothing could be used to perform this additional step. Furthermore, if adequate DWI spatial resolution is acquired, smoothing may be performed over the entire 3-dimensional spatial field (x, y, z). Adequate choice of is considered later in this article.

Noise Properties of σ ADC Maps
Maps of σ ADC generated by the methods proposed are affected by imaging noise, and it is of interest to understand the expected distribution for σ ADC given true underlying values for ADC, S 0 and σ 2 ν . Under the assumption of Gaussian distributed noise for signal intensities, the estimated data variance, σ 2 ν , will be chi-square distributed: we conclude that the ADC uncertainty is now chi-distributed (after appropriate scaling) with k degrees of freedom: where A 11 (first element of matrix A) is a function of the true values of ADC, S 0 , and the acquisition b-values. Given the expectation and variance of a chi-distribution are and Var ( z) = k − E ( z) 2 , respectively, and that E( z 2 ) = k and Var( z 2 ) = 2k ( z 2 is chi-square distributed), we have the results: An important result occurs from these formulae; consider an experiment in which M acquisitions are acquired at N different b-values, assumed to have data variance σ 2 ν , and weighting value A 11 . Then consider a second dataset for which only the average over all M acquisitions is provided for each of the N b-values. By the central limit theorem the data variance would be σ 2 ν(M) = σ 2 ν /M (the number in parentheses represent the number of signal averages), which when paired with the fact that A (M)11 = M · A (1)11 results in the relationships: If we inspect the case of N = 3 and M = 9 (as per the patient the examples explored in this article), then we would expect to see an 18-fold improvement in the SNR of σ ADC maps when keeping individually acquired images over using averaged data only: These equations also demonstrate that an unbiased estimator for σ ADC is given by: where α11 is derived from the DWLS or IWLS algorithms, and thus the estimated variance is modified to For the SIWLS modification, it is important to consider that the effective value for k will change following the smoothing process, resulting in a change in the bias correction factor c(k). For the box kernel used in this article we have k → k × 2 (i.e., for a box kernel of width = 3, the effective k would be increase by a factor of 9).

Monte Carlo Simulations
Monte Carlo simulations were performed to confirm the accuracy of the assumptions made in section 2 over a range of "true" SNR and ADC values (100 SNR values and 100 ADC values in the ranges 1 → 300 and 0 → 4 × 10 −3 mm 2 /s, respectively). For each combination of SNR and ADC we sampled nine noisy signals for each b-value in the set b i ∈ {50, 600, 900} from a Rician distribution with parameters µ = exp{−b i · ADC} and σ = 1/SNR (15). Using these simulated data we performed DWLS and IWLS fitting routines to obtain (unbiased) estimates ADC and σ ′ ADC . Furthermore, we performed IWLS fitting for data consisting of the arithmetic average of the nine values simulated at each b-value (denoted IWLS (9) ). Simulations were performed N s = 10 5 times for each combination of SNR and ADC. For each fitting method, at each SNR/ADC combination, we compared estimates ADC, σ ′ ADC , and Var( σ ′ ADC ) against simulated values for these parameters (ground truth) using the formulae: where A 11 can be calculated directly from the known ADC values, and k = 25 for IWLS and k = 3 for IWLS (9) . Simulations were performed on a 3.5 GHz, 16 GB personal computer running Python.

Phantom Study
In order to demonstrate the contrast-no-noise (CNR) characteristics of niceDWI, we performed a phantom study using an in-house developed test-object. Our phantom consists of five cylinders containing solutions of different proportions of water, manganese chloride and sucrose, in order to generate environments with similar T2 relaxivity and ADC ranges to those observed in tumor tissues (75-1,408 ms and 0.7-1.1 ×10 −3 mm 2 /s, respectively); cylinders were bathed in ice-water, which was then allowed to equilibrate at room temperature for ∼1 h in order to achieve an expected temperature of 0 • C within each vial (16). Diffusion-weighted images were acquired on a 1.5T system (Siemens Aera, Erlangen, Germany) with the following acquisition parameters: echo time TE = 93 ms, repetition time TR = 6 s, in-plane parallel imaging factor R = 2 (GRAPPA), image size = 128 × 128, pixel spacing 1.91 × 1.91 mm 2 , slice thickness 5 mm, pixel bandwidth 1,955 Hz/pixel, b-values = 50/600/900 s/mm 2 , three orthogonal diffusion encoding directions, number of signal averages NSA = 1. Imaging was repeated three times, and images acquired for each of the three orthogonal diffusion encoding directions were considered independent (under the assumption of isotropic diffusion), providing a total of nine images for each b-value. These data were used to calculate maps of ADC and σ ADC using the IWLS algorithm, which in turn were used to generate niceDWI images over a range of a c and b c values: a c ∈ (0, 1, 000, . . . , 50, 000) and b c ∈ (0, 100, . . . , 5, 000). This experiment was repeated twice, so that the noise statistics of niceDWI could be calculated as the standard deviation of the difference in voxel intensities from both experiments within each vial for each combination of a c and b c (scaled by a factor of 1 √ 2 ). From this, the CNR was calculated between cylinders 2-5 and cylinder 1 as the ratio of the difference between the average of niceDWI voxel intensities at each a c /b c combination, and the calculated standard deviation of signal noise.

Patient Study
We performed whole-body niceDWI experiments in a pilot population of 16 patients with suspect metastatic bone disease from primary prostate cancer. Axial images were acquired from skull base to mid-thigh over 4-5 imaging stations, with each station consisting of 40 slices. Images were acquired on a 1.5T scanner (Siemens Aera, Erlangen, Germany) with the following acquisition parameters: echo time TE = 79 ms, repetition time TR = 12.7 s, in-plane parallel imaging factor R = 2 (GRAPPA), image size = 128 × 104 (interpolated to 256 × 208), pixel spacing 3.36×3.36 mm 2 (interpolated to 1.68×1.68 mm 2 ), slice thickness 5 mm, pixel bandwidth 1,955 Hz/pixel, b-values = 50/600/900 s/mm 2 , three orthogonal diffusion encoding directions, number of signal averages NSA = 1. For each imaging station, image acquisition was repeated three times such that a total of nine images were available for each b-value at each slice location (three orthogonal diffusion encoding directions × three repeat acquisitions; we make the assumption of isotropic diffusion for metastatic disease).
We synthesized conventionally acquired data by determining a trace weighted image for each acquisition station as the geometric mean of signal intensities over all diffusion-encoding directions, followed by the arithmetic mean of the three trace-weighted datasets acquired at each b-value (equivalent to NSA = 3). For these data we calculated maps of σ ADC using the SIWLS algorithm using different smoothing widths ∈ (0, 1, . . . , 50) for a box smoothing function (Equation 9). We then calculated the root-mean-square-error (RMSE) between the σ ADC maps calculated using SIWLS and σ ADC maps calculated using conventional IWLS when retaining all nine independent acquisitions for each b-value (gold standard). The RMSE was calculated for each axial image within all patient datasets, resulting in a total of 2,680 measurements from which statistics could be calculated. Furthermore, metastatic bone disease was delineated in all patients by a dedicated radiologist using an in-house developed semi-automated segmentation pipeline (11). The RMSE was also calculated within regions-ofinterest for each axial image that contained disease (resulting in 1,122 measurements). Figure 2 illustrates the results from the simulation study. All algorithms are able to accurately measure ADC values in the regime SNR min > 5 (SNR min = minimum SNR across all bvalues) with negligible bias. Where SNR min < 5, ADC estimates demonstrate a negative bias, which is likely due to the effect of Rician noise in magnitude imaging data, as explored in previous studies (14) (IWLS appears to be slightly more robust to this effect). For ADC uncertainty measurements, σ ADC , the DWLS algorithm demonstrates a negative bias over all SNR and ADC values explored. Conversely, the IWLS and IWLS (9) algorithms FIGURE 2 | Each row demonstrates simulation results for ADC, σ ADC and Var ( σ ADC ) for each of the fitting methods explored (columns). In each case, red surface plots demonstrate the estimated value derived from the equations presented in section 2, whilst the surface mesh represents the ground truth derived from simulations. The difference between estimated values and ground truth (Estimate-Truth) is also shown as a color figure on the right of each surface plot. Results for σ ADC and Var [ σ ADC ] are log-transformed for visual clarity. The IWLS scheme results in the lowest bias across all statistical estimators, especially in the region where the minimum SNR over all b-values is >5 (bottom-right region from dashed curve on difference images). In regions where the minimum SNR is <5, the assumption of Gaussian noise no longer holds, and results in observed differences between estimated parameters and ground truth. This is lowest for IWLS (9) , where averaging of data at each b-value improves the effective SNR. However, the statistical noise present in σ ADC estimates is demonstrably poorer for IWLS (9) , as indicated by the generally higher Var ( σ ADC ) across all simulated ADC and SNR values (green arrows). demonstrate no bias in the region of SNR min > 5, confirming that the theoretical justifications provided in section 2 hold over these values. The same trend was observed for the noise properties of σ ADC , as parameterized by Var(σ ADC ). Importantly, simulations agreed with our observation that Var(σ ADC ) is larger when only the averaged data is available as opposed to retaining each individual acquisition at each b-value. Due to the bias observed in the DWLS method, we have excluded this technique from all further analysis.

Phantom Study
The results from our ice-water phantom study are presented in Figure 3. Our results demonstrate that contrast between the different vials can be generated in real-time by varying computed a c and b c values (representing increased σ ADC and ADC weighting, respectively). By allowing users to adjust a c and b c independently, we observe it is possible to have more flexibility in the desired contrast within the image when compared to cDWI, and optimize the CNR between two regions of interest within the image. This provides the ability to increase the CNR between signal within vials and background noise compared to using eDWI alone. Figure 4 illustrates the results of our optimization approach for the smoothing weight used in the SIWLS algorithm, along with exemplar maps of σ ADC for different values of in a single patient image. Results indicate that when considering FIGURE 3 | A comparison of cDWI, eDWI, and niceDWI performed for our ice-water test object over a range of computed a c and b c values. It is clear that a wide range of contrasts can be generated by adjusting the (a c , b c ) pair in real-time, compared with cDWI and eDWI where only b c can be adjusted. The CNR for each outer vial compared to the central vial can improved by increasing a c and/or b c , as evidenced by the CNR surface plots illustrated on the right (color coding for these plots is presented in the bottom-right figure). It is evident that both cDWI and eDWI present with poor signal-to-background contrast (green and red arrows, respectively), and that niceDWI provides a means to alleviate this issue. FIGURE 4 | Maps of σ ADC calculated using the SIWLS algorithm over different box-kernel widths (right), compared to the gold-standard IWLS algorithm (bottom-right) for a single axial image from one of our patient datasets. It is evident that as the kernel width is increased, the noise field within the σ ADC map is reduced and begins to more closely resemble the gold-standard map. However, as is increased even further, a "haloing" effect is observed (red arrow) suggesting that an optimum can be selected. On the left we present plots of RMSE calculated over the entire σ ADC map, and within disease only (an example of a region of interest around disease is illustrated as the red line on the bottom-right map). Over the entire map, a minimum RMSE can be observed at = 9, whilst for the disease regions only, a plateau is observed after approximately = 20.

Patient Studies
the RMSE between estimates of σ ADC from IWLS and SIWLS methods over the entire image, an optimum value for is ∼9 voxels (equivalent to 15.1 mm in this study). However, when considering voxels containing metastatic disease only, an optimum is minimum is reached once > 20. By inspecting the exemplar σ ADC maps we observe that for < 20, maps appear to be affected by background "mottle" effect, which may influence the poorer RMSE within diseased regions. Conversely, as reaches much higher values, σ ADC maps are affected by a "halo" around the patient surface. We therefore select a value of = 20 for all further exploration of the SIWLS algorithm. Figure 5 present results for four of the patient datasets, comparing total-body niceDWI images obtained from IWLS and SIWLS algorithms (a c = 20, 000 s/mm 2 , b c = 900 s/mm 2 ) with conventional cDWI images (b c = 900 s/mm 2 ). Maximum Intensity Projections (MIPs) of these datasets illustrate that images generated with IWLS and SIWLS algorithms produce similar image contrast for evaluating tumor load in these patients. Comparing these images with MIPs from cDWI demonstrates the ability of niceDWI to provide more uniform signal intensities over the entire patient field-of-view. Furthermore, by inspecting maps of σ ADC obtained by IWLS and SIWLS in these cases, it is apparent that the latter method can provide a good surrogate method where only averaged data is provided by the vendor. Regions of discrepancy between these two techniques include (i) areas of soft-tissue motion such as the bowel loop and liver, (ii) peripheral fat, and (iii) the spinal cord.
An overview of niceDWI volumes (calculated using IWLS and SIWLS) for all patients is presented in Figure 6. Comparing total-body MIPs, it is evident that niceDWI provides more uniform contrast for visualization of metastatic bone disease in these patients; good agreement is observed between IWLS and SIWLS estimations methods. Some differences in contrast may be seen between cDWI and niceDWI, which is to be expected given the different properties that give rise to the derived voxel intensities, but both provide good visualization of the extent of metastatic disease.
Calculation of σ ADC maps for entire WBDWI datasets using IWLS/SIWLS took ∼12 s on a personal computer with a 3.5 GHz processor running Python; vast increases in speed could FIGURE 5 | Exemplar results for four of the patients explored in this study. ADC maps for two slice locations are presented alongside σ ADC maps calculated directly through the IWLS scheme, or following signal averaging with SIWLS ( = 20). In addition, coronal total-body MIPs of niceDWI datasets are compared alongside conventional cDWI datasets (bottom row in each case); within MIPs dark regions represent voxels with low ADC and low σ ADC . From the σ ADC maps it is clear that the SIWLS algorithm is able to accurately measure estimates of ADC uncertainty within regions of disease (green arrows-zoomed regions also displayed) when (Continued) FIGURE 5 | compared with the IWLS algorithm (gold-standard). However, SIWLS is less able to recreate the appearance of σ ADC in regions of peripheral fat (yellow arrows), and may be less accurate in regions where motion occurs including the bowel loop (red arrows), or close to tissue boundaries (orange arrow). By comparing the MIPs, it is shown that niceDWI can remove some of the intensity non-uniformities observed between different imaging stations on cDWI (blue arrows). Furthermore, niceDWI may be able to enhance the visualization of bone disease (purple arrows). be expected once the algorithm has been optimized using parallel processing.

DISCUSSION
In this article we have developed niceDWI, a post-processing approach for generating a new contrast mechanism in WBDWI studies that combines estimates of ADC with voxel-wise measurement of the uncertainty in these ADC estimates. By producing computed images in this way, we generate purely quantitative signal intensities, where bright regions (or dark on volume rendered MIPs) reflect areas of low ADC and low uncertainty in these ADC measurements. Generating purely quantitative contrast through niceDWI could offer improved visualization of treatment effects in longitudinal studies; by mitigating the effects of coil sensitivity, T1/T2-weighting and variable gain settings between scanners, this technique could enhance comparability between diffusion-weighted imaging studies acquired at different institutions. This work would be further strengthened through multi-center studies, using data acquired from multiple MR-vendors to demonstrate how well the niceDWI approach can standardize signal intensity across different patients and scanners. The diagnostic sensitivity and specificity of niceDWI could be assessed in these larger cohort studies. We have evaluated three weighted least-squares fitting approaches (DWLS, IWLS and SIWLS) for joint estimation ADC and σ ADC from WBDWI data, and established a statistical framework for exploring noise properties of σ ADC maps derived from these approaches. Using this framework, we demonstrated a clear statistical advantage for retaining individual image acquisitions from WBDWI studies due to the considerable increase in the resulting SNR of calculated σ ADC maps (though negligible difference is observed in ADC estimation). Although we have set up our scanners to acquire individual image acquisitions for this study, such individual acquisitions are unfortunately not conventionally retained by scanner vendors, perhaps due to data storage limitations. We have therefore also suggested a minor alteration to our IWLS algorithm to provide a potential solution when only averaged data are available (SIWLS), which could improve uptake of niceDWI in the clinical setting. As this latter methodology can be easily adapted to most clinical WBDWI protocols, it would also facilitate retrospective evaluation of niceDWI.
Through Monte Carlo simulations, we have (i) demonstrated the validity of the assumptions made in our statistical derivations, and (ii) evaluated the performance of our DWLS and IWLS model fitting approaches. Results from our study indicate that the assumptions hold in the limit of minimum SNR >5 across all b-values; for our patient cohort the median SNR [calculated as S(900)/σ ν ] within disease voxels on b = 900 s/mm 2 images was 19.8 (5th-95th percentile range: 7.4-39.8), which is well above this limit. Furthermore, it is evident that DWLS suffers from a bias in the estimation of σ ADC , whilst the IWLS approach remains unbiased. We therefore conclude that IWLS should be used for future studies when estimating σ ADC . This has the minor drawback that this approach requires multiple iterations and can thus take longer to compute, although in our patient studies we found that 97.5% of the 143 × 10 6 voxels investigated had converged after just five iterations.
Our methodology could easily be adopted in clinical trials that utilize DWI for screening and assessing treatment response of other disease types, including in multiple myeloma (17), lymphoma, mesothelioma, and breast cancer metastasis. In future studies, repeatability measurements would be valuable to quantify what changes in niceDWI signal intensity could be considered significant following therapy; measurement of σ ADC also holds promising potential for providing repeatability of ADC measurements within individual patients, although this would need to be further explored. Whilst our methodology relies on manual selection of computed b-values to determine the ADC weighting within niceDWI images (the choice of these parameters can be highly operator dependant), approaches such as the one developed by Gatidis et al. (18) for automatic bvalue selection in conventional cDWI could be extended to niceDWI. Furthermore, due to low number of b-values typically acquired with WBDWI, we have focussed on the use of a monoexponential diffusion model in our experiments. Our approach could be extended to investigate the use of voxel-wise uncertainty estimation for parameters obtained from Intra-Voxel Incoherent Motion (IVIM) (19)(20)(21), kurtosis (22), and stretched exponential (23) diffusion imaging models.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Royal Marsden Hospital NHS Foundation Trust. The ethics committee waived the requirement of written informed consent for participation.