{'fraction_non_alphanumeric': 0.03433434650455927, 'fraction_numerical': 0.01052887537993921, 'mean_word_length': 4.721480244852532, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SNR-enhanced diffusion MRI with structure-preserving low-rank denoising in reproducing kernel Hilbert spaces

Purpose: To introduce, develop, and evaluate a novel denoising technique for diffusion MRI that leverages non-linear redundancy in the data to boost the SNR while preserving signal information. Methods: We exploit non-linear redundancy of the dMRI data by means of Kernel Principal Component Analysis (KPCA), a non-linear generalization of PCAto reproducing kernel Hilbert spaces. By mapping the signal to a high-dimensional space, better redundancy is achieved despite nonlinearities in the data thereby enabling better denoising than linear PCA. We implement KPCA with a Gaussian kernel, with parameters automatically selected from knowledge of the noise statistics, and validate it on realistic Monte-Carlo simulations as well as with in-vivo human brain submillimeter resolution dMRI data. We demonstrate KPCA denoising using multi-coil dMRI data also. Results: SNR improvements up to 2.7 X were obtained in real in-vivo datasets denoised with KPCA, in comparison to SNR gains of up to 1.8 X when using state-of-the-art PCA denoising, e.g., Marchenko- Pastur PCA (MPPCA). Compared to gold-standard dataset references created from averaged data, we showed that lower normalized root mean squared error (NRMSE) was achieved with KPCA compared to MPPCA. Statistical analysis of residuals shows that only noise is removed. Improvements in the estimation of diffusion model parameters such as fractional anisotropy, mean diffusivity, and fiber orientation distribution functions (fODFs)were demonstrated. Conclusion:Non-linear redundancy of the dMRI signal can be exploited with KPCA, which allows superior noise reduction/ SNR improvements than state-of-the-art PCA methods, without loss of signal information.


INTRODUCTION
Diffusion MRI (dMRI) is a noninvasive imaging modality that allows the characterization of tissue microstructure of biological tissues with an unrivaled level of quality and detail. When diffusion-encoding gradients are played out, the MR signal becomes sensitive to the diffusion of water molecules and their interaction with the surrounding microstructure 1 . Hence, the dMR signal, carrying unique information, can be used to probe the microstructural environment of tissues. Unfortunately, diffusion encoding gradients as well as the 2 spin-to-spin relaxation, both contribute to the characteristic attenuation of the dMRI signal during the diffusion probing time, which makes the signal-to-noise ratio (SNR) of the diffusion weighted MR images inherently low 2 . This is of particular concern for highresolution dMRI, as the SNR decreases even further due to a decrease in the voxel size. Low SNR not only complicates visual inspection but hampers quantitative analysis of informative tissue parameters, for example, by reducing the accuracy and precision of the parameter estimates 3 .
Increasing the SNR in dMRI is a major goal for the MRI community. Ultrahigh field dMRI 4,5,6 , advanced dMRI acquisition protocols 7,8,9,10,11,12,13,14 , or noise reduction techniques 15,16,17,3 , to name a few, are complementary approaches that have been shown to enhance the SNR. In this work, we focus on noise removal techniques, in particular, the thermal noise at the reception of the Radio-Frequency signal 15 , and which further propagates to the dMRI signal domain during image reconstruction 18 .
Noise reduction or denoising can be done by averaging 16 : the subject is imaged several times with identical image parameters, and the resulting images averaged. Under some statistical assumptions of the data, this simple technique can increase the SNR by a factor of √ scans with scans the number of repetitions. Evidently, this approach requires additional scan time, and as the acquisition time of a conventional dMRI protocol is already relatively long, averaging becomes impractical for routine use.
From a signal processing perspective, denoising can be seen as a post-processing approach, where an algorithm attempts to remove the noise, i.e., reduce the noise standard deviation, while maintaining the noise-free signal undistorted. Trivially formulated, denoising has been a longstanding problem in image processing, where many challenges need to be confronted, which are further aggravated in quantitative image modalities such as diffusion MRI. Indeed, early computer vision-based denoising algorithms applied to dMRI have shown to be detrimental to parameter estimation quality 3 . Some exemplar cases are the popular Total-Variationbased noise removal techniques 19,20 or more recently, nonlocal means algorithms 21,22 . The artefactual components that these algorithms introduce, though invisible to human eye, necessarily result in biased estimates in subsequent parameter quantification, thereby turning the improvement in noise reduction futile for quantitative analysis 3 . Other shortcomings are the loss of spatial resolution due to blurring in the presence of partial volume effects.
On the other hand, denoising algorithms that depart from exploiting spatial similarity but leverage "redundancy" of the dMRI signal along diffusion direction have been shown to suppress noise significantly while preserving the dMRI signal, with no apparent blurring or biases. Principal Component Analysis (PCA) based methods belong to this category of methods. PCA-based methods were first used in 23 , and ever since, have been thoroughly validated for diffusion MRI parameter quantification and substantially refined and improved. Most of the improvements that PCA algorithms have witnessed are based on the estimation of the number of signal-only principal components. By removing the rest of the components, attributed to noise, denoising is performed. While this threshold was heuristically set in 23 , the criterion was formalized in 3 , where an elegant approach was proposed relying on the theory of random covariance matrix theory, in particular, exploiting the universal Marchenko-Pastur law for eigenvalues. The MPPCA method of 3 has also been extended to general noise models other than the additive white Gaussian noise case 24 .
All dMRI PCA denoising methods are limited by the degree of redundancy in the signal with respect to the dimensionality of the data, i.e., the number of gradient directions. Several factors, including spatial resolution, number of b-values, number of gradient directions, determine the amount of redundancy in the data. The capability for noise reduction is hampered if the SNR is very low as in high-resolution diffusion MRI. We would like to emphasize here that this limitation is not attributable to the PCA method at hand but in the assumption that the covariance matrix of the signal is of sufficiently low rank.
Fortunately, there exist non-linear redundancies which are not captured by the linearity assumption implicitly adopted in PCA, but that can be exploited to enhance the SNR of dMRI substantially more than possible with current state-of-the-art approaches. The fundamental idea is to look for high dimensional non-linear spaces where the covariance matrix of the transformed dMRI signal turns out to be of low rank. Enforcing this prior knowledge, i.e., applying PCA in the transformed domain, we can denoise the signal in the high dimensional space and map it back to the original space. As a unified approach, this operation exploits the non-linear relationships within the data, and is referred to as a non-linear generalization of PCA. The whole process can be carried out by a technique called Kernel Principal Component Analysis (KPCA) 25 , where a kernel function implicitly determines the mapping to a high dimensional Hilbert space, where the dMRI signal is denoised 26 In this work, we introduce KPCA denoising to the dMRI community, and showcase it with a Gaussian kernel that maps the dMRI signal to an infinite dimensional space where redundancy is exploited. The parameters of the kernel as well as the rank of the covariance matrix are chosen in a data-driven manner, as those that provide the best signal representation according to the Mean Square Error (MSE) between the denoised and the underlying noise-free signal. The Stein Unbiased Risk Estimate (SURE) is used as a proxy of the MSE 27,28,29 , circumventing the need for the unobservable noise-free signal. Only knowledge of the noise statistics is required, which we input to the algorithm with state-of-the-art noise mapping techniques.
We thoroughly validate KPCA with realistic Monte-Carlo simulations as well as with several in-vivo human brain datasets acquired with submillimeter spatial resolution. In addition, KPCA was validated with an in vivo human brain multi-coil dMRI dataset, capturing non-linear redundancies in the coil and gradient directional domains simultaneously. In all cases, we confirmed superior noise reduction compared to the linear MPPCA method, which immediately translates to higher SNR enhancement, while preserving signal reliably, as confirmed by residual analysis. Finally, improved diffusion parameter estimation was invariably found when dMRI datasets were denoised with KPCA. A preliminary version of this work has been presented as an abstract at the ISMRM 2020 30 .

Redundancy of dMRI signals in canonical spaces: PCA denoising
In dMRI, it is often assumed that the diffusion signal carries redundant information between the gradient directions (also referred to as q-space). Redundancy can be elegantly captured by covariance matrix analysis. Let ∈ ℝ be the diffusion MR signal with the number gradient diffusion directions. Its second-order statistical characterization is given by its mean ∈ ℝ and covariance matrix ∈ ℝ × . The diffusion signal is said to be "redundant" if is rank-deficient, with rank substantially smaller than the dimensionality . In that case, , a low-rank matrix, can be written with eigenvalue decomposition = ∑

=1
. The low-rank diffusion signal with these statistical properties is given by the so-called "spike" model where are independent and identically distributed (IID) random variables with zero mean and unit variance. The diffusion signal is always corrupted by noise, ∈ ℝ , which is typically modeled as additive and with zero mean, i.e., = + . Denoising the signal can be formulated as an estimation problem, where the goal is to estimate the noise-free signal from the observed data . As is a full-rank matrix (noise is not redundant), it is precisely the low-rank nature of in comparison to ( << ) that allows to "separate" signal from noise, i.e., estimating reliably. Low-rank denoising is usually performed in image patches, each one containing diffusion signals , = 1, … , , with mean and covariance matrix . For the noisy spike model, the optimal estimate of in terms of the norm ∑ =1 || −̂ || 2 2 (assuming rank known) is given by 31 wherê ,̂ and̂ are obtained from the singular value decomposition (SVD) of the noisy data matrix by nullifying components with index > , and beinĝ the sample mean of . This denoising framework is also called Principal Component Analysis (PCA) denoising, witĥ as the principal components. As already mentioned in the introduction, all of the PCA-based dMRI denoising methods ( 22,3,32,24 ) fundamentally differs in the the way is estimated. No improvements are made in the intrinsic data model. In this work, we instead reconsider the low-rank model of Eq. 3. In particular, we are interested in the intrinsic redundancy in the noise-free image . If is comparable to , the optimal solution of Eq. 3 can hardly denoise, no matter which PCA algorithm we employ. MRI artifacts, partial volume effects, higher spatial resolution, non-conventional dMRI sequences, to name a few, all of these factors may increase the original rank of the signal. A possible way to increase the redundancy is to focus on , rather than , and increase the dimensionality of the signal by adding more gradient directions. Theoretically appealing, this modus operandi will necessarily necessitate additional scan time, which is often of concern in clinical settings.
In the quest for data redundancy, we look for information redundancy in data domains which are not necessarily the canonical space where conventional PCA is applied. The next section is devoted to motivate and formalize our approach, after which we present the novel denoising method in Section 3.1.1.

Redundancy of dMRI signals in high-dimensional Hilbert spaces
Our starting point is a function (⋅) that maps the diffusion signal from the native space ℝ to a "feature" space  , often high-dimensional, >> , being = ( ). We will deal with the definition of (⋅) later, but for now let us assume that the transformed data ( ) is redundant, i.e., the covariance matrix of ( ) is of rank  << , even if native rank is high. For example, any mapping which makes data "sparser" in the feature domain will "compress" the information more than in the native space. This redundancy translates into a lowrank covariance matrix ( ) . We can exploit this redundancy to denoise data in the feature space and eventually return to the native space to get the denoised signalŝ . Note that the entire process can be seen as a way to exploit non-linear redundancy that the dMRI signal can carry, and that could be otherwise difficult to capture with conventional PCA.
In  , similar to PCA, the optimal estimates of ( ) of rank  are those minimizing the error . They can be shown to be the projection of the mapped noisy signals ( ), = 1, … , , onto the feature space, ( ), In Eq. 5,̄ is the mean of ( ), = 1, … , ,̂ are the non-linear principal component directions, and the rest of parameters are obtained (we maintain the notation of Eq. 3) from the SVD of centered noisy projected data matrix While the low rank denoising is performed in  , we would like to come back to the native space. If we want to denoise the signal * at the center of the patch, we then look for that , which after being mapped to the feature space, ( ), turns out to be the closest to the projection ( * ) (Eq. 5), i.e., * = arg min || ( ) − ( * )|| 2 2 .
A fundamental result that is of high relevance for this work is the following. To apply PCA in the feature space  defined by the mapping (⋅), and to solve Eq. 7 in order to return to the native space, we do not need to know (⋅) explicitly, but just the inner product of the form ⟨ ( ), ( )⟩ for and in ℝ . Since ⟨ (⋅), (⋅)⟩ is a symmetric, positive definite function, it automatically defines a kernel function in ℝ × ℝ as ( , ) = ⟨ ( ), ( )⟩. Conversely, choosing a kernel function (⋅, ⋅) implicitly defines a mapping (⋅) 26 . Therefore, it is the kernel function that implicitly defines the feature space. Features spaces with this property are called reproducing kernel Hilbert spaces, and applying PCA in the feature space is termed Kernel Principal Component Analysis 26 . In the next section, we present our Kernel PCA denoising method in detail, elaborating on the selection of the kernel as well as the rank  .

Kernel Principal Component Analysis (KPCA) denoising 3.1.1 KPCA algorithm
Given a kernel (⋅, ⋅), the denoised signal̂ * in the feature space defined by (⋅, ⋅) can be written as 33 , an -dimensional column vector with all entries equal to one, and , the first  eigenvectors that solve the following eigenvalue problem where = − 1 is a "center" matrix and , the socalled kernel matrix ( × ), = ( , ). Finally, the coefficients , the components of projection of ( * ) onto the k-th non-linear principal component̂ , can be computed as with the -th coefficient of , and̃ ( * , ) equal to 34 The interested reader can find the mathematical proof of the derivation of KPCA in 33 and in Section 1.1 of the supplementary file of this submission

The choice of the kernel function
We showcase KPCA denoising for dMRI with a Gaussian kernel function, with ℎ the scale parameter. Gaussian kernels have shown excellent performance in machine learning tasks and are particularly interesting for dMRI denoising for the following reasons. The implicit feature space that the Gaussian kernel function generates can be shown to be infinite-dimensional 26 . As data tend to be sparser in high-dimensional spaces, higher redundancy is achieved by mapping the data with ℎ (⋅). As implied by the notation, we can control the shape of the mapping with the scale parameter ℎ, and, in fact, the components of ℎ ( ) decay with increasing ℎ. In that sense, by varying ℎ, we can adapt the level of redundancy of the dMRI signal in the feature space.
This aspect will be of high interest for the automatic selection of parameters. Finally, it is possible to demonstrate that, when ℎ → ∞, KPCA with a Gaussian kernel behaves as linear PCA in the canonical space 35 . Hence, being linear PCA then a particular case of KPCA with Gaussian kernel functions, it is expected that our KPCA denoising will perform typically better, as we confirm in this paper. More details about the implicit mapping related to the Gaussian kernel and the demonstration of the asymptotic equivalence of KPCA and PCA are given in the supplementary file, (Sections 1.2 and 1.3, respectively).
In addition, there are computational advantages in choosing the Gaussian kernel. The solution of Eq. 8 can be obtained in very short computational time with the approximation given with || ( * ) − ( )|| 2 2 calculated analytically. Details are given in Section 1.4 of supplementary file. Finally, in the discussion session, we elaborate on possible improvements of KPCA denoising by selecting more complex kernels.

Automatic parameter selection driven by noise statistics
Two parameters need to be selected for our KPCA method: the scale parameter ℎ and the rank  . Ideally, we would like to select those that best represent the noise-free signal * , for example, by quantifying the Mean Squared Error, (risk) 2 2 } for different choices of ℎ and  . Obviously, the ground-truth signal * is unobservable, and hence MSE is not computable. Instead, we use the Stein Unbiased Risk Estimate (SURE) 27 . Minimizing SURE can act as a surrogate for minimizing the MSE, with the critical difference that it does not require knowledge of * 28 . For an AWGN model like that of Eq. 2, SURE can be computed from the noisy signals , the denoised signal̂ * (ℎ,  ), and the standard deviation of the noise, . We estimate the noise maps of the DWI images using the method presented in 37 , with the assumption of Gaussian distributed data, which holds in our experiments as we show in the subsequent section. Then, for every voxel in the image patches, we fix and applied grid search minimization to get the optimal ℎ and  w.r.t. the SURE cost-function. We used the efficient implementation of the SURE method based on Monte-Carlo sampling 29 . We refer the reader to Section 1.5 of the supplementary file for more details about the SURE method for optimal parameter selection. In addition, a discussion is provided at the end of the paper about the extension of SURE to other noise models as well as different techniques to estimate ℎ and  that may be of interest.
An illustrative scheme of the KPCA denoising method used in this work is presented in Fig. 1 .

