Denoising High-Field Multi-Dimensional MRI With Local Complex PCA

Modern high field and ultra high field magnetic resonance imaging (MRI) experiments routinely collect multi-dimensional data with high spatial resolution, whether multi-parametric structural, diffusion or functional MRI. While diffusion and functional imaging have benefited from recent advances in multi-dimensional signal analysis and denoising, structural MRI has remained untouched. In this work, we propose a denoising technique for multi-parametric quantitative MRI, combining a highly popular denoising method from diffusion imaging, over-complete local PCA, with a reconstruction of the complex-valued MR signal in order to define stable estimates of the noise in the decomposition. With this approach, we show signal to noise ratio (SNR) improvements in high resolution MRI without compromising the spatial accuracy or generating spurious perceptual boundaries.


INTRODUCTION
Ultra-high field magnetic resonance (MR) imaging at 7 Tesla and beyond has enabled neuroscientists to probe the human brain in vivo beyond the macroscopic scale (Weiskopf et al., 2015). In particular, quantitative MRI techniques (Cercignani et al., 2018) have become more readily available and offer the promise of quantitative information about the underlying microstructure. Unfortunately, as the size of the imaging voxel decreases well below the cubic millimeter, so does the signal to noise ratio (SNR), and the achievable resolution even when imaging a small portion of the brain remains limited, and requires multiple averages at the highest resolutions (Stucht et al., 2015;Federau and Gallichan, 2016;Fracasso et al., 2016).
Multi-parametric quantitative methods acquire multiple images within a single sequence, in order to estimate the underlying quantity (Helms et al., 2008;Metere et al., 2017;Caan et al., 2018). This is comparable to diffusion-weighted MR imaging (DWI) where multiple images with different weighting are acquired. Recent work in diffusion analysis demonstrated that the intrinsic redundancy of the signal across these images can be employed to separate signal from noise with a principal component analysis (PCA) over small patches of the images (Manjón et al., 2013;Veraart et al., 2016). This principle can be transferred to multi-parametric quantitative MRI: in this work we present an extension of the PCA denoising to the recently described MP2RAGEME sequence, which provides estimates of R1 and R2 * relaxation rates as well as quantitative susceptibility maps (QSM) in a very compact imaging sequence (Caan et al., 2018).
The MP2RAGEME data comes with additional challenges for the classical PCA approaches. First, the images are complexvalued, comprising of both magnitude and phase needed for the estimation of quantitative MRI parameters. This complex nature of the data needs to be taken into account. Second, they include only a few images compared to the high number of directions acquired in DWI. Interestingly, we can make fairly simple assumptions about the dimensionality of the signal, as it contains primarily R1, R2 * , susceptibility and proton density (PD) weighting. Taking these features into account, we were able to effectively denoise high resolution MP2RAGEME data without averaging, and maintain spatial precision. We studied the impact of denoising on the computation of quantitative maps, and its practical impact for reconstructing thin vessels and delineating small, low contrast subcortical nuclei such as the habenula, a notoriously difficult-to-image, small structure in the subcortex.

Data Acquisition
Our denoising method focuses on the recently developed MP2RAGEME sequence (Caan et al., 2018) which combines multiple inversions and multiple echoes from a magnetizationprepared gradient echo sequence (MPRAGE) to simultaneously obtain an estimate of quantitative T1, T2 * and susceptibility. The specific protocol of interest is comprised of five different images: a first inversion with T1 weighting, followed by a second inversion with predominantly PD weighting and four echoes with increasing T2 * weighting (Figure 1).
For the experiments, we used a high-resolution imaging slab centered on the subcortex with the following parameters: resolution = 0.5 mm isotropic, TR = 8.33 s, TR 1 = 8 ms, TR 2 = 32 ms, TE 1 = 4.6 ms, TE 2A−D = 4.6/12.6/20.6/28.6 ms, TI 1 /TI 2 = 670/3738 ms, α 1 /α 2 = 7/8 • . Data was obtained on a 7T system (Philips Achieva, Netherlands) with a 32 channel rf-coil (Nova Medical Inc., United States) for five healthy human subjects as part of an ongoing atlasing study approved by the Ethics Review Board of the Faculty of Social and Behavioural Sciences, University of Amsterdam, Netherlands (approval number: 2016-DP-6897). All subjects provided written informed consent for the study.

