Image denoising and model-independent parameterization for IVIM MRI

Objective. To improve intravoxel incoherent motion imaging (IVIM) magnetic resonance Imaging quality using a new image denoising technique and model-independent parameterization of the signal versus b-value curve. Approach. IVIM images were acquired for 13 head-and-neck patients prior to radiotherapy. Post-radiotherapy scans were also acquired for five of these patients. Images were denoised prior to parameter fitting using neural blind deconvolution, a method of solving the ill-posed mathematical problem of blind deconvolution using neural networks. The signal decay curve was then quantified in terms of several area under the curve (AUC) parameters. Improvements in image quality were assessed using blind image quality metrics, total variation (TV), and the correlations between parameter changes in parotid glands with radiotherapy dose levels. The validity of blur kernel predictions was assessed by the testing the method's ability to recover artificial ‘pseudokernels’. AUC parameters were compared with monoexponential, biexponential, and triexponential model parameters in terms of their correlations with dose, contrast-to-noise (CNR) around parotid glands, and relative importance via principal component analysis. Main results. Image denoising improved blind image quality metrics, smoothed the signal versus b-value curve, and strengthened correlations between IVIM parameters and dose levels. Image TV was reduced and parameter CNRs generally increased following denoising. AUC parameters were more correlated with dose and had higher relative importance than exponential model parameters. Significance. IVIM parameters have high variability in the literature and perfusion-related parameters are difficult to interpret. Describing the signal versus b-value curve with model-independent parameters like the AUC and preprocessing images with denoising techniques could potentially benefit IVIM image parameterization in terms of reproducibility and functional utility.


Introduction
Intravoxel incoherent motion (IVIM) magnetic resonance imaging (MRI) is a diffusion imaging technique that quantifies the 'translational motion of molecules within imaging voxels' (Le Bihan 2019).This motion is primarily attributed to the random diffusion of water molecules within tissue as well as the randomly oriented microcirculation of blood in capillary beds.Motion is inferred by analyzing the signal in each voxel as a function of diffusion gradient strengths (b-values) used in the pulse sequence.Le Bihan first introduced a method of quantifying these in vivo motions using (Stejskal and Tanner 1965) diffusion gradient method in 1986 (Bihan et al 1986).This paper introduced the concept of an 'apparent diffusion coefficient' (ADC).The ADC quantifies the combined effects of diffusion and perfusion.By acquiring images with multiple b-values, the ADC can be fit using the equation This equation is analogous to Stejkal and Tanner's in vitro diffusion equation for a still liquid, S(b)/S(0) = e − D• b (Stejskal and Tanner 1965), where D is the self-diffusion coefficient of the liquid.The ADC is a simple metric for describing the signal decay curve in vivo, which in reality is a consequence of complicated biological mechanisms.While Stejskal and Tanner's diffusion equation was derived from the Bloch equations (Bloch 1946), Le Bihan's extension of this equation to in vivo circumstances is an approximation.It is theorised that diffusion and perfusion effects can be disentangled by fitting appropriate models to the signal decay curve, as originally posed by Le Bihan et al (1988).This is done by modelling the collective perfusion of blood through microcapillaries as a 'pseudodiffusion' process, with a pseudodiffusion coefficient, D * .D * is presumably on the order of 10 times D (Le Bihan 2019).As a consequence the signal decay curve is largely due to perfusion effects at low b-values, and is dominated by diffusion effects at high b-values (Zhou et al 2016).A common approach for separating diffusion and pseudodiffusion in practice is to apply Le Bihan's biexponential model to the decay curve (Le Bihan et al 1988) where f denotes the fraction of the signal decay attributed to perfusion.
In recent years, the literature has trended towards adopting a triexponential model (Cercueil et al 2014, Kuai et al 2017, Wurnig et al 2017, Chevallier et al 2019, Riexinger et al 2019, 2020) of the form 2. Methods

Dataset
This study was approved by an institutional review board, and written, informed consent was obtained from all patients.The inclusion criteria for recruitment were: (1) a diagnosis of nasopharynx, base-of-tongue, or tonsil cancer to be treated with external beam radiotherapy; (2) expected to receive dose in parotid glands during treatment; (3) scheduled to have MR images acquired prior to radiotherapy and at 3 months following the last treatment day.IVIM MR images encompassing the parotid glands were acquired during routine clinical appointments for 12 head-and-neck cancer patients prior to radiotherapy (Age: 35-78 Y, average: 60.4 ; sex: 11M, 1F).Of these, 5 had follow-up images acquired at 3 months post-RT (Age: 49-78 Y, average: 59.3 ; sex: 4M, 1F). 4 of these 5 patients received weekly chemotherapy concurrent with radiotherapy (three patients: weekly 40 mg m −2 cisplatin, one patient: weekly 20-40 mg m −2 per week and 240 mg carboplatin).
Images were acquired on a 1.5 T Magnetom Sola scanner (Siemens Healthineers) with an echo planar imaging sequence (EPI).Diffusion weighted images were acquired with a voxel size of 1.6 × 1.6 × (3.4-3.9)mm 3 .This slight inter-patient variation in slice thickness was imposed to encompass parotid glands within 32 slices while maximizing the signal-to-noise ratio (SNR).The repetition time (TR) and echo time (TE) were 4900 and 105 ms, with a flip angle of 90°.GRAPPA acceleration was used and the entire scan time was approximately 15 min.The EPI echo train length was 65, and multi-band excitation was used.Images were acquired for 16 b-values (0, 20, 30, 40, 50, 60, 70, 80, 90, 100, 120, 150, 250, 400, 800, 1000 s mm −2 ).These values were chosen to follow recommendations for effectively separating perfusion and diffusion parameters (Cohen et al 2014, Ye et al 2019), typically with a partition margin around 50 s mm −2 (Raju 2021).Signals were acquired in 6 directions, spanning the three principal Cartesian axes and diagonal vectors.Two signal averages were computed in each direction.Signals from various directions were reduced to a single quantity by taking the geometric mean (Mukherjee et al 2008).
Parotid glands were manually contoured by a graduate physics student on clinical T1 images.These images were acquired with a turbo spin-echo (TSE) sequence (TR: 564 ms, TE: 20 ms, flip angle: 148°).T1 images were acquired with a voxel size of 0.6 × 0.6 × 3 mm.Contours were used to calculate IVIM parameter statistics within parotid glands.DICOM radiotherapy dose and planning structure set files used for external beam radiotherapy treatment planning were exported from the ARIA Oncology Information System for computing dose statistics inside parotid glands.Whole-mean parotid gland dose levels were calculated within computed tomography (CT) contours, manually defined by a single, senior clinical radiation oncologist.

