{'fraction_non_alphanumeric': 0.037790843760993016, 'fraction_numerical': 0.010804563043369013, 'mean_word_length': 4.768115942028985, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian Image Quality Transfer with CNNs: Exploring Uncertainty in dMRI Super-Resolution

In this work, we investigate the value of uncertainty modeling in 3D super-resolution with convolutional neural networks (CNNs). Deep learning has shown success in a plethora of medical image transformation problems, such as super-resolution (SR) and image synthesis. However, the highly ill-posed nature of such problems results in inevitable ambiguity in the learning of networks. We propose to account for intrinsic uncertainty through a per-patch heteroscedastic noise model and for parameter uncertainty through approximate Bayesian inference in the form of variational dropout. We show that the combined benefits of both lead to the state-of-the-art performance SR of diffusion MR brain images in terms of errors compared to ground truth. We further show that the reduced error scores produce tangible benefits in downstream tractography. In addition, the probabilistic nature of the methods naturally confers a mechanism to quantify uncertainty over the super-resolved output. We demonstrate through experiments on both healthy and pathological brains the potential utility of such an uncertainty measure in the risk assessment of the super-resolved images for subsequent clinical use.


Introduction and Background
Algorithmic and hardware advancements of non-invasive imaging techniques, such as MRI, continue to push the envelope of quality and diversity of obtainable information of the underlying anatomy. However, their prohibitive cost and lengthy acquisition time often hinder the translation of such technological innovations into clinical practice. Poor image quality limits the accuracy of subsequent analysis, potentially leading to false clinical conclusions. Therefore, methods which can efficiently and reliably boost scan quality are in demand.
Numerous machine learning based methods have been proposed for various forms of image enhancement, generally via supervised regression of low quality (e.g., clinical) against high quality (e.g., experimental) image content. Alexander et al. [1] propose a general framework for supervised image quality enhancement, which they call image quality transfer (IQT). They demonstrated this with a random forest (RF) implementation of super-resolution (SR) of brain diffusion tensor images (DTIs) and estimation of advanced microstructure parameter maps from sparse measurements. More recently, deep learning has shown additional promise in this kind of task. For example, [2] proposed a CNN model to upsample a stack of 2D MRI cardiac volumes in the through-plane direction. The SR mapping is learnt from 3D cardiac volumes of nearly isotropic voxels, requiring a clinically impractical sequence due to its long scan time. Another application of CNNs is the prediction of 7T images from 3T MRI [3], where both contrast and resolution are enhanced. Current methods typically commit to a single prediction, leaving users with no measure of prediction reliability. One exception is Bayesian IQT [4], which uses a piece-wise linear Bayesian model, implemented with a RF to quantify predictive uncertainty over high-resolution (HR) DTIs and demonstrate its utility as a surrogate measure of accuracy.
This paper proposes a new implementation of Bayesian IQT via CNNs. This involves two key innovations in CNN-based models: 1) we extend the subpixel CNN of [5], previously limited to 2D images, to 3D volumes, outperforming previous models in accuracy and speed on a DTI SR task; 2) we devise new architectures enabling estimates of different components of the uncertainty in the SR mapping. The first enables us to bring the performance benefits of deep learning to this important problem, as well as reducing computation time to super-resolve the entire brain DTI in 1 s. For our second contribution, we describe two kinds of uncertainty which arise when tackling image enhancement problems. The first kind of uncertainty, which we call intrinsic uncertainty is defined as the irreducible variance of the statistical mapping from low-resolution (LR) to HR. This inherent ambiguity arises from the fact that the LR to HR problem is one-to-many, and is present independent of the amount of data we collect. We model the variation in intrinsic uncertainty over different structures within the anatomy through a per-patch heteroscedastic noise model [6]. The second kind of uncertainty, which we call parameter uncertainty, quantifies the degree of ambiguity in the model parameters that best explain the observed data, which arises from the finite training set. We account for it through approximate Bayesian inference in the form of variational dropout [7].
We first evaluate the performance of the proposed CNN methods and the benefits of uncertainty modelling by measuring the deviation from the ground truth on standard metrics. Human Connectome Project (HCP) dataset [8] and the Lifespan dataset (http://lifespan.humanconnectome.org/) are used for the quantitative comparison. We also test the benefits of the CNN-based methods in downstream tractography through SR of Mean Apparent Propagator (MAP)-MRI [9]. Lastly, we investigate the utility of uncertainty maps over the predicted HR image from the probabilistic CNN methods by testing on images of both healthy subjects and brain tumour patients.

Method
As in [1,2,3], we formulate the SR task as a patch-wise regression where an input LR image is split into smaller overlapping sub-volumes and the resolution of each is sequentially enhanced. For this, we propose a baseline CNN, on which we build by introducing two complementary ways of accounting for uncertainty to achieve a more robust model. Baseline network: Efficient subpixel-shifted convolutional network (ESPCN) [5] is a recently proposed method with the capacity to perform real-time perframe SR of videos while retaining cutting-edge performance. We extend this method to 3D and use this as our baseline model (3D-ESPCN). Most CNNbased SR techniques [2,10,11] first up-sample a low-resolution (LR) input image (e.g. through bilinear interpolation, deconvolution, fractional-strided convolution, etc.) and then refine the high-resolution (HR) estimate through a series of convolutions. These methods suffer from the fact that (1) the up-sampling can be a lossy process and (2) refinement in the HR-space has a higher computational cost than in the LR-space. ESPCN performs convolutions in the LR-space, upsampling afterwards. The reduced resolution of feature maps dramatically decreases the computational and memory costs, which is more pronounced in 3D. More specifically the ESPCN is a fully convolutional network, with a special shuffling operation on the output (see Fig. 1). The fully convolutional part of the network consists of 3 convolutional layers, each followed by a ReLU, where the final layer has cr 2 channels, r being the upsampling rate. The shuffling operation takes an input of shape h × w × cr 2 and remaps pixels from different channels into different spatial locations in the HR output, producing a rh×rw×c image, where h, w and c denote height, width and number of channels. The 3D version of this shufflling operation S is mathematically given by S(F ) i,j,k,c = F [i/r],[j/r],[k/r],(r 3 −1)c+mod(i,r)+r · mod(j,r)+r 3 · mod(k,r) where F is the pre-shuffled feature maps. [5] showed that the combined effects of the last convolution and shuffling is effectively a learned interpolation.
At test time, the network takes each subvolume x in a LR image, and predicts the corresponding HR sub-volume y. The network increases the resolution of the central voxel of each receptive field, e.g. the central 2 3 output voxels are estimated from the corresponding 5 3 receptive field in the input, coloured yellow in Fig. 1. By tessellating the predictions from shifted inputs x, the whole HR volume is reconstructed.
Given a training set , we optimize the network parameters by minimising the sum of per-pixel mean-squared-error (MSE) between the ground truth y and the predicted HR patch µ θ (x) over the training set. θ denotes all network parameters. This is equivalent to minimising the negative log likelihood (NLL) under the Gaussian noise model p(y|x, θ) = N (y; µ θ (x), σ 2 I).
Here, HR patches are modelled as a deterministic function of LR patches cor-rupted by isotropic noise with variance σ 2 . On the assumption that the model is correct, the variance σ 2 signifies the degree of irreducible uncertainty in the prediction of y given x, and thus the intrinsic uncertainty in the SR mapping defined in the introduction. However, the quality of this intrinsic uncertainty estimation is limited by the quality of likelihood model; the baseline network assumes constant uncertainty across all spatial locations and image channels, which is over-simplistic for most medical images. Heteroscedastic likelihood: Here we introduce a heteroscedastic noise model to approximate the variation in intrinsic uncertainty across the image. The likelihood becomes p(y|x, θ 1 , θ 2 ) = N (y; µ θ1 (x), Σ θ2 (x)) where both mean and covariance are estimated by two separate 3D-ESPCNs µ θ1 ( · ) and Σ θ2 ( · ) as a function of the input. The mean network makes predictions and the covariance network estimates the intrinsic uncertainty (see Fig. 2). We only use the diagonal of Σ θ2 (x), which quantifies estimated intrinsic uncertainty over individual components in µ θ1 (x). Adapting the NLL to heteroscedastic model gives (xi) i.e. mean squared Mahalanobis distance and H θ (D) = 1 N N i=1 log detΣ θ2 (x i ) i.e. mean differential entropy. Intuitively, M θ (D) seeks to minimise the weighted MSE under the covariance while H θ (D) keeps the 'spread' of Σ θ2 (x) under control. Bayesian inference through variational dropout: The baseline 3D-ESPCN and heteroscedastic model neglect parameter uncertainty, relying on a single estimate of the network parameters. In medical imaging where data size is commonly limited, this point-estimate approach potentially leads to overfitting. We combat this with a Bayesian approach, averaging over all possible models p(y|x, θ) weighted by the (posterior) probability of the parameters given the training data, p(θ|D). Mathematically this is p(y|x, D) = E p(θ|D) [p(y|x, θ)]. However, this expectation is intractable because: 1) evaluation of p(θ|D) is intractable, and 2) p(y|x, θ) is too complicated for the expectation to be computed in closed form. Variational dropout [7] addresses this problem for neural networks, using a form of variational inference where the posterior p(θ|D) is approximated by a factored Gaussian distribution q ϕ (θ) = ij N (θ ij ; m ij , s 2 ij ). The algorithm injects Gaussian noise into the weights during training, where the amount of noise is controlled by ϕ = {m ij , s 2 ij }, which are learnt. At test time, given a LR input x, we estimate the mean and covariance of the approximate predictive distribution q * are samples of the network parameters (i.e. the filters) from the approximate posterior q ϕ (θ). We use the sample mean as the final prediction of an HR patch and the diagonal of the sample variance as the corresponding uncertainty. When we use the baseline model, the first term in the sample variance reduces to σ 2 I. Implementation details: We employed a common protocol for the training of all networks. We minimized the loss using ADAM [12] for 200 epochs with learning rate 10 −3 . The best performing model was selected on a validation set.
As in [5], we use a minimal architecture for the baseline 3D-ESPCN, consisting of 3 convolutional layers with filters (3, 3, 3, 50) → (1, 1, 1, 100) → (3, 3, 3, r 3 c) where r is upsampling rate and c is the number of channels. The filter sizes are chosen so a (5, 5, 5) LR patch maps to a (r, r, r) HR patch, which mirrors competing random forest based methods [1,4] for a fair comparison. The heteroscedastic network of Section 2 is formed of two 3D-ESPCNs, separately estimating the mean and standard deviations. Positivity of the standard deviations is enforced by passing the output through a softplus function. For variational dropout we tried two flavours: Var.(I) optimises per-weight dropout rates, and Var.(II) optimises per-filter dropout rates. Variational dropout is applied to both the baseline and heteroscedastic models without changing the architectures.
All models are trained on datasets generated from 8 randomly selected HCP subjects [8], each consisting of 90 diffusion weighted images (DWIs) of voxel size 1.25 3 mm 3 with b = 1000 s/mm 2 . The training set is created by sampling HR subvolumes inside the brain region from the ground truth DTIs (or MAP-MRI coefficients) and then artificially downsampling to generate the LR counterparts. Downsampling is done in the raw DWI by a factor of r by taking a block-wise mean and then the DT or MAP coefficients are subsequently computed. Each network is trained on ∼ 4000 pairs of input/output patches of size 11 3 c and (7r) 3 c, amounting to ∼ 1.4 × 10 6 receptive field patch pairs of dimensions 5 3 c and r 3 c, which is roughly the same size as the maximal training set used in RF-IQT [1]. It takes under 30/120 mins to train a single network on DTI/MAP-MRI data on 1 TITAN X GPU.

