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

^{1}Integrative Model-Based Cognitive Neuroscience Research Unit, Department of Psychology, Universiteit van Amsterdam, Amsterdam, Netherlands^{2}Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany^{3}Spinoza Centre for Neuroimaging, Amsterdam, Netherlands^{4}Brain Imaging Centre, Amsterdam University Medical Center, Amsterdam, Netherlands^{5}Department of Psychology, Universiteit Utrecht, Utrecht, Netherlands

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 complex-valued, 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.

## Materials and Methods

### 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 magnetization-prepared 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).

**Figure 1.** The MP2RAGEME sequence: **(A)** First inversion, **(B–E)**: second inversion, first to fourth echo (top: magnitude, bottom: phase).

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_{2}_{A–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) and Eichner 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.

**Figure 2.** Phase pre-processing. **(A)** original phase map, **(B)** after unwrapping and global phase removal (estimating global phase variations with TV-smoothing and keeping only the residual phase), **(C)** reconstructed real signal combining magnitude and cosine of the phase, **(D)** reconstructed imaginary signal combining magnitude and sine of the phase. These reconstructed signals are the input for the local PCA algorithm.

### 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}=\frac{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.

**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.

### 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.

**Figure 4.** Noise properties of the PCA thresholding: estimation of the threshold to denoise a 4 × 4 × 4 voxel patch with 10 dimensions and 3 non-zero eigenvalues, with added Gaussian noise and a SNR of 10 **(left)** or 20 **(right)**, with no interpolation **(top)**, a single interpolation of all 10 dimensions **(middle)**, or multiple linear interpolations for each of the 10 dimensions **(bottom)**. All interpolations are linear interpolations shifting the 3-dimensional voxel grid by a uniformly distributed random offset of up to a half voxel. In each subfigure, the left side shows 20 eigen- and singular value examples with the corresponding threshold as a vertical bar, for the random matrix approach (in red) and the linear fitting approach (in blue). The right side shows the histogram of thresholds estimated over 10,000 simulated noise patterns. Dotted lines indicate the ideal number of eigenvalues to keep.

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 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).

**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.

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).

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 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 signal-to-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.

## Results

### 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.

**Figure 7.** Comparison of quantitative MRI reconstruction: **left**, from the original MP2RAGEME images; **right** from the denoised result.

### 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 obtained for the smaller structures such as substantia nigra, subthalamic nucleus and red nucleus. Results were similar across the five subjects of the experiment.

**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.

### Impact on Manual Delineations

Manual delineations of the habenula are challenging, and inter-rater 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.

**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.

### 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, however, it is worth noting that these fine structures are preserved in the denoising.

**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.

## 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 multi-echo 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).

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

We would like to thank Josephine Groot for her help in recruiting and scanning subjects and Jill van Dorp and Judith Blankenagel for their help in anatomical delineations. This manuscript has been released as a pre-print at BioRXiv (https://doi.org/10.1101/606582).

## Footnotes

**^**https://www.github.com/nighres/nighres/**^**https://www.github.com/imcn-uva/imcn-imaging/**^**https://uvaauas.figshare.com/projects/Denoising_High-field_Multi-dimensional_MRI_with_Local_Complex_PCA/61832**^**http://www.neuroimaging.at/pages/qsm.php

## References

Abdul-Rahman, H., Gdeisat, M., Burton, D., and Lalor, M. (2005). “Fast three-dimensional phase-unwrapping algorithm based on sorting by reliability following a non-continuous path,” in *Presented at the Optical Metrology*, eds W. Osten, C. Gorecki, and E. L. Novak, (Munich: International Society for Optics and Photonics), 32–40.

Avants, B. B., Epstein, C. L., Grossman, M., and Gee, J. (2008). Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain. *Med. Image Anal.* 12, 26–41. doi: 10.1016/j.media.2007.06.004

Baselice, F., Ferraioli, G., and Pascazio, V. (2009). Relaxation time estimation from complex magnetic resonance images. *Sensors* 10, 3611–3625. doi: 10.3390/s100403611

Bazin, P. L., Plessis, V., Fan, A. P., Villringer, A., Bazin, P. L., and Gauthier, C. J. (2016). “Vessel segmentation from quantitative susceptibility maps for local oxygenation venography,” in *Proceedings of the 2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI)*, (Piscataway, NJ: IEEE), 1135–1138.