Complex Signal Reconstruction
One of the main drawbacks of existing local PCA methods for denoising is that they work on magnitude images, which follow Rician distributions. A simple correction proposed in Baselice et al. (2009) andEichner et al. (2015) for removing the Rician bias in DWI averaging is to revert to the complex signal, taking into account the phase information. Methods exist to perform PCA directly in the complex domain (Horel, 1984) but those are focusing estimating variations in spatial dynamic processes. The MRI signal being static in nature, its real and imaginary parts can simply be treated as a two-dimensional vector for real-valued PCA.
Phase and magnitude of the MRI signal have very different noise properties, so it is not adequate to just stack them into a vector for the PCA. Instead, we follow the complex signal strategy by reconstructing the real and imaginary parts of the complex signal, which have the same noise characteristics. However, phase contains additional variations due to non-local effects of air cavities around the brain, which bring severe ringing artifacts in the reconstructed data (Figure 2A). In Eichner et al. (2015), the global phase information is removed with a total variation method, which generally respects the location of phase wraps. However, we found that this approach is not effective in regions where the phase wraps have high frequency, and there are residual phase artifacts in the local phase. While these have little consequences for their averaging application, they provide systematic variations for the PCA decomposition, which is undesirable here. We therefore start with a full phase unwrapping (Abdul-Rahman et al., 2005) followed by a total variation smoothing (Chambolle, 2004) of the unwrapped phase ( Figure 2B). The residual phase variations are used as local phase, and combined with the magnitude to reconstruct real and imaginary parts of the complex signal ( Figures 2C,D). Note that here, unlike (Eichner et al., 2015), we do not discard the imaginary part of the signal as it contains valuable information about the noise and residual anatomical information.

Local PCA of Complex Signal
The complex signal, now comprising ten image dimensions all similar in nature, is then processed following the local overcomplete PCA approach of Manjón et al. (2013). In short, the images are cut into small, overlapping patches of NxNxN voxels, and the M contrasts combined into a N 3 xM matrix. The average patch value per contrast is subtracted, and the matrix decomposed via singular value decomposition (SVD) to yield the eigenvectors and associated singular values (square roots of the eigenvalues) of the covariance matrix across the patch. In this work, we used patches of size N = 4, and M = 10.
The overlapping patches are combined following the technique of Manjón et al. (2013), weighting each patch by the number of kept eigenvectors W patch = 1 1+M kept . Once recombined, the eigenvectors rapidly change from a rich information content, encoding boundaries, to pure noise (Figure 3). However, the decision boundary between actual signal and noise is variable across the brain, due to the presence of different tissue types, multiple tissue boundaries, etc. In order to infer for each patch the number of eigenvectors to keep, we thus need to first quantify the expected noise distribution over the SVD.

Estimating the Noise
Noise estimation in advanced MRI is challenging: variations in coil sensitivity, non-local susceptibility effects and dependencies to the B0 field as well as various acquisition techniques will affect the signal and the noise differently in different regions. The original local PCA denoising methods of Manjón et al. (2013)   used elaborate estimates of the magnitude image noise, taking into account its Rician properties. A recent extension of the work uses random matrix theory to model the expected distribution of random noise eigenvalues and determine its threshold (Veraart et al., 2016). While this approach is theoretically grounded, it may be sensitive to disturbances when applied to data with a limited number of eigenvalues. In Figure 4, we simulated small 4 × 4 × 4 voxel patches with 10 dimensions (as in the MP2RAGEME sequence) and 3 non-zero eigenvalues. While the random matrix theory method performs well in non-interpolated data, its performance degrades strongly when the noise becomes correlated across voxels from image interpolation.
In this work, we propose a more robust approach based on two properties of the complex signal: first, that the noise is locally Gaussian, and second that the spread of singular values for Gaussian noise can be reasonably well approximated by a straight line, at least when the number of dimensions in the decomposition is small. This property is retained by interpolation, which makes it possible to perform denoising after image registration, for instance when motion correction is FIGURE 3 | Local PCA decomposition: first five singular value maps (top) and eigenvector maps (bottom) from highest to lowest. Eigenvectors were first reformatted into 4 × 4 × 4 patches and averaged using the same averaging weights as in the denoising, singular values were set as a constant over the patch and averaged in the same way. needed as in DWI or fMRI, or after fat navigator-based motion correction in structural MRI (Gallichan and Marques, 2017).
Our algorithm to estimate the noise level proceeds as follows: first, a line is fitted to the M/2 lowest singular values of the patch decomposition using linear least squares (here, M/2 = 5). Then, every singular value above a factor of α above the expected noise level given by the fitted line is kept, while the others are removed (in this work, we used α = 5%). Thus, the main requirements of our method are that: (1) local signal variations across contrasts in each individual patch are Gaussian-distributed, and (2) the intrinsic dimension of the data is lower than half of the number of acquired images. The patches are then reconstructed and averaged across the image, the complex images separated into magnitude and phase, and finally the discarded global phase variations are reintroduced and wrapped, to obtain data as similar as possible to the original input. In addition, a map of the number of kept eigenvectors and of the noise fitting residuals weighted by patch are computed for quality control (see Figure 5).
The entire algorithm, from phase pre-processing to noise estimation and recombination is summarized in Figure 6.