Denoising
Images were denoised using neural blind deconvolution (Ren et al 2020).Blind deconvolution is the ill-posed mathematical problem of estimating hypothetical denoised images, x, and their associated blur kernel, k, which convolve to yield the original image, y = x * k.This method is particularly useful in situations where spatial resolution and precision of signal localization is low, hence many studies have attempted to employ it for positron emission tomography (PET) (Alpert et al 2006, Lee et al 2008, Šroubek et al 2009, Erlandsson et al 2012, Guérit et al 2015).The use of a single kernel for entire images is an approximation, as it has been found (Levin et al 2009) that the shift-invariant assumption for the blur kernel is often incorrect.Neural blind deconvolution is a method of solving the blind deconvolution problem using neural networks which are simultaneously optimized to predict x and k.It was originally implemented by Ren et al (2020) in 2020 and applied to images from a two-dimensional natural image database.It was significantly improved upon by Kotera et al (2021) in 2021.Neural blind deconvolution is an unsupervised learning methodology which trains and predicts on a caseby-case basis, without the requirement of a separate training set with ground truth images.
Neural blind deconvolution was previously adapted for 3D prostate specific membrane antigen (PSMA) positron emission tomography (PET) in 2023, while incorporating simultaneous super-sampling into the methodology (Sample et al 2023).It was shown to improve blind image quality metrics, and strengthen correlations between PSMA PET uptake and sub-regional importance estimates in the parotid gland for predicting post-radiotherapy xerostomia (subjective dry mouth) (Sample et al 2024).Here we build off this methodology for suitability with IVIM MRI.Networks were adapted to handle four-dimensional input images (b, x, y, z), composed of all 3D images with different b-values stacked.The network architecture was maintained to predict a single kernel, such that a single blur kernel was predicted using all b-value images for each patient.Images were normalized as  , where μ body and σ body are the mean and standard deviation of signal in body voxels.
The formulation of the blind deconvolution problem used in previous neural blind deconvolution studies, y = x * k, accounts for bleeding of voxel signal into neighbouring voxels, but it does not account for background noise.As MRI signals contain rician noise (Gudbjartsson and Patz 1995), and IVIM MRI using EPI is particularly prone to noise, we modified the formulation of the problem to include a noise term: This required the introduction of an additional symmetric convolutional auto-encoder network, G δ , for predicting the noise term δ.The networks used are thus, where θ x , θ δ and θ k represent trainable model parameters of G x , G δ and G k , respectively.The following implicit constraints exist on network outputs which are satisfied automatically The optimization problem can be written in terms of the model parameters θ x , θ δ and θ k as arg min , .7

*
Optimization then proceeds as a type of 'zero-shot' self-supervised learning method (Shocher et al 2018, Ren et al 2020), such that no separate ground-truth image is present, and the models train on the input data that they are ultimately to predict denoised images for.A joint-optimization strategy is used (Ren et al 2020), where θ x , θ k , and θ δ are updated simultaneously in each optimization iteration.Each optimization iteration proceeds by first computing the loss function and its gradient, then backpropagates to compute gradients for all network weights, which are then used to update these network weights according to the given learning rate using the ADAM optimization algorithm (Kingma and Ba 2014).After a set number of iterations, N, the final denoised image, x, blur kernel, k, and noise term, δ, are predicted using ( ) q G x x N , ( ) q G k k N , and ( ) G N .The networks, G x , G k and G δ all require arbitrary input images to be fed through network layers to produce final output predictions.In the original implementation of neural blind deconvolution (Ren et al 2020), input images were sampled from uniform distributions.The full architecture is summarized visually in figure 1.
Optimization was performed using a single NVIDIA GeForce GTX 1060 GPU and a 2.8 GHz Intel Core TM i5-8400 CPU.As neural blind deconvolution is a computationally expensive procedure, further exacerbated by stacking multiple b-value images, images were cropped to extend at least five voxels outside of parotid gland borders in all 3 cartesian directions.This was performed separately for the left and right parotid glands, such that cropped images around each gland were denoised separately.This cropping procedure was necessary to limit GPU memory usage.
A multi-scale optimization procedure was implemented, as recommended for neural blind deconvolution by Kotera et al (2021) and recently implemented for PET denoising (Sample et al 2023(Sample et al , 2024)).This involved pretraining networks using images that were first down-sampled by a factor of two, then up-sampling final outputs by 2 before performing a second round of pre-training.After two pre-training rounds, the outputs, G x , G k , and G δ , are used as initial inputs for the regular training procedure.An 11 × 11 × 11 kernel was used, which was down-sampled to 5 × 5 × 5, and then 7 × 7 × 7 during pre-training stages.The algorithm for updating network weights to predict denoised diffusion images is summarized in table 1 The training procedure employed 5000 iterations.This was split into two stages, with two different loss functions.These stages and their corresponding loss terms are summarized in table 2. It is necessary to include regularization terms in the loss function to avoid convergence to a trivial solution (k becomes a unit impulse).Therefore, the mean squared error (MSE) of kernel values above 0.7 were penalized.We only penalized values above this cutoff to avoid a trivial solution while still allowing small contributions from neighbouring voxels to go unpenalized.We furthermore penalized the MSE of noise voxels above μ δ + σ δ , where μ δ and σ δ are the mean and standard deviation of the normalized original image voxels located outside of the body.
A total variation (TV) loss term for x was included after the 3500th iteration.Including TV loss in later stages of training was shown by Kotera et al (2021) to improve convergence.Kotera et al. introduced a TV loss term after the 2000th iteration due to its tendency to drive optimization towards the trivial solution during initial iterations.Potentially due to the high total variation seen in EPI compared to natural images, we found during initial testing that it was advantageous to introduce TV loss after 3500 iterations.The fidelity loss term switched from MSE to the structural similarity index metric (SSIM) after the 3500th iteration.This was similar to the method recommended by Kotera et al (2021) for improving final image quality.(Ren et al 2020) and implements modifications suggested by Kotera et al (2021).
( ) ( ) 11. Upsample k by 7/5 and upsample spatial dimensions of x and δ by a factor of 2 12. Downsample spatial dimensions of y by 2 to match new convolution size 13.z x = x, z δ = δ, z k = k 14.Repeat steps 4 through 10. 15.Upsample k by 11/7 and upsample spatial dimensions of x and δ by a factor of 2 16.
Compute loss and back-propagate gradients 22. Update G i x , d G i and G i k using the ADAM optimizer (Kingma and Ba 2014) 23.