Bazin, P. L., Weiss, M., Dinse, J., Schäfer, A., Trampel, R., and Turner, R. (2014). A computational framework for ultra-high resolution cortical segmentation at 7Tesla. *NeuroImage* 93, 201–209. doi: 10.1016/j.neuroimage.2013.03.077

Caan, M. W. A., Bazin, P., Marques, J. P., Hollander, G., Dumoulin, S., and Zwaag, W. (2018). MP2RAGEME: T 1, T 2 ^{∗}, and QSM mapping in one sequence at 7 tesla. *Hum. Brain Mapp.* 40, 1786–1798. doi: 10.1002/hbm.24490

Cercignani, M., Dowell, N. G., and Tofts, P. S. (eds) (2018). *Quantitative MRI of the Brain: Principles of Physical Measurement, Second Edition.* Boca Raton, FL: CRC Press.

Chambolle, A. (2004). An algorithm for total variation minimization and applications. *J. Math. Imaging Vis.* 20, 89–97. doi: 10.1023/b:jmiv.0000011321.19549.88

Eichner, C., Cauley, S. F., Cohen-Adad, J., Möller, H. E., Turner, R., Setsompop, K., et al. (2015). Real diffusion-weighted MRI enabling true signal averaging and increased diffusion contrast. *NeuroImage* 122, 373–384. doi: 10.1016/j.neuroimage.2015.07.074

Federau, C., and Gallichan, D. (2016). Motion-correction enabled ultra-high resolution in-vivo 7T-MRI of the brain. *PLoS One* 11:e0154974. doi: 10.1371/journal.pone.0154974

Fracasso, A., van Veluw, S. J., Visser, F., Luijten, P. R., Spliet, W., Zwanenburg, J. J. M., et al. (2016). Lines of Baillarger in vivo and ex vivo: myelin contrast across lamina at 7 T MRI and histology. *NeuroImage* 133, 163–175. doi: 10.1016/j.neuroimage.2016.02.072

Gallichan, D., and Marques, J. P. (2017). Optimizing the acceleration and resolution of three-dimensional fat image navigators for high-resolution motion correction at 7T. *Magn. Reson. Med.* 77, 547–558. doi: 10.1002/mrm.26127

Helms, G., Dathe, H., and Dechent, P. (2008). Quantitative FLASH MRI at 3T using a rational approximation of the Ernst equation. *Magn. Reson. Med.* 59, 667–672. doi: 10.1002/mrm.21542

Hikosaka, O. (2010). The habenula: from stress evasion to value-based decision-making. *Nat. Rev. Neurosci.* 11, 503–513. doi: 10.1038/nrn2866

Horel, J. D. (1984). Complex principal component analysis: theory and examples. *J. Clim. Appl. Meteorol.* 23, 1660–1673. doi: 10.1186/s13742-016-0117-6

Huck, J., Wanner, Y., Fan, A. P., Jäger, A. T., Grahl, S., Schneider, U., et al. (2019). High resolution atlas of the venous brain vasculature from 7 T quantitative susceptibility maps. *Brain Struct. Funct.* 224, 2467–2485. doi: 10.1007/s00429-019-01919-4

Huntenburg, J. M., Steele, C. J., and Bazin, P.-L. (2018). Nighres: processing tools for high-resolution neuroimaging. *GigaScience* 7:giy082. doi: 10.1093/gigascience/giy082

Kundu, P., Voon, V., Balchandani, P., Lombardo, M. V., Poser, B. A., and Bandettini, P. A. (2017). Multi-echo fMRI: a review of applications in fMRI denoising and analysis of BOLD signals. *NeuroImage* 154, 59–80. doi: 10.1016/j.neuroimage.2017.03.033