Experimental validation
We validated KPCA denoising using simulated and in-vivo human brain dMRI data, both quantitatively and qualitatively, and compared our results with MPPCA denoising 3 . Both algorithms were implemented in a sliding-window fashion, where for each patch, only the signal at the center was denoised. The selection of the parameters for KPCA was done as follows. The standard deviation of the noise was estimated with the method of 37 . The set of possible values of  to minimize the SURE was chosen to be in the range [1,30], and the scale parameter of the Gaussian kernel ℎ, was parameterized by ℎ = min-class where min-class is the average minimum distance between all pairs of signals in the patch 36 , and the values for was chosen from ten equidistant points in the interval [0.6, 6] (see Section 1.4 of the supplementary file).

Simulations
A Monte-Carlo based experiment was conducted to assess the benefits of KPCA denoising in subsequent diffusion parameter estimation. Similar to the patch-based simulation experiment in 3 , we generated 5 × 5 × 5 signals based on a diffusion tensor (axially symmetric) model and a total of gradient directions uniformly spread on the sphere with a given b-value, . The underlying fractional anisotropy (FA) and mean diffusivity (MD) for each tensor in the pach was sampled from a distribution with fixed mean FA GT and MD GT , and a standard deviation of 10% with respect to the mean. MC = 5000 zeromean uncorrelated Gaussian noise realizations were added to each of the noise-free = 125 signals. The standard deviation of the noise was parameterized by a nominal SNR value, i.e., = 1∕SNR.
The Normalized Root Mean Square Error (NRMSE) was used to compare the denoised signals with respect to the ground-truth signal. Diffusion tensor parameters were estimated from the log-linearized signals with a Linear Least . Solving the eigenvalue problem (2), we obtain the eigenvectors and eigenvalueŝ associated to the mapping ℎ (⋅), which defines the feature space where linear PCA is performed. The best feature space, i.e, optimal ℎ, and optimal rank  are selected as those that, after reprojection of the low-rank denoised signal to the native space, give the best signal representation in terms of the Stein Unbiased Risk Estimate (SURE) (3). Finally, the denoising signal at the center of the patch (4) is obtained by applying low rank denoising (with optimal  ) in the optimal feature space and reprojecting. Squares (LLS) estimator. Next, the FA and MD were estimated and compared to the ground-truth FA and MD, FA GT and MD GT .