Evaluating denoised images
As neural blind deconvolution is an unsupervised training procedure, ground truth images were not available to evaluate denoised images.Image quality was first assessed in terms of blind image quality metrics from the PyTorch image quality (PIQ) library (Kastryulin et al 2019).The blind/referenceless image spatial quality evaluator (BRISQUE) (Mittal et al 2012) and contrastive language image pre-training (CLIP) (Wang et al 2022) score were evaluated.BRISQUE ranks image quality on a scale between 0 and 100 based on a series of features derived from various signal intensity and distribution statistics.These features are used to predict the deviation of the input image from a natural, undistorted image.CLIP uses a language/vision neural network trained using millions of natural image/caption pairs to assess either image quality or caption quality.We also determined blind image quality metrics of IVIM parameter maps.The total variation (TV) of diffusion images and parameter maps was also compared before and after denoising.BRISQUE, CLIP, and TV were the only available referenceless quality metrics available with the PIQ library (Kastryulin et al 2019).
To evaluate the model's ability to predict accurate blur kernels, we employed the same strategy used to evaluate PSMA PET image deblurring (Sample et al 2023).This involves generating artificial kernels, convolving them with previously denoised images, then restarting neural blind deconvolution.The inner product between generated and predicted pseudokernels was then assessed.Four types of pseudokernels were generated: a regular Gaussian with a standard deviation of 1 voxel in each direction, and three oblong Gaussians, each having a standard deviation truncated by a factor of three in one of the three Cartesian directions.Pseudokernels were normalized to have unit sums.
We furthermore assessed the similarity between predicted kernels amongst patients.As all patient images were acquired on the same scanner, predicted blur kernels were expected to be similar, albeit different.This is because the blur kernel contains components due to both scanner-related imperfections as well as patient motion and other physiological processes.We compared kernels by computing the inner product.As denoising was performed separately on images cropped around each parotid gland, intra-patient blur kernels were compared as well, which were expected to have higher similarity than inter-patient blur kernels.
Original and denoised images were compared in terms of their practical ability for model-building.The ability of radiotherapy dose levels to predict changes in IVIM parameters following treatment was assessed.Spearman's rank correlation coefficient (r s ) was compared between the whole-mean dose and relative changes in mean IVIM parameters within each parotid gland ( ) , where P i and P f are the parameters before and after radiotherapy).Parameter dose responses were also analysed using a multiple regression analysis.For this, postradiotherapy parameters were dependent variables while pre-radiotherapy parameters and mean dose levels served as two independent variables.The best-fit slope with respect to dose was then compared between parameters.Images were first normalized to statistical Z-values.

Extracting exponential model parameters
To serve as a basis of comparison for newly derived, model-independent IVIM parameters, we extracted ADC, biexponential, and triexponential parameter maps.Parameters were fit in each voxel using the IVIM 3 -NET (Kaandorp et al 2021) methodology.IVIM 3 -NET trains neural networks on a given dataset to optimally predict IVIM parameters.A major benefit of IVIM 3 -NET is that it simultaneously predicts S 0 , the signal at b = 0, which is used to reduce the influence of noisy S 0 values on parameter estimates.
Table 2. Optimization consisted of 5000 iterations, split into 3 stages, each having a different loss function.Up to the 3500th iteration, the fidelity loss function was mean squared error (MSE).The structural similarity index metric (SSIM) was used in remaining iterations.In all iterations, a kernel and noise regularization term were included.The MSE of kernel values were penalized to avoid convergence to the trivial solution.The MSE of noise values above μ δ + σ δ , where μ δ and σ δ are the mean and standard deviation of the normalized original image voxels located outside of the body.Finally, the total variation (TV) between denoised image voxels was penalized after the 3500th iteration, as it has been shown to improve results when employed in later stages of optimization (Kotera et al 2021).