Software Package
To make our method more widely usable, we released the denoising algorithm in Open Source in the Nighres python library (Huntenburg et al., 2018) 1 . The algorithm can be called as the nighres.intensity.lcpca_denoising module in python (3.0+) scripts and notebooks. The algorithm itself is implemented in Java 8 for increased efficiency in the IMCN Toolkit 2 and wrapped in python with JCC (see Huntenburg et al., 2018 for details on the integration). 1 https://www.github.com/nighres/nighres/ 2 https://www.github.com/imcn-uva/imcn-imaging/ A sample data set along with a complete processing script as well as the statistics from our experiments are available on FigShare 3 .

Quantitative Mapping
The main interest of quantitative MR mapping techniques such as the MP2RAGEME sequence is to obtain estimates of the MR parameters of T1, T2 * relaxation times and susceptibility χ from the measured images. In order to do so, we used the look-up table method of Marques et al. (2010) to generate T1 estimates, a simple regression in log domain to obtain T2 * (Miller and Joseph, 1993) and the TGV-QSM reconstruction algorithm of Langkammer et al. (2015) to create quantitative susceptibility maps (QSM), following the approach of Caan et al. (2018) to estimate QSM for each of the three last echoes of the second inversion and take the median. As a byproduct, the T1 estimation also generates bias-corrected T1-weighted images, and the T2 * fitting produces an S0 baseline image devoid of T2 * effects, which will be used in some of the processing. T1 and T2 * mapping methods are included in the Nighres library, and the TGV-QSM algorithm is also freely available 4 .

Manual Delineations
In order to obtain structure-specific measures of quality and to evaluate the practical benefits of the method for segmentation purposes, we performed a manual segmentation study. To evaluate the benefits of the denoising, we selected the habenula, a small structure ventral to the thalamus, which together with the pineal gland and stria medullaris forms the epithalamus (Hikosaka, 2010;Strotmann et al., 2013;Mai et al., 2016). Since we cannot distinguish between the medial and lateral habenula on MRI, both were included in our delineations. Given its size, complex shape, and heterogeneous composition, the habenula requires high resolution images and detailed anatomical knowledge to be reliably delineated on MRI.
In this study, the habenula was delineated by two raters on the reconstructed quantitative T1 maps obtained with or without denoising. The consensus mask (regions labeled by all raters on all images) were used to measure noise properties, while the overlap of the masks was used to measure inter-rater reliability.

Co-registration to Standard Space
Previously proposed methods for local PCA denoising have shown decreased performance when handling interpolated data due to the changes in the noise distribution (Manjón et al., 2013;Veraart et al., 2016). To test the performance of our noise estimation approach, we co-registered our data to standard space by aligning it first with a whole brain image of the same subject and then to the MNI152 template at 0.5 mm resolution. Both registrations used the first inversion of the MP2RAGEME sequence, and were performed with linear registration in ANTS (Avants et al., 2008). Phase images were unwrapped before transformation to avoid interpolating across phase wraps, and quantitative maps were computed before and after denoising in standard space.

Region Labeling
To evaluate the signal improvement across a variety of structures, we used manual delineations of the striatum, globus pallidus internal and external segment, subthalamic nucleus, substantia nigra and red nucleus obtained from lower FIGURE 5 | Maps of estimated noise threshold from actual data: number of kept eigenvectors (left) and R 2 goodness of fit of the linear fit of the lowest singular values of the PCA (right) for data in original acquisition space (A) and linearly co-registered into MNI space (B). As in Figure 3 full-sized images were obtained by setting the value of each patch to either the corresponding number of kept eigenvectors or the R 2 score and averaging the patches over the image. resolution MP2RAGEME images (0.64 × 0.64 × 0.7 mm resolution) with higher SNR. The delineations were performed independently by two raters and the conjunction between them was used as the structure mask. The masks were co-registered linearly to the high resolution slab based on their scanner coordinates followed by a rigid registration in ANTS. The intensities of the quantitative maps inside each region were averaged and used to compute signalto-noise ratio statistics (i.e., mean intensity over standard deviation inside each structure) of the original images and reconstructed quantitative maps before and after denoising in original space.