In vivo human brain submillimeter resolution dMRI data
Whole human brain in-vivo submillimeter dMRI data were acquired and reconstructed with the generalized slice dithered enhanced resolution (gSlider) technique 8 , and denoised with KPCA and MPPCA. Two datasets with different spatial resolutions were considered.
Data was acquired 10 with a single-shot EPI sequence: Muti-Band = 2, partial Fourier = 6/8, phase-encoding (superiorinferior axis) under-sampling factor in-plane = 2, TR/TE = 4400/80 ms, five radio-frequency encoding pulses. The total acquisition time was about 25 min. Three repetitions were acquired to construct a gold-standard reference. Conventional gSlider 8 was used to reconstruct the data and obtain wholbrain isotropic 660 resolution. Prior to gSlider reconstruction, slice and in-plane GRAPPA was used for k-space and SMS reconstruction, and real-valued data was obtained with background phase correction 16 . Eddy-current and motion were corrected between all acquisitions using the FSL technique FLIRT. The three datasets were then denoised with KPCA and MPPCA, and compared to the averaged dataset, which is considered here as the gold-standard reference. Both algorithms were implemented in a voxel-wise fashion, with a sliding window of [5 × 5 × 5] voxels.

isotropic gSlider data
Whole human brain gSlider-SMS data were collected from a healthy male volunteer on a Siemens 3T Prisma scanner. Four scans of the full brain (FOV = 220 × 220 × 163 mm 3 ) were obtained. A total of 38 thick axial slices were acquired with in-plane resolution of 860 and matrix size 256 × 256. The diffusion protocol consisted of 64 diffusion-weighted images (diffusion directions uniformly distributed along sphere) at = 2000 ∕mm 2 and 8 b0-images. Data was acquired 11 with a single-shot EPI sequence: Muti-Band = 2, partial Fourier = 6/8, phase-encoding (posterior-anterior axis) under-sampling factor in-plane = 3, TR/TE = 3500/81 ms, five radio-frequency encoding pulses. The total acquisition time was about 20 min. Four repetitions were acquired to construct a gold-standard reference. Data was preprocessed and reconstructed as described for the 660 case. After affine registration, one of the datasets was denoised with KPCA and MPPCA (identical window size as before), and compared to the averaged dataset, the gold-standard reference.