Fidelity term MSE SSIM
Regularization terms 1. Kernel 1. Kernel 2. Noise 2. Noise 3. TV Separate IVIM 3 -NET models were trained to predict ADC, biexponential, and triexponential parameters.The dataset was filtered to exclude voxels outside the body from training.This resulted in approximately 3.6 million training voxels.Voxels for each b-value were normalized to signals acquired with b = 0 s mm −2 .IVIM 3 -NET hyper-parameters were tuned according to previously derived optimal values (Kaandorp et al 2021).Models are fully-connected neural networks with two hidden layers, which take normalized signal versus b-value arrays as input, and output the parameters corresponding to the given model.Each layer consists of a linear function, an exponential linear unit (ELU) activation function (Clevert et al 2015), and batch normalization.Sigmoid functions were used to constrain parameter predictions within boundaries chosen to over-encompass anticipated parameter ranges (table 3).Network weights were then optimized via backpropagation using a batch size of 256 and an ADAM optimizer (Kingma and Ba 2014).Refer to the original IVIM 3 -NET (Kaandorp et al 2021) article for more information.
To mediate the effect of image noise at S(b = 0) affecting all normalized measurements acquired at other bvalues, IVIM 3 -NET models simultaneously predicted S(b = 0) along with other model parameters.This was shown to improve fit results in a previous study (Kaandorp et al 2021).Hence, the loss term is of the form where S(b) is the actual signal in a given voxel, and Ŝ0 and ˆ( ) S b are model predictions.For example, a biexponential model outputs Ŝ0 , f, D * and D, for which ˆ( ) ( )

Model-independent IVIM parameters
The AUC was calculated for each voxel using a middle Riemann sum.The AUC was furthermore calculated in 3 subsets of b-values, defining a low-, middle-, and high-range AUC (AUC L , AUC M , AUC H , respectively).The ranges were defined as: The method of calculating AUC, AUC l , AUC m , and AUC h using middle Riemann sums is shown graphically in figure 2. The choice of boundaries between the AUC subregions was chosen to separate the AUC into regions of perfusion and diffusion effects.Small and large b values are thought to correspond to perfusion and diffusion effects, respectively (Le Bihan 2019).The 'threshold' region separating perfusion and diffusion effects is thought to be approximately 200 s mm −2 , so the AUC m was included as a threshold region of the signal versus b-value curve.
During AUC calculation, S 0 values predicted from the ADC IVIM 3 -Net model were used in place of actual S 0 values.This was done to reduce the impact of noise associated with a single voxel on the determination of the AUC.
The anticipated relationship between AUC and exponential parameters can be found by expressing the AUC in terms of the area under monoexponential, biexponential, and triexponential functions.This is done by integrating the fitted curves from b = 0 to b = 1000 s mm −2 .For example, the area under a monoexponential curve is:  From this it is seen that AUC and ADC are negatively correlated.For the biexponential model, which shows that AUC is negatively correlated with D, and is roughly independent of D * and negatively correlated with f assuming D * ≈ 10D (Le Bihan 2019), and f is small.For the triexponential model which shows that AUC is roughly independent of D 1 * and D 2 * and negatively correlated with f 1 and f 2 , assuming and f i is small, i ä {1, 2}.Therefore, while AUC depends on all parameters in a given exponential model, it is mostly affected by (apparent) diffusion coefficients and perfusion fractions.

Comparing parameters
Whole-mean parameters were computed within all parotid glands.To compare functional utility of AUC parameters with exponential model parameters, AUC parameters were included in the dose-response comparison of IVIM parameters in parotid glands.This allowed the practical utility of parameters for use in dose response studies to be analysed, based on the hypothesis that parameter statistics are related to tissue composition and characteristics at a histological level.The significance of correlation differences between AUC and exponential parameters was determined by applying a Pearson-corrected Fisher transformation to correlations.For each of the four AUC parameters separately, we also tested the overall significance of AUC correlations being higher than exponential parameter correlations.For this, the significance of correlation differences between the AUC parameter and each exponential parameter was then determined using a 1-sample t-test with all Z differences (9 differences between a given AUC parameter and all exponential parameters).
The contrast-to-noise ratio (CNR) of IVIM parameters inside parotid glands and in a ring just outside the border of the glands was determined before and after denoising.This ring had a thickness of two voxels, generated by dilating parotid gland masks by two voxels before subtracting original masks.Voxels located outside the body were not included in the surrounding ring mask.All images before and after radiotherapy were Lastly, the relative importance of IVIM parameters for capturing the variance in the dataset was determined by evaluating parameter contributions to the primary principal component obtained through singular value decomposition.The primary principal component is the optimal linear combination of all input features for describing the maximum variation in the data.This is also the eigenvector with the greatest eigenvalue in the eigendecomposition of the covariance matrix.This was computed by first creating a 2D array, where each row represents a single voxel in an image and each column represents a given feature.Voxels within the body in all cropped images were included as rows in this array.All biexponential, triexponential, ADC, and AUC features were included as columns.This yielded an m × 13 matrix for which the singular value decomposition was computed.The squared projection of the primary principal component on each of the 13 parameters was then calculated and plotted.