Vascular Reconstructions
Finally, we tested the capabilities of the denoising to maintain fine details by segmenting the vasculature with the method of Huck et al. (2019) on R2 * quantitative maps. The algorithm uses a spatial vessel filter followed by global diffusion, and is particularly sensitive to small vessels (Bazin et al., 2016). The filter was run with standard parameters on both the original and denoised R2 * maps after skull stripping using the S0 image and T1 map obtained during quantitative mapping (Bazin et al., 2014), and the result was overlaid on the bias-corrected T1-weighted image for orientation.

Impact on Quantitative Maps
As shown in Figure 7, the LCPCA denoising strongly impacts the estimation of quantitative signals, as the SNR gains of multiple images are combined. Visually, we found a strong gain in particular for R2 * mapping, which is quite noise-sensitive with the low number of echo times of the MP2RAGEME sequence. While the improvements are subtle at the whole brain scale, they are very clear when focusing on small regions, particularly in the subcortex where the original SNR is low.

SNR Measures
The denoising systematically improved the SNR in all structures, although not identically across regions and contrasts (Figure 8). On individual MP2RAGEME images, the improvement was modest, while denoising had a stronger cumulative impact on quantitative map estimation, especially for R1. While R2 * appears visually improved, the SNR measures were more similar partly due to the heterogeneity of R2 * values in many anatomical structures. QSM is the least affected, both visually and quantitatively, which is expected as the quantitative susceptibility reconstruction algorithm includes its own spatial regularization method [in this case, total generalized variation (TGV)]. Still, some more heterogeneous regions of the subcortex appear smoother after denoising. Interestingly, the highest gains were FIGURE 7 | Comparison of quantitative MRI reconstruction: left, from the original MP2RAGEME images; right from the denoised result.
FIGURE 8 | SNR improvement for the delineated structures and the different contrasts for the five subjects of the experiment (from left to right: first inversion, second inversion echo 1 to 4, estimated R1 map, estimated R2 * map, estimated QSM; mean and standard deviation across subjects). SNR is computed as the mean signal over its standard deviation in each structure, and the improvement is the difference between original and denoised data.
obtained for the smaller structures such as substantia nigra, subthalamic nucleus and red nucleus. Results were similar across the five subjects of the experiment.

Impact on Manual Delineations
Manual delineations of the habenula are challenging, and interrater agreement is low. Our denoising improved the consensus (Figure 9), but not to the point of an acceptable level of reproducibility for measuring anatomical quantities such as volume or shape. Again, there were no notable variations of the improvement across subjects.

Effects of Interpolation
When performing the denoising step after image co-registration to a standard space and interpolation, we observe only small differences in the LCPCA algorithm behavior (Figure 5). The interpolation procedure tends to slightly increase the dimension of the kept signal, probably due to the mixing of signals across voxels. However, the proposed dimensionality estimation procedure appears largely robust to interpolation.