Quantitative validation
We assessed the performance of KPCA denoising in signal preservation and parameter estimation, quantitatively. The NRMSE was used to compare the signal of the denoised diffusion-weighted (DW) images with the signal of the averaged data set. To assess the ability of KPCA denoising for SNR enhancement, we estimated the noise maps (noise standard deviation) of the denoised datasets with the homomorphic approach 37 . The SNR gain was defined as the ratio between the standard deviation of the noise in the original dataset and that of the denoised datasets. To demonstrate that KPCA preserves the underlying diffusion signal reliably, we calculated the normalized residuals between the noisy datasets and the denoised versions, and checked if any anatomical structure was present 3 .
We conducted DTI analysis and High Angular Resolution Diffusion Imaging (HARDI) validation. FA and MD maps were estimated with dtifit from FSL, and compared to the maps from the reference set. The NRMSE was used to assess the improved quality in parameter estimation. HARDI analysis was carried out with MRtrix3 38 . Fiber Orientation Distribution Functions (fODFs) were calculated in white matter area only, with the single-shell single-tissue Constrained Spherical Deconvolution (CSD) technique 39 . For each voxel, main fiber peaks were extracted, and the angular error compared to those from the reference set were calculated. The variability in the estimation of the ODF peaks was probed with the coherence metric, , ( ∈ [0, 1]), which was originally proposed in 40 and used in 3 . A high value of indicates low angular variability, that is, high angular precision.