Results
Denoising effectively suppressed random noise while enhancing fine detail within images.This improvement is supported quantitatively by significant (p < 0.001) reductions in image TV and higher blind image quality scores (table 4).
Predicted blur kernels were similar among patients (mean inter-patient inner product of 0.88 ± 0.10) and intra-patient blur kernels were nearly identical (mean inner product was 0.97 ± 0.02).Projections of the mean predicted kernel are shown in figure 3. Kernel voxels were mainly confined to axial planes.Highest variance was seen in the anterior-posterior direction.The general shape of pseudokernels were predicted; however, kernel voxels were mostly confined within one voxel of the centre (figure 4  The average blur kernel predicted for patients is projected onto the three primary imaging planes.The kernel was mostly confined to axial planes, which was expected due to elongation of voxels in the superior-inferior direction (slice thickness >2× axial slice pixel spacing).Axial pixel spacing size was 1.6 × 1.6 mm 2 .The average slice thickness was 3.8 mm (range 3.4-3.9)Image denoising strengthened correlations between whole-mean parotid dose and relative IVIM parameter changes following radiotherapy (figure 7).Changes in the AUC, AUC M , AUC H and ADC had the strongest correlations with dose, especially after denoising.As expected with such a small dataset, individual correlations and their changes following denoising were statistically insignificant.However, the overall occurrence of AUC correlations being higher than exponential parameter correlations was statistically significant (p < 0.001) for all four AUC parameters.This 'overall' significance suggests that the higher correlations of AUC parameters with dose are unlikely to have occurred randomly.Similarly, while denoising did not significantly increase any individual parameter's correlation with dose, the overall increase in correlations seen after denoising was significant (p < 0.02).
Regression slopes for changes in parameters versus dose are also shown in figure 7. AUC parameters were found to decrease following radiotherapy.Changes in ADC had the highest slope with respect to dose in both denoised and original images.AUC, AUC M , AUC H , D 2 *, f, and D biexp had similarly high dose slopes, while D biexp * , D triexp , D 1 *, f 1 , f 2 , and AUC L had lower slopes.Figure 7 also shows the relative contribution of each parameter to the first principal component obtained via the SVD analysis.AUC, AUC H and ADC captured the highest proportion of variance in the data, with f 1 also capturing a relatively high proportion of the variance.Pairwise correlations between all parameter combinations are displayed in figure 7. AUC, AUC M , and AUC H were highly correlated with one another, and also with ADC, D (biexponential and triexponential), f, and f 1 , while AUC L was most correlated with f 2 and D 1 *.
Denoising increased the CNR of IVIM parameters, as shown in figure 8. Denoising improved the CNR of all parameters.AUC parameters had slightly higher CNRs than exponential parameters.Figure 9   The neural blind deconvolution denoising process was tested by convolving fixed pseudokernels with previously denoised diffusion images before restarting the denoising process.Four pseudokernels, including a regular Gaussian, as well as a Gaussian elongated in each of the 3 primary image axes, were applied.Axial pixel spacing size was 1.6 × 1.6 mm 2 .The average slice thickness was 3.8 mm (range 3.4-3.9)

Discussion
The proposed model-independent parameter, the AUC, was effective for characterizing the signal decay versus b-value curve.AUC parameters may be more effective than traditional IVIM parameters for predicting various clinical outcomes, based on their high correlations with radiotherapy dose levels, high relative CNRs, and high contributions to the primary principal component of the total dataset.Changes in the AUC, AUC M and AUC H following radiotherapy were more correlated with dose than diffusion coefficients derived from exponential models.It is notable that changes in the ADC following radiotherapy were more correlated with dose than diffusion coefficients derived from biexponential and triexponential models.Perfusion fractions for biexponential and triexponential models, f and f 1 , were more correlated with dose than their corresponding model's diffusion coefficient, D.
Neural blind deconvolution for denoising IVIM MRI enhanced fine detail in images while suppressing noise.This was reflected by higher blind image quality metrics, BRISQUE and CLIP (Mittal et al 2012, Wang et al 2022), higher CNR, and reduced TV.While both BRISQUE and CLIP scores significantly increased following denoising, it should be noted that the CLIP score was only raised from 0.37 to 0.35 (on average).This is likely a cause of CLIP being ill-suited for non-natural images.While BRISQUE is also trained using natural images, it analyzes voxel distributions based on correlations between certain pixel distributions and image quality, and is likely a favourable metric for IVIM MRI.The BRISQUE score was more than doubled by denoising.As image post-processing algorithms become more prominent, a useful follow-up work could develop specific blind image quality metrics for various medical image modalities.Table 5.All IVIM parameter averages and standard deviations inside parotid glands are listed.Averages were calculated over parotid glands from pre-radiotherapy (pre-RT) and post-RT scans separately, along with the significance of their difference.Results are calculated using both denoised and original images.Lastly, the final column lists the significance of differences between denoised and original values before radiotherapy.Significance values were determined using a paired t-test.

Denoised images
Original images  This analysis was performed with a small dataset of 5 patients having scans acquired before and after radiotherapy.While no individual correlation was statistically significant (as expected with such a small dataset), the overall increase of correlations after denoising was significant (p < 0.02).Similarly, the overall higher observed correlations between AUC parameters over exponential parameters was statistically significant (p < 0.001).The singular value decomposition (SVD) of the entire dataset (all patients before and after radiotherapy) was computed, and the relative variance captured by each parameter in the first principal component was plotted (bottom).
Denoising resulted in overall higher correlations between dose and IVIM parameter values, especially with AUC, AUC M , and AUC H .This suggests that denoising may be beneficial for future studies seeking to use IVIM imaging to assess radiotherapy-induced tissue damage.Overall AUC parameters had significantly higher correlations with dose than other exponential parameters.Denoising was found to significantly increase the correlations of parameter changes with dose, when considering all parameters.It should be noted that no individual correlations for a given parameter were significant on their own, nor were there changes with dose.This was expected due to our small dataset, and is a weakness of this study.It would be beneficial for these relationships to be further evaluated on a larger dataset of patients with scans before and after radiotherapy.While the correlation of dose with IVIM parameters is an interesting function to consider, and higher correlations may indicate better model-building potential, there are relevant issues to consider when using it to assess parameter quality.The dose response is expected to vary from patient-to-patient, so perfect correlations are not expected.However, when comparing multiple parameters, relatively higher correlations do seem desirable for feature selection purposes.
In both denoised and original images, AUC and AUC H had the highest relative importance as determined by principal component analysis.This is not surprising, as the AUC is a more condensed parameter for describing the signal decay curve (versus, for instance, three for the biexponential model and five for the triexponential model).Another interesting result was that the ADC had a larger contribution to the first principal component than the diffusion coefficients of biexponential and triexponential models.Perfusion fractions, f (biexponential) and f 1 (triexponential) had a higher contribution to the first principal component than diffusion or pseudodiffusion coefficients.It could be argued that comparing AUC parameters with any single parameter from the biexponential or triexponential model is not relevant, as it is the combination of parameters that constitutes their model building potential.The takeaway point is that when choosing a single parameter to describe the signal decay curve, which is often the case anyways with the popularity of the monoexponential model, the AUC may be more effective than any exponential parameter, and at least has potential to improve model building when used in addition to monoexponential parameters.We do not recommend abandoning exponential models entirely and substituting in the AUC.Instead, the AUC could be included and compared with exponential parameters in future studies with relative ease, as it is easy to compute.
The results of this study suggest that AUC parameters and the ADC are effective IVIM parameters for modelling dose-related outcomes.However, this study included a small number of patients having post- radiotherapy data, rendering our dose response quantities prone to uncertainty.Further validation is needed to assess the efficacy of AUC parameters and neural blind deconvolution for denoising IVIM MRI.It appears that denoising images could improve parameter reproducibility, but validation is needed.
Predicted blur kernels were highly similar between patients, and kernels predicted for images encompassing the left and right parotid gland of each patient were nearly identical (mean inner product: 0.97 ± 0.02).Predicted kernels were mainly confined to axial planes, which was expected, due to the axial elongation of the voxel geometry used in this study.To avoid convergence towards a trivial solution, where the kernel becomes a delta function, a regularization term penalizing the MSE of kernel values above 0.7 was included.Previous neural blind deconvolution studies have simply applied MSE to all kernel voxels (Kotera et al 2021, Ren et al 2020), but we only penalized voxels values above a given cutoff to allow small contributions from neighbouring voxels to be unpenalized.The choice of cutoff threshold to use for this loss term was somewhat arbitrary, but it did not interfere with kernel predictions as maximum kernel values were much less than the threshold after optimization.
Denoising was effective for correcting cross-talk between neighbouring voxels, as seen in the pseudokernel analysis.Predictions were mostly confined to within one voxel of the kernel centres.Similar truncation was reported when deblurring PSMA PET images (Sample et al 2023); however, voxels were more truncated in the present study.This may be partially explained by the smaller kernel size employed in this study.Considering that a relatively small kernel for IVIM MRI images was anticipated (due to generally sharper MR images than PET) while also considering memory constraints, we decreased the kernel size to 11 × 11 × 11 from the 15 × 15 × 15 size used for PSMA PET (Sample et al 2023(Sample et al , 2024)).Furthermore, the relatively high noise present in EPI-based  IVIM MRI is expected to complicate blur kernel estimation.Neighbouring voxels had high total variation, making it difficult for 'true' voxel signal to be distinguished from surrounding 'blur' signal and noise contributions.It is also possible that the uniform blur kernel assumption used in this work was unsatisfactory and contributed to the noted truncation of predicted kernels.Modifying the current approach to accommodate spatially varying kernels could improve the overall quality of neural blind deconvolution for denoising IVIM MRI.
Parameters obtained from denoised and original images had similar population means and standard deviations (table 5), but paired t-tests between values found all parameters to be significantly different, except for the ADC.AUC values tended to be slightly higher in denoised images than original images, although differences did not exceed one standard deviation.Mean ADC values agreed with those previously reported in the literature (Bruvo and Mahmood 2021, Nada et al 2017, Kimura et al 2022).D and D * were lower than previously reported values (Becker et al 2016, Zhou et al 2016, Marzi et al 2018, Kimura et al 2022), both before and after radiotherapy.However, D and D * were within a standard deviation of the mean values reported by (Kimura et al 2022).ADC, and biexponential parameters, D and f, increased after radiotherapy with similar proportions to reported values (Zhou et al 2016).The disagreement between D * changes from previous literature values is not surprising given literature variability seen in pseudodiffusion coefficient estimates.Furthermore, our dataset of post-radiotherapy patients was relatively small, and accumulating more data will likely shift means.f within parotid glands was lower than values reported by Zhou et al (2016) and Becker et al (2016) but smaller than values reported by Kimura et al (2022).
It is not uncommon for IVIM parameters reported in the literature to be in conflict (Lemke et al 2011), which was the motivation for this work.For example, the significant variability of reported IVIM parameters in the liver has created difficulties establishing normal, baseline perfusion parameters (Cieszanowski et al 2018).The low SNR associated with echo-planar imaging is likely a significant contributor to these issues.Disregarding time constraints, a spin echo sequence may lead to much greater reliability in parameter estimates.Variations in SNR between scanners used in different studies has been reported to affect parameter estimates (Cieszanowski et  The intrinsically low SNR of echo planar imaging sequences makes denoising an important task for IVIM imaging.Incorporating a noise term into the neural blind deconvolution methodology allowed scanner related noise to be corrected.A notable result is the increased monotonicity of the signal decay curve found after image denoising (figure 6).It is expected that increasing diffusion gradient strength should result in strictly decreasing signal acquisition; however, this expectation was not enforced in the loss function used for optimization.This perceivable suppression of noise in the decay curve seems to validate the denoising process.
Characterizing the signal decay curve non-parametrically can be advantageous, as it has the potential to improve the utility and reproducibility of IVIM parameters without presuming the underlying mathematical model.While it is more desirable in theory to have physically interpretable parameters, it has not proven to be practical for IVIM imaging.The AUC and its sub-quantities are only one method of describing the signal decay curve.While this study affirms their utility, there remains room for further exploration and validation of other model-independent parameters.
Many methods of denoising 3D MRI images have been suggested (Wiest-Daesslé et al 2007, Coupé et al 2008, Manjón et al 2010, Manjón 2013, Lam et al 2014, Huang and Lin 2019), but neural blind deconvolution offers several notable advantages.The added autocoder network in this study allowed the noise term to be predicted, resulting in denoising effects, but neural blind deconvolution is fundamentally a deblurring technique as well.The denoising algorithm developed by Huang and Lin (2019) demonstrated performance improvements of kernel-based denoising over other traditional approaches.The present work advances this effort by adding deblurring and the unique advantage of using neural networks for improving image quality, as demonstrated by Ren et al (2020), Kotera et al (2021).A major benefit of neural blind deconvolution over supervised learning approaches with neural networks is the lack of a need for large training sets.Instead, neural blind deconvolution is a self-supervised approach that trains networks on each predicted image, as over-fitting is not an issue.The neural blind deconvolution architecture can be further optimized in future studies.The code containing models and training procedures is available online (Sample 2023).
Poor reproducibility of perfusion-related parameters could be partially explained by over-simplification of physiological processes inherent in simple exponential models, as it has been shown that the optimal choice of model is dependent on tissue type (Liao et al 2021).In reality, each image voxel containing biological tissue will encompass a complex arrangement of micro-structures, all having various perfusion fractions and molecular motion.Kuai et al (2017) used simulated data with 2-5 perfusion components to show that the variance of f and D * tend to increase as the number of perfusion components increases, or the difference between pseudodiffusion components increases.
This study was done on MR images acquired in the head and neck, as future goals involve investigating salivary glands using IVIM MRI.The present methodology can also be applied and optimized in other anatomic sites, MRI sequences, and other medical imaging modalities.The goal of this study was to improve overall image and parameter map quality so that inter-parotid gland parameter heterogeneity can be assessed with higher stability.An interesting follow-up study could assess improvements in tumor volume delineation after denoising with neural blind deconvolution.
This study was performed in a clinical centre with limited scanner access.Institutions with greater access to scanners could continue investigating this methodology with further patient and phantom studies.For example, a useful study could acquire images with both low and high SNRs to estimate how denoising improves low SNR diffusion image quality.Phantom images would be particularly useful to assess this methodology with known ground truth values and without patient motion effects.

Conclusion
Model-independent parameterization of the IVIM signal decay curve and denoising with neural blind deconvolution could improve the reproducibility of parameters and the practical utility of IVIM imaging for developing outcome-predictive models.IVIM parameters have shown substantial variability throughout the literature.A shift in the methodology for characterizing the signal decay curve offers a potential solution to mitigate this inconsistency.

Figure 1 .
Figure 1.The blind deconvolution architecture used for denoising IVIM images is illustrated.G x and G δ are symmetric, convolutional auto-encoder networks for predicting the denoised image, x, and its additive noise contributions, δ.G k is a fully-connected network for predicting the blur kernel, k.The fidelity loss term used in optimization compares the original image with x * k + δ by using backpropagated model gradients from the loss function for iterative optimization.

Figure 2 .
Figure2.AUC parameters were calculated using middle Riemann sums.On the left, the AUC within a single voxel is calculated graphically.In a separate voxel on the right, AUC L , AUC M , and AUC H are calculated (yellow, rose, blue), which sum to the regular AUC.
). Axial slices of denoised and original images with several b-values are shown in figure 5. Image denoising significantly (p < 0.001) increased the monotonicity of signal versus b-value curves, as shown in figure 6. Spearman's rank correlation coefficient, r s , between signal and b-values was −0.80 ± 0.31 in denoised images and −0.60 ± 0.29 in original images.

Figure 3 .
Figure3.The average blur kernel predicted for patients is projected onto the three primary imaging planes.The kernel was mostly confined to axial planes, which was expected due to elongation of voxels in the superior-inferior direction (slice thickness >2× axial slice pixel spacing).Axial pixel spacing size was 1.6 × 1.6 mm 2 .The average slice thickness was 3.8 mm (range 3.4-3.9) ADC and biexponential parameters, D and f, increased following radiotherapy with similar proportions to values noted in the literature within parotid glands (Zhou et al 2016, Marzi et al 2018).Changes in D * following radiotherapy did not agree with previously reported values (Zhou et al 2016, Marzi et al 2018).
shows reduction in TV upon denoising.TV was normalized by standard deviations of parameters for comparison's sake.Parameter maps of biexponential parameters corresponding to the same slice are shown in figure 10.AUC parameter maps for denoised and original diffusion images are shown in figure 11.All slices in figures 10 and 11 correspond to the same location, as shown with b = 0 s mm −2 in figure 10.Denoising suppressed image noise, and clearly reduced TV within images.Anatomical borders and details had better definition after denoising.The perfusion fraction parameter map, f, had particularly improved anatomical detail after denoising (figure 10).AUC parameter maps in figure 11 had similar patterns, with different intensity emphases in the AUC, AUC L , AUC M , and AUC H .

Figure
Figure4.The neural blind deconvolution denoising process was tested by convolving fixed pseudokernels with previously denoised diffusion images before restarting the denoising process.Four pseudokernels, including a regular Gaussian, as well as a Gaussian elongated in each of the 3 primary image axes, were applied.Axial pixel spacing size was 1.6 × 1.6 mm 2 .The average slice thickness was 3.8 mm (range 3.4-3.9)

Figure 5 .
Figure 5. Axial slices through parotid glands of denoised (left) and original (right) diffusion images at 4 different b-values are shown.Denoising appears to effectively suppress noise while enhancing fine detail within images.Here in particular, the posterior region of the parotid gland becomes discernible only upon denoising.The parotid gland mask as contoured using anatomical T1 images is shown overlayed with the b = 0 s mm −2 image at the top.

Figure 6 .
Figure 6.Axial slices through the parotid gland of denoised and original diffusion images with b = 1000 s mm −2 are shown, along with the signal decay curve of each in a single representative voxel (indicated in the denoised image).Signal values were normalized to statistical Z-values within the body.Denoising appears to suppress random noise while illuminating fine detail within images.The loss function used during image denoising did not impose a decrease in signal value with increasing b-value, yet signal values in denoised images appear to naturally decrease more monotonically and more smoothly than in original images.

Figure 7 .
Figure7.Pairwise Spearman's rank correlation coefficients, r s , between dose and relative changes in biexponential, triexponential, ADC, and AUC parameters derived from IVIM images are displayed (top) for denoised and original images (left and right).Best-fit regression slopes of parameters versus dose are shown underneath in green.This analysis was performed with a small dataset of 5 patients having scans acquired before and after radiotherapy.While no individual correlation was statistically significant (as expected with such a small dataset), the overall increase of correlations after denoising was significant (p < 0.02).Similarly, the overall higher observed correlations between AUC parameters over exponential parameters was statistically significant (p < 0.001).The singular value decomposition (SVD) of the entire dataset (all patients before and after radiotherapy) was computed, and the relative variance captured by each parameter in the first principal component was plotted (bottom).

Figure 8 .
Figure8.The contrast to noise ratio (CNR) between IVIM parameters in the parotid glands and parameters in the periphery of the parotid inside the body (within two voxels).Denoising tended to increase the CNR, and was highest for AUC parameters.CNR was calculated as the difference between means divided by the RMS standard deviation.

Figure 9 .
Figure 9.The total variation (TV) of diffusion images and derived parameter maps decreased following denoising.TV in each voxel is defined as the average difference in the voxel value from neighbouring voxel values in the  x ,  y , and  z directions.Top: The TV for all IVIM parameters in the parotid gland before and after denoising.Bottom: the TV for diffusion images at all b values before and after denoising.Error bars denote standard deviations over patients.

Figure 10 .
Figure 10.Axial slices through the parotid gland of biexponential model parameter maps are shown, derived using denoised (top) and original (bottom) diffusion MR images.

Figure 11 .
Figure 11.Axial slices through the parotid gland of AUC parameter maps are shown, derived using denoised (top) and original (bottom) diffusion MR images.Axial b = 0, AUC, AUC L , AUC M , and AUC H image slices derived from denoised (top) and original (bottom) are shown.
al 2018, Liao et al 2021).Variability of b-value distributions (Dyvorne et al 2014) and voxel sizes also contribute to variability among parameter estimates in the literature.It has been reported that fitting a monoexponential function, S(b)/S(0) = (1 − f ) e − b• D is more reliable than a biexponential fit for differentiating pathological grades of esophageal squamous cell carcinoma (ESCC) (Liu et al 2021).Heightened measurement error of signal at low b-values (Bihan et al 1991, Pekar et al 1992) leads to sub-optimal conditions for estimating perfusion effects, and estimation of ADC in low b-value regions has shown poor reproducibility (Koh et al 2009).Noise-levels have also been shown to greatly impact parameter estimates (Iima et al 2015).Reproducibility of IVIM parameters is impacted by a lack of standardization of voxel sizes and b-value distributions, which several studies have attempted to optimize (Zhu et al 2019, Perucho et al 2020, Riexinger et al 2020, Lemke et al 2011, Raju 2021, Paganelli et al 2023).One interesting follow-up study could compare the impact of voxel size on AUC parameter estimates.Due to limited scanner resources, we were unable to acquire images with multiple voxel sizes or fields of view.Another follow-up study could optimize the smallest b-value distribution needed to accurately estimate the AUC.In the present study, b-values were selected to separate perfusion and diffusion effects in the small b-value regime.For this reason, values were biased towards b-values less than 100 s mm −2 .More b-values between 100 and 1000 s mm −2 would likely improve AUC parameter quality.Despite the suboptimal b-value distribution used for estimating the AUC in this study, the AUC appeared to have higher quality/utility than exponential parameters.

Table 1 .
The optimization algorithm for updating network weights to predict deblurred PSMA PET images.This algorithm builds off of Ren et al.ʼs proposed joint optimization algorithm Compute loss and back-propagate gradients 9. Update G i x , d G i and G i k using the ADAM optimizer (Kingma and Ba 2014) 10.

Table 3 .
(Barbieri et al 2019)for IVIM 3 -NET predictions of exponential model parameters.Bounds were set to over-encompass anticipated parameter ranges, in order to compensate for diminishing gradients at sigmoid asymptotes(Barbieri et al 2019).

Table 4 .
Blind image quality metrics, CLIP and BRISQUE are compared for original (y) and denoised (x) diffusion MR images.Results were averaged over all slices and b-values.The significance of the difference was tested using a paired t-test.
ADC, biexponential, triexponential, and all AUC parameter statistics in parotid glands are summarized in table 5.