Impact on Vascular Reconstruction
Brain vasculature is very difficult to image, as most of the vessels have sub-millimeter resolution even at the surface of the brain. Because most methods use local information to identify vessels, they are very sensitive to noise, and prone to large amounts of false detections. Figure 10 illustrates the potential benefits of the proposed denoising in separating vessels from noise, while preserving the fine details of the vasculature. A full test of the consistency of detected vessels across multiple scans would be needed to assess that the method does enhance vessels, FIGURE 9 | Comparison of manual delineations of the lateral habenula with or without denoising for the five subjects of the experiment. (A) Inter-rater agreement between delineations, (B) 3D rendering of the lateral habenula in one subject, as delineated by each rater (red and green outline, respectively), (C) Zoomed-in rendering.
FIGURE 10 | Comparison of vasculature reconstruction with or without denoising (maximum intensity projection over 15 mm of the estimated vessel probability map (overlaid on the T1-weighted image for the subject, for orientation). Insets show regions of the cerebellum (left) and subcortex (right) where lower SNR seems to artificially increase the detected vascular density in the original image.
however, it is worth noting that these fine structures are preserved in the denoising.

DISCUSSION
In this work we presented a new local PCA denoising technique, using the full complex signal acquired to correct acquisitions with multiple image contrasts such as DWI or quantitative MRI sequences. Even with a sequence like MP2RAGEME, where only five different images are used in order to recover T1, T2 * , and QSM, the proposed approach was able to effectively reduce the noise without introducing blurring artifacts. The homogeneity of structures compared to their boundary was improved over the basal ganglia structures, and delineation of small brain structures such as the habenula was shown to be more reproducible in high resolution images with high levels of imaging noise.
By using both magnitude and phase information, the proposed method captures the entire noise distribution, and stays in the statistically simpler domain of complex Gaussian perturbations. However, obtaining high quality phase images is not trivial for all scanners and imaging sequences and can be a limitation also in retrospective processing, as the phase is commonly discarded. Note also that certain acceleration techniques such as partial Fourier encoding may sacrifice the phase signal estimation. In the case of quantitative MRI sequences based on the MP2RAGE method or focusing on QSM, phase data is required to obtain the quantitative maps, so it is important to take it into account and denoise it. While there is no conceptual requirement to use both phase and magnitude in the denoising, lowering the number of dimensions reduces the applicability of the method: in the case of the MP2RAGEME using only five dimensions in order to separate four-dimensional signals from noise is challenging, and the linear fit of the noise to the last half of the local PCA eigenvalues becomes unreliable. Yet, preliminary experiments with other relaxometry methods such as multi-parameter mapping (Weiskopf et al., 2013), which acquire between 14 and 20 images, indicated that the separation of signal and noise based on magnitude alone was possible (unpublished data).
The main requirements of the noise estimation method are that: (1) local signal variations across contrasts in each individual patch are Gaussian-distributed, and (2) the intrinsic dimension of the data is generally lower than half of the number of acquired images. The first requirement is easily met in small patches, regardless of the type of data acquired. The second requirement depends on the type of MR sequence and tissue properties under consideration, but can be met simply by running twice the same sequence, as is commonly done for increasing SNR by classical averaging. Note also that in the case of a single contrast PCA approach reduces to simple averaging with no added value.
The main interest of the PCA-based approaches for denoising high resolution images is that unlike most other methods, they do not impose spatial regularity but rather enforce regularity across contrasts, preserving small anatomical details without blurring or inducing artificial boundaries. Recent advances in spatial filtering techniques have demonstrated significant improvements, but still introduce some form of spatial averaging (see Milanfar, 2013 for a review). Note that a thorough comparison with advanced spatial denoising techniques was performed by Manjón et al. (2013) in the context of DWI data, showing comparable or better performance for LPCA. In addition, while the impact of noise removal on conventional MR images remains subtle, it offers the option to push MR sequences beyond the usually accepted limits of thermal noise as long as enough signal remains to be reliably detected. As ultra-high field advances toward higher and higher resolutions, such denoising methods may become essential part of the imaging protocol for multi-contrast anatomical imaging in the same way they have become a standard tool for advanced DWI pre-processing (Veraart et al., 2016).
Finally, the proposed method is generally applicable to other relaxometry sequences such as multi-echo GRE, multi-echo MPRAGE or multi-parametric maps (MPM), as well as DWI or (multi-echo) fMRI. In some cases, e.g., T2 * relaxometry or multiecho fMRI, it is important to note that the known relationship between echoes could also be used to further discriminate signal from noise (Kundu et al., 2017). The robustness of the proposed noise threshold estimation technique to interpolation makes it also interesting for other MR imaging protocols that include multiple images with different contrasts, even if acquired sequentially or even at different resolutions, as long as the ratio of total number of images to actual number of measured contrast mechanisms is sufficient to estimate noise properties.

CONCLUSION
Here we presented a new local PCA method to denoise high resolution multi-parametric quantitative MRI data. Combining magnitude and phase data, we could differentiate between signal-and noise-induced variations with a simple model that is robust to interpolation. The resulting quantitative images are automatically regularized and additional anatomical detail is visible in low contrast regions. The denoising software is openly available as part of the IMCN Toolkit (see text footnote 2) and the Nighres library (see text footnote 1). The proposed method can be extended to denoise other MR imaging sequences with similar properties, namely partially redundant contrasts and low intrinsic dimensionality. While the method is already efficient to denoise MR images acquired with cutting edge methods at the lower limits of SNR, we hope they may help further to push MR imaging toward acquiring even more challenging data where the noise may visually dominate but significant amounts of signal are still available.

DATA AVAILABILITY STATEMENT
All of the data from this study cannot be shared publicly because of current GDPR privacy regulations. However, statistical tables and a representative data set for a single subject can be accessed openly here: https://uvaauas.figshare.com/projects/ Denoising_High-field_Multi-dimensional_MRI_with_Local_ Complex_PCA/61832.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethics Review Board of the Faculty of Social and Behavioural Sciences, University of Amsterdam, Netherlands. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
P-LB developed the study. MM, WZ, and MC acquired and pre-processed the data. P-LB and AA performed the validation experiments. All authors prepared and edited the manuscript.

FUNDING
This work was partially supported by a NWO-Vici grant (BF) and a NWO-STW grant (AA and BF).