Langkammer, C., Bredies, K., Poser, B. A., Barth, M., Reishofer, G., Fan, A. P., et al. (2015). Fast quantitative susceptibility mapping using 3D EPI and total generalized variation. *NeuroImage* 111, 622–630. doi: 10.1016/j.neuroimage.2015.02.041

Mai, J. K., Majtanik, M., and Paxinos, G. (2016). *Atlas of the Human Brain*, 4 edn, San Diego, CA: Academic Press/Elsevier.

Manjón, J. V., Coupé, P., Concha, L., Buades, A., Collins, D. L., and Robles, M. (2013). Diffusion weighted image denoising using overcomplete local PCA. *PLoS One* 8:e73021. doi: 10.1371/journal.pone.0073021

Marques, J. P., Kober, T., Krueger, G., van der Zwaag, W., Van de Moortele, P.-F., and Gruetter, R. (2010). MP2RAGE, a self bias-field corrected sequence for improved segmentation and T_{1}-mapping at high field. *NeuroImage* 49, 1271–1281. doi: 10.1016/j.neuroimage.2009.10.002

Metere, R., Kober, T., Möller, H. E., and Schäfer, A. (2017). Simultaneous quantitative MRI mapping of T1, T2^{∗} and magnetic susceptibility with multi-echo MP2RAGE. *PLoS One* 12:e0169265. doi: 10.1371/journal.pone.0169265

Milanfar, P. (2013). A tour of modern image filtering: new insights and methods, both practical and theoretical. *IEEE Signal Process. Mag.* 30, 106–128. doi: 10.1109/MSP.2011.2179329

Miller, A. J., and Joseph, P. M. (1993). The use of power images to perform quantitative analysis on low SNR MR images. *Magn. Reson. Imaging* 11, 1051–1056. doi: 10.1016/0730-725x(93)90225-3

Strotmann, B., Kögler, C., Bazin, P.-L., Weiss, M., Villringer, A., and Turner, R. (2013). Mapping of the internal structure of human habenula with ex vivo MRI at 7T. *Front. Hum. Neurosci.* 7:878. doi: 10.3389/fnhum.2013.00878

Stucht, D., Danishad, K. A., Schulze, P., Godenschweger, F., Zaitsev, M., and Speck, O. (2015). Highest resolution in vivo human brain MRI using prospective motion correction. *PLoS One* 10:e0133921. doi: 10.1371/journal.pone.0133921

Veraart, J., Novikov, D. S., Christiaens, D., Ades-aron, B., Sijbers, J., and Fieremans, E. (2016). Denoising of diffusion MRI using random matrix theory. *NeuroImage* 142, 394–406. doi: 10.1016/j.neuroimage.2016.08.016

Weiskopf, N., Mohammadi, S., Lutti, A., and Callaghan, M. F. (2015). Advances in MRI-based computational neuroanatomy: from morphometry to *in-vivo* histology. *Curr. Opin. Neurol.* 28, 313–322. doi: 10.1097/WCO.0000000000000222

Keywords: denoising, ultra-high field MRI, quantitative MRI, local PCA, complex MRI signal

Citation: Bazin P-L, Alkemade A, van der Zwaag W, Caan M, Mulder M and Forstmann BU (2019) Denoising High-Field Multi-Dimensional MRI With Local Complex PCA. *Front. Neurosci.* 13:1066. doi: 10.3389/fnins.2019.01066

Received: 15 July 2019; Accepted: 24 September 2019;

Published: 09 October 2019.

Edited by:

John Ashburner, University College London, United KingdomReviewed by:

Karsten Tabelow, Weierstrass Institute for Applied Analysis and Stochastics (LG), GermanySuyash P. Awate, Indian Institute of Technology Bombay, India

Copyright © 2019 Bazin, Alkemade, van der Zwaag, Caan, Mulder and Forstmann. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Pierre-Louis Bazin, pilou.bazin@uva.nl; bazin@cbs.mpg.de