Capturing non-linear coil and diffusion redundancy simultaneously
We investigated whether KPCA denoising can work at the reconstruction level, for example, by denoising multi-coil and diffusion data simultaneously. To that end, we used an in-vivo human brain DW image data set comprising of one b0-image and 15 diffusion gradient directions that were uniformly spread over the sphere (b = 1200 ∕mm 2 ). The acquisition protocol was as follows. With a 3T Philips scanner, an axial slice was acquired with a single-shot EPI sequence, matrix size = 70 × 91, in-plane resolution of 2 mm, multi-coil system with eight channels and no-undersampling factor. To create a goldstandard reference, 20 repetitions of the same axial slice were obtained. Prior to denoising, k-space data was transformed to image space with an inverse Fourier transform. Phase estimation was obtained by taking the complex argument of the image resulting from the inverse Fourier transform of low-pass filtered k-space data (center of the k-space). Complex conjugate phase correction was applied, and the real part was retained. No statistical correlation was assumed between the different coils.
To apply denoising, the coil and the diffusion dimension were merged into a single dimension, = 8 channels × 15 diffusion directions = 120. The size of the patch was [7 × 7]. As in the previous experiment, the NRMSE, the noise maps, and the normalized residuals were calculated both from the denoised dataset with MPPCA and KPCA. NRMSE maps for diffusion derived metrics, FA and MD were also computed.