Experiments and Results
Performance comparison for DTI SR: We evaluate the performance of our models for DTI SR on two HCP datasets. The first contains 8 unseen subjects from the same HCP cohort used for training. The second consists of 10 subjects from the HCP Lifespan dataset. The latter tests generalisability, as they are acquired with different protocols, at lower resolution (1.5 mm isotropic), and on  subjects of a different age range (45-75) to the original HCP data (22-36). We perform ×2 upsampling in each direction, measuring reconstruction accuracy with RMSE, PSNR and MSSIM on the interior and the exterior separately as shown in Fig. 3(b). This is important, as the estimation problem is quite different in boundary regions, but remains valuable particularly for applications like tractography where seed or target regions are often in the cortical surface of the brain. We only present the RMSE results, but the derived conclusions remain the same for the other two metrics. Fig. 3(a) shows our baseline achieves 8.5%/39.8% reduction in RMSE on the HCP dataset on the interior/exterior regions with respect to the best published method, BIQT-RF [4]. Note that IQT-RF and BIQT-RF are only trained on interior patches, and SR on boundary patches requires a separate ad-hoc procedure. Despite including exterior patches in training our model, which complicates the learning task, the baseline CNN out-performs the RF methods on both regions-this goes for the Lifespan dataset too. The 3D-ESPCN estimates whole HR volumes <10 s on a CPU and ∼ 1 s on a GPU, while BIQT-RF takes ∼ 10 mins with 8 trees. Faster reconstruction is achieved using input patches of size 22 3 , twice as large as the training input size 11 3 .
Heteroscedastic network further improves on the performance of 3D-ESPCN with high statistical significance on the interior region for both HCP and Lifespan data (p < 10 −3 ). However, poorer performance is observed on the exterior than the baseline. Using 200 weight samples, we see Var.(I)-CNN performs best on both datasets on the exterior region. Combination of heteroscedastic model and variational dropout (i.e. Hetero+Var.(I) or (II)) leads to the top 2 performance on both datasets on the interior region and reduces errors on the exterior to the level comparable or better than the baseline.
The performance difference of heteroscedastic network between the interior and the exterior region roots from the loss function. The term M θ (D) imposes a larger penalty on the regions with smaller intrinsic uncertainty. The network therefore allocates more of its resources towards the lower noise regions where the statistical mapping from the LR to HR space is less ambiguous. The dramatic reduction on the exterior error from variational dropout indicates its regularisation effect against such overfitting, and as a result also improves the robustness of prediction on the interior. Tractography with MAP-MRI SR: Reconstruction accuracy does not reflect real world utility. We thus further assessed SR quality with a tractography experiment on the Prisma dataset, which contains two DWIs of the same subject from a Siemens Prisma 3T scanner, with 1.35 mm and 2.5 mm resolution. The b-values and gradient directions match the HCP protocol. An ensemble of 8 het-ero+var.(I) CNNs super-resolves the MAP-MRI coefficients [9] derived from the LR DWIs (2.5 mm), then the HR MAP volume is used to predict the HR DWIs (1.25 mm). The final prediction is computed as the average estimate weighted by the inverse covariance as in RF-IQT. We also generate HR datasets by using IQT-RF and linear interpolation. Fig. 4 shows streamline maps of the probabilistic tractography [13] for the original LR/HR data and various upsampled images, and focuses on examples that highlight the benefits of reduced SR reconstruction errors. In the top row, tractography on the LR data produces a false-positive tract under the corpus callosum (yellow arrow in the 1st row), which tractography at HR avoids. Reconstructured HR images from IQT-RF and CNN avoid the false positive better than linear intepolation. Note that we do not expect to reproduce the HR tractography map exactly, as the HR and LR images do not aligned exactly. The bottom row shows shaper recovery of small gyral white matter pathways (green arrow) at HR than LR resulting from reduced partial volume effect. CNN reconstruction produces a sharper pathway than RF-IQT and linear interpolation. Visualisation of predictive uncertainty: We measure the expectation and variance of mean diffusivity (MD) and fractional anisotropy (FA) with respect to the predictive distribution q * ϕ (y|x) of Hetero+Var. (I) by MC sampling. Comparative results are shown in Fig. 5(a), where we drew 200 samples of HR DTIs from the predictive distribution of the Hetero+Var.(I). The uncertainty map is highly correlated with the error maps. In particular, the MD uncertainty map captures subtle variations within the white matter and central CSF, which demonstrates the potential utility of the uncertainty map as a surrogate measure of accuracy.    5(b) shows the best-performing SR model (Hetero+Var (II)) trained on a healthy HCP cohort applied to the DTI of a brain tumour patient. The raw data (DWI) with b = 700 s/mm 2 is processed as before with input voxel size 2 3 mm 3 . We show the input, SR image and uncertainty map. The ground truth is unavailable but the estimated image sharpens the input without introducing noticeable artifacts. The uncertainty map shows high uncertainty on the tumour, not represented in the training data, again illustrating the potential of the uncertainty maps to flag potential low accuracy areas.

Discussion
We present a super-resolution algorithm based on 3D subpixel-CNNs with stateof-the-art accuracy and reconstruction efficiency on diffusion MRI datasets. An application to the MAP-MRI coefficients indicates benefits to tractography in comparison with the previous methods. We also demonstrate that assimilation of intrinsic and parameter uncertainty in the model leads to best predictive performance. The uncertainty map highly correlates with reconstruction errors and is able to highlight pathologies. This can be used to gain insight into the 'black-box' prediction from CNNs. Understanding the behaviours of these uncertainty measures in unfamiliar test environments (e.g. pathologies) and their relations to predictive performance is an important future work for designing a more generalisable method. The presented ideas extend to many other quality enhancement problems in medical image analysis and beyond.