Simulations
NMRSE results for the case b = 1200 ∕mm 2 and = 64 direction, are shown as bar plots in Figure 2 , whereas the rest of results are shown in Table format in the supplementary file  for brevity (Table S1, S2 and S3). In general, KPCA achieves the lowest NRMSE results in signal quality, FA and MD, and shows particularly drastic improvement in FA for low nominal SNR values, e.g., SNR = 5 and SNR = 8, which correspond to real SNR values (defined as the the noise-free dMRI signal divided by sigma) of 2 and 3, typically encountered in real data noisy scenarios as those like submillimeter resolution data shown in this paper. Interestingly, the superiority of KPCA over MPPCA denoising becomes more notorious for reduced number of diffusion directions = 32 or = 64. This is attributed to the lack of enough "linear" redundancy of the dMRI signal compared to low values of dimensionality . KPCA, however, as it performs low-rank denoising in a highdimensional space ( >> ), can achieved superior noise reduction, and, hence, improved parameter estimation.

In vivo human brain submillimeter resolution dMRI data
Denoised images with MPPCA and KPCA from the 660 and 860 m resolution datasets are shown in Fig.3 and Fig.4 , respectively. The original dataset (no denoising) as well as the gold-standard reference, three averages for the 660 m case and four averages for 860 m case are also shown.
Visually, KPCA denoising achieves a higher noise suppression than MPPCA without signal loss. No anatomical structure can be seen in the 'residual' images (original -denoised DW image) that are shown in Fig. 5 and Fig. S2 of the supplementary file (860 m dataset).
Statistical analysis of residuals confirms signal preservation in KPCA denoising. Any anatomical structure in the residual dataset will make the standard deviation higher than the noise standard deviation, 3 . By analyzing the -normalized residuals 3 , we found that in both cases, 660 and 860 resolution data, approximately follows a zero-mean Gaussian distribution with standard deviation 0.82 and 0.79, respectively, see blue dotted-blue line graphs representing the estimated pdf, ( ), of normalized residuals (logarithmic plot on the left, linear plot on the right). As the standard deviation is lower than the unit (see solid black-line representing zero-mean standard Gaussian distribution), we then can conclude that no anatomical structure is lost in KPCA denoising. The pdf ( ) was estimated with a kernel density estimator. The solid blue-line represents the analytical Gaussian pdf that best fits the data in a Maximum Likelihood sense.
The estimated noise maps after denoising are presented in Fig. 6 and Fig. S3 of the supplementary material. Note that the noise mapping method 37 we use to estimate assumes either a Gaussian or Rician distribution. Since, by assumption, the original data is Gaussian distributed (real-value phase corrected images 16 ) and the residuals are shown to be Gaussian, our assumptions are well founded.
KPCA achieves higher noise suppression while still reliably preserving signal, supported by the previous experiment with residuals. The lowest levels of noise are found when the original data is denoised with KPCA, indicating that KPCA enhances the SNR to a greater extent than what is achievable with MPPCA denoising. The SNR gain is more than 60% higher than those obtained with MPPCA, see Table. 1 . Superior noise removal performance as well as reliable signal preservation make the NMRSE (compared to the averaged data case) substantially lower than MPPCA, both in white and gray matter (Table. 1 ).
In good agreement with the findings from the simulation experiment, denoising improves diffusion parameter estimation, and in particular, KPCA denoising helps estimate quantitative parameters with lower statistical error. NRMSE of both estimated FA and MD are considerably lower (Table 1 ) when DTI is applied after KPCA denoising compared to MPPCA. The improvement is significantly noticeable in gray matter. Cortical gray matter seems better delineated in the FA maps that are obtained after denoising data with KPCA (see colorencoded FA maps in Fig. 7 ). This is highly relevant since mapping cortical gray matter is one of the main motivations of ultra-high resolution dMRI protocols 41 . In fact, error maps also suggest that the estimation of FA is improved to a greater extent in cortical areas. Color-encoded FA maps as well as errors map for the 860 case are shown in Fig. S4 of the supplementary material.   TABLE 1 Quantitative results from experiment with in-vivo human brain submillimeter resolution dMRI data. Note that in all cases i.e., brain, white matter (WM), and gray matter (GM), KPCA denoising achieves better results than MPPCA.
fODF estimation becomes more robust after denoising, and more accurate and precise angular directions can be achieved if data is first denoised with KPCA. As shown in Table 1 , lower angular errors (mean of the errors for the first, second, and third

FIGURE 5
Residual maps from the 660 resolution datasets after being denoised with KPCA. On top of the figure, the residual map from a given DW image, which shows no anatomical information. On the bottom, the probability density function of the residuals ( ) normalized by the level of noise, . For the statistics, the normalized residuals are taken for all diffusion directions and number of repetitions. Note that the residuals for KPCA approximately follows a Gaussian distribution (blue dotted line on both plots representing the estimated pdf). On blue solid-line the optimal analytical zeromean Gaussian distribution that best fits the data (Maximum Likelihood sense). Note that the standard deviation of the normalized residual, 0.82, is lower than 1 (black-line represents a zero-mean standard Gaussian distribution). peak) are achieved with KPCA. Graphs of the prevalence/probability of angular errors in the 660 data are shown in Fig.  8 and Fig. S5 (860 ). Clearly, the distribution of the KPCA angular errors is skewed to the left more than that of MPPCA and the original data, demonstrating lower angular errors in the white matter map obtained from the dataset denoised with KPCA rather than that obtained from the MPPCA denoised or original data. Graphs of the fODFs plotted in the three-fiber crossing area of Fig. 8 shows lower peak variability with KPCA denoising, a direct consequence of higher noise suppression. The coherence metric, , proposed in 40 is in agreement with this observation. As shown in the plot, both in Fig. 8 and Fig. S5, the prevalence graphs of for the KPCA are more right skewed than that of MPPCA or the original data. As a result, overall coherence metric values are higher for KPCA (Table 1 ), indicating higher angular precision could be achieved if data is denoised first with KPCA.

Capturing non-linear coil and diffusion redundancy simultaneously
Coil DWI images as well as coil-combined DW images are presented in Fig. 9 . The Sum of Squares (SoS) method was employed for coil combination.
As in the previous experiments, KPCA achieves higher noise suppression than MPPCA, and the result is comparable to the twenty average case. As expected, differences in noise reduction are less notorious in the SoS images, as this technique already denoises the data due to averaging. However, higher noise reduction and good structure preservation is still observed by inspecting the DW images that are denoised with KPCA. Similar to the experiment with high-resolution dMRI data, -normalized residual maps of the multi-coil data show no anatomical features, and signal preservation is confirmed by statistical residual analysis (Fig. S6 of the supplementary material). Furthermore, normalized residuals follow a Gaussian distribution with standard deviation less than one. Noise maps presented in Fig. S7 show higher SNR enhancement when KPCA denoising is applied, compared to MPPCA, i.e., 2.48 X and 1.73 X, respectively. NRMSE values (maps in Fig.  S8) in the whole brain were also lower for KPCA than MPPCA. Improved estimation of FA and MD compared to MPPCA is achieved as well (see Table S4. in the supplementary file).

DISCUSSION
We have shown using realistic simulations and in-vivo dMRI experiments that it is possible to achieve superior noise suppression than state-of-the-art linear PCA denoising while preserving the dMRI image structure, if non-linear redundancies in the data are exploited. No signal structure is removed, as confirmed by the residual analysis. The KPCA denoising methodology can be used to enhance the typically low SNR of dMRI protocols, without compromising signal integrity.
To exploit non-linear redundancy of the dMRI signal the key point is to apply PCA in a reproducing Kernel Hilbert space where the low-rank assumption of the covariance matrix holds to a greater extent than in the canonical linear PCA space. As the kernel determines the feature space, the choice of the kernel is an interesting problem that deserves to be discussed. We motivated the selection of the Gaussian kernel in Section 3.1.2. Though it provided excellent results, other kernels that are specifically tailored to certain features of the dMRI signal, for example, the angular information can be used. This could be accomplished by defining a corresponding spherical

FIGURE 6
Maps of the NRMSE (hot colormap) and noise level (gray colormap) for the denoised DW images at 660 isotropic resolution and b-value of = 1500 ∕mm 2 . Observe that KPCA denoising obtains the lowest level of noise (highest SNR gains) and lowest NRMSE.
covariance function for the diffusion directions, as done in 42,43 , and incorporating this covariance matrix into the conventional Gaussian kernel.
The selection of the rank  and the kernel parameters clearly affects the performance. The SURE method allows us to rely on the statistical distribution model of the data, providing the optimal representation in the MSE sense. Originally, the SURE approach was conceived for additive noise models where the covariance of the noise is diagonal and parameterized by a single noise level . This case gives excellent results in all of our experiments. However, it can be extended to other statistical distributions 44 , including Gaussian noise models with arbitrarily complex covariance matrices. This could be of interest for denoising multi-coil data where correlation between channels do exist, and different noise levels can be found for different channels, very often happening with parallel FIGURE 8 Angular error as well as angular precision, probed by coherence metric , for the peaks of the fODFS estimated with CSD after denoising the 660 isotropic resolution DW images ( = 1500 ∕mm 2 ). Further, corresponding fODFs maps in a representative crossing-fibers area are displayed. Observe the lower variability in the fODFs of KPCA denoising compared to MPPCA.
imaging reconstructions methods. Note that in the experiment with in-vivo multi-coil data, no under-sampling was carried out, and hence no statistical correlation existed in our dataset. Finally, it is worth noting that random matrix theory in kernel matrices is of much less interest for optimal rank selection than in the conventional PCA case. The difficulty of tracking noise statistics over the kernel transformation, and the asymptotic approximations necessary to obtain meaningful theoretical results 45 makes this line of action impractical. This is one of the main benefits of the SURE method. Indeed, though KPCA denoising performs the low-rank decomposition in the feature transformed space, where noise statistics are difficult to model, the optimal selection of the rank and the scale parameter is done in the native space after reconstruction, where the noise distribution model is well defined.
We would like to emphasize the broad applicability of KPCA denoising beyond conventional dMRI pulse sequences.
We envisage even further benefits of using KPCA denoising compared to conventional PCA in situations where the complexity of the dMRI signal increases. New developments in diffusion sequences such as multidimensional dMRI 46,47,48,49 are highly attractive applications for KPCA denoising. It is part of our future work to extensively evaluate KPCA denoising in tensor-encoding dMRI data 50 and extend our preliminary experiments on this kind of data. Combination of relaxometry and diffusion MRI data may be another application where the non-trivial redundancy between different modalities could be better exploited with KPCA 51,52,53,54 .

CONCLUSION
We introduce to the dMRI community a novel denoising technique, Kernel PCA, which goes beyond the linear compressibility assumption of PCA-based methods and exploits

FIGURE 9
Mid-axial, coronal and sagittal slices of multi-coil denoised DW images. DW images coil-combined with the SoS method are also shown.
the non-linear redundancies that is intrinsic to dMRI data. Substantially superior SNR-enhanced dMRI data can be obtained compared to PCA, without compromising signal integrity, in a short-computation time, and with no manual parameter tunning. We showcase the power of KPCA denoising with several in-vivo whole human brain submillimeter resolution datasets as well as conventional spatial resolution multi-coil dMRI data. Improved diffusion parameter estimation was observed in all cases compared to state-of-the-art PCA denoising, e.g., Marchenko-Pastur PCA (MPPCA). We believe KPCA denoising could be beneficial in any diffusion MRI processing pipeline and particularly critical when processing very low SNR data, as in high-resolution dMRI.